home *** CD-ROM | disk | FTP | other *** search
/ C/C++ Interactive Guide / c-cplusplus-interactive-guide.iso / c_ref / csource4 / 247_01 / bnflsh3.c < prev    next >
Encoding:
C/C++ Source or Header  |  1989-04-19  |  5.9 KB  |  289 lines
  1. /*
  2.  *  MIRACL flash trig.
  3.  *  bnflsh3.c
  4.  */
  5.  
  6. #include <stdio.h>
  7. #include <math.h>
  8. #include "miracl.h"
  9. #define TAN 1
  10. #define SIN 2
  11. #define COS 3
  12.  
  13. /* Access Global variables */
  14.  
  15. extern int  depth;    /* error tracing... */
  16. extern int  trace[];  /* .. mechanism     */
  17. extern int  nib;      /* length of bigs   */
  18. extern int  workprec;
  19. extern int  stprec;   /* start precision  */
  20. extern int  lg2b;     /* no. bits in base */
  21.  
  22. extern big  w8,w9,w10,w11; /* workspace variables */
  23. extern big  w12,w13,w14,w15;
  24. extern big  w5,w6;
  25.  
  26. int norm(type,y)
  27. flash y;
  28. int type;
  29. { /* convert y to first quadrant angle *
  30.    * and return sign of result         */
  31.     int s;
  32.     if (ERNUM) return 0;
  33.     s=PLUS;
  34.     if (size(y)<0)
  35.     {
  36.         negate(y,y);
  37.         if (type!=COS) s=(-s);
  38.     }
  39.     fpi(w15);
  40.     fpmul(w15,1,2,w8);
  41.     if (fcomp(y,w8)<=0) return s;
  42.     fpmul(w15,2,1,w8);
  43.     if (fcomp(y,w8)>0)
  44.     { /* reduce mod 2.pi */
  45.         fdiv(y,w8,w9);
  46.         ftrunc(w9,w9,w9);
  47.         fmul(w9,w8,w9);
  48.         fsub(y,w9,y);
  49.     }
  50.     if (fcomp(y,w15)>0)
  51.     { /* if greater than pi */
  52.         fsub(y,w15,y);
  53.         if (type!=TAN) s=(-s);
  54.     }
  55.     fpmul(w15,1,2,w8);
  56.     if (fcomp(y,w8)>0)
  57.     { /* if greater than pi/2 */
  58.         fsub(w15,y,y);
  59.         if (type!=SIN) s=(-s);
  60.     }
  61.     return s;
  62. }
  63.  
  64. int tan1(w,n)
  65. big w;
  66. int n;
  67. {  /* generator for C.F. of tan(1) */ 
  68.     if (n==0) return 1;
  69.     if (n%2==1) return 2*(n/2)+1;
  70.     else return 1;
  71. }
  72.  
  73. void ftan(x,y)
  74. flash x,y;
  75. { /* calculates y=tan(x) */
  76.     int i,n,nsq,m,sqrn,sgn,op[5];
  77.     copy(x,y);
  78.     if (ERNUM || size(y)==0) return;
  79.     depth++;
  80.     trace[depth]=57;
  81.     if (TRACER) track();
  82.     sgn=norm(TAN,y);
  83.     ftrunc(y,y,w10);
  84.     n=size(y);
  85.     if (n!=0) build(y,tan1);
  86.     if (size(w10)==0)
  87.     {
  88.         insign(sgn,y);
  89.         depth--;
  90.         return;
  91.     }
  92.     sqrn=isqrt(lg2b*workprec,lg2b);
  93.     nsq=0;
  94.     copy(w10,w8);
  95.     frecip(w10,w10);
  96.     ftrunc(w10,w10,w10);
  97.     m=logb2(w10);
  98.     if (m<sqrn)
  99.     { /* scale fraction down */
  100.         nsq=sqrn-m;
  101.         expb2(w10,nsq);
  102.         fdiv(w8,w10,w8);
  103.     }
  104.     zero(w10);
  105.     fmul(w8,w8,w9);
  106.     negate(w9,w9);
  107.     op[0]=0x4B;    /* set up for x/(y+C) */
  108.     op[1]=op[3]=1;
  109.     op[2]=0;
  110.     for (m=sqrn;m>1;m--)
  111.     { /* Unwind C.F for tan(x)=z/(1-z^2/(3-z^2/(5-...)))))) */
  112.         op[4]=2*m-1;
  113.         flop(w9,w10,op,w10);
  114.     }
  115.     op[4]=1;
  116.     flop(w8,w10,op,w10);
  117.     op[0]=0x6C;     /* set up for tan(x+y)=tan(x)+tan(y)/(1-tan(x).tan(y)) */
  118.     op[1]=op[2]=op[3]=1;
  119.     op[4]=(-1);
  120.     for (i=0;i<nsq;i++)
  121.     { /* now square it back up again */
  122.         flop(w10,w10,op,w10);
  123.     }
  124.     flop(y,w10,op,y);
  125.     insign(sgn,y);
  126.     depth--;
  127. }
  128.  
  129. void fatan(x,y)
  130. flash x,y;
  131. { /* calculate y=atan(x) */
  132.     int s,c,op[5];
  133.     bool inv,hack;
  134.     copy(x,y);
  135.     if (ERNUM || size(y)==0) return;
  136.     depth++;
  137.     trace[depth]=58;
  138.     if (TRACER) track();
  139.     s=exsign(y);
  140.     insign(PLUS,y);
  141.     inv=FALSE;
  142.     fpi(w15);
  143.     fconv(1,1,w11);
  144.     c=fcomp(y,w11);
  145.     if (c==0)
  146.     { /* atan(1)=pi/4 */
  147.         fpmul(w15,1,4,y);
  148.         insign(s,y);
  149.         depth--;
  150.         return;
  151.     }
  152.     if (c>0)
  153.     { /* atan(x)=pi/2 - atan(1/x) for x>1 */
  154.         inv=TRUE;
  155.         frecip(y,y);
  156.     }
  157.     hack=FALSE;
  158.     if (lent(y)<=2)
  159.     { /* for 'simple' values of y */
  160.         hack=TRUE;
  161.         fconv(3,1,w11);
  162.         froot(w11,2,w11);
  163.         op[0]=0xC6;
  164.         op[2]=op[3]=op[4]=1;
  165.         op[1]=(-1);
  166.         flop(y,w11,op,y);
  167.     }
  168.     op[0]=0x4B;
  169.     op[1]=op[3]=op[4]=1;
  170.     op[2]=0;
  171.     workprec=stprec;
  172.     dconv(atan(fdsize(y)),w11);
  173.     while (workprec!=nib)
  174.     { /* Newtons iteration w11=w11+(y-tan(w11))/(tan(w11)^2+1) */
  175.         if (workprec<nib) workprec*=2;
  176.         if (workprec>=nib) workprec=nib;
  177.         else if (workprec*2>nib) workprec=(nib+1)/2;
  178.         ftan(w11,w12);
  179.         fsub(y,w12,w8);
  180.         fmul(w12,w12,w12);
  181.         flop(w8,w12,op,w12);  /* w12=w8/(w12+1) */
  182.         fadd(w12,w11,w11);
  183.     }
  184.     copy(w11,y);
  185.     op[0]=0x6C;
  186.     op[1]=op[3]=6;
  187.     op[2]=1;
  188.     op[4]=0;
  189.     if (hack) flop(y,w15,op,y);
  190.     op[1]=(-2);
  191.     op[3]=2;
  192.     if (inv) flop(y,w15,op,y);
  193.     insign(s,y);
  194.     depth--;
  195. }
  196.  
  197. void fsin(x,y)
  198. flash x,y;
  199. { /*  calculate y=sin(x) */
  200.     int sgn,op[5];
  201.     copy(x,y);
  202.     if (ERNUM || size(y)==0) return;
  203.     depth++;
  204.     trace[depth]=59;
  205.     if (TRACER) track();
  206.     sgn=norm(SIN,y);
  207.     fpmul(y,1,2,y);
  208.     ftan(y,y);
  209.     op[0]=0x6C;
  210.     op[1]=op[2]=op[3]=op[4]=1;
  211.     flop(y,y,op,y);
  212.     insign(sgn,y);
  213.     depth--;
  214. }
  215.  
  216. void fasin(x,y)
  217. flash x,y;
  218. { /* calculate y=asin(x) */
  219.     int s,op[5];
  220.     copy(x,y);
  221.     if (ERNUM || size(y)==0) return;
  222.     depth++;
  223.     trace[depth]=60;
  224.     if (TRACER) track();
  225.     s=exsign(y);
  226.     insign(PLUS,y);
  227.     op[0]=0x3C;
  228.     op[1]=(-1);
  229.     op[2]=op[3]=1;
  230.     op[4]=0;
  231.     flop(y,y,op,w11);  /* w11 = -(y.y-1) */
  232.     froot(w11,2,w11);
  233.     if (size(w11)==0)
  234.     {
  235.         fpi(w15);
  236.         fpmul(w15,1,2,y);
  237.     }
  238.     else
  239.     {
  240.         fdiv(y,w11,y);
  241.         fatan(y,y);
  242.     }
  243.     insign(s,y);    
  244.     depth--;
  245. }
  246.  
  247. void fcos(x,y)
  248. flash x,y;
  249. { /* calculate y=cos(x) */
  250.     int sgn,op[5];
  251.     copy(x,y);
  252.     if (ERNUM || size(y)==0)
  253.     {
  254.         convert(1,y);
  255.         return;
  256.     }
  257.     depth++;
  258.     trace[depth]=61;
  259.     if (TRACER) track();
  260.     sgn=norm(COS,y);
  261.     fpmul(y,1,2,y);
  262.     ftan(y,y);
  263.     op[0]=0x33;
  264.     op[1]=op[3]=op[4]=1;
  265.     op[2]=(-1);
  266.     flop(y,y,op,y);
  267.     insign(sgn,y);
  268.     depth--;
  269. }
  270.  
  271. void facos(x,y)
  272. flash x,y;
  273. { /* calculate y=acos(x) */
  274.     int op[5];
  275.     if (ERNUM) return;
  276.     depth++;
  277.     trace[depth]=62;
  278.     if (TRACER) track();
  279.     fasin(x,y);
  280.     fpi(w15);
  281.     op[0]=0x6C;
  282.     op[1]=(-2);
  283.     op[2]=1;
  284.     op[3]=2;
  285.     op[4]=0;
  286.     flop(y,w15,op,y);    /* y = pi/2 - y */
  287.     depth--;
  288. }
  289.