summaryrefslogtreecommitdiff
path: root/usr/src/lib/libmp/common/gcd.c
blob: 18b6b58bead8cef4066ae556f8f5369d41bc3f5c (plain)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
/*	Copyright (c) 1984, 1986, 1987, 1988, 1989 AT&T	*/
/*	  All Rights Reserved  	*/


/*
 * Copyright (c) 1980 Regents of the University of California.
 * All rights reserved.  The Berkeley software License Agreement
 * specifies the terms and conditions for redistribution.
 */
/* 	Portions Copyright(c) 1988, Sun Microsystems Inc.	*/
/*	All Rights Reserved					*/

/*
 * Copyright (c) 1997, by Sun Microsystems, Inc.
 * All rights reserved.
 */

#ident	"%Z%%M%	%I%	%E% SMI"	/* SVr4.0 1.1	*/

/* LINTLIBRARY */

#include <mp.h>
#include "libmp.h"
#include <sys/types.h>

void
mp_gcd(MINT *a, MINT *b, MINT *c)
{
	MINT x, y, z, w;

	x.len = y.len = z.len = w.len = 0;
	_mp_move(a, &x);
	_mp_move(b, &y);
	while (y.len != 0) {
		mp_mdiv(&x, &y, &w, &z);
		_mp_move(&y, &x);
		_mp_move(&z, &y);
	}
	_mp_move(&x, c);
	_mp_xfree(&x);
	_mp_xfree(&y);
	_mp_xfree(&z);
	_mp_xfree(&w);
}

void
mp_invert(MINT *x1, MINT *x0, MINT *c)
{
	MINT u2, u3;
	MINT v2, v3;
	MINT zero;
	MINT q, r;
	MINT t;
	MINT x0_prime;
	static MINT *one = NULL;

	/*
	 * Minimize calls to allocators.  Don't use pointers for local
	 * variables, for the one "initialized" multiple precision
	 * variable, do it just once.
	 */
	if (one == NULL)
		one = mp_itom(1);

	zero.len = q.len = r.len = t.len = 0;

	x0_prime.len = u2.len = u3.len = 0;
	_mp_move(x0, &u3);
	_mp_move(x0, &x0_prime);

	v2.len = v3.len = 0;
	_mp_move(one, &v2);
	_mp_move(x1, &v3);

	while (mp_mcmp(&v3, &zero) != 0) {
		/* invariant: x0*u1 + x1*u2 = u3 */
		/* invariant: x0*v1 + x2*v2 = v3 */
		/* invariant: x(n+1) = x(n-1) % x(n) */
		mp_mdiv(&u3, &v3, &q, &r);
		_mp_move(&v3, &u3);
		_mp_move(&r, &v3);

		mp_mult(&q, &v2, &t);
		mp_msub(&u2, &t, &t);
		_mp_move(&v2, &u2);
		_mp_move(&t, &v2);
	}
	/* now x0*u1 + x1*u2 == 1, therefore,  (u2*x1) % x0  == 1 */
	_mp_move(&u2, c);
	if (mp_mcmp(c, &zero) < 0) {
		mp_madd(&x0_prime, c, c);
	}
	_mp_xfree(&zero);
	_mp_xfree(&v2);
	_mp_xfree(&v3);
	_mp_xfree(&u2);
	_mp_xfree(&u3);
	_mp_xfree(&q);
	_mp_xfree(&r);
	_mp_xfree(&t);
}