home *** CD-ROM | disk | FTP | other *** search
/ The C Users' Group Library 1994 August / wc-cdrom-cusersgrouplibrary-1994-08.iso / vol_200 / 247_01 / bnflsh2.c < prev    next >
Text File  |  1989-04-19  |  4KB  |  188 lines
  1. /*
  2.  *  MIRACL flash exponential and logs
  3.  *  bnflsh2.c
  4.  */
  5.  
  6. #include <stdio.h>
  7. #include <math.h>
  8. #include "miracl.h"
  9. #define abs(x)  ((x)<0? (-(x)) : (x))
  10.  
  11. /* Access Global variables */
  12.  
  13. extern int  depth;    /* error tracing... */
  14. extern int  trace[];  /* .. mechanism     */
  15. extern int  nib;      /* length of bigs   */
  16. extern int  workprec;
  17. extern int  stprec;   /* start precision  */
  18. extern int  lg2b;     /* no. bits in base */
  19.  
  20. extern big  w8,w9,w10,w11; /* workspace variables */
  21. extern big  w12,w13;
  22. extern big  w5,w6;
  23.  
  24. int expon(w,n)
  25. big w;
  26. int n;
  27. {  /* generator for C.F. of e */ 
  28.     if (n==0) return 2;
  29.     if (n%3==2) return 2*(n/3)+2;
  30.     else return 1;
  31. }
  32.  
  33. void fexp(x,y)
  34. flash x,y;
  35. { /* calculates y=exp(x) */
  36.     int i,n,nsq,m,sqrn,op[5];
  37.     bool minus;
  38.     if (ERNUM) return;
  39.     if (size(x)==0)
  40.     {
  41.         convert(1,y);
  42.         return;
  43.     }
  44.     copy(x,y);
  45.     depth++;
  46.     trace[depth]=54;
  47.     if (TRACER) track();
  48.     minus=FALSE;
  49.     if (size(y)<0)
  50.     {
  51.         minus=TRUE;
  52.         negate(y,y);
  53.     }
  54.     ftrunc(y,y,w9);
  55.     n=size(y);
  56.     if (n==TOOBIG)
  57.     {
  58.         berror(13);
  59.         depth--;
  60.         return;
  61.     }
  62.     if (n==0) convert(1,y);
  63.     else
  64.     {
  65.         build(y,expon);
  66.         fpower(y,n,y);
  67.     }
  68.     if (size(w9)==0)
  69.     {
  70.         if (minus) frecip(y,y);
  71.         depth--;
  72.         return;
  73.     }
  74.     sqrn=isqrt(lg2b*workprec,lg2b);
  75.     nsq=0;
  76.     copy(w9,w8);
  77.     frecip(w9,w9);
  78.     ftrunc(w9,w9,w9);
  79.     m=logb2(w9)+1;
  80.     if (m<sqrn)
  81.     { /* scale fraction down */
  82.         nsq=sqrn-m;
  83.         expb2(w9,nsq);
  84.         fdiv(w8,w9,w8);
  85.     }
  86.     zero(w10);
  87.     op[0]=0x4B;  /* set up for x/(C+y) */
  88.     op[1]=1;
  89.     op[2]=0;
  90.     for (m=sqrn;m>0;m--)
  91.     { /* Unwind C.F. expansion for exp(x)-1 */
  92.         if (m%2==0) op[4]=2,op[3]=1;
  93.         else        op[4]=m,op[3]=(-1);
  94.         flop(w8,w10,op,w10);
  95.     }
  96.     op[0]=0x2C;  /* set up for (x+2).y */
  97.     op[1]=op[3]=1;
  98.     op[2]=2;
  99.     op[4]=0;
  100.     for (i=0;i<nsq;i++)
  101.     { /* now square it back up again */
  102.         flop(w10,w10,op,w10);
  103.     }
  104.     op[2]=1;
  105.     flop(w10,y,op,y);
  106.     if (minus) frecip(y,y);
  107.     depth--;
  108. }
  109.  
  110. void flog(x,y)
  111. flash x,y;
  112. { /* calculate y=log(x) to base e */
  113.     bool hack;
  114.     int op[5];
  115.     copy(x,y);
  116.     if (ERNUM) return;
  117.     if (size(y)==1)
  118.     {
  119.         zero(y);
  120.         return;
  121.     }
  122.     depth++;
  123.     trace[depth]=55;
  124.     if (TRACER) track();
  125.     if (size(y)<=0)
  126.     {
  127.         berror(15);
  128.         depth--;
  129.         return;
  130.     }
  131.     hack=FALSE;
  132.     if (lent(y)<=2)
  133.     { /* for 'simple' values of y */
  134.         hack=TRUE;
  135.         build(w11,expon);
  136.         fdiv(y,w11,y);
  137.     }
  138.     op[0]=0x68;
  139.     op[1]=op[3]=1;
  140.     op[2]=(-1);
  141.     op[4]=0;
  142.     workprec=stprec;
  143.     dconv(log(fdsize(y)),w11);
  144.     while (workprec!=nib)
  145.     { /* Newtons iteration w11=w11+(y-exp(w11))/exp(w11) */
  146.         if (workprec<nib) workprec*=2;
  147.         if (workprec>=nib) workprec=nib;
  148.         else if (workprec*2>nib) workprec=(nib+1)/2;
  149.         fexp(w11,w12);
  150.         flop(y,w12,op,w12);
  151.         fadd(w12,w11,w11);
  152.     }
  153.     copy(w11,y);
  154.     if (hack) fincr(y,1,1,y);
  155.     depth--;
  156. }
  157.     
  158. void fpowf(x,y,z)
  159. flash x,y,z;
  160. { /* raise flash number to flash power *
  161.    *     z=x^y  -> z=exp(y.log(x))     */
  162.     int n;
  163.     if (ERNUM) return;
  164.     depth++;
  165.     trace[depth]=56;
  166.     if (TRACER) track();
  167.     n=size(y);
  168.     if (abs(n)<TOOBIG)
  169.     { /* y is small int */
  170.         fpower(x,n,z);
  171.         depth--;
  172.         return;
  173.     }
  174.     frecip(y,w13);
  175.     n=size(w13);
  176.     if (abs(n)<TOOBIG)
  177.     { /* 1/y is small int */
  178.         froot(x,n,z);
  179.         depth--;
  180.         return;
  181.     }
  182.     copy(x,z);
  183.     flog(z,z);
  184.     fdiv(z,w13,z);
  185.     fexp(z,z);
  186.     depth--;
  187. }
  188.