home *** CD-ROM | disk | FTP | other *** search
/ The C Users' Group Library 1994 August / wc-cdrom-cusersgrouplibrary-1994-08.iso / vol_200 / 247_01 / bnxgcd.c < prev    next >
Text File  |  1989-04-19  |  4KB  |  153 lines
  1. /*
  2.  *   MIRACL Extended Greatest Common Divisor module.
  3.  *   bnxgcd.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 small base;     /* number base     */
  15.  
  16. extern big w1;       /* workspace variables  */
  17. extern big w2;
  18. extern big w3;
  19. extern big w4;
  20. extern big w5;
  21. extern big w0;
  22.  
  23. int xgcd(x,y,xd,yd,z)
  24. big x,y,xd,yd,z;
  25. { /* greatest common divisor by Euclids method  *
  26.    * extended to also calculate xd and yd where *
  27.    *      z = x.xd + y.yd = gcd(x,y)            *
  28.    * if xd, yd not distinct, only xd calculated *
  29.    * z only returned if distinct from xd and yd *
  30.    * xd will always be positive, yd negative    */
  31.     small q,r,a,b,c,d,m;
  32.     long u,v,lq,lr;
  33.     int s,n;
  34.     bool last;
  35.     big t;
  36.     if (ERNUM) return 0;
  37.     depth++;
  38.     trace[depth]=30;
  39.     if (TRACER) track();
  40.     copy(x,w1);
  41.     copy(y,w2);
  42.     s=exsign(w1);
  43.     insign(PLUS,w1);
  44.     insign(PLUS,w2);
  45.     copy(w1,w3);
  46.     copy(w2,w4);
  47.     convert((small)1,w3);
  48.     zero(w4);
  49.     last=FALSE;
  50.     a=b=c=d=0;
  51.     while (size(w2)!=0)
  52.     {
  53.         if (b==0)
  54.         { /* update w1 and w2 */
  55.             divide(w1,w2,w5);
  56.             t=w1,w1=w2,w2=t;    /* swap(w1,w2) */
  57.             multiply(w4,w5,w0);
  58.             subtract(w3,w0,w3);
  59.             t=w3,w3=w4,w4=t;    /* swap(xd,yd) */
  60.         }
  61.         else
  62.         {
  63.             premult(w1,c,w5);
  64.             premult(w1,a,w1);
  65.             premult(w2,b,w0);
  66.             premult(w2,d,w2);
  67.             add(w1,w0,w1);
  68.             add(w2,w5,w2);
  69.             premult(w3,c,w5);
  70.             premult(w3,a,w3);
  71.             premult(w4,b,w0);
  72.             premult(w4,d,w4);
  73.             add(w3,w0,w3);
  74.             add(w4,w5,w4);
  75.         }
  76.         if (ERNUM || size(w2)==0) break;
  77.         n=w1[0];
  78.         a=1;
  79.         b=0;
  80.         c=0;
  81.         d=1;
  82.         if (n==1)
  83.         {
  84.             last=TRUE;
  85.             u=w1[1];
  86.             v=w2[1];
  87.         }
  88.         else
  89.         {
  90.             m=w1[n]+1;
  91.             if (sizeof(long)==2*sizeof(small))
  92.             { /* use longs if they are double length */
  93.                 if (n>2)
  94.                 { /* squeeze out as much significance as possible */
  95.                     u=muldiv(w1[n],base,w1[n-1],m,&r);
  96.                     u=u*base+muldiv(r,base,w1[n-2],m,&r);
  97.                     v=muldiv(w2[n],base,w2[n-1],m,&r);
  98.                     v=v*base+muldiv(r,base,w2[n-2],m,&r);
  99.                 }
  100.                 else
  101.                 {
  102.                     u=(long)base*w1[n]+w1[n-1];
  103.                     v=(long)base*w2[n]+w2[n-1];
  104.                     last=TRUE;
  105.                 }
  106.             }
  107.             else
  108.             {
  109.                 u=muldiv(w1[n],base,w1[n-1],m,&r);
  110.                 v=muldiv(w2[n],base,w2[n-1],m,&r);
  111.             }
  112.         }
  113.         forever
  114.         { /* work only with most significant piece */
  115.             if (last)
  116.             {
  117.                 if (v==0) break;
  118.                 lq=u/v;
  119.             }
  120.             else
  121.             {
  122.                 if (((v+c)==0) || ((v+d)==0)) break;
  123.                 lq=(u+a)/(v+c);
  124.                 if (lq!=(u+b)/(v+d)) break;
  125.             }
  126.             if (lq>=MAXBASE/abs(d)) break;
  127.             q=lq;
  128.             r=a-q*c;
  129.             a=c;
  130.             c=r;
  131.             r=b-q*d;
  132.             b=d;
  133.             d=r;
  134.             lr=u-lq*v;
  135.             u=v;
  136.             v=lr;
  137.         }
  138.     }
  139.     if (s==MINUS) negate(w3,w3);
  140.     if (exsign(w3)==MINUS) add(w3,y,w3);
  141.     if (xd!=yd)
  142.     {
  143.         negate(x,w2);
  144.         mad(w2,w3,w1,y,w4,w4);
  145.         copy(w4,yd);
  146.     }
  147.     copy(w3,xd);
  148.     if (z!=xd && z!=yd) copy(w1,z);
  149.     depth--;
  150.     return (size(w1));
  151. }
  152.  
  153.