home *** CD-ROM | disk | FTP | other *** search
/ The C Users' Group Library 1994 August / wc-cdrom-cusersgrouplibrary-1994-08.iso / vol_200 / 247_01 / bnbuild.c < prev    next >
Text File  |  1989-04-19  |  3KB  |  130 lines
  1. /*
  2.  *   MIRACL flash number builder: uses generator of
  3.  *   regular continued fraction expansion to create
  4.  *   a flash number, rounded if necessary.
  5.  *   bnbuild.c
  6.  */
  7.  
  8. #include <stdio.h>
  9. #include "miracl.h"
  10. #define abs(x)  ((x)<0? (-(x)) : (x))
  11.  
  12. /* Access global variables */
  13.  
  14. extern int  depth;    /* error tracing .. */
  15. extern int  trace[];  /* .. mechanism     */
  16. extern int  check;    /* error checking   */
  17. extern small base;    /* number base      */
  18. extern int  nib;      /* no. ints in big  */
  19. extern int  workprec;
  20.  
  21. extern big w0;
  22. extern big w1;       /* workspace variables  */
  23. extern big w2;
  24. extern big w3;
  25. extern big w4;
  26. extern big w7;
  27.  
  28. void build(x,gen)
  29. flash x;
  30. int (*gen)();
  31. { /* Build x from its regular c.f. *
  32.    * generated by gen()            */
  33.     small a,b,c,d,rm,ex1,ex2,ex;
  34.     int q,n,prc,lw2,lw4,lz;
  35.     bool finoff,last;
  36.     big t;
  37.     if (ERNUM) return;
  38.     depth++;
  39.     trace[depth]=48;
  40.     if (TRACER) track();
  41.     zero(w1);
  42.     convert((small)1,w2);
  43.     convert((small)1,w3);
  44.     zero(w4);
  45.     finoff=FALSE;
  46.     last=FALSE;
  47.     n=0;
  48.     q=(*gen)(x,n);   /* Note - first quotient may be zero */
  49.     ex=base;
  50.     if (nib==workprec) prc=nib;
  51.     else prc=workprec+1;
  52.     while (!ERNUM && q>=0)
  53.     {
  54.         if (q==TOOBIG || n==0 || finoff)
  55.         {
  56.             if (q!=TOOBIG) convert((small)q,x);
  57.             else last=FALSE;
  58.             check=OFF;
  59.             multiply(w2,x,w0);
  60.             subtract(w1,w0,w7);
  61.             check=ON;
  62.             if (abs(w7[0])>nib) break;
  63.             copy(w7,w1);
  64.             t=w1,w1=w2,w2=t;   /* swap(w1,w2) */
  65.             check=OFF;
  66.             multiply(w4,x,w0);
  67.             subtract(w3,w0,w7);
  68.             check=ON;
  69.             if (abs(w7[0])>nib) break;
  70.             copy(w7,w3);
  71.             t=w3,w3=w4,w4=t;   /* swap(w3,w4) */
  72.             n++;
  73.         }
  74.         lw2=abs(w2[0]);
  75.         lw4=abs(w4[0]);
  76.         lz=lw2+lw4;
  77.         if (lz > prc) break;  /* too big - exit */
  78.         if (last)
  79.         {
  80.             if (finoff) break;
  81.             finoff=TRUE;
  82.             q=(*gen)(x,n);
  83.             continue;
  84.         }
  85.         if (lz>=prc-1)
  86.         { /* nearly finished - so be careful not to overshoot */
  87.             ex1=base/(w2[lw2]+1);
  88.             ex2=base/(w4[lw4]+1);
  89.             if (ex2>ex1) ex=ex1,ex1=ex2,ex2=ex;
  90.             if (lz==prc) ex=ex2;
  91.             else         ex=ex1;
  92.             last=TRUE;
  93.         }
  94.         a=1;
  95.         b=0;
  96.         c=0;
  97.         d=1;
  98.         forever
  99.         {
  100.             q=(*gen)(x,n);
  101.             if (q<0 || q==TOOBIG || q>=MAXBASE/abs(d)) break;
  102.             rm=b-q*d;
  103.             b=d;
  104.             d=rm;
  105.             rm=a-q*c;
  106.             a=c;
  107.             c=rm;
  108.             n++;
  109.             if (abs(c-d)>ex) break;
  110.         }
  111.         premult(w1,c,w7);
  112.         premult(w1,a,w1);
  113.         premult(w2,b,w0);
  114.         premult(w2,d,w2);
  115.         add(w1,w0,w1);
  116.         add(w2,w7,w2);
  117.         premult(w3,c,w7);
  118.         premult(w3,a,w3);
  119.         premult(w4,b,w0);
  120.         premult(w4,d,w4);
  121.         add(w3,w0,w3);
  122.         add(w4,w7,w4);
  123.     }
  124.     if (abs(w2[0]-w4[0]) <= nib) fpack(w2,w4,x);
  125.     else                         fpack(w1,w3,x);
  126.     negate (x,x);
  127.     if (q!=(-1)) EXACT=FALSE;
  128.     depth--;
  129. }
  130.