tmpeuclid.c - plan9port - [fork] Plan 9 from user space
 (HTM) git clone git://src.adamsgaard.dk/plan9port
 (DIR) Log
 (DIR) Files
 (DIR) Refs
 (DIR) README
 (DIR) LICENSE
       ---
       tmpeuclid.c (1326B)
       ---
            1 #include "os.h"
            2 #include <mp.h>
            3 
            4 /* extended euclid */
            5 /* */
            6 /* For a and b it solves, d = gcd(a,b) and finds x and y s.t. */
            7 /* ax + by = d */
            8 /* */
            9 /* Handbook of Applied Cryptography, Menezes et al, 1997, pg 67 */
           10 
           11 void
           12 mpeuclid(mpint *a, mpint *b, mpint *d, mpint *x, mpint *y)
           13 {
           14         mpint *tmp, *x0, *x1, *x2, *y0, *y1, *y2, *q, *r;
           15 
           16         if(a->sign<0 || b->sign<0)
           17                 sysfatal("mpeuclid: negative arg");
           18 
           19         if(mpcmp(a, b) < 0){
           20                 tmp = a;
           21                 a = b;
           22                 b = tmp;
           23                 tmp = x;
           24                 x = y;
           25                 y = tmp;
           26         }
           27 
           28         if(b->top == 0){
           29                 mpassign(a, d);
           30                 mpassign(mpone, x);
           31                 mpassign(mpzero, y);
           32                 return;
           33         }
           34 
           35         a = mpcopy(a);
           36         b = mpcopy(b);
           37         x0 = mpnew(0);
           38         x1 = mpcopy(mpzero);
           39         x2 = mpcopy(mpone);
           40         y0 = mpnew(0);
           41         y1 = mpcopy(mpone);
           42         y2 = mpcopy(mpzero);
           43         q = mpnew(0);
           44         r = mpnew(0);
           45 
           46         while(b->top != 0 && b->sign > 0){
           47                 /* q = a/b */
           48                 /* r = a mod b */
           49                 mpdiv(a, b, q, r);
           50                 /* x0 = x2 - qx1 */
           51                 mpmul(q, x1, x0);
           52                 mpsub(x2, x0, x0);
           53                 /* y0 = y2 - qy1 */
           54                 mpmul(q, y1, y0);
           55                 mpsub(y2, y0, y0);
           56                 /* rotate values */
           57                 tmp = a;
           58                 a = b;
           59                 b = r;
           60                 r = tmp;
           61                 tmp = x2;
           62                 x2 = x1;
           63                 x1 = x0;
           64                 x0 = tmp;
           65                 tmp = y2;
           66                 y2 = y1;
           67                 y1 = y0;
           68                 y0 = tmp;
           69         }
           70 
           71         mpassign(a, d);
           72         mpassign(x2, x);
           73         mpassign(y2, y);
           74 
           75         mpfree(x0);
           76         mpfree(x1);
           77         mpfree(x2);
           78         mpfree(y0);
           79         mpfree(y1);
           80         mpfree(y2);
           81         mpfree(q);
           82         mpfree(r);
           83         mpfree(a);
           84         mpfree(b);
           85 }