xref: /illumos-gate/usr/src/lib/libmp/common/gcd.c (revision 598f4ceed9327d2d6c2325dd67cae3aa06f7fea6)
1 /*	Copyright (c) 1984, 1986, 1987, 1988, 1989 AT&T	*/
2 /*	  All Rights Reserved  	*/
3 
4 
5 /*
6  * Copyright (c) 1980 Regents of the University of California.
7  * All rights reserved.  The Berkeley software License Agreement
8  * specifies the terms and conditions for redistribution.
9  */
10 /* 	Portions Copyright(c) 1988, Sun Microsystems Inc.	*/
11 /*	All Rights Reserved					*/
12 
13 /*
14  * Copyright (c) 1997, by Sun Microsystems, Inc.
15  * All rights reserved.
16  */
17 
18 #ident	"%Z%%M%	%I%	%E% SMI"	/* SVr4.0 1.1	*/
19 
20 /* LINTLIBRARY */
21 
22 #include <mp.h>
23 #include "libmp.h"
24 #include <sys/types.h>
25 
26 void
27 mp_gcd(MINT *a, MINT *b, MINT *c)
28 {
29 	MINT x, y, z, w;
30 
31 	x.len = y.len = z.len = w.len = 0;
32 	_mp_move(a, &x);
33 	_mp_move(b, &y);
34 	while (y.len != 0) {
35 		mp_mdiv(&x, &y, &w, &z);
36 		_mp_move(&y, &x);
37 		_mp_move(&z, &y);
38 	}
39 	_mp_move(&x, c);
40 	_mp_xfree(&x);
41 	_mp_xfree(&y);
42 	_mp_xfree(&z);
43 	_mp_xfree(&w);
44 }
45 
46 void
47 mp_invert(MINT *x1, MINT *x0, MINT *c)
48 {
49 	MINT u2, u3;
50 	MINT v2, v3;
51 	MINT zero;
52 	MINT q, r;
53 	MINT t;
54 	MINT x0_prime;
55 	static MINT *one = NULL;
56 
57 	/*
58 	 * Minimize calls to allocators.  Don't use pointers for local
59 	 * variables, for the one "initialized" multiple precision
60 	 * variable, do it just once.
61 	 */
62 	if (one == NULL)
63 		one = mp_itom(1);
64 
65 	zero.len = q.len = r.len = t.len = 0;
66 
67 	x0_prime.len = u2.len = u3.len = 0;
68 	_mp_move(x0, &u3);
69 	_mp_move(x0, &x0_prime);
70 
71 	v2.len = v3.len = 0;
72 	_mp_move(one, &v2);
73 	_mp_move(x1, &v3);
74 
75 	while (mp_mcmp(&v3, &zero) != 0) {
76 		/* invariant: x0*u1 + x1*u2 = u3 */
77 		/* invariant: x0*v1 + x2*v2 = v3 */
78 		/* invariant: x(n+1) = x(n-1) % x(n) */
79 		mp_mdiv(&u3, &v3, &q, &r);
80 		_mp_move(&v3, &u3);
81 		_mp_move(&r, &v3);
82 
83 		mp_mult(&q, &v2, &t);
84 		mp_msub(&u2, &t, &t);
85 		_mp_move(&v2, &u2);
86 		_mp_move(&t, &v2);
87 	}
88 	/* now x0*u1 + x1*u2 == 1, therefore,  (u2*x1) % x0  == 1 */
89 	_mp_move(&u2, c);
90 	if (mp_mcmp(c, &zero) < 0) {
91 		mp_madd(&x0_prime, c, c);
92 	}
93 	_mp_xfree(&zero);
94 	_mp_xfree(&v2);
95 	_mp_xfree(&v3);
96 	_mp_xfree(&u2);
97 	_mp_xfree(&u3);
98 	_mp_xfree(&q);
99 	_mp_xfree(&r);
100 	_mp_xfree(&t);
101 }
102