home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 9 / 09.iso / l / l042 / 1.ddi / CHAP2.ARC / ROOTSEQU.INC < prev    next >
Encoding:
Text File  |  1987-12-30  |  63.0 KB  |  1,584 lines
  1. {----------------------------------------------------------------------------}
  2. {-                                                                          -}
  3. {-     Turbo Pascal Numerical Methods Toolbox                               -}
  4. {-     Copyright (c) 1986, 87 by Borland International, Inc.                -}
  5. {-                                                                          -}
  6. {----------------------------------------------------------------------------}
  7.  
  8. procedure Bisect{(LeftEnd : Float;
  9.                  RightEnd : Float;
  10.                  Tol      : Float;
  11.                  MaxIter  : integer;
  12.              var Answer   : Float;
  13.              var fAnswer  : Float;
  14.              var Iter     : integer;
  15.              var Error    : byte;
  16.                  FuncPtr  : Pointer)};
  17.  
  18. var
  19.   Found : boolean;
  20.  
  21. procedure TestInput(LeftEndpoint  : Float;
  22.                     RightEndpoint : Float;
  23.                     Tol           : Float;
  24.                     MaxIter       : integer;
  25.                 var Answer        : Float;
  26.                 var fAnswer       : Float;
  27.                 var Error         : byte;
  28.                 var Found         : boolean);
  29.  
  30. {--------------------------------------------------------------}
  31. {- Input:  LeftEndpoint, RightEndpoint, Tol, MaxIter          -}
  32. {- Output: Answer, fAnswer, Error, Found                      -}
  33. {-                                                            -}
  34. {- This procedure tests the input data for errors.  If        -}
  35. {- LeftEndpoint > RightEndpoint, Tol <= 0, or MaxIter < 0,    -}
  36. {- then an error is returned.  If one the of the endpoints    -}
  37. {- (LeftEndpoint, RightEndpoint) is a root, then Found = TRUE -}
  38. {- and Answer and fAnswer are returned.                       -}
  39. {--------------------------------------------------------------}
  40.  
  41. var
  42.   yLeft, yRight : Float; { The values of the function at the endpoints. }
  43.  
  44. begin
  45.   yLeft := UserFunction(LeftEndpoint, FuncPtr);
  46.   yRight := UserFunction(RightEndpoint, FuncPtr);
  47.  
  48.   if ABS(yLeft) <= TNNearlyZero then
  49.   begin
  50.     Answer := LeftEndpoint;
  51.     fAnswer := UserFunction(Answer, FuncPtr);
  52.     Found := true;
  53.   end;
  54.  
  55.   if ABS(yRight) <= TNNearlyZero then
  56.   begin
  57.     Answer := RightEndpoint;
  58.     fAnswer := UserFunction(Answer, FuncPtr);
  59.     Found := true;
  60.   end;
  61.  
  62.   if not Found then   { Test for errors }
  63.   begin
  64.     if yLeft * yRight > 0 then
  65.       Error := 2;
  66.     if Tol <= 0 then
  67.       Error := 3;
  68.     if MaxIter < 0 then
  69.       Error := 4;
  70.   end;
  71. end; { procedure Tests }
  72.  
  73. procedure Converge(var LeftEndpoint  : Float;
  74.                    var RightEndpoint : Float;
  75.                        Tol           : Float;
  76.                    var Found         : boolean;
  77.                        MaxIter       : integer;
  78.                    var Answer        : Float;
  79.                    var fAnswer       : Float;
  80.                    var Iter          : integer;
  81.                    var Error         : byte);
  82.  
  83. {------------------------------------------------------------------}
  84. {- Input:  LeftEndpoint, RightEndpoint, Tol, MaxIter              -}
  85. {- Output: Found, Answer, fAnswer, Iter, Error                    -}
  86. {-                                                                -}
  87. {- This procedure applies the bisection method to find a          -}
  88. {- root to TNTargetF(x) on the interval [LeftEndpoint,            -}
  89. {- RightEndpoint]. The root must be found within MaxIter          -}
  90. {- iterations to a tolerance of Tol.  If a root is found, then it -}
  91. {- is returned in Answer, the value of the function at the        -}
  92. {- approximated root is returned in fAnswer (should be close to   -}
  93. {- zero), and the number of iterations is returned in Iter.       -}
  94. {------------------------------------------------------------------}
  95.  
  96. var
  97.   yLeft : Float;
  98.   MidPoint, yMidPoint : Float;
  99.  
  100. procedure Initial(LeftEndpoint : Float;
  101.               var Iter         : integer;
  102.               var yLeft        : Float);
  103. { Initialize variables. }
  104. begin
  105.   Iter := 0;
  106.   yLeft := UserFunction(LeftEndpoint, FuncPtr);
  107. end; { procedure Initial }
  108.  
  109. function TestForRoot(X, OldX, Y, Tol : Float) : boolean;
  110. {--------------------------------------------------------------------}
  111. {-  These are the stopping criteria.  Four different ones are       -}
  112. {-  provided.  If you wish to change the active criteria, simply    -}
  113. {-  comment off the current criteria (including the appropriate OR) -}
  114. {-  and remove the comment brackets from the criteria (including    -}
  115. {-  the appropriate OR) you wish to be active.                      -}
  116. {--------------------------------------------------------------------}
  117. begin
  118.   TestForRoot :=                      {---------------------------}
  119.     (ABS(Y) <= TNNearlyZero)          {- Y = 0                   -}
  120.                                       {-                         -}
  121.            or                         {-                         -}
  122.                                       {-                         -}
  123.     (ABS(X - OldX) < ABS(OldX * Tol)) {- Relative change in X    -}
  124.                                       {-                         -}
  125.                                       {-                         -}
  126. (*         or                      *) {-                         -}
  127. (*                                 *) {-                         -}
  128. (*  (ABS(X - OldX) < Tol)          *) {- Absolute change in X    -}
  129. (*                                 *) {-                         -}
  130. (*         or                      *) {-                         -}
  131. (*                                 *) {-                         -}
  132. (*  (ABS(Y) <= Tol)                *) {- Absolute change in Y    -}
  133.                                       {---------------------------}
  134.  
  135. {-----------------------------------------------------------------------}
  136. {- The first criteria simply checks to see if the value of the         -}
  137. {- function is zero.  You should probably always keep this criteria    -}
  138. {- active.                                                             -}
  139. {-                                                                     -}
  140. {- The second criteria checks the relative error in x. This criteria   -}
  141. {- evaluates the fractional change in x between interations. Note      -}
  142. {- that x has been multiplied through the inequality to avoid divide   -}
  143. {- by zero errors.                                                     -}
  144. {-                                                                     -}
  145. {- The third criteria checks the absolute difference in x between      -}
  146. {- iterations.                                                         -}
  147. {-                                                                     -}
  148. {- The fourth criteria checks the absolute difference between          -}
  149. {- the value of the function and zero.                                 -}
  150. {-----------------------------------------------------------------------}
  151.  
  152. end; { function TestForRoot }
  153.  
  154. begin { procedure Converge }
  155.   Initial(LeftEndpoint, Iter, yLeft);
  156.   while not(Found) and (Iter < MaxIter) do
  157.   begin
  158.     Iter := Succ(Iter);
  159.     MidPoint := (LeftEndpoint + RightEndpoint) / 2;
  160.     yMidPoint := UserFunction(MidPoint, FuncPtr);
  161.     Found := TestForRoot(MidPoint, LeftEndpoint, yMidPoint, Tol);
  162.     if (yLeft * yMidPoint) < 0 then
  163.        RightEndpoint := MidPoint
  164.     else
  165.       begin
  166.         LeftEndpoint := MidPoint;
  167.         yLeft := yMidPoint;
  168.       end;
  169.   end;
  170.   Answer := MidPoint;
  171.   fAnswer := yMidPoint;
  172.   if Iter >= MaxIter then
  173.     Error := 1;
  174. end; { procedure Converge }
  175.  
  176. begin { procedure Bisect }
  177.   Found := false;
  178.   TestInput(LeftEnd, RightEnd, Tol, MaxIter, Answer, fAnswer, Error, Found);
  179.   if (Error = 0) and (Found = false) then  { i.e. no error }
  180.     Converge(LeftEnd, RightEnd, Tol, Found, MaxIter,
  181.              Answer, fAnswer, Iter, Error);
  182. end; { procedure Bisect }
  183.  
  184. procedure Newton_Raphson{(Guess   : Float;
  185.                          Tol      : Float;
  186.                          MaxIter  : integer;
  187.                      var Root     : Float;
  188.                      var Value    : Float;
  189.                      var Deriv    : Float;
  190.                      var Iter     : integer;
  191.                      var Error    : byte;
  192.                          FuncPtr1 : Pointer;
  193.                          FuncPtr2 : Pointer)};
  194.  
  195. var
  196.   Found : boolean;              { Flags that a root has been Found }
  197.   OldX, OldY, OldDeriv,
  198.   NewX, NewY, NewDeriv : Float;  { Iteration variables }
  199.  
  200. procedure CheckSlope(Slope : Float;
  201.                  var Error : byte);
  202.  
  203. {---------------------------------------------------}
  204. {- Input:  Slope                                   -}
  205. {- Output: Error                                   -}
  206. {-                                                 -}
  207. {- This procedure checks the slope to see if it is -}
  208. {- zero.  The Newton Raphson algorithm may not be  -}
  209. {- applied at a point where the slope is zero.     -}
  210. {---------------------------------------------------}
  211.  
  212. begin
  213.   if ABS(Slope) <= TNNearlyZero then
  214.     Error := 2;     { Slope is zero }
  215. end; { procedure CheckSlope }
  216.  
  217. procedure Initial(Guess    : Float;
  218.                   Tol      : Float;
  219.                   MaxIter  : integer;
  220.               var OldX     : Float;
  221.               var OldY     : Float;
  222.               var OldDeriv : Float;
  223.               var Found    : boolean;
  224.               var Iter     : integer;
  225.               var Error    : byte);
  226.  
  227.  
  228. {-------------------------------------------------------------}
  229. {- Input:  Guess, Tol, MaxIter                               -}
  230. {- Output: OldX, OldY, OldDeriv, Found, Iter, Error          -}
  231. {-                                                           -}
  232. {- This procedure sets the initial values of the above       -}
  233. {- variables. If OldY is zero, then a root has been          -}
  234. {- Found and Found = TRUE.  This procedure also checks       -}
  235. {- the tolerance (Tol) and the maximum number of iterations  -}
  236. {- (MaxIter) for errors.                                     -}
  237. {-------------------------------------------------------------}
  238.  
  239. begin
  240.   Found := false;
  241.   Iter := 0;
  242.   Error := 0;
  243.   OldX := Guess;
  244.   OldY := UserFunction(OldX, FuncPtr1);
  245.   OldDeriv := UserFunction(OldX, FuncPtr2);
  246.   if ABS(OldY) <= TNNearlyZero then
  247.     Found := true
  248.   else
  249.     CheckSlope(OldDeriv, Error);
  250.   if Tol <= 0 then
  251.     Error := 3;
  252.   if MaxIter < 0 then
  253.     Error := 4;
  254. end; { procedure Initial }
  255.  
  256. function TestForRoot(X, OldX, Y, Tol : Float) : boolean;
  257.  
  258. {-----------------------------------------------------------------}
  259. {- These are the stopping criteria.  Four different ones are     -}
  260. {- provided.  If you wish to change the active criteria, simply  -}
  261. {- comment off the current criteria (including the preceding OR) -}
  262. {- and remove the comment brackets from the criteria (including  -}
  263. {- the following OR) you wish to be active.                      -}
  264. {-----------------------------------------------------------------}
  265.  
  266.     begin
  267.       TestForRoot :=                      {---------------------------}
  268.         (ABS(Y) <= TNNearlyZero)          {- Y = 0                   -}
  269.                                           {-                         -}
  270.                or                         {-                         -}
  271.                                           {-                         -}
  272.         (ABS(X - OldX) < ABS(OldX*Tol))   {- Relative change in X    -}
  273.                                           {-                         -}
  274. (*             or                      *) {-                         -}
  275. (*                                     *) {-                         -}
  276. (*      (ABS(OldX - X) < Tol)          *) {- Absolute change in X    -}
  277. (*                                     *) {-                         -}
  278. (*             or                      *) {-                         -}
  279. (*                                     *) {-                         -}
  280. (*      (ABS(Y) <= Tol)                *) {- Absolute change in Y    -}
  281.                                           {---------------------------}
  282.  
  283. {-----------------------------------------------------------------------}
  284. {- The first criteria simply checks to see if the value of the         -}
  285. {- function is zero.  You should probably always keep this criteria    -}
  286. {- active.                                                             -}
  287. {-                                                                     -}
  288. {- The second criteria checks the relative error in X. This criteria   -}
  289. {- evaluates the fractional change in X between interations. Note      -}
  290. {- that X has been multiplied throught the inequality to avoid divide  -}
  291. {- by zero errors.                                                     -}
  292. {-                                                                     -}
  293. {- The third criteria checks the absolute difference in X between      -}
  294. {- iterations.                                                         -}
  295. {-                                                                     -}
  296. {- The fourth criteria checks the absolute difference between          -}
  297. {- the value of the function and zero.                                 -}
  298. {-----------------------------------------------------------------------}
  299.  
  300. end; { procedure TestForRoot }
  301.  
  302. begin { procedure Newton_Raphson }
  303.   Initial(Guess, Tol, MaxIter, OldX, OldY, OldDeriv,
  304.           Found, Iter, Error);
  305.  
  306.   while not(Found) and (Error = 0) and (Iter<MaxIter) do
  307.   begin
  308.     Iter := Succ(Iter);
  309.     NewX := OldX - OldY / OldDeriv;
  310.     NewY := UserFunction(NewX, FuncPtr1);
  311.     NewDeriv := UserFunction(NewX, FuncPtr2);
  312.     Found := TestForRoot(NewX, OldX, NewY, Tol);
  313.     OldX := NewX;
  314.     OldY := NewY;
  315.     OldDeriv := NewDeriv;
  316.     if not(Found) then
  317.       CheckSlope(OldDeriv, Error);
  318.   end;
  319.   Root := OldX;
  320.   Value := OldY;
  321.   Deriv := OldDeriv;
  322.   if not(Found) and (Error = 0) and (Iter >= MaxIter) then
  323.     Error := 1;
  324. end; { procedure Newton_Raphson }
  325.  
  326. procedure Secant{(Guess1 : Float;
  327.                  Guess2  : Float;
  328.                  Tol     : Float;
  329.                  MaxIter : integer;
  330.              var Root    : Float;
  331.              var Value   : Float;
  332.              var Iter    : integer;
  333.              var Error   : byte;
  334.                  FuncPtr : Pointer)};
  335.  
  336. var
  337.   Found : boolean;            { Flags that a root has been Found }
  338.   OldX, OldY, X, Y,
  339.   NewX, NewY  : Float;         { Iteration variables }
  340.  
  341. procedure Initial(Guess1  : Float;
  342.                   Guess2  : Float;
  343.                   Tol     : Float;
  344.                   MaxIter : integer;
  345.               var OldX    : Float;
  346.               var X       : Float;
  347.               var OldY    : Float;
  348.               var Y       : Float;
  349.               var Found   : boolean;
  350.               var Iter    : integer;
  351.               var Error   : byte);
  352.  
  353. {-------------------------------------------------------------}
  354. {- Input:  Guess1, Guess2, Tol, MaxIter                      -}
  355. {- Output: OldX, X, OldY, Y, Found, Iter, Error              -}
  356. {-                                                           -}
  357. {- This procedure sets the initial values of the above       -}
  358. {- variables. If OldY or Y is zero, then a root has been     -}
  359. {- Found and Found = TRUE.  This procedure also checks       -}
  360. {- the tolerance (Tol) and the maximum number of iterations  -}
  361. {- (MaxIter) for errors.                                     -}
  362. {-------------------------------------------------------------}
  363.  
  364. begin
  365.   Found := false;
  366.   Iter := 0;
  367.   Error := 0;
  368.   OldX := Guess1;
  369.   X := Guess2;
  370.   OldY := UserFunction(OldX, FuncPtr);
  371.   Y := UserFunction(X, FuncPtr);
  372.   if ABS(OldY) <= TNNearlyZero then         { OldX is the root }
  373.     begin
  374.       X := OldX;
  375.       Y := OldY;
  376.       Found := true;
  377.     end
  378.   else
  379.     if ABS(Y) <= TNNearlyZero then
  380.       Found := true                         {  X is the root  }
  381.     else
  382.       if ABS(OldY - Y) <= TNNearlyZero then
  383.         Error := 2;                  { Slope of line is zero; no intercept }
  384.   if Tol <= 0 then
  385.     Error := 3;
  386.   if MaxIter < 0 then
  387.     Error := 4;
  388. end; { procedure Initial }
  389.  
  390. function TestForRoot(X, OldX, y, Tol : Float) : boolean;
  391.  
  392. {-----------------------------------------------------------------}
  393. {- These are the stopping criteria.  Four different ones are     -}
  394. {- provided.  If you wish to change the active criteria, simply  -}
  395. {- comment off the current criteria (including the preceding OR) -}
  396. {- and remove the comment brackets from the criteria (including  -}
  397. {- the following OR) you wish to be active.                      -}
  398. {-----------------------------------------------------------------}
  399.  
  400. begin
  401.   TestForRoot :=                      {---------------------------}
  402.     (ABS(y) <= TNNearlyZero)          {- Y = 0                   -}
  403.                                       {-                         -}
  404.            or                         {-                         -}
  405.                                       {-                         -}
  406.     (ABS(X - OldX) < ABS(OldX*Tol))   {- Relative change in X    -}
  407.                                       {-                         -}
  408. (*         or                      *) {-                         -}
  409. (*                                 *) {-                         -}
  410. (*  (ABS(OldX - X) < Tol)          *) {- Absolute change in X    -}
  411. (*                                 *) {-                         -}
  412. (*         or                      *) {-                         -}
  413. (*                                 *) {-                         -}
  414. (*  (ABS(y) <= Tol)                *) {- Absolute change in Y    -}
  415.                                       {---------------------------}
  416.  
  417. {-----------------------------------------------------------------------}
  418. {- The first criteria simply checks to see if the value of the         -}
  419. {- function is zero.  You should probably always keep this criteria    -}
  420. {- active.                                                             -}
  421. {-                                                                     -}
  422. {- The second criteria checks the relative error in X. This criteria   -}
  423. {- evaluates the fractional change in X between interations. Note      -}
  424. {- that X has been multiplied through the inequality to avoid divide   -}
  425. {- by zero errors.                                                     -}
  426. {-                                                                     -}
  427. {- The third criteria checks the absolute difference in X between      -}
  428. {- iterations.                                                         -}
  429. {-                                                                     -}
  430. {- The fourth criteria checks the absolute difference between          -}
  431. {- the value of the function and zero.                                 -}
  432. {-----------------------------------------------------------------------}
  433.  
  434. end; { procedure TestForRoot }
  435.  
  436. begin { procedure Secant }
  437.   Initial(Guess1, Guess2, Tol, MaxIter, OldX, X, OldY, Y,
  438.           Found, Iter, Error);
  439.   while not(Found) and (Error = 0) and (Iter<MaxIter) do
  440.   begin
  441.     Iter := Succ(Iter);
  442.     NewX := X - Y * (X - OldX) / (Y - OldY);    { The secant formula }
  443.     NewY := UserFunction(NewX, FuncPtr);
  444.     Found := TestForRoot(NewX, OldX, NewY, Tol);
  445.     OldX := X;
  446.     OldY := Y;
  447.     X := NewX;
  448.     Y := NewY;
  449.   end;
  450.   Root := X;
  451.   Value := Y;
  452.   if not(Found) and (Error = 0) and (Iter >= MaxIter) then
  453.     Error := 1;
  454. end; { procedure Secant }
  455.  
  456. procedure Newt_Horn_Defl{(InitDegree : integer;
  457.                           InitPoly   : TNvector;
  458.                           Guess      : Float;
  459.                           Tol        : Float;
  460.                           MaxIter    : integer;
  461.                       var Degree     : integer;
  462.                       var NumRoots   : integer;
  463.                       var Poly       : TNvector;
  464.                       var Root       : TNvector;
  465.                       var Imag       : TNvector;
  466.                       var Value      : TNvector;
  467.                       var Deriv      : TNvector;
  468.                       var Iter       : TNIntVector;
  469.                       var Error      : byte)};
  470.  
  471. var
  472.   Iter1 : integer;
  473.  
  474. procedure SynDiv(Degree  : integer;
  475.                  Poly    : TNvector;
  476.                  X       : Float;
  477.              var NewPoly : TNvector);
  478.  
  479. {--------------------------------------------------------------------------}
  480. {- Input:  Degree, Poly                                                   -}
  481. {- Output: NewPoly                                                        -}
  482. {-                                                                        -}
  483. {- This procedure applies the technique of synthetic division             -}
  484. {- to a polynomial, Poly, at the value X.  The kth element of NewPoly     -}
  485. {- is the (k-1)th element of the new polynomial. The 0th element          -}
  486. {- of the resulting polynomial, NewPoly, is the value of the polynomial   -}
  487. {- at X                                                                   -}
  488. {--------------------------------------------------------------------------}
  489.  
  490. var
  491.   Term : integer;
  492.  
  493. begin
  494.   NewPoly[Degree] := Poly[Degree];
  495.   for Term := Degree - 1 downto 0 do
  496.     NewPoly[Term] := NewPoly[Term + 1] * X + Poly[Term];
  497. end; { procedure SynDiv }
  498.  
  499. procedure QuadraticFormula(A       : Float;
  500.                            B       : Float;
  501.                            C       : Float;
  502.                        var ReRoot1 : Float;
  503.                        var ImRoot1 : Float;
  504.                        var ReRoot2 : Float;
  505.                        var ImRoot2 : Float);
  506.  
  507. {----------------------------------------------------------------}
  508. {- Input:  A, B, C                                              -}
  509. {- Output: ReRoot1, ImRoot1, ReRoot2, ImRoot2                   -}
  510. {-                                                              -}
  511. {- This procedure applies the quadratic formula to the equation -}
  512. {- AX^2 + BX + C, where A,B,C are real.  It returns the real or -}
  513. {- complex roots of the equation                                -}
  514. {----------------------------------------------------------------}
  515.  
  516. var
  517.   Denominator, Discrim : Float;
  518.  
  519. begin
  520.   Denominator := 2 * A;
  521.   ReRoot1 := -B / Denominator;
  522.   ReRoot2 := ReRoot1;
  523.   Discrim := B * B - 4 * A * C;
  524.   if Discrim < 0 then
  525.     begin
  526.       ImRoot1 := -Sqrt(-Discrim) / Denominator;
  527.       ImRoot2 :=  Sqrt(-Discrim) / Denominator;
  528.     end
  529.   else
  530.     begin
  531.       if B < 0 then  { Choose ReRoot1 to have the greatest absolute value }
  532.         ReRoot1 := ReRoot1 + Sqrt(Discrim) / Denominator
  533.       else
  534.         ReRoot1 := ReRoot1 - Sqrt(Discrim) / Denominator;
  535.       ReRoot2 := C / (A * ReRoot1);   { The product of the 2 roots is C/A }
  536.       ImRoot1 := 0;
  537.       ImRoot2 := 0;
  538.     end;
  539. end;    { procedure QuadraticFormula }
  540.  
  541. procedure TestData(InitDegree : integer;
  542.                    InitPoly   : TNvector;
  543.                    Tol        : Float;
  544.                    MaxIter    : integer;
  545.                var Degree     : integer;
  546.                var Poly       : TNvector;
  547.                var NumRoots   : integer;
  548.                var Roots      : TNvector;
  549.                var yRoots     : TNvector;
  550.                var Iter       : TNIntVector;
  551.                var Error      : byte);
  552.  
  553. {----------------------------------------------------------}
  554. {- Input:  InitDegree, InitPoly, Tol, MaxIter             -}
  555. {- Output: Degree, Poly, NumRoots, Roots, yRoots,         -}
  556. {-         Iter, Error                                    -}
  557. {-                                                        -}
  558. {- This procedure sets the initial value of the above     -}
  559. {- variables.  This procedure also tests the tolerance    -}
  560. {- (Tol), maximum number of iterations (MaxIter), and     -}
  561. {- Degree for errors and returns the appropriate error    -}
  562. {- code.  Finally, it examines the coefficients of Poly.  -}
  563. {- If the constant term is zero, then zero is one of the  -}
  564. {- roots and the polynomial is deflated accordingly. Also -}
  565. {- if the leading coefficient is zero, then Degree is     -}
  566. {- reduced until the leading coefficient is non-zero.     -}
  567. {----------------------------------------------------------}
  568.  
  569. var
  570.   Term : integer;
  571.  
  572. begin
  573.   Error := 0;
  574.   NumRoots := 0;
  575.   Degree := InitDegree;
  576.   Poly := InitPoly;
  577.   if Tol <= 0 then
  578.     Error := 4;
  579.   if MaxIter < 0 then
  580.     Error :=5;
  581.   { Reduce Degree until leading coefficient <> zero }
  582.   while (Degree > 0) and (ABS(Poly[Degree]) < TNNearlyZero) do
  583.     { Reduce Degree until leading coefficient <> zero }
  584.     Degree  :=  Pred(Degree);
  585.   if Degree <= 0 then
  586.     Error := 3;
  587.   { Deflate polynomial until the constant term <> zero }
  588.   while (ABS(Poly[0]) < TNNearlyZero) and (Degree > 0) do
  589.   begin
  590.     NumRoots := Succ(NumRoots);
  591.     Roots[NumRoots] := 0;
  592.     yRoots[NumRoots] := 0;
  593.     Iter[NumRoots] := 0;
  594.     Degree := Pred(Degree);
  595.     for Term := 0 to Degree do
  596.       Poly[Term] := Poly[Term + 1];
  597.   end;
  598. end; { procedure TestData }
  599.  
  600. procedure FindValueAndDeriv(Degree : integer;
  601.                             Poly   : TNvector;
  602.                             X      : Float;
  603.                         var Value  : Float;
  604.                         var Deriv  : Float);
  605.  
  606. {--------------------------------------------------------------------}
  607. {- Input:  Degree, Poly, X                                          -}
  608. {- Output: Value, Deriv                                             -}
  609. {-                                                                  -}
  610. {- This procedure applies the technique of synthetic division to    -}
  611. {- determine both the Value and derivative of the polynomial, Poly, -}
  612. {- at X.  The 0th element of the first synthetic division is the    -}
  613. {- value of the polynomial at X, and the 1st element of the second  -}
  614. {- synthetic division is the derivative of the polynomial at X.     -}
  615. {--------------------------------------------------------------------}
  616.  
  617. var
  618.   Poly1, Poly2 : TNvector;
  619. begin
  620.   SynDiv(Degree, Poly, X, Poly1);
  621.   Value := Poly1[0];
  622.   SynDiv(Degree, Poly1, X, Poly2);
  623.   Deriv := Poly2[1];
  624. end; { procedure FindValueAndDeriv }
  625.  
  626. procedure FindOneRoot(Degree  : integer;
  627.                       Poly    : TNvector;
  628.                       Guess   : Float;
  629.                       Tol     : Float;
  630.                       MaxIter : integer;
  631.                   var Root    : Float;
  632.                   var Value   : Float;
  633.                   var Deriv   : Float;
  634.                   var Iter    : integer;
  635.                   var Error   : byte);
  636.  
  637. {-------------------------------------------------------------------}
  638. {- Input:  Degree, Poly, Guess, Tol, MaxIter                       -}
  639. {- Output: Root, Value, Deriv, Iter, Error                         -}
  640. {-                                                                 -}
  641. {- A single root of the polynomial Poly.  The root must be         -}
  642. {- approximated within MaxIter iterations to a tolerance of Tol.   -}
  643. {- The root, value of the polynomial at the root (Value), and the  -}
  644. {- value of the derivative of the polynomial at the root (Deriv),  -}
  645. {- and the number of iterations (Iter) are returned. If no root    -}
  646. {- is found, the appropriate error code (Error) is returned.       -}
  647. {-------------------------------------------------------------------}
  648.  
  649. var
  650.   Found : boolean;
  651.   OldX, OldY, OldDeriv,
  652.   NewX, NewY, NewDeriv : Float;
  653.  
  654. procedure CheckSlope(Slope : Float;
  655.                  var Error : byte);
  656.  
  657. {---------------------------------------------------}
  658. {- Input:  Slope                                   -}
  659. {- Output: Error                                   -}
  660. {-                                                 -}
  661. {- This procedure checks the slope to see if it is -}
  662. {- zero.  The Newton Raphson algorithm may not be  -}
  663. {- applied at a point where the slope is zero.     -}
  664. {---------------------------------------------------}
  665.  
  666. begin
  667.   if ABS(Slope) <= TNNearlyZero then
  668.     Error := 2;
  669. end; { procedure CheckSlope }
  670.  
  671. procedure Initial(Degree   : integer;
  672.                   Poly     : TNvector;
  673.                   Guess    : Float;
  674.               var OldX     : Float;
  675.               var OldY     : Float;
  676.               var OldDeriv : Float;
  677.               var Found    : boolean;
  678.               var Iter     : integer;
  679.               var Error    : byte);
  680.  
  681. {-------------------------------------------------------------}
  682. {- Input:  Degree, Poly, Guess                               -}
  683. {- Output: OldX, OldY, OldDeriv, Found, Iter, Error          -}
  684. {-                                                           -}
  685. {- This procedure sets the initial values of the above       -}
  686. {- variables. If OldY is zero, then a root has been          -}
  687. {- found and Found = TRUE.                                   -}
  688. {-------------------------------------------------------------}
  689.  
  690. begin
  691.   Found := false;
  692.   Iter := 0;
  693.   Error := 0;
  694.   OldX := Guess;
  695.   FindValueAndDeriv(Degree, Poly, OldX, OldY, OldDeriv);
  696.   if ABS(OldY) <= TNNearlyZero then
  697.     Found := true
  698.   else
  699.     CheckSlope(OldDeriv, Error);
  700. end; { procedure Initial }
  701.  
  702. function TestForRoot(X, OldX, Y, Tol : Float) : boolean;
  703.  
  704. {----------------------------------------------------------------}
  705. { These are the stopping criteria.  Four different ones are     -}
  706. { provided.  If you wish to change the active criteria, simply  -}
  707. { comment off the current criteria (including the preceding OR) -}
  708. { and remove the comment brackets from the criteria (including  -}
  709. { the following OR) you wish to be active.                      -}
  710. {----------------------------------------------------------------}
  711.  
  712. begin
  713.   TestForRoot :=                      {---------------------------}
  714.     (ABS(Y) <= TNNearlyZero)          {- Y=0                     -}
  715.                                       {-                         -}
  716.            or                         {-                         -}
  717.                                       {-                         -}
  718.     (ABS(X - OldX) < ABS(OldX*Tol))   {- Relative change in X    -}
  719.                                       {-                         -}
  720.  (*        or                      *) {-                         -}
  721.  (*                                *) {-                         -}
  722.  (* (ABS(OldX - X) < Tol)          *) {- Absolute change in X    -}
  723.  (*                                *) {-                         -}
  724.  (*        or                      *) {-                         -}
  725.  (*                                *) {-                         -}
  726.  (* (ABS(Y) <= Tol)                *) {- Absolute change in Y    -}
  727.                                       {---------------------------}
  728.  
  729. {-----------------------------------------------------------------------}
  730. {- The first criteria simply checks to see if the value of the         -}
  731. {- function is zero.  You should probably always keep this criteria    -}
  732. {- active.                                                             -}
  733. {-                                                                     -}
  734. {- The second criteria checks the relative error in X. This criteria   -}
  735. {- evaluates the fractional change in X between interations. Note      -}
  736. {- that X has been multiplied through the inequality to avoid divide   -}
  737. {- by zero errors.                                                     -}
  738. {-                                                                     -}
  739. {- The third criteria checks the absolute difference in X between      -}
  740. {- iterations.                                                         -}
  741. {-                                                                     -}
  742. {- The fourth criteria checks the absolute difference between          -}
  743. {- the value of the function and zero.                                 -}
  744. {-----------------------------------------------------------------------}
  745.  
  746. end; { procedure TestForRoot }
  747.  
  748. begin { procedure FindOneRoot }
  749.   Initial(Degree, Poly, Guess, OldX, OldY, OldDeriv, Found, Iter, Error);
  750.   while not(Found) and (Error = 0) and (Iter<MaxIter) do
  751.   begin
  752.     Iter := Succ(Iter);
  753.     NewX := OldX - OldY / OldDeriv;
  754.     FindValueAndDeriv(Degree, Poly, NewX, NewY, NewDeriv);
  755.     Found := TestForRoot(NewX, OldX, NewY, Tol);
  756.     OldX := NewX;
  757.     OldY := NewY;
  758.     OldDeriv := NewDeriv;
  759.     if not(Found) then
  760.       CheckSlope(OldDeriv, Error);
  761.   end;
  762.   Root := OldX;
  763.   Value := OldY;
  764.   Deriv := OldDeriv;
  765.   if not(Found) and (Error = 0) and (Iter >= MaxIter) then
  766.     Error := 1;
  767.   if Found then
  768.     Error := 0;
  769. end; { procedure FindOneRoot }
  770.  
  771. procedure ReducePolynomial(var Degree : integer;
  772.                            var Poly   : TNvector;
  773.                                Root   : Float);
  774.  
  775. {------------------------------------------------------}
  776. {- Input:  Degree, Poly, Root                         -}
  777. {- Output: Degree, Poly                               -}
  778. {-                                                    -}
  779. {- This procedure deflates the polynomial Poly by     -}
  780. {- factoring out the Root.  Degree is reduced by one. -}
  781. {------------------------------------------------------}
  782.  
  783. var
  784.   NewPoly : TNvector;
  785.   Term : integer;
  786. begin
  787.   SynDiv(Degree, Poly, Root, NewPoly);
  788.   Degree := Pred(Degree);
  789.   for Term := 0 to Degree do
  790.     Poly[Term] := NewPoly[Term+1];
  791. end; { procedure ReducePolynomial }
  792.  
  793. begin  { procedure Newt_Horn_Defl }
  794.   TestData(InitDegree, InitPoly, Tol, MaxIter,
  795.            Degree, Poly, NumRoots, Root, Value, Iter, Error);
  796.   while (Error=0) and (Degree>2) do
  797.   begin
  798.     FindOneRoot(Degree, Poly, Guess, Tol, MaxIter, Root[NumRoots+1],
  799.                 Value[NumRoots+1], Deriv[NumRoots+1], Iter[NumRoots+1], Error);
  800.     if Error = 0 then
  801.     begin
  802.       NumRoots := Succ(NumRoots);
  803.       {------------------------------------------------------}
  804.       {- The next statement refines the approximate root by -}
  805.       {- plugging it into the original polynomial.  This    -}
  806.       {- eliminates a lot of the round-off error            -}
  807.       {- accumulated through many iterations                -}
  808.       {------------------------------------------------------}
  809.       if NumRoots > 1 then
  810.       begin
  811.         Iter1 := 0;
  812.         FindOneRoot(InitDegree, InitPoly, Root[NumRoots],
  813.                     Tol, MaxIter, Root[NumRoots], Value[NumRoots],
  814.                     Deriv[NumRoots], Iter1, Error);
  815.         Iter[NumRoots] := Iter[NumRoots] + Iter1;
  816.       end;
  817.       ReducePolynomial(Degree, Poly, Root[NumRoots]);
  818.       Guess := Root[NumRoots];
  819.     end;
  820.   end;
  821.   case Degree of
  822.     1 : begin           { Solve this linear }
  823.           Degree := 0;
  824.           NumRoots := Succ(NumRoots);
  825.           Root[NumRoots] := -Poly[0] / Poly[1];
  826.           FindOneRoot(InitDegree, InitPoly, Root[NumRoots], Tol,
  827.                       MaxIter, Root[NumRoots], Value[NumRoots],
  828.                       Deriv[NumRoots], Iter[NumRoots], Error);
  829.         end;
  830.  
  831.     2 : begin           { Solve this quadratic }
  832.           Degree := 0;
  833.           NumRoots := Succ(Succ(NumRoots));
  834.           QuadraticFormula(Poly[2], Poly[1], Poly[0],
  835.                            Root[NumRoots - 1], Imag[NumRoots - 1],
  836.                            Root[NumRoots], Imag[NumRoots]);
  837.           if ABS(Imag[NumRoots]) < TNNearlyZero then
  838.             { if the roots are real, they can be     }
  839.             { made more accurate using Newton-Horner }
  840.             begin
  841.               FindOneRoot(InitDegree, InitPoly, Root[NumRoots-1], Tol,
  842.                           MaxIter, Root[NumRoots-1], Value[NumRoots-1],
  843.                           Deriv[NumRoots-1], Iter[NumRoots-1], Error);
  844.  
  845.               FindOneRoot(InitDegree, InitPoly, Root[NumRoots], Tol,
  846.                           MaxIter, Root[NumRoots], Value[NumRoots],
  847.                           Deriv[NumRoots], Iter[NumRoots], Error);
  848.             end
  849.           else
  850.             { If the roots are complex, then assign    }
  851.             { the value to be zero (which is true,     }
  852.             { except for some roundoff error) and the  }
  853.             { derivative to be zero (which is usually  }
  854.             { FALSE; the derivative is usually complex }
  855.             begin
  856.               Value[NumRoots-1] := 0;    Value[NumRoots] := 0;
  857.               Deriv[NumRoots-1] := 0;    Deriv[NumRoots] := 0;
  858.               Iter[NumRoots-1] := 0;     Iter[NumRoots] := 0;
  859.             end;
  860.         end;
  861.   end; { case }
  862. end; { procedure Newt_Horn_Defl }
  863.  
  864. procedure Muller{(Guess  : TNcomplex;
  865.                  Tol     : Float;
  866.                  MaxIter : integer;
  867.              var Answer  : TNcomplex;
  868.              var yAnswer : TNcomplex;
  869.              var Iter    : integer;
  870.              var Error   : byte;
  871.                  FuncPtr : Pointer)};
  872.  
  873. type
  874.   TNquadratic = record
  875.                   A, B, C : TNcomplex;
  876.                 end;
  877.  
  878. var
  879.   X0, X1, OldApprox,
  880.   NewApprox, yNewApprox : TNcomplex;    { Iteration variables }
  881.   Factor : TNquadratic;                 { Factor of polynomial }
  882.   Found : boolean;                      { Flags that a factor }
  883.                                         { has been found }
  884.  
  885. {----------- Here are a few complex operations ------------}
  886.  
  887. procedure Conjugate(C1 : TNcomplex; var C2 : TNcomplex);
  888. begin
  889.   C2.Re := C1.Re;
  890.   C2.Im := -C1.Im;
  891. end; { procedure Conjugate }
  892.  
  893. function Modulus(var C1 : TNcomplex) : Float;
  894. begin
  895.   Modulus := Sqrt(Sqr(C1.Re) + Sqr(C1.Im));
  896. end; { function Modulus }
  897.  
  898. procedure Add(C1, C2 : TNcomplex; var C3 : TNcomplex);
  899. begin
  900.   C3.Re := C1.Re + C2.Re;
  901.   C3.Im := C1.Im + C2.Im;
  902. end; { procedure Add }
  903.  
  904. procedure Sub(C1, C2 : TNcomplex; var C3 : TNcomplex);
  905. begin
  906.   C3.Re := C1.Re - C2.Re;
  907.   C3.Im := C1.Im - C2.Im;
  908. end; { procedure Sub }
  909.  
  910. procedure Mult(C1, C2 : TNcomplex; var C3 : TNcomplex);
  911. begin
  912.   C3.Re := C1.Re * C2.Re - C1.Im * C2.Im;
  913.   C3.Im := C1.Im * C2.Re + C1.Re * C2.Im;
  914. end; { procedure Mult }
  915.  
  916. procedure Divide(C1, C2 : TNcomplex; var C3 : TNcomplex);
  917. var
  918.   Dum1, Dum2 : TNcomplex;
  919.   E : Float;
  920. begin
  921.   Conjugate(C2, Dum1);
  922.   Mult(C1, Dum1, Dum2);
  923.   E := Sqr(Modulus(C2));
  924.   C3.Re := Dum2.Re / E;
  925.   C3.Im := Dum2.Im / E;
  926. end; { procedure Divide }
  927.  
  928. procedure SquareRoot(C1 : TNcomplex; var C2 : TNcomplex);
  929. var
  930.   R, Theta : Float;
  931. begin
  932.   R := Sqrt(Sqr(C1.Re) + Sqr(C1.Im));
  933.   if ABS(C1.Re) < TNNearlyZero then
  934.     begin
  935.       if C1.Im < 0 then
  936.         Theta := Pi / 2
  937.       else
  938.         Theta := -Pi / 2;
  939.     end
  940.   else
  941.     if C1.Re < 0 then
  942.       Theta := ArcTan(C1.Im / C1.Re) + Pi
  943.     else
  944.       Theta := ArcTan(C1.Im / C1.Re);
  945.   C2.Re := Sqrt(R) * Cos(Theta / 2);
  946.   C2.Im := Sqrt(R) * Sin(Theta / 2);
  947. end; { procedure SquareRoot }
  948.  
  949.  
  950. procedure MakeParabola(X0     : TNcomplex;
  951.                        X1     : TNcomplex;
  952.                        Approx : TNcomplex;
  953.                    var Factor : TNquadratic;
  954.                    var Error  : byte);
  955.  
  956. {--------------------------------------------------------}
  957. {- Input:  X0, X1, Approx                               -}
  958. {- Output: Factor, Error                                -}
  959. {-                                                      -}
  960. {- This procedure constructs a parabola to fit the      -}
  961. {- three points X0, X1, Approx.  The intersection of    -}
  962. {- this parabola with the x-axis will yield the next    -}
  963. {- approximation.  If the parabola is a horizontal line -}
  964. {- then Error = 2 since a horizontal line will not      -}
  965. {- intersect the x-axis.                                -}
  966. {--------------------------------------------------------}
  967.  
  968. var
  969.   H1, H2, H3, H : TNcomplex;
  970.   Delta1, Delta2 : TNcomplex;
  971.   Dum1, Dum2, Dum3, Dum4 : TNcomplex;     { Dummy variables }
  972.  
  973. begin
  974.   with Factor do
  975.   begin
  976.     Sub(X0, Approx, H1);
  977.     Sub(X1, Approx, H2);
  978.     Sub(X0, X1, H3);
  979.     Mult(H2, H3, Dum1);
  980.     Mult(H1, Dum1, H);
  981.     if Modulus(H) < TNNearlyZero then
  982.       Error := 2;                     { Can't fit a quadratic to these points }
  983.     UserProcedure(X1, Dum1, FuncPtr);
  984.     Sub(Dum1, C, Delta1);             { C was passed in }
  985.     UserProcedure(X0, Dum1, FuncPtr);
  986.     Sub(Dum1, C, Delta2);
  987.  
  988.     if Error = 0 then                 { Calculate coefficients of quadratic }
  989.     begin
  990.       Mult(H1, H1, Dum1);             { Calculate B }
  991.       Mult(Dum1, Delta1, Dum2);
  992.       Mult(H2, H2, Dum1);
  993.       Mult(Dum1, Delta2, Dum3);
  994.       Sub(Dum2, Dum3, Dum4);
  995.       Divide(Dum4, H, B);
  996.  
  997.       Mult(H2, Delta2, Dum1);         { Calculate A }
  998.       Mult(H1, Delta1, Dum2);
  999.       Sub(Dum1, Dum2, Dum3);
  1000.       Divide(Dum3, H, A);
  1001.     end;
  1002.  
  1003.     if (Modulus(A) <= TNNearlyZero) and (Modulus(B) <= TNNearlyZero) then
  1004.       Error := 2;            { Test if parabola is actually a constant }
  1005.   end; { with }
  1006. end; { procedure MakeParabola }
  1007.  
  1008. procedure Initial(Guess      : TNcomplex;
  1009.                   Tol        : Float;
  1010.                   MaxIter    : integer;
  1011.               var Error      : byte;
  1012.               var Iter       : integer;
  1013.               var Found      : boolean;
  1014.               var NewApprox  : TNcomplex;
  1015.               var yNewApprox : TNcomplex;
  1016.               var X0         : TNcomplex;
  1017.               var X1         : TNcomplex;
  1018.               var OldApprox  : TNcomplex;
  1019.               var Factor     : TNquadratic);
  1020.  
  1021. {------------------------------------------------------------------}
  1022. {- Input:  Guess, Tol, MaxIter                                    -}
  1023. {- Output: Error, Iter, Found, NewApprox, yNewApprox,             -}
  1024. {-         X0, X1, OldApprox, Factor                              -}
  1025. {-                                                                -}
  1026. {- This procedure initializes all the above variables. It         -}
  1027. {- sets OldApprox equal to Guess, X0 and X1 are set close to      -}
  1028. {- Guess.  The procedure also checks the tolerance (Tol) and      -}
  1029. {- maximum number of iterations (MaxIter) for errors.             -}
  1030. {------------------------------------------------------------------}
  1031.  
  1032. var
  1033.   cZero, yOldApprox : TNcomplex;
  1034.  
  1035. begin
  1036.   cZero.Re := 0;    { Complex zero }
  1037.   cZero.Im := 0;    { Complex zero }
  1038.   Error := 0;
  1039.   Found := false;
  1040.   Iter := 0;
  1041.   NewApprox := cZero;
  1042.   yNewApprox := cZero;
  1043.           { X0 and X1 are points which are close to Guess }
  1044.   X0.Re := -0.75 * Guess.Re - 1;     X0.Im := Guess.Im;
  1045.   X1.Re := 0.75 * Guess.Re;          X1.Im := 1.2 * Guess.Im + 1;
  1046.   OldApprox := Guess;
  1047.   UserProcedure(OldApprox, yOldApprox, FuncPtr);  { Evaluate the function at OldApprox }
  1048.   Factor.A := cZero;
  1049.   Factor.B := cZero;
  1050.   Factor.C := yOldApprox;
  1051.   MakeParabola(X0, X1, OldApprox, Factor, Error);
  1052.   if Tol <= 0 then
  1053.     Error := 3;
  1054.   if MaxIter < 0 then
  1055.     Error := 4;
  1056. end; { procedure Initial }
  1057.  
  1058. procedure QuadraticFormula(Factor    : TNquadratic;
  1059.                            OldApprox : TNcomplex;
  1060.                        var NewApprox : TNcomplex);
  1061.  
  1062. {-------------------------------------------------------------}
  1063. {- Input:  Factor, OldApprox                                 -}
  1064. {- Output: NewApprox                                         -}
  1065. {-                                                           -}
  1066. {- This procedure applies the complex quadratic formula      -}
  1067. {- to the quadratic Factor to determine where the parabola   -}
  1068. {- represented by Factor intersects the x-axis. The solution -}
  1069. {- of the quadratic formula is subtracted from OldApprox to  -}
  1070. {- yield NewApprox.                                          -}
  1071. {-------------------------------------------------------------}
  1072.  
  1073. var
  1074.   Discrim, Difference, Dum1, Dum2, Dum3 : TNcomplex;
  1075.  
  1076. begin
  1077.   with Factor do
  1078.   begin
  1079.     Mult(B, B, Dum1);                     { B^2 }
  1080.     Mult(A, C, Dum2);
  1081.     Discrim.Re := Dum1.Re - 4 * Dum2.Re;  { B^2 - 4AC }
  1082.     Discrim.Im := Dum1.Im - 4 * Dum2.Im;
  1083.     SquareRoot(Discrim, Discrim);
  1084.     Sub(B, Discrim, Dum1);                { B +/- sqrt(B^2 - 4AC) }
  1085.     Add(B, Discrim, Dum2);
  1086.     { Choose the root with B +/- Discrim greatest }
  1087.     if Modulus(Dum1) < Modulus(Dum2) then
  1088.       Dum1 := Dum2;
  1089.     Add(C, C, Dum3);
  1090.     if Modulus(Dum1) < TNNearlyZero then  { if B +/- sqrt(B^2 - 4AC) = 0 }
  1091.       begin
  1092.         NewApprox.Re := 0;
  1093.         NewApprox.Im := 0;
  1094.       end
  1095.     else
  1096.       begin                               { 2C/[B +/- sqrt(B^2 - 4AC)] }
  1097.         Divide(Dum3, Dum1, Difference);
  1098.         { Calculate NewApprox }
  1099.         Sub(OldApprox, Difference, NewApprox);
  1100.       end;
  1101.   end; { with }
  1102. end; { procedure QuadraticFormula }
  1103.  
  1104. function TestForRoot(X    : TNcomplex;
  1105.                      OldX : TNcomplex;
  1106.                      Y    : TNcomplex;
  1107.                      Tol  : Float) : boolean;
  1108.  
  1109. var
  1110.   Dif, FracDif : TNcomplex;
  1111.  
  1112. {-----------------------------------------------------------------}
  1113. {- These are the stopping criteria.  Four different ones are     -}
  1114. {- provided.  If you wish to change the active criteria, simply  -}
  1115. {- comment off the current criteria (including the preceding OR) -}
  1116. {- and remove the comment brackets from the criteria (including  -}
  1117. {- the following OR) you wish to be active.                      -}
  1118. {-----------------------------------------------------------------}
  1119.  
  1120. begin
  1121.   Sub(X, OldX, Dif);
  1122.   FracDif.Re := X.Re * Tol;
  1123.   FracDif.Im := X.Im * Tol;
  1124.  
  1125.   TestForRoot :=                      {---------------------------}
  1126.     (Modulus(Y) <= TNNearlyZero)      {- Y = 0                   -}
  1127.                                       {-                         -}
  1128.            or                         {-                         -}
  1129.                                       {-                         -}
  1130.     (Modulus(Dif) < Modulus(FracDif)) {- Relative change in X    -}
  1131.                                       {-                         -}
  1132.                                       {-                         -}
  1133. (*         or                      *) {-                         -}
  1134. (*                                 *) {-                         -}
  1135. (*  (Modulus(Dif) < Tol)           *) {- Absolute change in X    -}
  1136. (*                                 *) {-                         -}
  1137. (*         or                      *) {-                         -}
  1138. (*                                 *) {-                         -}
  1139. (*  (Modulus(Y) <= Tol)            *) {- Absolute change in y   -}
  1140.                                       {---------------------------}
  1141.  
  1142. {-----------------------------------------------------------------------}
  1143. {- The first criteria simply checks to see if the value of the         -}
  1144. {- function is zero.  You should probably always keep this criteria    -}
  1145. {- active.                                                             -}
  1146. {-                                                                     -}
  1147. {- The second criteria checks the relative error in x. This criteria   -}
  1148. {- evaluates the fractional change in x between interations. Note      -}
  1149. {- that x has been multiplied through the inequality to avoid Divide   -}
  1150. {- by zero errors.                                                     -}
  1151. {-                                                                     -}
  1152. {- The third criteria checks the absolute difference in x between      -}
  1153. {- iterations.                                                         -}
  1154. {-                                                                     -}
  1155. {- The fourth criteria checks the absolute difference between          -}
  1156. {- the value of the function and zero.                                 -}
  1157. {-----------------------------------------------------------------------}
  1158.  
  1159. end; { function TestForRoot }
  1160.  
  1161. begin { procedure  Muller }
  1162.   Initial(Guess, Tol, MaxIter, Error, Iter, Found, NewApprox, yNewApprox,
  1163.           X0, X1, OldApprox, Factor);
  1164.  
  1165.   while not Found and (Error = 0) and (Iter < MaxIter) do
  1166.   begin
  1167.     Iter := Succ(Iter);
  1168.     QuadraticFormula(Factor, OldApprox, NewApprox);
  1169.     UserProcedure(NewApprox, yNewApprox, FuncPtr);    { Calculate a new yNewApprox }
  1170.     Found := TestForRoot(NewApprox, OldApprox, yNewApprox, Tol);
  1171.     X0 := X1;
  1172.     X1 := OldApprox;
  1173.     OldApprox := NewApprox;
  1174.     Factor.C := yNewApprox;
  1175.     MakeParabola(X0, X1, OldApprox, Factor, Error);
  1176.   end;
  1177.   Answer := NewApprox;
  1178.   yAnswer := yNewApprox;
  1179.   if Found then
  1180.     Error := 0
  1181.   else
  1182.     if (Error = 0) and (Iter >= MaxIter) then
  1183.       Error := 1;
  1184. end; { procedure Muller }
  1185.  
  1186. procedure Laguerre{(var Degree   : integer;
  1187.                    var Poly      : TNCompVector;
  1188.                        InitGuess : TNcomplex;
  1189.                        Tol       : Float;
  1190.                        MaxIter   : integer;
  1191.                    var NumRoots  : integer;
  1192.                    var Roots     : TNCompVector;
  1193.                    var yRoots    : TNCompVector;
  1194.                    var Iter      : TNIntVector;
  1195.                    var Error     : byte)};
  1196.  
  1197. type
  1198.   TNquadratic = record
  1199.                   A, B, C : Float;
  1200.                 end;
  1201.  
  1202. var
  1203.   AddIter    : integer;
  1204.   InitDegree : integer;
  1205.   InitPoly   : TNCompVector;
  1206.   GuessRoot  : TNcomplex;
  1207.  
  1208. {----------- Here are a few complex operations ------------}
  1209.  
  1210. procedure Conjugate(var C1, C2 : TNcomplex);
  1211. begin
  1212.   C2.Re := C1.Re;
  1213.   C2.Im := -C1.Im;
  1214. end; { procedure Conjugate }
  1215.  
  1216. function Modulus(var C1 : TNcomplex) : Float;
  1217. begin
  1218.   Modulus := Sqrt(Sqr(C1.Re) + Sqr(C1.Im));
  1219. end; { function Modulus }
  1220.  
  1221. procedure Add(var C1, C2, C3 : TNcomplex);
  1222. begin
  1223.   C3.Re := C1.Re + C2.Re;
  1224.   C3.Im := C1.Im + C2.Im;
  1225. end; { procedure Add }
  1226.  
  1227. procedure Sub(var C1, C2, C3 : TNcomplex);
  1228. begin
  1229.   C3.Re := C1.Re - C2.Re;
  1230.   C3.Im := C1.Im - C2.Im;
  1231. end; { procedure Sub }
  1232.  
  1233. procedure Mult(var C1, C2, C3 : TNcomplex);
  1234. begin
  1235.   C3.Re := C1.Re * C2.Re - C1.Im * C2.Im;
  1236.   C3.Im := C1.Im * C2.Re + C1.Re * C2.Im;
  1237. end; { procedure Mult }
  1238.  
  1239. procedure Divide(var C1, C2, C3 : TNcomplex);
  1240. var
  1241.   Dum1, Dum2 : TNcomplex;
  1242.   E : Float;
  1243. begin
  1244.   Conjugate(C2, Dum1);
  1245.   Mult(C1, Dum1, Dum2);
  1246.   E := Sqr(Modulus(C2));
  1247.   C3.Re := Dum2.Re / E;
  1248.   C3.Im := Dum2.Im / E;
  1249. end;  { procedure Divide }
  1250.  
  1251. procedure SquareRoot(var C1, C2 : TNcomplex);
  1252. const
  1253.   NearlyZero = 1E-015;
  1254. var
  1255.   R, Theta : Float;
  1256. begin
  1257.   R := Sqrt(Sqr(C1.Re) + Sqr(C1.Im));
  1258.   if ABS(C1.Re) < NearlyZero then
  1259.     begin
  1260.       if C1.Im < 0 then
  1261.         Theta := Pi / 2
  1262.       else
  1263.         Theta := -Pi / 2;
  1264.     end
  1265.   else
  1266.     if C1.Re < 0 then
  1267.       Theta := ArcTan(C1.Im / C1.Re) + Pi
  1268.     else
  1269.       Theta := ArcTan(C1.Im / C1.Re);
  1270.   C2.Re := Sqrt(R) * Cos(Theta / 2);
  1271.   C2.Im := Sqrt(R) * Sin(Theta / 2);
  1272. end; { procedure SquareRoot }
  1273.  
  1274. procedure InitAndTest(var Degree     : integer;
  1275.                       var Poly       : TNCompVector;
  1276.                           Tol        : Float;
  1277.                           MaxIter    : integer;
  1278.                           InitGuess  : TNcomplex;
  1279.                       var NumRoots   : integer;
  1280.                       var Roots      : TNCompVector;
  1281.                       var yRoots     : TNCompVector;
  1282.                       var Iter       : TNIntVector;
  1283.                       var GuessRoot  : TNcomplex;
  1284.                       var InitDegree : integer;
  1285.                       var InitPoly   : TNCompVector;
  1286.                       var Error      : byte);
  1287.  
  1288. {----------------------------------------------------------}
  1289. {- Input:  Degree, Poly, Tol, MaxIter, InitGuess          -}
  1290. {- Output: InitDegree, InitPoly, Degree, Poly, NumRoots,  -}
  1291. {-         Roots, yRoots, Iter, GuessRoot, Error          -}
  1292. {-                                                        -}
  1293. {- This procedure sets the initial value of the above     -}
  1294. {- variables.  This procedure also tests the tolerance    -}
  1295. {- (Tol), maximum number of iterations (MaxIter), and     -}
  1296. {- code.  Finally, it examines the coefficients of Poly.  -}
  1297. {- If the constant term is zero, then zero is one of the  -}
  1298. {- roots and the polynomial is deflated accordingly. Also -}
  1299. {- if the leading coefficient is zero, then Degree is     -}
  1300. {- reduced until the leading coefficient is non-zero.     -}
  1301. {----------------------------------------------------------}
  1302.  
  1303. var
  1304.   Term : integer;
  1305.  
  1306. begin
  1307.   Error := 0;
  1308.   if Degree <= 0 then
  1309.     Error := 2;      { degree is less than 2 }
  1310.   if Tol <= 0 then
  1311.     Error := 3;
  1312.   if MaxIter < 0 then
  1313.     Error := 4;
  1314.  
  1315.   if Error = 0 then
  1316.   begin
  1317.     NumRoots := 0;
  1318.     GuessRoot := InitGuess;
  1319.     InitDegree := Degree;
  1320.     InitPoly := Poly;
  1321.     { Reduce degree until leading coefficient <> zero }
  1322.     while (Degree > 0) and (Modulus(Poly[Degree]) < TNNearlyZero) do
  1323.       Degree := Pred(Degree);
  1324.     { Deflate polynomial until the constant term <> zero }
  1325.     while (Modulus(Poly[0]) = 0) and (Degree > 0) do
  1326.     begin
  1327.       { Zero is a root }
  1328.       NumRoots := Succ(NumRoots);
  1329.       Roots[NumRoots].Re := 0;
  1330.       Roots[NumRoots].Im := 0;
  1331.       yRoots[NumRoots].Re := 0;
  1332.       yRoots[NumRoots].Im := 0;
  1333.       Iter[NumRoots] := 0;
  1334.       Degree := Pred(Degree);
  1335.       for Term := 0 to Degree do
  1336.         Poly[Term] := Poly[Term + 1];
  1337.     end;
  1338.   end;
  1339. end; { procedure InitAndTest }
  1340.  
  1341. procedure FindOneRoot(Degree    : integer;
  1342.                       Poly      : TNCompVector;
  1343.                       GuessRoot : TNcomplex;
  1344.                       Tol       : Float;
  1345.                       MaxIter   : integer;
  1346.                   var Root      : TNcomplex;
  1347.                   var yValue    : TNcomplex;
  1348.                   var Iter      : integer;
  1349.                   var Error     : byte);
  1350.  
  1351. {-------------------------------------------------------------------}
  1352. {- Input:  Degree, Poly, GuessRoot, Tol, MaxIter                   -}
  1353. {- Output: Root, yValue, Iter, Error                               -}
  1354. {-                                                                 -}
  1355. {- This procedure approximates a single root of the polynomial     -}
  1356. {- Poly.  The root must be approximated within MaxIter             -}
  1357. {- iterations to a tolerance of Tol.  The root, value of the       -}
  1358. {- polynomial at the root (yValue), and the number of iterations   -}
  1359. {- (Iter) are returned. If no root is found, the appropriate error -}
  1360. {- code (Error) is returned.                                       -}
  1361. {-------------------------------------------------------------------}
  1362.  
  1363. var
  1364.   Found : boolean;
  1365.   Dif : TNcomplex;
  1366.   yPrime, yDoublePrime : TNcomplex;
  1367.  
  1368. procedure EvaluatePoly(Degree       : integer;
  1369.                        Poly         : TNCompVector;
  1370.                        X            : TNcomplex;
  1371.                    var yValue       : TNcomplex;
  1372.                    var yPrime       : TNcomplex;
  1373.                    var yDoublePrime : TNcomplex);
  1374.  
  1375. {--------------------------------------------------------------------}
  1376. {- Input:  Degree, Poly, X                                          -}
  1377. {- Output: yValue, yPrime, yDoublePrime                             -}
  1378. {-                                                                  -}
  1379. {- This procedure applies the technique of synthetic division to    -}
  1380. {- determine value (yValue), first derivative (yPrime) and second   -}
  1381. {- derivative (yDoublePrime) of the  polynomial, Poly, at X.        -}
  1382. {- The 0th element of the first synthetic division is the           -}
  1383. {- value of Poly at X, the 1st element of the second synthetic      -}
  1384. {- division is the first derivative of Poly at X, and twice the     -}
  1385. {- 2nd element of the third synthetic division is the second        -}
  1386. {- derivative of Poly at X.                                         -}
  1387. {--------------------------------------------------------------------}
  1388.  
  1389. var
  1390.   Loop : integer;
  1391.   Dummy, yDPdummy : TNcomplex;
  1392.   Deriv, Deriv2 : TNCompVector;
  1393.  
  1394. begin
  1395.   Deriv[Degree] := Poly[Degree];
  1396.   for Loop := Degree - 1 downto 0 do
  1397.   begin
  1398.     Mult(Deriv[Loop + 1], X, Dummy);
  1399.     Add(Dummy, Poly[Loop], Deriv[Loop]);
  1400.   end;
  1401.   yValue := Deriv[0];    { Value of Poly at X }
  1402.  
  1403.   Deriv2[Degree] := Deriv[Degree];
  1404.   for Loop := Degree - 1 downto 1 do
  1405.   begin
  1406.     Mult(Deriv2[Loop + 1], X, Dummy);
  1407.     Add(Dummy, Deriv[Loop], Deriv2[Loop]);
  1408.   end;
  1409.   yPrime := Deriv2[1];   { 1st deriv. of Poly at X }
  1410.  
  1411.   yDPdummy := Deriv2[Degree];
  1412.   for Loop := Degree - 1 downto 2 do
  1413.   begin
  1414.     Mult(yDPdummy, X, Dummy);
  1415.     Add(Dummy, Deriv2[Loop], yDPdummy);
  1416.   end;
  1417.   yDoublePrime.Re := 2 * yDPdummy.Re;    { 2nd derivative of Poly at X }
  1418.   yDoublePrime.Im := 2 * yDPdummy.Im;
  1419. end; { procedure EvaluatePoly }
  1420.  
  1421. procedure ConstructDifference(Degree       : integer;
  1422.                               yValue       : TNcomplex;
  1423.                               yPrime       : TNcomplex;
  1424.                               yDoublePrime : TNcomplex;
  1425.                           var Dif          : TNcomplex);
  1426.  
  1427. {------------------------------------------------------------------}
  1428. {- Input:  Degree, yValue, yPrime, yDoublePrime                   -}
  1429. {- Output: Dif                                                    -}
  1430. {-                                                                -}
  1431. {- This procedure computes the difference between approximations; -}
  1432. {- given information about the function and its first two         -}
  1433. {- derivatives.                                                   -}
  1434. {-----------------------------------------------------------------}
  1435.  
  1436. var
  1437.   yPrimeSQR, yTimesyDPrime, Sum, SRoot,
  1438.   Numer1, Numer2, Numer, Denom : TNcomplex;
  1439.  
  1440. begin
  1441.   Mult(yPrime, yPrime, yPrimeSQR);
  1442.   yPrimeSQR.Re := Sqr(Degree - 1) * yPrimeSQR.Re;
  1443.   yPrimeSQR.Im := Sqr(Degree - 1) * yPrimeSQR.Im;
  1444.   Mult(yValue, yDoublePrime, yTimesyDPrime);
  1445.   yTimesyDPrime.Re := (Degree - 1) * Degree * yTimesyDPrime.Re;
  1446.   yTimesyDPrime.Im := (Degree - 1) * Degree * yTimesyDPrime.Im;
  1447.   Sub(yPrimeSQR, yTimesyDPrime, Sum);
  1448.   SquareRoot(Sum, SRoot);
  1449.   Add(yPrime, SRoot, Numer1);
  1450.   Sub(yPrime, SRoot, Numer2);
  1451.   if Modulus(Numer1) > Modulus(Numer2) then
  1452.     Numer := Numer1
  1453.   else
  1454.     Numer := Numer2;
  1455.   Denom.Re := Degree * yValue.Re;
  1456.   Denom.Im := Degree * yValue.Im;
  1457.   if Modulus(Numer) < TNNearlyZero then
  1458.     begin
  1459.       Dif.Re := 0;
  1460.       Dif.Im := 0;
  1461.     end
  1462.   else
  1463.     Divide(Denom, Numer, Dif);  { The difference is the   }
  1464.                                 { inverse of the fraction }
  1465. end; { procedure ConstructDifference }
  1466.  
  1467. function TestForRoot(X, Dif, Y, Tol : Float) : boolean;
  1468.  
  1469. {--------------------------------------------------------------------}
  1470. {-  These are the stopping criteria.  Four different ones are       -}
  1471. {-  provided.  If you wish to change the active criteria, simply    -}
  1472. {-  comment off the current criteria (including the appropriate OR) -}
  1473. {-  and remove the comment brackets from the criteria (including    -}
  1474. {-  the appropriate OR) you wish to be active.                      -}
  1475. {--------------------------------------------------------------------}
  1476.  
  1477. begin
  1478.   TestForRoot :=                      {---------------------------}
  1479.     (ABS(Y) <= TNNearlyZero)          {- Y=0                     -}
  1480.                                       {-                         -}
  1481.            or                         {-                         -}
  1482.                                       {-                         -}
  1483.     (ABS(Dif) < ABS(X * Tol))         {- Relative change in X    -}
  1484.                                       {-                         -}
  1485.                                       {-                         -}
  1486.  (*       or                      *)  {-                         -}
  1487.  (*                               *)  {-                         -}
  1488.  (* (ABS(Dif) < Tol)              *)  {- Absolute change in X    -}
  1489.  (*                               *)  {-                         -}
  1490.  (*       or                      *)  {-                         -}
  1491.  (*                               *)  {-                         -}
  1492.  (* (ABS(Y) <= Tol)               *)  {- Absolute change in Y    -}
  1493.                                       {---------------------------}
  1494.  
  1495. {-----------------------------------------------------------------------}
  1496. {- The first criteria simply checks to see if the value of the         -}
  1497. {- function is zero.  You should probably always keep this criteria    -}
  1498. {- active.                                                             -}
  1499. {-                                                                     -}
  1500. {- The second criteria checks the relative error in X. This criteria   -}
  1501. {- evaluates the fractional change in X between interations. Note      -}
  1502. {- that X has been multiplied throught the inequality to avoid divide  -}
  1503. {- by zero errors.                                                     -}
  1504. {-                                                                     -}
  1505. {- The third criteria checks the absolute difference in X between      -}
  1506. {- iterations.                                                         -}
  1507. {-                                                                     -}
  1508. {- The fourth criteria checks the absolute difference between          -}
  1509. {- the value of the function and zero.                                 -}
  1510. {-----------------------------------------------------------------------}
  1511.  
  1512. end; { procedure TestForRoot }
  1513.  
  1514. begin { procedure FindOneRoot }
  1515.   Root := GuessRoot;
  1516.   Found := false;
  1517.   Iter := 0;
  1518.   EvaluatePoly(Degree, Poly, Root, yValue, yPrime, yDoublePrime);
  1519.   while (Iter < MaxIter) and not(Found) do
  1520.   begin
  1521.     Iter := Succ(Iter);
  1522.     ConstructDifference(Degree, yValue, yPrime, yDoublePrime, Dif);
  1523.     Sub(Root, Dif, Root);
  1524.     EvaluatePoly(Degree, Poly, Root, yValue, yPrime, yDoublePrime);
  1525.     Found := TestForRoot(Modulus(Root), Modulus(Dif), Modulus(yValue), Tol);
  1526.   end;
  1527.   if not(Found) then Error := 1;   { Iterations execeeded MaxIter }
  1528. end; { procedure FindOneRoot }
  1529.  
  1530. procedure ReducePoly(var Degree : integer;
  1531.                      var Poly   : TNCompVector;
  1532.                      Root       : TNcomplex);
  1533.  
  1534. {------------------------------------------------------}
  1535. {- Input: Degree, Poly, Root                          -}
  1536. {- Output: Degree, Poly                               -}
  1537. {-                                                    -}
  1538. {- This procedure deflates the polynomial Poly by     -}
  1539. {- factoring out the Root.  Degree is reduced by one. -}
  1540. {------------------------------------------------------}
  1541.  
  1542. var
  1543.   Term : integer;
  1544.   NewPoly : TNCompVector;
  1545.   Dummy : TNcomplex;
  1546.  
  1547. begin
  1548.   NewPoly[Degree - 1] := Poly[Degree];
  1549.   for Term := Degree - 1 downto 1 do
  1550.   begin
  1551.     Mult(NewPoly[Term], Root, Dummy);
  1552.     Add(Dummy, Poly[Term], NewPoly[Term - 1]);
  1553.   end;
  1554.   Degree := Pred(Degree);
  1555.   Poly := NewPoly;
  1556. end; { procedure ReducePoly }
  1557.  
  1558. begin  { procedure Laguerre }
  1559.   InitAndTest(Degree, Poly, Tol, MaxIter, InitGuess, NumRoots, Roots,
  1560.               yRoots, Iter, GuessRoot, InitDegree, InitPoly, Error);
  1561.   while (Degree > 0) and (Error = 0) do
  1562.   begin
  1563.     FindOneRoot(Degree, Poly, GuessRoot, Tol, MaxIter,
  1564.                 Roots[NumRoots + 1], yRoots[NumRoots + 1],
  1565.                 Iter[NumRoots + 1], Error);
  1566.     if Error = 0 then
  1567.     begin
  1568.       {------------------------------------------------------}
  1569.       {- The next statement refines the approximate root by -}
  1570.       {- plugging it into the original polynomial.  This    -}
  1571.       {- eliminates a lot of the round-off error            -}
  1572.       {- accumulated through many iterations                -}
  1573.       {------------------------------------------------------}
  1574.       FindOneRoot(InitDegree, InitPoly, Roots[NumRoots + 1],
  1575.                   Tol, MaxIter, Roots[NumRoots + 1],
  1576.                   yRoots[NumRoots + 1], AddIter, Error);
  1577.       Iter[NumRoots + 1] := Iter[NumRoots + 1] + AddIter;
  1578.       NumRoots := Succ(NumRoots);
  1579.       ReducePoly(Degree, Poly, Roots[NumRoots]); { Reduce polynomial }
  1580.     end;
  1581.     GuessRoot := Roots[NumRoots];
  1582.   end;
  1583. end; { procedure Laguerre }
  1584.