home *** CD-ROM | disk | FTP | other *** search
/ The C Users' Group Library 1994 August / wc-cdrom-cusersgrouplibrary-1994-08.iso / vol_200 / 247_01 / bncore.c < prev    next >
Text File  |  1989-04-19  |  14KB  |  698 lines
  1. /*
  2.  *   MIRACL Core module - contains initialisation code and general purpose 
  3.  *   utilities.
  4.  *   bncore.c
  5.  */
  6.  
  7. #include <stdio.h>
  8. #include "mirdef.h"
  9. #define abs(x)  ((x)<0? (-(x)) : (x))
  10. #define sign(x) ((x)<0? (-1) : 1)
  11.  
  12. #define NK      55
  13. #define NJ      24
  14. #define NV      21
  15.                    /* Constants for random number   *
  16.                     * generator. See Knuth Chapt. 3 */
  17. typedef small *big;
  18. typedef small *flash;
  19.  
  20. /* Declare global variables */
  21.  
  22. small base;        /* number base     */
  23. small apbase;      /* apparent base   */
  24. int pack;          /* packing density */
  25. int lg2b;          /* bits in base    */
  26. small base2;       /* 2^lg2b          */
  27. int nib;           /* length of bigs  */
  28. int workprec;
  29. int stprec;        /* start precision */
  30. int depth;         /* error tracing ..*/
  31. int trace[20];     /* .. mechanism    */
  32. bool check;        /* overflow check  */
  33. long ira[NK];      /* random number   */
  34. int rndptr;        /* array & pointer */
  35. char *names[] =
  36. {"your program","innum","otnum","cinbin","cotbin",
  37. "multiply","divide","incr","decr","premult",
  38. "subdiv","fdsize","gcd","cbase",
  39. "cinnum","cotnum","root","power",
  40. "powmod","bigdig","bigrand","nxprime","prime",
  41. "mirvar","mad","mirsys","putdig",
  42. "add","subtract","select","xgcd",
  43. "fpack","dconv","shift","round","fmul",
  44. "fdiv","fadd","fsub","fcomp","fconv",
  45. "frecip","fpmul","fincr","fsize","ftrunc",
  46. "frand","strongp","build","log2","exp2",
  47. "fpower","froot","fpi","fexp","flog","fpowf",
  48. "ftan","fatan","fsin","fasin","fcos","facos",
  49. "ftanh","fatanh","fsinh","fasinh","fcosh",
  50. "facosh","flop","gprime"}; /* 70 */
  51.  
  52. big w0;            /* workspace bigs  */
  53. big w1,w2,w3,w4;
  54. big w5,w6,w7;
  55. big w8,w9,w10,w11;
  56. big w12,w13,w14,w15;
  57.  
  58. unsigned char IBUFF[IBSIZ+1]; /* input buffer    */
  59. unsigned char OBUFF[OBSIZ+1]; /* output buffer   */
  60. bool ERCON;        /* error control   */
  61. int  ERNUM;        /* last error code */
  62. bool STROUT;       /* redirect output */
  63. int  LINE;         /* Output line length */
  64. int  IOBASE;       /* base for input and output */
  65. bool WRAP;         /* =ON for wrapround on output, otherwise =OFF       */
  66. bool EXACT;        /* exact flag      */
  67. bool POINT;        /* =ON for radix point, =OFF for fractions in output */
  68. int  NTRY;         /* no. of tries for probablistic primality testing   */
  69. bool TRACER;       /* turns trace tracker on/off */
  70. int  INPLEN;       /* input length               */
  71.  
  72. extern char *calloc();
  73.  
  74. void berror(nerr)
  75. int nerr;
  76. {  /*  Big number error routine  */
  77. int i;
  78. if (ERCON)
  79. {
  80.     ERNUM=nerr;
  81.     return;
  82. }
  83. printf("MIRACL error from routine %s\n",
  84. names[trace[depth]]);
  85. for (i=depth-1;i>=0;i--)
  86. printf("              called from %s\n",
  87. names[trace[i]]);
  88. switch (nerr)
  89. {
  90. case 1 :
  91. printf("Number base too big for representation\n");
  92. exit(0);
  93. case 2 :
  94. printf("Division by zero attempted\n");
  95. exit(0);
  96. case 3 : 
  97. printf("Overflow\n");
  98. exit(0);
  99. case 4 :
  100. printf("Internal result is negative\n");
  101. exit(0);
  102. case 5 : 
  103. printf("Input format error\n");
  104. exit(0);
  105. case 6 :
  106. printf("Illegal number base\n");
  107. exit(0);
  108. case 7 : 
  109. printf("Illegal parameter usuage\n");
  110. exit(0);
  111. case 8 :
  112. printf("Out of space\n");
  113. exit(0);
  114. case 9 :
  115. printf("Even root of a negative number\n");
  116. exit(0);
  117. case 10:
  118. printf("Raising integer to negative power\n");
  119. exit(0);
  120. case 11:
  121. printf("Attempt to take illegal root\n");
  122. exit(0);
  123. case 12:
  124. printf("Integer operation attempted on Flash number\n");
  125. exit(0);
  126. case 13:
  127. printf("Flash overflow\n");
  128. exit(0);
  129. case 14:
  130. printf("Numbers too big\n");
  131. exit(0);
  132. case 15:
  133. printf("Log of a non-positive number\n");
  134. exit(0);
  135. case 16:
  136. printf("Double to flash conversation failure\n");
  137. exit(0);
  138. case 17:
  139. printf("I/O buffer overflow\n");
  140. exit(0);
  141. }
  142. }
  143.  
  144. void track()
  145. { /* track course of program execution *
  146.    * through the MIRACL routines       */
  147.     int i;
  148.     for (i=0;i<depth;i++) putchar('-');
  149.     putchar('>');
  150.     printf("%s\n",names[trace[depth]]);
  151. }
  152.  
  153. small brand()
  154. { /* random number generator - uses Knuth's *
  155.    * subtractive method                     */
  156.     int i;
  157.     small t;
  158.     long k;   
  159.     rndptr++;
  160.     t=(small)ira[rndptr];
  161.     if (rndptr<NK) return abs(t);
  162.     rndptr=0;
  163.     for (i=0;i<NJ;i++)
  164.     {
  165.         k=ira[i]-ira[i+NK-NJ];
  166.         if (k<0) k+=MAXNUM;
  167.         ira[i]=k;
  168.     }
  169.     for (i=NJ;i<NK;i++)
  170.     {
  171.         k=ira[i]-ira[i-NJ];
  172.         if (k<0) k+=MAXNUM;
  173.         ira[i]=k;
  174.     }
  175.     t=(small)ira[0];
  176.     return abs(t);
  177. }
  178.  
  179. void irand(seed)
  180. long seed;
  181. { /* initialise random number system */
  182.     int i,in;
  183.     long k;
  184.     k=1L;
  185.     ira[0]=seed;
  186.     for (i=1;i<NK;i++)
  187.     {
  188.         in=(NV*i)%NK;
  189.         ira[in]=k;
  190.         k=seed-k;
  191.         if (k<0) k+=MAXNUM;
  192.         seed=ira[in];
  193.     }
  194.     for (i=1;i<4;i++)
  195.     { /* warm-up the generator */
  196.         rndptr=NK;
  197.         brand();
  198.     }
  199. }
  200.  
  201. void setbase(nb)
  202. small nb;
  203. {  /* set base. Pack as many digits as  *
  204.     * possible into each computer word  */
  205.     small temp;
  206.     apbase=nb;
  207.     pack=1;
  208.     base=nb;
  209.     temp=MAXBASE/nb;
  210.     while (temp>=nb)
  211.     {
  212.         temp/=nb;
  213.         base*=nb;
  214.         pack++;
  215.     }
  216. }
  217.  
  218. bool fit(x,y,f)
  219. big x,y;
  220. int f;
  221. { /* returns TRUE if x/y would fit flash format of length f */
  222.     int n,d;
  223.     n=abs(x[0]);
  224.     d=abs(y[0]);
  225.     if (n==1 && x[1]==1) n=0;
  226.     if (d==1 && y[1]==1) d=0;
  227.     if (n+d<=f) return TRUE;
  228.     return FALSE;
  229. }
  230.  
  231. int lent(x)
  232. flash x;
  233. { /* return length of big or flash in words */
  234.     small lx;
  235.     lx=abs(x[0]);
  236.     return (int)((lx&MSK)+((lx>>BTS)&MSK));
  237. }
  238.  
  239. void zero(x)
  240. flash x;
  241. { /* set big/flash number to zero */
  242.     int i,n;
  243.     if (ERNUM) return;
  244.     n=lent(x);
  245.     for (i=0;i<=n;i++)
  246.         x[i]=0;
  247. }
  248.  
  249. void convert(n,x)
  250. big x;
  251. small n;
  252. {  /*  convert integer n to big number format  */
  253.     int m,s;
  254.     zero(x);
  255.     if (n==0) return;
  256.     s=sign(n);
  257.     n=abs(n);
  258.     m=0;
  259.     while (n>0)
  260.     {
  261.         m++;
  262.         x[m]=n%base;
  263.         n/=base;
  264.     }
  265.     x[0]=m*s;
  266. }
  267.  
  268. void lconv(n,x)
  269. big x;
  270. long n;
  271. { /* convert long integer to big number format */
  272.     int m,s;
  273.     zero(x);
  274.     if (n==0) return;
  275.     s=sign(n);
  276.     n=abs(n);
  277.     m=0;
  278.     while (n>0)
  279.     {
  280.         m++;
  281.         x[m]=n%base;
  282.         n/=base;
  283.     }
  284.     x[0]=m*s;
  285. }
  286.  
  287. flash mirvar(iv)
  288. small iv;
  289. { /* initialize big/flash number */
  290.     flash x;
  291.     if (ERNUM) return x;
  292.     depth++;
  293.     trace[depth]=23;
  294.     if (TRACER) track();
  295.     x=(small *)calloc(nib+1,sizeof(small));
  296.     if (x==NULL)
  297.     {
  298.         berror(8);
  299.         depth--;
  300.         return x;
  301.     }
  302.     convert(iv,x);
  303.     depth--;
  304.     return x;
  305. }
  306.  
  307. void mirsys(nd,nb)
  308. int nd;
  309. small nb;
  310. {  /*  Initialize MIRACL system to   *
  311.     *  use numbers to base nb, and   *
  312.     *  nd digits or (-nd) bytes long */
  313.     small b;
  314.     depth=0;
  315.     trace[0]=0;
  316.     depth++;
  317.     trace[depth]=25;
  318.     if (nb<2 || nb>MAXBASE)
  319.     {
  320.         berror(6);
  321.         depth--;
  322.         return;
  323.     }
  324.     setbase(nb);
  325.     b=base;
  326.     lg2b=0;
  327.     base2=1;
  328.     while (b>1)
  329.     {
  330.         b/=2;
  331.         lg2b++;
  332.         base2*=2;
  333.     }
  334.     if (nd>0)
  335.         nib=(nd-1)/pack+1;
  336.     else
  337.         nib=(lg2b-8*nd-1)/lg2b;
  338.     if (nib<2) nib=2;
  339.     workprec=nib;
  340.     stprec=nib;
  341.     while (stprec>WPERD) stprec=(stprec+1)/2;
  342.     check=ON;
  343.     IOBASE=10;   /* defaults */
  344.     WRAP=ON;
  345.     STROUT=FALSE;
  346.     ERCON=FALSE;
  347.     ERNUM=0;
  348.     POINT=OFF;
  349.     NTRY=6;           /* <=27 */
  350.     EXACT=TRUE;
  351.     IBUFF[0]='\0';
  352.     TRACER=OFF;
  353.     LINE=80;
  354.     INPLEN=0;
  355.     irand(0L);     /* start rand number generator */
  356.     nib*=2;
  357.     if (nib!=(nib&MSK) || nib>TOOBIG)
  358.     {
  359.         berror(14);
  360.         nib/=2;
  361.         depth--;
  362.         return;
  363.     }
  364.     w0=mirvar((small)0); /* w0 is double length  */
  365.     nib/=2;
  366.     w1=mirvar((small)0); /* initialize workspace */
  367.     w2=mirvar((small)0);
  368.     w3=mirvar((small)0);
  369.     w4=mirvar((small)0);
  370.     nib*=2;              /* double length */
  371.     w5=mirvar((small)0);
  372.     w6=mirvar((small)0);
  373.     w7=mirvar((small)0);
  374.     nib/=2;
  375.     w8=mirvar((small)0);
  376.     w9=mirvar((small)0);
  377.     w10=mirvar((small)0);
  378.     w11=mirvar((small)0);
  379.     w12=mirvar((small)0);
  380.     w13=mirvar((small)0);
  381.     w14=mirvar((small)0);
  382.     w15=mirvar((small)0); /* used to store pi */
  383.     depth--;
  384. } 
  385.  
  386. int exsign(x)
  387. flash x;
  388. { /* extract sign of big/flash number */
  389.     return sign(x[0]);
  390. }
  391.  
  392. void insign(s,x)
  393. flash x;
  394. int s;
  395. {  /* assert sign of big/flash number */
  396.     x[0]=s*abs(x[0]);
  397. }   
  398.  
  399. void lzero(x)
  400. big x;
  401. {  /*  strip leading zeros from big number  */
  402.     int m;
  403.     bool neg;
  404.     if (x[0]<0)
  405.     {
  406.         m=(-x[0]);
  407.         neg=TRUE;
  408.     }
  409.     else
  410.     {
  411.         m=x[0];
  412.         neg=FALSE;
  413.     }
  414.     while (m>0 && x[m]==0)
  415.         m--;
  416.     if (neg) x[0]=(-m);
  417.     else     x[0]=m;
  418. }
  419.  
  420. int getdig(x,i)
  421. big x;
  422. int i;
  423. {  /* extract a packed digit */
  424.     int k;
  425.     small n;
  426.     i--;
  427.     n=x[i/pack+1];
  428.     k=i%pack;
  429.     for (i=1;i<=k;i++)
  430.         n/=apbase;
  431.     return n%apbase;
  432. }
  433.  
  434. int numdig(x)
  435. big x;
  436. {  /* returns number of digits in x */
  437.     int nd;
  438.     nd=abs(x[0])*pack;
  439.     while (getdig(x,nd)==0)
  440.         nd--;
  441.     return nd;
  442. } 
  443.  
  444. void putdig(n,x,i)  
  445. big x;
  446. int i;
  447. int n;
  448. {  /* insert a digit into a packed word */
  449.     int j,k,lx;
  450.     small m,p;
  451.     if (ERNUM) return;
  452.     depth++;
  453.     trace[depth]=26;
  454.     if (TRACER) track();
  455.     lx=abs(x[0]);
  456.     m=getdig(x,i);
  457.     p=n;
  458.     i--;
  459.     j=i/pack+1;
  460.     k=i%pack;
  461.     for (i=1;i<=k;i++)
  462.     {
  463.         m*=apbase;
  464.         p*=apbase;
  465.     }
  466.     if (check && j>nib)
  467.     {
  468.         berror(3);
  469.         depth--;
  470.         return;
  471.     }
  472.     x[j]=x[j]-m+p;
  473.     if (j>lx) x[0]=j*exsign(x);
  474.     lzero(x);
  475.     depth--;
  476. }
  477.  
  478. void copy(x,y)
  479. flash x;
  480. flash y;
  481. {  /* copy x to y: y=x  */
  482.     int i,nx,ny;
  483.     if (ERNUM || x==y) return;
  484.     ny=lent(y);
  485.     nx=lent(x);
  486.     for (i=nx+1;i<=ny;i++)
  487.         y[i]=0;
  488.     for (i=0;i<=nx;i++)
  489.         y[i]=x[i];
  490. }
  491.  
  492. void negate(x,y)
  493. flash x;
  494. flash y;
  495. { /* negate a big/flash variable: y=-x */
  496.     copy(x,y);
  497.     y[0]=(-y[0]);
  498. }
  499.  
  500. void absol(x,y)
  501. flash x;
  502. flash y;
  503. { /* y=abs(x) */
  504.     copy(x,y);
  505.     y[0]=abs(y[0]);
  506. }
  507.  
  508. bool notint(x)
  509. flash x;
  510. { /* returns TRUE if x is Flash */
  511.     if (((abs(x[0])>>BTS)&MSK)!=0) return TRUE;
  512.     return FALSE;
  513. }
  514.  
  515. void shift(x,n,w)
  516. big x,w;
  517. int n;
  518. { /* set w=x.(base^n) by shifting */
  519.     int i,bl;
  520.     copy(x,w);
  521.     if (ERNUM || w[0]==0 || n==0) return;
  522.     depth++;
  523.     trace[depth]=33;
  524.     if (TRACER) track();
  525.     if (notint(w)) berror(12);
  526.     bl=abs(w[0])+n;
  527.     if (check && bl>nib) berror(3);
  528.     if (ERNUM)
  529.     {
  530.        depth--;
  531.        return;
  532.     }
  533.     if (n>0)
  534.     {
  535.         for (i=bl;i>n;i--)
  536.             w[i]=w[i-n];
  537.         for (i=1;i<=n;i++)
  538.             w[i]=0;
  539.     }
  540.     else
  541.     {
  542.         n=(-n);
  543.         for (i=1;i<=bl;i++)
  544.             w[i]=w[i+n];
  545.         for (i=1;i<=n;i++)
  546.             w[bl+i]=0;
  547.     }
  548.     w[0]=bl*exsign(w);
  549.     depth--;
  550. }
  551.  
  552. int size(x)
  553. big x;
  554. {  /*  get size of big number;  convert to *
  555.     *  integer - if possible               */
  556.     int s;
  557.     small m;
  558.     s=sign(x[0]);
  559.     m=abs(x[0]);
  560.     if (m==0) return 0;
  561.     if (m==1 && x[1]<TOOBIG) return s*x[1];
  562.     return s*TOOBIG;
  563. }
  564.  
  565. int compare(x,y)
  566. big x;
  567. big y;
  568. {  /* compare x and y: =1 if x>y  =-1 if x<y *
  569.     *  =0 if x=y                             */
  570.     int m,s;
  571.     if (x==y) return 0;
  572.     if (x[0]>y[0]) return 1;
  573.     if (x[0]<y[0]) return -1;
  574.     s=sign(x[0]);
  575.     m=abs(x[0]);
  576.     while (m>0)
  577.     { /* check digit by digit */
  578.         if (x[m]>y[m]) return s;
  579.         if (x[m]<y[m]) return -s;
  580.         m--;
  581.     }
  582.     return 0;
  583. }
  584.  
  585. void fpack(n,d,x)
  586. big n,d;
  587. flash x;
  588. { /* create floating-slash number x=n/d from *
  589.    * big integer numerator and denominator   */
  590.     int i,s,ld,ln;
  591.     if (ERNUM) return;
  592.     depth++;
  593.     trace[depth]=31;
  594.     if (TRACER) track();
  595.     ld=abs(d[0]);
  596.     if (ld==0) berror(13);
  597.     if (x==d) berror(7);
  598.     if (notint(n) || notint(d)) berror(12);
  599.     s=exsign(n);
  600.     ln=abs(n[0]);
  601.     if (abs(size(n))==1) ln=0;
  602.     if (ld+ln>nib) berror(13);
  603.     if (ERNUM)
  604.     {
  605.        depth--;
  606.        return;
  607.     }
  608.     copy(n,x);
  609.     if (size(n)==0)
  610.     {
  611.         depth--;
  612.         return;
  613.     }
  614.     s*=exsign(d);
  615.     if (abs(size(d))==1)
  616.     {
  617.         insign(s,x);
  618.         depth--;
  619.         return;
  620.     }
  621.     for (i=1;i<=ld;i++)
  622.         x[ln+i]=d[i];
  623.     x[0]=s*(ln+((small)ld<<BTS));
  624.     depth--;
  625. }
  626.  
  627. void numer(x,y)
  628. flash x;
  629. big y;
  630. { /* extract numerator of x */
  631.     int i,s,ln,ld;
  632.     small ly;
  633.     if (ERNUM) return;
  634.     copy(x,y);
  635.     if (notint(y))
  636.     {
  637.         s=exsign(y);
  638.         ly=abs(y[0]);
  639.         ln=(ly&MSK);
  640.         if (ln==0)
  641.         {
  642.             convert((small)s,y);
  643.             return;
  644.         }
  645.         ld=((ly>>BTS)&MSK);
  646.         for (i=1;i<=ld;i++)
  647.             y[ln+i]=0;
  648.         y[0]=ln*s;
  649.     }
  650. }
  651.  
  652. void denom(x,y)
  653. flash x;
  654. big y;
  655. { /* extract denominator of x */
  656.     int i,ln,ld;
  657.     small ly;
  658.     if (ERNUM) return;
  659.     if (!notint(x))
  660.     {
  661.         convert((small)1,y);
  662.         return;
  663.     }
  664.     copy(x,y);
  665.     ly=abs(y[0]);
  666.     ln=(ly&MSK);
  667.     ld=((ly>>BTS)&MSK);
  668.     y[0]=ld;
  669.     if (ln==0) return;
  670.     for (i=1;i<=ld;i++)
  671.         y[i]=y[ln+i];
  672.     for (i=1;i<=ln;i++)
  673.         y[ld+i]=0;
  674. }
  675.  
  676. small igcd(x,y)
  677. small x,y;
  678. { /* integer GCD, returns GCD of x and y */
  679.     small r;
  680.     if (y==0) return x;
  681.     while ((r=x%y)!=0)
  682.         x=y,y=r;
  683.     return y;
  684. }
  685.  
  686. small isqrt(num,guess)
  687. small num,guess;
  688. { /* square root of an integer */
  689.     small sqr;
  690.     forever
  691.     { /* Newtons iteration */
  692.         sqr=guess+(((num/guess)-guess)/2);
  693.         if (sqr==guess) return sqr;
  694.         guess=sqr;
  695.     }
  696. }
  697.  
  698.