Skip to content

Commit d548c21

Browse files
committed
Bring in libtommath bigint library. Build it and link it all into an nqp_bigint.ops file (which is just a stub at the moment).
1 parent 7bc5fcd commit d548c21

File tree

126 files changed

+11277
-1
lines changed

Some content is hidden

Large Commits have some content hidden by default. Use the searchbox below for content that may be hidden.

126 files changed

+11277
-1
lines changed

3rdparty/libtommath/LICENSE

Lines changed: 29 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,29 @@
1+
LibTomMath is licensed under DUAL licensing terms.
2+
3+
Choose and use the license of your needs.
4+
5+
[LICENSE #1]
6+
7+
LibTomMath is public domain. As should all quality software be.
8+
9+
Tom St Denis
10+
11+
[/LICENSE #1]
12+
13+
[LICENSE #2]
14+
15+
DO WHAT THE FUCK YOU WANT TO PUBLIC LICENSE
16+
Version 2, December 2004
17+
18+
Copyright (C) 2004 Sam Hocevar <sam@hocevar.net>
19+
20+
Everyone is permitted to copy and distribute verbatim or modified
21+
copies of this license document, and changing it is allowed as long
22+
as the name is changed.
23+
24+
DO WHAT THE FUCK YOU WANT TO PUBLIC LICENSE
25+
TERMS AND CONDITIONS FOR COPYING, DISTRIBUTION AND MODIFICATION
26+
27+
0. You just DO WHAT THE FUCK YOU WANT TO.
28+
29+
[/LICENSE #2]

3rdparty/libtommath/bn_error.c

Lines changed: 47 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,47 @@
1+
#include <tommath.h>
2+
#ifdef BN_ERROR_C
3+
/* LibTomMath, multiple-precision integer library -- Tom St Denis
4+
*
5+
* LibTomMath is a library that provides multiple-precision
6+
* integer arithmetic as well as number theoretic functionality.
7+
*
8+
* The library was designed directly after the MPI library by
9+
* Michael Fromberger but has been written from scratch with
10+
* additional optimizations in place.
11+
*
12+
* The library is free for all purposes without any express
13+
* guarantee it works.
14+
*
15+
* Tom St Denis, tomstdenis@gmail.com, http://libtom.org
16+
*/
17+
18+
static const struct {
19+
int code;
20+
const char *msg;
21+
} msgs[] = {
22+
{ MP_OKAY, "Successful" },
23+
{ MP_MEM, "Out of heap" },
24+
{ MP_VAL, "Value out of range" }
25+
};
26+
27+
/* return a char * string for a given code */
28+
char *mp_error_to_string(int code)
29+
{
30+
int x;
31+
32+
/* scan the lookup table for the given message */
33+
for (x = 0; x < (int)(sizeof(msgs) / sizeof(msgs[0])); x++) {
34+
if (msgs[x].code == code) {
35+
return msgs[x].msg;
36+
}
37+
}
38+
39+
/* generic reply for invalid code */
40+
return "Invalid error code";
41+
}
42+
43+
#endif
44+
45+
/* $Source$ */
46+
/* $Revision$ */
47+
/* $Date$ */
Lines changed: 148 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,148 @@
1+
#include <tommath.h>
2+
#ifdef BN_FAST_MP_INVMOD_C
3+
/* LibTomMath, multiple-precision integer library -- Tom St Denis
4+
*
5+
* LibTomMath is a library that provides multiple-precision
6+
* integer arithmetic as well as number theoretic functionality.
7+
*
8+
* The library was designed directly after the MPI library by
9+
* Michael Fromberger but has been written from scratch with
10+
* additional optimizations in place.
11+
*
12+
* The library is free for all purposes without any express
13+
* guarantee it works.
14+
*
15+
* Tom St Denis, tomstdenis@gmail.com, http://libtom.org
16+
*/
17+
18+
/* computes the modular inverse via binary extended euclidean algorithm,
19+
* that is c = 1/a mod b
20+
*
21+
* Based on slow invmod except this is optimized for the case where b is
22+
* odd as per HAC Note 14.64 on pp. 610
23+
*/
24+
int fast_mp_invmod (mp_int * a, mp_int * b, mp_int * c)
25+
{
26+
mp_int x, y, u, v, B, D;
27+
int res, neg;
28+
29+
/* 2. [modified] b must be odd */
30+
if (mp_iseven (b) == 1) {
31+
return MP_VAL;
32+
}
33+
34+
/* init all our temps */
35+
if ((res = mp_init_multi(&x, &y, &u, &v, &B, &D, NULL)) != MP_OKAY) {
36+
return res;
37+
}
38+
39+
/* x == modulus, y == value to invert */
40+
if ((res = mp_copy (b, &x)) != MP_OKAY) {
41+
goto LBL_ERR;
42+
}
43+
44+
/* we need y = |a| */
45+
if ((res = mp_mod (a, b, &y)) != MP_OKAY) {
46+
goto LBL_ERR;
47+
}
48+
49+
/* 3. u=x, v=y, A=1, B=0, C=0,D=1 */
50+
if ((res = mp_copy (&x, &u)) != MP_OKAY) {
51+
goto LBL_ERR;
52+
}
53+
if ((res = mp_copy (&y, &v)) != MP_OKAY) {
54+
goto LBL_ERR;
55+
}
56+
mp_set (&D, 1);
57+
58+
top:
59+
/* 4. while u is even do */
60+
while (mp_iseven (&u) == 1) {
61+
/* 4.1 u = u/2 */
62+
if ((res = mp_div_2 (&u, &u)) != MP_OKAY) {
63+
goto LBL_ERR;
64+
}
65+
/* 4.2 if B is odd then */
66+
if (mp_isodd (&B) == 1) {
67+
if ((res = mp_sub (&B, &x, &B)) != MP_OKAY) {
68+
goto LBL_ERR;
69+
}
70+
}
71+
/* B = B/2 */
72+
if ((res = mp_div_2 (&B, &B)) != MP_OKAY) {
73+
goto LBL_ERR;
74+
}
75+
}
76+
77+
/* 5. while v is even do */
78+
while (mp_iseven (&v) == 1) {
79+
/* 5.1 v = v/2 */
80+
if ((res = mp_div_2 (&v, &v)) != MP_OKAY) {
81+
goto LBL_ERR;
82+
}
83+
/* 5.2 if D is odd then */
84+
if (mp_isodd (&D) == 1) {
85+
/* D = (D-x)/2 */
86+
if ((res = mp_sub (&D, &x, &D)) != MP_OKAY) {
87+
goto LBL_ERR;
88+
}
89+
}
90+
/* D = D/2 */
91+
if ((res = mp_div_2 (&D, &D)) != MP_OKAY) {
92+
goto LBL_ERR;
93+
}
94+
}
95+
96+
/* 6. if u >= v then */
97+
if (mp_cmp (&u, &v) != MP_LT) {
98+
/* u = u - v, B = B - D */
99+
if ((res = mp_sub (&u, &v, &u)) != MP_OKAY) {
100+
goto LBL_ERR;
101+
}
102+
103+
if ((res = mp_sub (&B, &D, &B)) != MP_OKAY) {
104+
goto LBL_ERR;
105+
}
106+
} else {
107+
/* v - v - u, D = D - B */
108+
if ((res = mp_sub (&v, &u, &v)) != MP_OKAY) {
109+
goto LBL_ERR;
110+
}
111+
112+
if ((res = mp_sub (&D, &B, &D)) != MP_OKAY) {
113+
goto LBL_ERR;
114+
}
115+
}
116+
117+
/* if not zero goto step 4 */
118+
if (mp_iszero (&u) == 0) {
119+
goto top;
120+
}
121+
122+
/* now a = C, b = D, gcd == g*v */
123+
124+
/* if v != 1 then there is no inverse */
125+
if (mp_cmp_d (&v, 1) != MP_EQ) {
126+
res = MP_VAL;
127+
goto LBL_ERR;
128+
}
129+
130+
/* b is now the inverse */
131+
neg = a->sign;
132+
while (D.sign == MP_NEG) {
133+
if ((res = mp_add (&D, b, &D)) != MP_OKAY) {
134+
goto LBL_ERR;
135+
}
136+
}
137+
mp_exch (&D, c);
138+
c->sign = neg;
139+
res = MP_OKAY;
140+
141+
LBL_ERR:mp_clear_multi (&x, &y, &u, &v, &B, &D, NULL);
142+
return res;
143+
}
144+
#endif
145+
146+
/* $Source$ */
147+
/* $Revision$ */
148+
/* $Date$ */
Lines changed: 172 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,172 @@
1+
#include <tommath.h>
2+
#ifdef BN_FAST_MP_MONTGOMERY_REDUCE_C
3+
/* LibTomMath, multiple-precision integer library -- Tom St Denis
4+
*
5+
* LibTomMath is a library that provides multiple-precision
6+
* integer arithmetic as well as number theoretic functionality.
7+
*
8+
* The library was designed directly after the MPI library by
9+
* Michael Fromberger but has been written from scratch with
10+
* additional optimizations in place.
11+
*
12+
* The library is free for all purposes without any express
13+
* guarantee it works.
14+
*
15+
* Tom St Denis, tomstdenis@gmail.com, http://libtom.org
16+
*/
17+
18+
/* computes xR**-1 == x (mod N) via Montgomery Reduction
19+
*
20+
* This is an optimized implementation of montgomery_reduce
21+
* which uses the comba method to quickly calculate the columns of the
22+
* reduction.
23+
*
24+
* Based on Algorithm 14.32 on pp.601 of HAC.
25+
*/
26+
int fast_mp_montgomery_reduce (mp_int * x, mp_int * n, mp_digit rho)
27+
{
28+
int ix, res, olduse;
29+
mp_word W[MP_WARRAY];
30+
31+
/* get old used count */
32+
olduse = x->used;
33+
34+
/* grow a as required */
35+
if (x->alloc < n->used + 1) {
36+
if ((res = mp_grow (x, n->used + 1)) != MP_OKAY) {
37+
return res;
38+
}
39+
}
40+
41+
/* first we have to get the digits of the input into
42+
* an array of double precision words W[...]
43+
*/
44+
{
45+
register mp_word *_W;
46+
register mp_digit *tmpx;
47+
48+
/* alias for the W[] array */
49+
_W = W;
50+
51+
/* alias for the digits of x*/
52+
tmpx = x->dp;
53+
54+
/* copy the digits of a into W[0..a->used-1] */
55+
for (ix = 0; ix < x->used; ix++) {
56+
*_W++ = *tmpx++;
57+
}
58+
59+
/* zero the high words of W[a->used..m->used*2] */
60+
for (; ix < n->used * 2 + 1; ix++) {
61+
*_W++ = 0;
62+
}
63+
}
64+
65+
/* now we proceed to zero successive digits
66+
* from the least significant upwards
67+
*/
68+
for (ix = 0; ix < n->used; ix++) {
69+
/* mu = ai * m' mod b
70+
*
71+
* We avoid a double precision multiplication (which isn't required)
72+
* by casting the value down to a mp_digit. Note this requires
73+
* that W[ix-1] have the carry cleared (see after the inner loop)
74+
*/
75+
register mp_digit mu;
76+
mu = (mp_digit) (((W[ix] & MP_MASK) * rho) & MP_MASK);
77+
78+
/* a = a + mu * m * b**i
79+
*
80+
* This is computed in place and on the fly. The multiplication
81+
* by b**i is handled by offseting which columns the results
82+
* are added to.
83+
*
84+
* Note the comba method normally doesn't handle carries in the
85+
* inner loop In this case we fix the carry from the previous
86+
* column since the Montgomery reduction requires digits of the
87+
* result (so far) [see above] to work. This is
88+
* handled by fixing up one carry after the inner loop. The
89+
* carry fixups are done in order so after these loops the
90+
* first m->used words of W[] have the carries fixed
91+
*/
92+
{
93+
register int iy;
94+
register mp_digit *tmpn;
95+
register mp_word *_W;
96+
97+
/* alias for the digits of the modulus */
98+
tmpn = n->dp;
99+
100+
/* Alias for the columns set by an offset of ix */
101+
_W = W + ix;
102+
103+
/* inner loop */
104+
for (iy = 0; iy < n->used; iy++) {
105+
*_W++ += ((mp_word)mu) * ((mp_word)*tmpn++);
106+
}
107+
}
108+
109+
/* now fix carry for next digit, W[ix+1] */
110+
W[ix + 1] += W[ix] >> ((mp_word) DIGIT_BIT);
111+
}
112+
113+
/* now we have to propagate the carries and
114+
* shift the words downward [all those least
115+
* significant digits we zeroed].
116+
*/
117+
{
118+
register mp_digit *tmpx;
119+
register mp_word *_W, *_W1;
120+
121+
/* nox fix rest of carries */
122+
123+
/* alias for current word */
124+
_W1 = W + ix;
125+
126+
/* alias for next word, where the carry goes */
127+
_W = W + ++ix;
128+
129+
for (; ix <= n->used * 2 + 1; ix++) {
130+
*_W++ += *_W1++ >> ((mp_word) DIGIT_BIT);
131+
}
132+
133+
/* copy out, A = A/b**n
134+
*
135+
* The result is A/b**n but instead of converting from an
136+
* array of mp_word to mp_digit than calling mp_rshd
137+
* we just copy them in the right order
138+
*/
139+
140+
/* alias for destination word */
141+
tmpx = x->dp;
142+
143+
/* alias for shifted double precision result */
144+
_W = W + n->used;
145+
146+
for (ix = 0; ix < n->used + 1; ix++) {
147+
*tmpx++ = (mp_digit)(*_W++ & ((mp_word) MP_MASK));
148+
}
149+
150+
/* zero oldused digits, if the input a was larger than
151+
* m->used+1 we'll have to clear the digits
152+
*/
153+
for (; ix < olduse; ix++) {
154+
*tmpx++ = 0;
155+
}
156+
}
157+
158+
/* set the max used and clamp */
159+
x->used = n->used + 1;
160+
mp_clamp (x);
161+
162+
/* if A >= m then A = A - m */
163+
if (mp_cmp_mag (x, n) != MP_LT) {
164+
return s_mp_sub (x, n, x);
165+
}
166+
return MP_OKAY;
167+
}
168+
#endif
169+
170+
/* $Source$ */
171+
/* $Revision$ */
172+
/* $Date$ */

0 commit comments

Comments
 (0)