tmpextendedgcd.c - plan9port - [fork] Plan 9 from user space
  HTML git clone git://src.adamsgaard.dk/plan9port
   DIR Log
   DIR Files
   DIR Refs
   DIR README
   DIR LICENSE
       ---
       tmpextendedgcd.c (1764B)
       ---
            1 #include "os.h"
            2 #include <mp.h>
            3 
            4 #define iseven(a)        (((a)->p[0] & 1) == 0)
            5 
            6 /* extended binary gcd */
            7 /* */
            8 /* For a anv b it solves, v = gcd(a,b) and finds x and y s.t. */
            9 /* ax + by = v */
           10 /* */
           11 /* Handbook of Applied Cryptography, Menezes et al, 1997, pg 608.   */
           12 void
           13 mpextendedgcd(mpint *a, mpint *b, mpint *v, mpint *x, mpint *y)
           14 {
           15         mpint *u, *A, *B, *C, *D;
           16         int g;
           17 
           18         if(a->sign < 0 || b->sign < 0){
           19                 mpassign(mpzero, v);
           20                 mpassign(mpzero, y);
           21                 mpassign(mpzero, x);
           22                 return;
           23         }
           24 
           25         if(a->top == 0){
           26                 mpassign(b, v);
           27                 mpassign(mpone, y);
           28                 mpassign(mpzero, x);
           29                 return;
           30         }
           31         if(b->top == 0){
           32                 mpassign(a, v);
           33                 mpassign(mpone, x);
           34                 mpassign(mpzero, y);
           35                 return;
           36         }
           37 
           38         g = 0;
           39         a = mpcopy(a);
           40         b = mpcopy(b);
           41 
           42         while(iseven(a) && iseven(b)){
           43                 mpright(a, 1, a);
           44                 mpright(b, 1, b);
           45                 g++;
           46         }
           47 
           48         u = mpcopy(a);
           49         mpassign(b, v);
           50         A = mpcopy(mpone);
           51         B = mpcopy(mpzero);
           52         C = mpcopy(mpzero);
           53         D = mpcopy(mpone);
           54 
           55         for(;;) {
           56 /*                print("%B %B %B %B %B %B\n", u, v, A, B, C, D); */
           57                 while(iseven(u)){
           58                         mpright(u, 1, u);
           59                         if(!iseven(A) || !iseven(B)) {
           60                                 mpadd(A, b, A);
           61                                 mpsub(B, a, B);
           62                         }
           63                         mpright(A, 1, A);
           64                         mpright(B, 1, B);
           65                 }
           66 
           67 /*                print("%B %B %B %B %B %B\n", u, v, A, B, C, D); */
           68                 while(iseven(v)){
           69                         mpright(v, 1, v);
           70                         if(!iseven(C) || !iseven(D)) {
           71                                 mpadd(C, b, C);
           72                                 mpsub(D, a, D);
           73                         }
           74                         mpright(C, 1, C);
           75                         mpright(D, 1, D);
           76                 }
           77 
           78 /*                print("%B %B %B %B %B %B\n", u, v, A, B, C, D); */
           79                 if(mpcmp(u, v) >= 0){
           80                         mpsub(u, v, u);
           81                         mpsub(A, C, A);
           82                         mpsub(B, D, B);
           83                 } else {
           84                         mpsub(v, u, v);
           85                         mpsub(C, A, C);
           86                         mpsub(D, B, D);
           87                 }
           88 
           89                 if(u->top == 0)
           90                         break;
           91 
           92         }
           93         mpassign(C, x);
           94         mpassign(D, y);
           95         mpleft(v, g, v);
           96 
           97         mpfree(A);
           98         mpfree(B);
           99         mpfree(C);
          100         mpfree(D);
          101         mpfree(u);
          102         mpfree(a);
          103         mpfree(b);
          104 
          105         return;
          106 }