1 #ifndef __MP_H__
2 #define __MP_H__ 1
3 #ifdef __cplusplus
4 extern "C" {
5 #endif
6 
7 /*
8 #pragma	src	"/sys/src/libmp"
9 #pragma	lib	"libmp.a"
10 */
11 
12 #define _MPINT 1
13 
14 typedef long mpdigit;
15 
16 // the code assumes mpdigit to be at least an int
17 // mpdigit must be an atomic type.  mpdigit is defined
18 // in the architecture specific u.h
19 
20 typedef struct mpint mpint;
21 
22 struct mpint
23 {
24 	int	sign;	// +1 or -1
25 	int	size;	// allocated digits
26 	int	top;	// significant digits
27 	mpdigit	*p;
28 	char	flags;
29 };
30 
31 enum
32 {
33 	MPstatic=	0x01,
34 	Dbytes=		sizeof(mpdigit),	// bytes per digit
35 	Dbits=		Dbytes*8		// bits per digit
36 };
37 
38 // allocation
39 void	mpsetminbits(int n);	// newly created mpint's get at least n bits
40 mpint*	mpnew(int n);		// create a new mpint with at least n bits
41 void	mpfree(mpint *b);
42 void	mpbits(mpint *b, int n);	// ensure that b has at least n bits
43 void	mpnorm(mpint *b);		// dump leading zeros
44 mpint*	mpcopy(mpint *b);
45 void	mpassign(mpint *old, mpint *new);
46 
47 // random bits
48 mpint*	mprand(int bits, void (*gen)(uchar*, int), mpint *b);
49 
50 // conversion
51 mpint*	strtomp(char*, char**, int, mpint*);	// ascii
52 int	mpfmt(Fmt*);
53 char*	mptoa(mpint*, int, char*, int);
54 mpint*	letomp(uchar*, uint, mpint*);	// byte array, little-endian
55 int	mptole(mpint*, uchar*, uint, uchar**);
56 mpint*	betomp(uchar*, uint, mpint*);	// byte array, little-endian
57 int	mptobe(mpint*, uchar*, uint, uchar**);
58 uint	mptoui(mpint*);			// unsigned int
59 mpint*	uitomp(uint, mpint*);
60 int	mptoi(mpint*);			// int
61 mpint*	itomp(int, mpint*);
62 uvlong	mptouv(mpint*);			// unsigned vlong
63 mpint*	uvtomp(uvlong, mpint*);
64 vlong	mptov(mpint*);			// vlong
65 mpint*	vtomp(vlong, mpint*);
66 
67 // divide 2 digits by one
68 void	mpdigdiv(mpdigit *dividend, mpdigit divisor, mpdigit *quotient);
69 
70 // in the following, the result mpint may be
71 // the same as one of the inputs.
72 void	mpadd(mpint *b1, mpint *b2, mpint *sum);	// sum = b1+b2
73 void	mpsub(mpint *b1, mpint *b2, mpint *diff);	// diff = b1-b2
74 void	mpleft(mpint *b, int shift, mpint *res);	// res = b<<shift
75 void	mpright(mpint *b, int shift, mpint *res);	// res = b>>shift
76 void	mpmul(mpint *b1, mpint *b2, mpint *prod);	// prod = b1*b2
77 void	mpexp(mpint *b, mpint *e, mpint *m, mpint *res);	// res = b**e mod m
78 void	mpmod(mpint *b, mpint *m, mpint *remainder);	// remainder = b mod m
79 
80 // quotient = dividend/divisor, remainder = dividend % divisor
81 void	mpdiv(mpint *dividend, mpint *divisor,  mpint *quotient, mpint *remainder);
82 
83 // return neg, 0, pos as b1-b2 is neg, 0, pos
84 int	mpcmp(mpint *b1, mpint *b2);
85 
86 // extended gcd return d, x, and y, s.t. d = gcd(a,b) and ax+by = d
87 void	mpextendedgcd(mpint *a, mpint *b, mpint *d, mpint *x, mpint *y);
88 
89 // res = b**-1 mod m
90 void	mpinvert(mpint *b, mpint *m, mpint *res);
91 
92 // bit counting
93 int	mpsignif(mpint*);	// number of sigificant bits in mantissa
94 int	mplowbits0(mpint*);	// k, where n = 2**k * q for odd q
95 
96 // well known constants
97 extern mpint	*mpzero, *mpone, *mptwo;
98 
99 // sum[0:alen] = a[0:alen-1] + b[0:blen-1]
100 // prereq: alen >= blen, sum has room for alen+1 digits
101 void	mpvecadd(mpdigit *a, int alen, mpdigit *b, int blen, mpdigit *sum);
102 
103 // diff[0:alen-1] = a[0:alen-1] - b[0:blen-1]
104 // prereq: alen >= blen, diff has room for alen digits
105 void	mpvecsub(mpdigit *a, int alen, mpdigit *b, int blen, mpdigit *diff);
106 
107 // p[0:n] += m * b[0:n-1]
108 // prereq: p has room for n+1 digits
109 void	mpvecdigmuladd(mpdigit *b, int n, mpdigit m, mpdigit *p);
110 
111 // p[0:n] -= m * b[0:n-1]
112 // prereq: p has room for n+1 digits
113 int	mpvecdigmulsub(mpdigit *b, int n, mpdigit m, mpdigit *p);
114 
115 // p[0:alen*blen-1] = a[0:alen-1] * b[0:blen-1]
116 // prereq: alen >= blen, p has room for m*n digits
117 void	mpvecmul(mpdigit *a, int alen, mpdigit *b, int blen, mpdigit *p);
118 
119 // sign of a - b or zero if the same
120 int	mpveccmp(mpdigit *a, int alen, mpdigit *b, int blen);
121 
122 // divide the 2 digit dividend by the one digit divisor and stick in quotient
123 // we assume that the result is one digit - overflow is all 1's
124 void	mpdigdiv(mpdigit *dividend, mpdigit divisor, mpdigit *quotient);
125 
126 // playing with magnitudes
127 int	mpmagcmp(mpint *b1, mpint *b2);
128 void	mpmagadd(mpint *b1, mpint *b2, mpint *sum);	// sum = b1+b2
129 void	mpmagsub(mpint *b1, mpint *b2, mpint *sum);	// sum = b1+b2
130 
131 // chinese remainder theorem
132 typedef struct CRTpre	CRTpre;		// precomputed values for converting
133 					//  twixt residues and mpint
134 typedef struct CRTres	CRTres;		// residue form of an mpint
135 
136 struct CRTres
137 {
138 	int	n;		// number of residues
139 	mpint	*r[1];		// residues
140 };
141 
142 CRTpre*	crtpre(int, mpint**);			// precompute conversion values
143 CRTres*	crtin(CRTpre*, mpint*);			// convert mpint to residues
144 void	crtout(CRTpre*, CRTres*, mpint*);	// convert residues to mpint
145 void	crtprefree(CRTpre*);
146 void	crtresfree(CRTres*);
147 
148 
149 /* #pragma	varargck	type	"B"	mpint* */
150 #ifdef __cplusplus
151 }
152 #endif
153 #endif