home *** CD-ROM | disk | FTP | other *** search
/ The C Users' Group Library 1994 August / wc-cdrom-cusersgrouplibrary-1994-08.iso / vol_200 / 247_01 / bnflash.c < prev    next >
Text File  |  1989-04-19  |  6KB  |  300 lines
  1. /*
  2.  *   MIRACL floating-Slash arithmetic
  3.  *   bnflash.c
  4.  */
  5.  
  6. #include <stdio.h>
  7. #include "miracl.h"
  8. #define abs(x)  ((x)<0? (-(x)) : (x))
  9.  
  10. /* Access global variables */
  11.  
  12. extern int  depth;    /* error tracing ..*/
  13. extern int  trace[];  /* .. mechanism    */
  14. extern bool check;    /* overflow check  */
  15. extern int  nib;      /* length of bigs  */
  16.  
  17. extern big w1;        /* workspace variables  */
  18. extern big w2;
  19. extern big w3;
  20. extern big w4;
  21. extern big w5;
  22. extern big w6;
  23. extern big w0;
  24.  
  25. void flop(x,y,op,z)
  26. flash x,y,z;
  27. int *op;
  28. { /* Do basic flash operation - depending on op[] *
  29.    * Performs operations of the form
  30.  
  31.                  A.f1(x,y) + B.f2(x,y)
  32.                  -------------------
  33.                  C.f3(x,y) + D.f4(x,y)
  34.  
  35.    * Four functions f(x,y) are supported and      *
  36.    * coded thus                                   *
  37.    *              00 - Nx.Ny                      *
  38.    *              01 - Nx.Dy                      *
  39.    *              10 - Dx.Ny                      *
  40.    *              11 - Dx.Dy                      *
  41.    * where Nx is numerator of x, Dx denominator   *
  42.    * of x, etc.                                   *
  43.    * op[0] contains the codes for f1-f4 in bottom *
  44.    * eight bits =  00000000f1f2f3f4               *
  45.    * op[1], op[2], op[3], op[4] contain the       *
  46.    * single precision multipliers A, B, C, D      */
  47.     int i,code;
  48.     if (ERNUM) return;
  49.     depth++;
  50.     trace[depth]=69;
  51.     if (TRACER) track();
  52.     numer(x,w1);
  53.     denom(x,w2);
  54.     numer(y,w3);
  55.     denom(y,w4);
  56.     check=OFF;
  57.     for (i=1;i<=4;i++)
  58.     {
  59.         zero(w0);
  60.         if (op[i]==0) continue;
  61.         code=(op[0]>>(2*(4-i)))&3;
  62.         switch (code)
  63.         {
  64.         case 0 : if (x==y) multiply(w1,w1,w0);
  65.                  else      multiply(w1,w3,w0);
  66.                  break;
  67.         case 1 : multiply(w1,w4,w0);
  68.                  break;
  69.         case 2 : multiply(w2,w3,w0);
  70.                  break;
  71.         case 3 : if(x==y) multiply(w2,w2,w0);
  72.                  else     multiply(w2,w4,w0);
  73.                  break;
  74.         }
  75.         premult(w0,(small)op[i],w0);
  76.         switch (i)
  77.         {
  78.         case 1 : copy(w0,w5);
  79.                  break;
  80.         case 2 : add(w5,w0,w5);
  81.                  break;
  82.         case 3 : copy(w0,w6);
  83.                  break;
  84.         case 4 : add(w6,w0,w6);
  85.                  break;
  86.         }
  87.     }
  88.     check=ON;
  89.     round(w5,w6,z);
  90.     depth--;
  91. }
  92.  
  93. void fmul(x,y,z)
  94. flash x,y,z;
  95. { /* Flash multiplication - z=x*y */
  96.     int op[5];
  97.     if (ERNUM) return;
  98.     depth++;
  99.     trace[depth]=35;
  100.     if (TRACER) track();
  101.     op[0]=0x0C;
  102.     op[1]=op[3]=1;
  103.     op[2]=op[4]=0;
  104.     flop(x,y,op,z);
  105.     depth--;
  106. }
  107.  
  108. void fdiv(x,y,z)
  109. flash x,y,z;
  110. { /* Flash divide - z=x/y */
  111.     int op[5];
  112.     if (ERNUM) return;
  113.     depth++;
  114.     trace[depth]=36;
  115.     if (TRACER) track();
  116.     op[0]=0x48;
  117.     op[1]=op[3]=1;
  118.     op[2]=op[4]=0;
  119.     flop(x,y,op,z);
  120.     depth--;
  121. }
  122.  
  123. void fadd(x,y,z)
  124. flash x,y,z;
  125. { /* Flash add - z=x+y */
  126.     int op[5];
  127.     if (ERNUM) return;
  128.     depth++;
  129.     trace[depth]=37;
  130.     if (TRACER) track();
  131.     op[0]=0x6C;
  132.     op[1]=op[2]=op[3]=1;
  133.     op[4]=0;   
  134.     flop(x,y,op,z);
  135.     depth--;
  136. }
  137.  
  138. void fsub(x,y,z)
  139. flash x,y,z;
  140. { /* Flash subtract - z=x-y */
  141.     int op[5];
  142.     if (ERNUM) return;
  143.     depth++;
  144.     trace[depth]=38;
  145.     if (TRACER) track();
  146.     op[0]=0x6C;
  147.     op[1]=op[3]=1;
  148.     op[2]=(-1);
  149.     op[4]=0;
  150.     flop(x,y,op,z);
  151.     depth--;
  152. }
  153.  
  154. int fcomp(x,y)
  155. flash x,y;
  156. { /* compares two Flash numbers             *
  157.    * returns -1 if y>x; +1 if x>y; 0 if x=y */
  158.     if (ERNUM) return 0;
  159.     depth++;
  160.     trace[depth]=39;
  161.     if (TRACER) track();
  162.     numer(x,w1);
  163.     denom(y,w2);
  164.     check=OFF;
  165.     multiply(w1,w2,w5);
  166.     numer(y,w1);
  167.     denom(x,w2);
  168.     multiply(w1,w2,w0);
  169.     check=ON;
  170.     depth--;
  171.     return (compare(w5,w0));
  172. }
  173.  
  174. void ftrunc(x,y,z)
  175. flash x,z;
  176. big y;
  177. { /* sets y=int(x), z=rem(x) - returns *
  178.    * y only for ftrunc(x,y,y)          */
  179.     if (ERNUM) return;
  180.     depth++;
  181.     trace[depth]=45;
  182.     if (TRACER) track();
  183.     numer(x,w1);
  184.     denom(x,w2);
  185.     divide(w1,w2,w3);
  186.     copy(w3,y);
  187.     if (y!=z) fpack(w1,w2,z);
  188.     depth--;
  189. }
  190.  
  191. void fconv(n,d,x)
  192. small n,d;
  193. flash x;
  194. { /* convert simple fraction n/d to Flash form */
  195.     small g;
  196.     if (ERNUM) return;
  197.     depth++;
  198.     trace[depth]=40;
  199.     if (TRACER) track();
  200.     if (d<0)
  201.     {
  202.         d=(-d);
  203.         n=(-n);
  204.     }
  205.     g=igcd(n,d);
  206.     n/=g;
  207.     d/=g;
  208.     convert(n,w5);
  209.     convert(d,w6);
  210.     fpack(w5,w6,x);
  211.     depth--;
  212. }
  213.  
  214. void frecip(x,y)
  215. flash x,y;
  216. { /* set y=1/x */
  217.     if (ERNUM) return;
  218.     depth++;
  219.     trace[depth]=41;
  220.     if (TRACER) track();
  221.     numer(x,w1);
  222.     denom(x,w2);
  223.     fpack(w2,w1,y);
  224.     depth--;
  225. }
  226.  
  227. void fpmul(x,n,d,y)
  228. flash x,y;
  229. small n,d;
  230. { /* multiply x by small fraction n/d - y=x*n/d */
  231.     small r,g;
  232.     if (ERNUM) return;
  233.     if (n==0 || size(x)==0)
  234.     {
  235.         zero(y);
  236.         return;
  237.     }
  238.     if (n==d)
  239.     {
  240.         copy(x,y);
  241.         return;
  242.     }
  243.     depth++;
  244.     trace[depth]=42;
  245.     if (TRACER) track();
  246.     if (d<0)
  247.     {
  248.         d=(-d);
  249.         n=(-n);
  250.     }
  251.     g=igcd(n,d);
  252.     n/=g;
  253.     d/=g;
  254.     numer(x,w1);
  255.     denom(x,w2);
  256.     r=subdiv(w1,d,w3);
  257.     g=igcd(d,r);
  258.     r=subdiv(w2,n,w3);
  259.     g*=igcd(n,r);
  260.     check=OFF;
  261.     premult(w1,n,w5);
  262.     premult(w2,d,w6);
  263.     subdiv(w5,g,w5);
  264.     subdiv(w6,g,w6);
  265.     check=ON;
  266.     if (abs(w5[0])+abs(w6[0])<=nib)
  267.         fpack(w5,w6,y);
  268.     else
  269.         round(w5,w6,y);
  270.     depth--;
  271. }
  272.  
  273. void fincr(x,n,d,y)
  274. flash x,y;
  275. small n,d;
  276. { /* increment x by small fraction n/d - y=x+(n/d) */
  277.     if (ERNUM) return;
  278.     depth++;
  279.     trace[depth]=43;
  280.     if (TRACER) track();
  281.     if (d<0)
  282.     {
  283.         d=(-d);
  284.         n=(-n);
  285.     }
  286.     numer(x,w1);
  287.     denom(x,w2);
  288.     check=OFF;
  289.     premult(w1,d,w5);
  290.     premult(w2,d,w6);
  291.     premult(w2,n,w0);
  292.     add(w5,w0,w5);
  293.     check=ON;
  294.     if (d==1 && (abs(w5[0])+abs(w6[0])<=nib))
  295.         fpack(w5,w6,y);
  296.     else
  297.         round(w5,w6,y);
  298.     depth--;
  299. }
  300.