home *** CD-ROM | disk | FTP | other *** search
/ C/C++ Interactive Guide / c-cplusplus-interactive-guide.iso / c_ref / csource4 / 247_01 / bnflsh1.c < prev    next >
Encoding:
C/C++ Source or Header  |  1989-04-19  |  4.1 KB  |  200 lines
  1. /*
  2.  *    MIRACL flash roots and powers
  3.  *    bnflsh1.c
  4.  */
  5.  
  6. #include <stdio.h>
  7. #include <math.h>
  8. #include "miracl.h"
  9.  
  10. /* Access global variables */
  11.  
  12. extern int depth;     /* error tracing .. */
  13. extern int trace[];   /* ... mechanism    */
  14. extern int nib;       /* length of bigs   */
  15. extern int lg2b;      /* no. bits in base */
  16. extern small base;    /* number base      */
  17. extern int workprec;
  18. extern int stprec;    /* start precision  */
  19.  
  20. extern big w7,w8,w9,w10,w15; /* workspace variables */
  21.  
  22. static int RS,RD;
  23.  
  24. int quad(w,n)
  25. big w;
  26. int n;
  27. { /* generator for C.F. of square root of small integer */
  28.     static int v,u,ud,vd,a,oldn;
  29.     int t;
  30.     if (n==0)
  31.     {
  32.         u=2*RD;
  33.         ud=u;
  34.         v=1;
  35.         vd=RS;
  36.         a=RD;
  37.         if (a>=TOOBIG)
  38.         {
  39.             convert(a,w);
  40.             return TOOBIG;
  41.         }
  42.         return a;
  43.     }
  44.     else if (n==oldn) return a;
  45.     t=v;
  46.     v=a*(ud-u)+vd;
  47.     vd=t;
  48.     a=u/v;
  49.     ud=u;
  50.     u=2*RD-u%v;
  51.     oldn=n;
  52.     if (a>=TOOBIG)
  53.     {
  54.         convert(a,w);
  55.         return TOOBIG;
  56.     }
  57.     return a;
  58. }
  59.  
  60. void fpower(x,n,w)
  61. flash x,w;
  62. int n;
  63. { /* raise floating-slash number to integer power w=x^n */
  64.     copy(x,w8);
  65.     zero(w);
  66.     if (ERNUM || size(w8)==0) return;
  67.     convert(1,w);
  68.     if (n==0) return;
  69.     depth++;
  70.     trace[depth]=51;
  71.     if (TRACER) track();
  72.     if (n<0)
  73.     {
  74.         n=(-n);
  75.         frecip(w8,w8);
  76.     }
  77.     if (n==1)
  78.     {
  79.         copy(w8,w);
  80.         depth--;
  81.         return;
  82.     }
  83.     forever
  84.     {
  85.         if (n%2!=0) fmul(w,w8,w);
  86.         n/=2;
  87.         if (ERNUM || n==0) break;
  88.         fmul(w8,w8,w8);
  89.     }
  90.     depth--;
  91. }
  92.  
  93. bool froot(x,n,w)
  94. flash w,x;
  95. int n;
  96. { /* extract nth root of x  - w=x^(1/n) using Newtons method */
  97.     bool minus,rn,rm,hack;
  98.     int nm,dn,s,op[5];
  99.     copy(x,w);
  100.     if (ERNUM || n==1) return TRUE;
  101.     if (n==(-1))
  102.     {
  103.         frecip(w,w);
  104.         return TRUE;
  105.     }
  106.     depth++;
  107.     trace[depth]=52;
  108.     if (TRACER) track();
  109.     minus=FALSE;
  110.     if (n<0)
  111.     {
  112.         minus=TRUE;
  113.         n=(-n);
  114.     }
  115.     s=exsign(w);
  116.     if (n%2==0 && s==MINUS)
  117.     {
  118.         berror(9);
  119.         depth--;
  120.         return FALSE;
  121.     }
  122.     insign(PLUS,w);
  123.     numer(w,w8);
  124.     denom(w,w9);
  125.     rn=root(w8,n,w8);
  126.     rm=root(w9,n,w9);
  127.     if (rn && rm)
  128.     {
  129.         fpack(w8,w9,w);
  130.         if (minus) frecip(w,w);
  131.         insign(s,w);
  132.         depth--;
  133.         return TRUE;
  134.     }
  135.     nm=size(w8);
  136.     dn=size(w9);
  137.     if (n==2 && ((nm<TOOBIG) || rn) && ((dn<TOOBIG) || rm))
  138.     {
  139.         if (!rn && nm<TOOBIG)
  140.         {
  141.             multiply(w8,w8,w8);
  142.             numer(w,w7);
  143.             subtract(w7,w8,w8);
  144.             RS=w8[1]+base*w8[2];
  145.             RD=nm;
  146.             build(w8,quad);
  147.         }
  148.         if (!rm && dn<TOOBIG)
  149.         {
  150.             multiply(w9,w9,w9);
  151.             denom(w,w7);
  152.             subtract(w7,w9,w9);
  153.             RS=w9[1]+base*w9[2];
  154.             RD=dn;
  155.             build(w9,quad);
  156.         }
  157.         if (size(w9)==1) copy(w8,w);
  158.         else             fdiv(w8,w9,w);
  159.         if (minus) frecip(w,w);
  160.         insign(s,w);
  161.         depth--;
  162.         return FALSE;
  163.     }
  164.     hack=FALSE;
  165.     if (lent(w)<=2)
  166.     { /* for 'simple' w only */
  167.         hack=TRUE;
  168.         fpi(w15);
  169.         fpmul(w15,1,3,w10);
  170.         fpower(w10,n,w10);
  171.         fmul(w,w10,w);
  172.     }
  173.     op[0]=0x6C;  /* set up for [(n-1).x+y]/n */
  174.     op[1]=n-1;
  175.     op[2]=1;
  176.     op[3]=n;
  177.     op[4]=0;
  178.     workprec=stprec;
  179.     dconv(pow(fdsize(w),1.0/(double)n),w10);
  180.     while (workprec!=nib)
  181.     { /* Newtons iteration w10=(w/w10^(n-1)+(n-1)*w10)/n */
  182.         if (workprec<nib) workprec*=2;
  183.         if (workprec>=nib) workprec=nib;
  184.         else if (workprec*2>nib) workprec=(nib+1)/2;
  185.         fpower(w10,n-1,w9);
  186.         fdiv(w,w9,w9);
  187.         flop(w10,w9,op,w10);
  188.     }
  189.     copy(w10,w);
  190.     op[0]=0x48;
  191.     op[1]=3;
  192.     op[3]=1;
  193.     op[2]=op[4]=0;
  194.     if (hack) flop(w,w15,op,w);
  195.     if (minus) frecip(w,w);
  196.     insign(s,w);
  197.     depth--;
  198.     return FALSE;
  199. }
  200.