Télécharger laguer.eso

Retour à la liste

Numérotation des lignes :

  1. C LAGUER SOURCE PV 07/11/23 21:17:36 5978
  2. subroutine laguer(a,m,x,its,kerre)
  3. IMPLICIT INTEGER(I-N)
  4. IMPLICIT REAL*8 (A-H,O-Z)
  5. dimension frac(8)
  6. real*8 err,abp,abm,abx,frac
  7. complex dx,x1,b,d,f,g,h,sq,gp,gm,g2
  8. complex a(m+1),x
  9. data frac/.5,.25,.75,.13,.38,.62,.88,1./
  10. * save frac
  11. * write(6,*) ' entree dans laguer m ',m
  12. epss = 2.e-7
  13. mr = 8
  14. mt=10
  15. maxit = mt*mr*3
  16. do 1 iter = 1,maxit
  17. its=iter
  18. b=a(m+1)
  19. err = abs(b)
  20. d=cmplx(0.,0.)
  21. f=cmplx(0.,0.)
  22. abx = abs(x)
  23. do 2 j=m,1,-1
  24. f = x*f+d
  25. d=x*d+b
  26. b=x*b+a(j)
  27. err=abs(b)+abx*err
  28. 2 continue
  29. err=epss*err
  30. if(abs(b).le.err) then
  31. return
  32. else
  33. g=d/b
  34. g2=g*g
  35. h=g2-2.*f/b
  36. sq=sqrt((m-1)*(m*h-g2))
  37. gp=g+sq
  38. gm=g-sq
  39. abp=abs(gp)
  40. abm=abs(gm)
  41. if(abp.lt.abm) gp=gm
  42. if(max(abp,abm).gt.0.) then
  43. dx=m/gp
  44. else
  45. dx=exp(cmplx(log(1.+abx),float(iter)))
  46. endif
  47. endif
  48. x1=x-dx
  49. if(x.eq.x1) return
  50. if(mod(iter,MT).ne.0) then
  51. x=x1
  52. else
  53. x=x-dx*frac(iter/mt)
  54. endif
  55. 1 continue
  56. kerre =944
  57. return
  58. end
  59.  
  60.  
  61.  
  62.  
  63.  
  64.  
  65.  

© Cast3M 2003 - Tous droits réservés.
Mentions légales