home *** CD-ROM | disk | FTP | other *** search
/ Otherware / Otherware_1_SB_Development.iso / amiga / misc / icalc.lha / ICalc / Scripts / zroots.ic < prev   
Encoding:
Text File  |  1992-08-03  |  2.9 KB  |  135 lines
  1. # This stuff is based on algorithms presented in
  2. # "Numerical Recipes - The Art of Scientific Computing", by Press et al.
  3. #
  4. # Finds roots (real and complex) of a polynomial.
  5. # For usage see note preceeding definition of zroots().
  6. #
  7. # MWS, July 1992.
  8.  
  9. silent
  10.  
  11. #
  12. # Given coeffs of poly of degree (m+1) in a, and estimate x of root,
  13. # this function improves x to desired accuracy; if polish != 0, this is
  14. # to machine's roundoff-limit, otherwise it is to eps.
  15. # NB: assumes arraybase is zero, as set in zroots() which calls laguer()
  16. #
  17. LEPSS = 6.e-8
  18. LMAXIT = 100    # limit on iterations for convergence
  19. #
  20. func laguer(@a,m,*x,eps,polish) = {
  21.     local j,iter,err,dxold,cdx,abx,sq,h,gp,gm,g2,g,b,d,dx,f,x1
  22.  
  23.     dxold=abs(x)
  24.     for (iter=1;iter<=LMAXIT;iter+=1) {
  25.         b=a[m]
  26.         err=abs(b)
  27.         d=f=0
  28.         abx=abs(x)
  29.         for (j=m-1;j>=0;j-=1) {
  30.             # compute value of poly and 1st & 2nd derivatives at x
  31.             f=x*f+d
  32.             d=x*d+b
  33.             b=x*b+a[j]
  34.             err=abs(b)+abx*err;
  35.         }
  36.         err *= LEPSS            # estimate roundoff error
  37.         if (abs(b) <= err) return 1    # we are at root
  38.         g=d/b                # use Laguerre's formula
  39.         g2=g*g
  40.         h=g2-2*f/b
  41.         sq=sqrt((m-1)*(m*h-g2))
  42.         gp=g+sq
  43.         gm=g-sq
  44.         if (abs(gp) < abs(gm)) gp=gm
  45.         dx=m/gp
  46.         x1=x-dx
  47.         if (x == x1) return 1        # convergence achieved
  48.         x=x1
  49.         cdx=abs(dx);
  50.         dxold=cdx;
  51.         if (!polish)
  52.             if (cdx <= eps*abs(x)) return 1
  53.     }
  54.     error(0)    # too many iterations
  55. }
  56.  
  57.  
  58. #
  59. # This is the main interface. Call zroots with two arrays:
  60. #     a    contains the coefficients of poly a[0] + a[1]*x +...+ a[m]*x^m
  61. #    roots    storage area for roots of poly
  62. # The polish parameter determines whether highest accuracy is required. Set
  63. # to 1 if this is desired, otherwise set it to zero.
  64. #
  65. ZEPS=2.0e-6        # desired accuracy
  66. array zrwkspc[1]    # resize as required
  67. #
  68. func zroots(@a,@roots,polish) = {
  69.     local m,jj,j,ii,x,b,c,finis,oldbase
  70.  
  71.     oldbase = arraybase(0)    # we want a[0] as first element
  72.     m = sizeof(a)-1        # m is degree of poly
  73.     resize(zrwkspc,m+1)    # make workspace big enough
  74.  
  75.     for (j=0;j<=m;j+=1) zrwkspc[j]=a[j]    # copy coefficients
  76.     for (j=m;j>=1;j-=1) {            # loop over each root to be found
  77.         x=0
  78.         laguer(zrwkspc,j,x,ZEPS,0)    # find root
  79.         if (abs(Im(x)) <= (2.0*ZEPS*abs(Re(x)))) x=Re(x)
  80.         roots[j-1]=x
  81.         b=zrwkspc[j]            # forward deflation
  82.  
  83.         for (jj=j-1;jj>=0;jj-=1) {
  84.             c=zrwkspc[jj]
  85.             zrwkspc[jj]=b
  86.             b=x*b+c
  87.         }
  88.     }
  89.     if (polish)        # polish roots
  90.         for (j=0;j<m;j+=1)
  91.             laguer(a,m,roots[j],ZEPS,1);
  92.  
  93.     arraybase(oldbase)
  94.     1            # standard default return-value
  95. }
  96.  
  97.  
  98. #
  99. # from poly.ic - evaluates poly in c at x
  100. #
  101. func epoly(@c, x) = {
  102.     local n,j,ss,oldbase
  103.  
  104.     oldbase = arraybase(1)
  105.     n = sizeof(c)
  106.  
  107.     ss = c[n]
  108.     for (j = n-1; j >= 1; j -= 1)
  109.         ss = ss*x+c[j]
  110.  
  111.     arraybase(oldbase)
  112.     ss
  113. }
  114.  
  115.  
  116. #
  117. #an example of usage
  118. #
  119.  
  120. N=7    #degree
  121. array a[N+1] = {1,2,3,4,5,6,7,8}    # polynomial
  122. array r[N]                # where to store roots
  123.  
  124. echo "Routines defined. Performing example:"
  125. echo "\tPoly is 1+2x+3x^2+...+8x^7 - be patient..."
  126. zroots(a,r,1)
  127. echo "done. Roots are:"
  128.  
  129. display(r)
  130. echo "\nPoly evaluated at each root:"
  131. for (j=1; j <= N; j+=1)
  132.     print(epoly(a,r[j]))
  133.  
  134. verbose
  135.