plan9port

[fork] Plan 9 from user space
git clone git://src.adamsgaard.dk/plan9port # fast
git clone https://src.adamsgaard.dk/plan9port.git # slow
Log | Files | Refs | README | LICENSE Back to index

mpextendedgcd.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 }