share/contrib/fresnel/fresnelCStaylor.mac

   1 fresnelC(z):=block([nn,BFLOAT_PRECISION,term,az:dfloat(abs(z)),
   2    BL:ceiling(abs(log10(DFLOAT_PRECISION))),FLOAT2BF:true],
   3    BFLOAT_PRECISION:BL,
   4    for n:0 step 1 do (
   5       term:dfloat(((%pi/2)^(2*n)/(((2*n)!)*(4*n+1)))*(az^(4*n+1))),
   6       lterm:ceiling(log10(term)),
   7       if (lterm+BL) > BFLOAT_PRECISION then BFLOAT_PRECISION:lterm+BL,
   8       if term < DFLOAT_EPSILON then (nn:n,return(nn))
   9    ),
  10    BFLOAT_PRECISION:BFLOAT_PRECISION+2,
  11    z:bfloat(z),
  12    sum(bfloat((-1)^n*((%pi/2)^(2*n)/(((2*n)!)*(4*n+1)))*(z^(4*n+1))),n,0,nn)
  13 )$
  14 fresnelS(z):=block([nn,BFLOAT_PRECISION,term,az:dfloat(abs(z)),
  15    BL:ceiling(abs(log10(DFLOAT_PRECISION))),FLOAT2BF:true],
  16    BFLOAT_PRECISION:BL,
  17    for n:0 step 1 do (
  18       term:dfloat(((%pi/2)^(2*n+1)/(((2*n+1)!)*(4*n+3)))*(az^(4*n+3))),
  19       lterm:ceiling(log10(term)),
  20       if (lterm+BL) > BFLOAT_PRECISION then BFLOAT_PRECISION:lterm+BL,
  21       if term < DFLOAT_EPSILON then (nn:n,return(nn))
  22    ),
  23    BFLOAT_PRECISION:BFLOAT_PRECISION+2,
  24    z:bfloat(z),
  25    sum(bfloat((-1)^n*((%pi/2)^(2*n+1)/(((2*n+1)!)*(4*n+3)))*(z^(4*n+3))),n,0,nn)
  26 )$