1*f9fbec18Smcpowers /*
2*f9fbec18Smcpowers * ***** BEGIN LICENSE BLOCK *****
3*f9fbec18Smcpowers * Version: MPL 1.1/GPL 2.0/LGPL 2.1
4*f9fbec18Smcpowers *
5*f9fbec18Smcpowers * The contents of this file are subject to the Mozilla Public License Version
6*f9fbec18Smcpowers * 1.1 (the "License"); you may not use this file except in compliance with
7*f9fbec18Smcpowers * the License. You may obtain a copy of the License at
8*f9fbec18Smcpowers * http://www.mozilla.org/MPL/
9*f9fbec18Smcpowers *
10*f9fbec18Smcpowers * Software distributed under the License is distributed on an "AS IS" basis,
11*f9fbec18Smcpowers * WITHOUT WARRANTY OF ANY KIND, either express or implied. See the License
12*f9fbec18Smcpowers * for the specific language governing rights and limitations under the
13*f9fbec18Smcpowers * License.
14*f9fbec18Smcpowers *
15*f9fbec18Smcpowers * The Original Code is the Multi-precision Binary Polynomial Arithmetic Library.
16*f9fbec18Smcpowers *
17*f9fbec18Smcpowers * The Initial Developer of the Original Code is
18*f9fbec18Smcpowers * Sun Microsystems, Inc.
19*f9fbec18Smcpowers * Portions created by the Initial Developer are Copyright (C) 2003
20*f9fbec18Smcpowers * the Initial Developer. All Rights Reserved.
21*f9fbec18Smcpowers *
22*f9fbec18Smcpowers * Contributor(s):
23*f9fbec18Smcpowers * Sheueling Chang Shantz <sheueling.chang@sun.com> and
24*f9fbec18Smcpowers * Douglas Stebila <douglas@stebila.ca> of Sun Laboratories.
25*f9fbec18Smcpowers *
26*f9fbec18Smcpowers * Alternatively, the contents of this file may be used under the terms of
27*f9fbec18Smcpowers * either the GNU General Public License Version 2 or later (the "GPL"), or
28*f9fbec18Smcpowers * the GNU Lesser General Public License Version 2.1 or later (the "LGPL"),
29*f9fbec18Smcpowers * in which case the provisions of the GPL or the LGPL are applicable instead
30*f9fbec18Smcpowers * of those above. If you wish to allow use of your version of this file only
31*f9fbec18Smcpowers * under the terms of either the GPL or the LGPL, and not to allow others to
32*f9fbec18Smcpowers * use your version of this file under the terms of the MPL, indicate your
33*f9fbec18Smcpowers * decision by deleting the provisions above and replace them with the notice
34*f9fbec18Smcpowers * and other provisions required by the GPL or the LGPL. If you do not delete
35*f9fbec18Smcpowers * the provisions above, a recipient may use your version of this file under
36*f9fbec18Smcpowers * the terms of any one of the MPL, the GPL or the LGPL.
37*f9fbec18Smcpowers *
38*f9fbec18Smcpowers * ***** END LICENSE BLOCK ***** */
39*f9fbec18Smcpowers /*
40*f9fbec18Smcpowers * Copyright 2007 Sun Microsystems, Inc. All rights reserved.
41*f9fbec18Smcpowers * Use is subject to license terms.
42*f9fbec18Smcpowers *
43*f9fbec18Smcpowers * Sun elects to use this software under the MPL license.
44*f9fbec18Smcpowers */
45*f9fbec18Smcpowers
46*f9fbec18Smcpowers #include "mp_gf2m.h"
47*f9fbec18Smcpowers #include "mp_gf2m-priv.h"
48*f9fbec18Smcpowers #include "mplogic.h"
49*f9fbec18Smcpowers #include "mpi-priv.h"
50*f9fbec18Smcpowers
51*f9fbec18Smcpowers const mp_digit mp_gf2m_sqr_tb[16] =
52*f9fbec18Smcpowers {
53*f9fbec18Smcpowers 0, 1, 4, 5, 16, 17, 20, 21,
54*f9fbec18Smcpowers 64, 65, 68, 69, 80, 81, 84, 85
55*f9fbec18Smcpowers };
56*f9fbec18Smcpowers
57*f9fbec18Smcpowers /* Multiply two binary polynomials mp_digits a, b.
58*f9fbec18Smcpowers * Result is a polynomial with degree < 2 * MP_DIGIT_BITS - 1.
59*f9fbec18Smcpowers * Output in two mp_digits rh, rl.
60*f9fbec18Smcpowers */
61*f9fbec18Smcpowers #if MP_DIGIT_BITS == 32
62*f9fbec18Smcpowers void
s_bmul_1x1(mp_digit * rh,mp_digit * rl,const mp_digit a,const mp_digit b)63*f9fbec18Smcpowers s_bmul_1x1(mp_digit *rh, mp_digit *rl, const mp_digit a, const mp_digit b)
64*f9fbec18Smcpowers {
65*f9fbec18Smcpowers register mp_digit h, l, s;
66*f9fbec18Smcpowers mp_digit tab[8], top2b = a >> 30;
67*f9fbec18Smcpowers register mp_digit a1, a2, a4;
68*f9fbec18Smcpowers
69*f9fbec18Smcpowers a1 = a & (0x3FFFFFFF); a2 = a1 << 1; a4 = a2 << 1;
70*f9fbec18Smcpowers
71*f9fbec18Smcpowers tab[0] = 0; tab[1] = a1; tab[2] = a2; tab[3] = a1^a2;
72*f9fbec18Smcpowers tab[4] = a4; tab[5] = a1^a4; tab[6] = a2^a4; tab[7] = a1^a2^a4;
73*f9fbec18Smcpowers
74*f9fbec18Smcpowers s = tab[b & 0x7]; l = s;
75*f9fbec18Smcpowers s = tab[b >> 3 & 0x7]; l ^= s << 3; h = s >> 29;
76*f9fbec18Smcpowers s = tab[b >> 6 & 0x7]; l ^= s << 6; h ^= s >> 26;
77*f9fbec18Smcpowers s = tab[b >> 9 & 0x7]; l ^= s << 9; h ^= s >> 23;
78*f9fbec18Smcpowers s = tab[b >> 12 & 0x7]; l ^= s << 12; h ^= s >> 20;
79*f9fbec18Smcpowers s = tab[b >> 15 & 0x7]; l ^= s << 15; h ^= s >> 17;
80*f9fbec18Smcpowers s = tab[b >> 18 & 0x7]; l ^= s << 18; h ^= s >> 14;
81*f9fbec18Smcpowers s = tab[b >> 21 & 0x7]; l ^= s << 21; h ^= s >> 11;
82*f9fbec18Smcpowers s = tab[b >> 24 & 0x7]; l ^= s << 24; h ^= s >> 8;
83*f9fbec18Smcpowers s = tab[b >> 27 & 0x7]; l ^= s << 27; h ^= s >> 5;
84*f9fbec18Smcpowers s = tab[b >> 30 ]; l ^= s << 30; h ^= s >> 2;
85*f9fbec18Smcpowers
86*f9fbec18Smcpowers /* compensate for the top two bits of a */
87*f9fbec18Smcpowers
88*f9fbec18Smcpowers if (top2b & 01) { l ^= b << 30; h ^= b >> 2; }
89*f9fbec18Smcpowers if (top2b & 02) { l ^= b << 31; h ^= b >> 1; }
90*f9fbec18Smcpowers
91*f9fbec18Smcpowers *rh = h; *rl = l;
92*f9fbec18Smcpowers }
93*f9fbec18Smcpowers #else
94*f9fbec18Smcpowers void
s_bmul_1x1(mp_digit * rh,mp_digit * rl,const mp_digit a,const mp_digit b)95*f9fbec18Smcpowers s_bmul_1x1(mp_digit *rh, mp_digit *rl, const mp_digit a, const mp_digit b)
96*f9fbec18Smcpowers {
97*f9fbec18Smcpowers register mp_digit h, l, s;
98*f9fbec18Smcpowers mp_digit tab[16], top3b = a >> 61;
99*f9fbec18Smcpowers register mp_digit a1, a2, a4, a8;
100*f9fbec18Smcpowers
101*f9fbec18Smcpowers a1 = a & (0x1FFFFFFFFFFFFFFFULL); a2 = a1 << 1;
102*f9fbec18Smcpowers a4 = a2 << 1; a8 = a4 << 1;
103*f9fbec18Smcpowers tab[ 0] = 0; tab[ 1] = a1; tab[ 2] = a2; tab[ 3] = a1^a2;
104*f9fbec18Smcpowers tab[ 4] = a4; tab[ 5] = a1^a4; tab[ 6] = a2^a4; tab[ 7] = a1^a2^a4;
105*f9fbec18Smcpowers tab[ 8] = a8; tab[ 9] = a1^a8; tab[10] = a2^a8; tab[11] = a1^a2^a8;
106*f9fbec18Smcpowers tab[12] = a4^a8; tab[13] = a1^a4^a8; tab[14] = a2^a4^a8; tab[15] = a1^a2^a4^a8;
107*f9fbec18Smcpowers
108*f9fbec18Smcpowers s = tab[b & 0xF]; l = s;
109*f9fbec18Smcpowers s = tab[b >> 4 & 0xF]; l ^= s << 4; h = s >> 60;
110*f9fbec18Smcpowers s = tab[b >> 8 & 0xF]; l ^= s << 8; h ^= s >> 56;
111*f9fbec18Smcpowers s = tab[b >> 12 & 0xF]; l ^= s << 12; h ^= s >> 52;
112*f9fbec18Smcpowers s = tab[b >> 16 & 0xF]; l ^= s << 16; h ^= s >> 48;
113*f9fbec18Smcpowers s = tab[b >> 20 & 0xF]; l ^= s << 20; h ^= s >> 44;
114*f9fbec18Smcpowers s = tab[b >> 24 & 0xF]; l ^= s << 24; h ^= s >> 40;
115*f9fbec18Smcpowers s = tab[b >> 28 & 0xF]; l ^= s << 28; h ^= s >> 36;
116*f9fbec18Smcpowers s = tab[b >> 32 & 0xF]; l ^= s << 32; h ^= s >> 32;
117*f9fbec18Smcpowers s = tab[b >> 36 & 0xF]; l ^= s << 36; h ^= s >> 28;
118*f9fbec18Smcpowers s = tab[b >> 40 & 0xF]; l ^= s << 40; h ^= s >> 24;
119*f9fbec18Smcpowers s = tab[b >> 44 & 0xF]; l ^= s << 44; h ^= s >> 20;
120*f9fbec18Smcpowers s = tab[b >> 48 & 0xF]; l ^= s << 48; h ^= s >> 16;
121*f9fbec18Smcpowers s = tab[b >> 52 & 0xF]; l ^= s << 52; h ^= s >> 12;
122*f9fbec18Smcpowers s = tab[b >> 56 & 0xF]; l ^= s << 56; h ^= s >> 8;
123*f9fbec18Smcpowers s = tab[b >> 60 ]; l ^= s << 60; h ^= s >> 4;
124*f9fbec18Smcpowers
125*f9fbec18Smcpowers /* compensate for the top three bits of a */
126*f9fbec18Smcpowers
127*f9fbec18Smcpowers if (top3b & 01) { l ^= b << 61; h ^= b >> 3; }
128*f9fbec18Smcpowers if (top3b & 02) { l ^= b << 62; h ^= b >> 2; }
129*f9fbec18Smcpowers if (top3b & 04) { l ^= b << 63; h ^= b >> 1; }
130*f9fbec18Smcpowers
131*f9fbec18Smcpowers *rh = h; *rl = l;
132*f9fbec18Smcpowers }
133*f9fbec18Smcpowers #endif
134*f9fbec18Smcpowers
135*f9fbec18Smcpowers /* Compute xor-multiply of two binary polynomials (a1, a0) x (b1, b0)
136*f9fbec18Smcpowers * result is a binary polynomial in 4 mp_digits r[4].
137*f9fbec18Smcpowers * The caller MUST ensure that r has the right amount of space allocated.
138*f9fbec18Smcpowers */
139*f9fbec18Smcpowers void
s_bmul_2x2(mp_digit * r,const mp_digit a1,const mp_digit a0,const mp_digit b1,const mp_digit b0)140*f9fbec18Smcpowers s_bmul_2x2(mp_digit *r, const mp_digit a1, const mp_digit a0, const mp_digit b1,
141*f9fbec18Smcpowers const mp_digit b0)
142*f9fbec18Smcpowers {
143*f9fbec18Smcpowers mp_digit m1, m0;
144*f9fbec18Smcpowers /* r[3] = h1, r[2] = h0; r[1] = l1; r[0] = l0 */
145*f9fbec18Smcpowers s_bmul_1x1(r+3, r+2, a1, b1);
146*f9fbec18Smcpowers s_bmul_1x1(r+1, r, a0, b0);
147*f9fbec18Smcpowers s_bmul_1x1(&m1, &m0, a0 ^ a1, b0 ^ b1);
148*f9fbec18Smcpowers /* Correction on m1 ^= l1 ^ h1; m0 ^= l0 ^ h0; */
149*f9fbec18Smcpowers r[2] ^= m1 ^ r[1] ^ r[3]; /* h0 ^= m1 ^ l1 ^ h1; */
150*f9fbec18Smcpowers r[1] = r[3] ^ r[2] ^ r[0] ^ m1 ^ m0; /* l1 ^= l0 ^ h0 ^ m0; */
151*f9fbec18Smcpowers }
152*f9fbec18Smcpowers
153*f9fbec18Smcpowers /* Compute xor-multiply of two binary polynomials (a2, a1, a0) x (b2, b1, b0)
154*f9fbec18Smcpowers * result is a binary polynomial in 6 mp_digits r[6].
155*f9fbec18Smcpowers * The caller MUST ensure that r has the right amount of space allocated.
156*f9fbec18Smcpowers */
157*f9fbec18Smcpowers void
s_bmul_3x3(mp_digit * r,const mp_digit a2,const mp_digit a1,const mp_digit a0,const mp_digit b2,const mp_digit b1,const mp_digit b0)158*f9fbec18Smcpowers s_bmul_3x3(mp_digit *r, const mp_digit a2, const mp_digit a1, const mp_digit a0,
159*f9fbec18Smcpowers const mp_digit b2, const mp_digit b1, const mp_digit b0)
160*f9fbec18Smcpowers {
161*f9fbec18Smcpowers mp_digit zm[4];
162*f9fbec18Smcpowers
163*f9fbec18Smcpowers s_bmul_1x1(r+5, r+4, a2, b2); /* fill top 2 words */
164*f9fbec18Smcpowers s_bmul_2x2(zm, a1, a2^a0, b1, b2^b0); /* fill middle 4 words */
165*f9fbec18Smcpowers s_bmul_2x2(r, a1, a0, b1, b0); /* fill bottom 4 words */
166*f9fbec18Smcpowers
167*f9fbec18Smcpowers zm[3] ^= r[3];
168*f9fbec18Smcpowers zm[2] ^= r[2];
169*f9fbec18Smcpowers zm[1] ^= r[1] ^ r[5];
170*f9fbec18Smcpowers zm[0] ^= r[0] ^ r[4];
171*f9fbec18Smcpowers
172*f9fbec18Smcpowers r[5] ^= zm[3];
173*f9fbec18Smcpowers r[4] ^= zm[2];
174*f9fbec18Smcpowers r[3] ^= zm[1];
175*f9fbec18Smcpowers r[2] ^= zm[0];
176*f9fbec18Smcpowers }
177*f9fbec18Smcpowers
178*f9fbec18Smcpowers /* Compute xor-multiply of two binary polynomials (a3, a2, a1, a0) x (b3, b2, b1, b0)
179*f9fbec18Smcpowers * result is a binary polynomial in 8 mp_digits r[8].
180*f9fbec18Smcpowers * The caller MUST ensure that r has the right amount of space allocated.
181*f9fbec18Smcpowers */
s_bmul_4x4(mp_digit * r,const mp_digit a3,const mp_digit a2,const mp_digit a1,const mp_digit a0,const mp_digit b3,const mp_digit b2,const mp_digit b1,const mp_digit b0)182*f9fbec18Smcpowers void s_bmul_4x4(mp_digit *r, const mp_digit a3, const mp_digit a2, const mp_digit a1,
183*f9fbec18Smcpowers const mp_digit a0, const mp_digit b3, const mp_digit b2, const mp_digit b1,
184*f9fbec18Smcpowers const mp_digit b0)
185*f9fbec18Smcpowers {
186*f9fbec18Smcpowers mp_digit zm[4];
187*f9fbec18Smcpowers
188*f9fbec18Smcpowers s_bmul_2x2(r+4, a3, a2, b3, b2); /* fill top 4 words */
189*f9fbec18Smcpowers s_bmul_2x2(zm, a3^a1, a2^a0, b3^b1, b2^b0); /* fill middle 4 words */
190*f9fbec18Smcpowers s_bmul_2x2(r, a1, a0, b1, b0); /* fill bottom 4 words */
191*f9fbec18Smcpowers
192*f9fbec18Smcpowers zm[3] ^= r[3] ^ r[7];
193*f9fbec18Smcpowers zm[2] ^= r[2] ^ r[6];
194*f9fbec18Smcpowers zm[1] ^= r[1] ^ r[5];
195*f9fbec18Smcpowers zm[0] ^= r[0] ^ r[4];
196*f9fbec18Smcpowers
197*f9fbec18Smcpowers r[5] ^= zm[3];
198*f9fbec18Smcpowers r[4] ^= zm[2];
199*f9fbec18Smcpowers r[3] ^= zm[1];
200*f9fbec18Smcpowers r[2] ^= zm[0];
201*f9fbec18Smcpowers }
202*f9fbec18Smcpowers
203*f9fbec18Smcpowers /* Compute addition of two binary polynomials a and b,
204*f9fbec18Smcpowers * store result in c; c could be a or b, a and b could be equal;
205*f9fbec18Smcpowers * c is the bitwise XOR of a and b.
206*f9fbec18Smcpowers */
207*f9fbec18Smcpowers mp_err
mp_badd(const mp_int * a,const mp_int * b,mp_int * c)208*f9fbec18Smcpowers mp_badd(const mp_int *a, const mp_int *b, mp_int *c)
209*f9fbec18Smcpowers {
210*f9fbec18Smcpowers mp_digit *pa, *pb, *pc;
211*f9fbec18Smcpowers mp_size ix;
212*f9fbec18Smcpowers mp_size used_pa, used_pb;
213*f9fbec18Smcpowers mp_err res = MP_OKAY;
214*f9fbec18Smcpowers
215*f9fbec18Smcpowers /* Add all digits up to the precision of b. If b had more
216*f9fbec18Smcpowers * precision than a initially, swap a, b first
217*f9fbec18Smcpowers */
218*f9fbec18Smcpowers if (MP_USED(a) >= MP_USED(b)) {
219*f9fbec18Smcpowers pa = MP_DIGITS(a);
220*f9fbec18Smcpowers pb = MP_DIGITS(b);
221*f9fbec18Smcpowers used_pa = MP_USED(a);
222*f9fbec18Smcpowers used_pb = MP_USED(b);
223*f9fbec18Smcpowers } else {
224*f9fbec18Smcpowers pa = MP_DIGITS(b);
225*f9fbec18Smcpowers pb = MP_DIGITS(a);
226*f9fbec18Smcpowers used_pa = MP_USED(b);
227*f9fbec18Smcpowers used_pb = MP_USED(a);
228*f9fbec18Smcpowers }
229*f9fbec18Smcpowers
230*f9fbec18Smcpowers /* Make sure c has enough precision for the output value */
231*f9fbec18Smcpowers MP_CHECKOK( s_mp_pad(c, used_pa) );
232*f9fbec18Smcpowers
233*f9fbec18Smcpowers /* Do word-by-word xor */
234*f9fbec18Smcpowers pc = MP_DIGITS(c);
235*f9fbec18Smcpowers for (ix = 0; ix < used_pb; ix++) {
236*f9fbec18Smcpowers (*pc++) = (*pa++) ^ (*pb++);
237*f9fbec18Smcpowers }
238*f9fbec18Smcpowers
239*f9fbec18Smcpowers /* Finish the rest of digits until we're actually done */
240*f9fbec18Smcpowers for (; ix < used_pa; ++ix) {
241*f9fbec18Smcpowers *pc++ = *pa++;
242*f9fbec18Smcpowers }
243*f9fbec18Smcpowers
244*f9fbec18Smcpowers MP_USED(c) = used_pa;
245*f9fbec18Smcpowers MP_SIGN(c) = ZPOS;
246*f9fbec18Smcpowers s_mp_clamp(c);
247*f9fbec18Smcpowers
248*f9fbec18Smcpowers CLEANUP:
249*f9fbec18Smcpowers return res;
250*f9fbec18Smcpowers }
251*f9fbec18Smcpowers
252*f9fbec18Smcpowers #define s_mp_div2(a) MP_CHECKOK( mpl_rsh((a), (a), 1) );
253*f9fbec18Smcpowers
254*f9fbec18Smcpowers /* Compute binary polynomial multiply d = a * b */
255*f9fbec18Smcpowers static void
s_bmul_d(const mp_digit * a,mp_size a_len,mp_digit b,mp_digit * d)256*f9fbec18Smcpowers s_bmul_d(const mp_digit *a, mp_size a_len, mp_digit b, mp_digit *d)
257*f9fbec18Smcpowers {
258*f9fbec18Smcpowers mp_digit a_i, a0b0, a1b1, carry = 0;
259*f9fbec18Smcpowers while (a_len--) {
260*f9fbec18Smcpowers a_i = *a++;
261*f9fbec18Smcpowers s_bmul_1x1(&a1b1, &a0b0, a_i, b);
262*f9fbec18Smcpowers *d++ = a0b0 ^ carry;
263*f9fbec18Smcpowers carry = a1b1;
264*f9fbec18Smcpowers }
265*f9fbec18Smcpowers *d = carry;
266*f9fbec18Smcpowers }
267*f9fbec18Smcpowers
268*f9fbec18Smcpowers /* Compute binary polynomial xor multiply accumulate d ^= a * b */
269*f9fbec18Smcpowers static void
s_bmul_d_add(const mp_digit * a,mp_size a_len,mp_digit b,mp_digit * d)270*f9fbec18Smcpowers s_bmul_d_add(const mp_digit *a, mp_size a_len, mp_digit b, mp_digit *d)
271*f9fbec18Smcpowers {
272*f9fbec18Smcpowers mp_digit a_i, a0b0, a1b1, carry = 0;
273*f9fbec18Smcpowers while (a_len--) {
274*f9fbec18Smcpowers a_i = *a++;
275*f9fbec18Smcpowers s_bmul_1x1(&a1b1, &a0b0, a_i, b);
276*f9fbec18Smcpowers *d++ ^= a0b0 ^ carry;
277*f9fbec18Smcpowers carry = a1b1;
278*f9fbec18Smcpowers }
279*f9fbec18Smcpowers *d ^= carry;
280*f9fbec18Smcpowers }
281*f9fbec18Smcpowers
282*f9fbec18Smcpowers /* Compute binary polynomial xor multiply c = a * b.
283*f9fbec18Smcpowers * All parameters may be identical.
284*f9fbec18Smcpowers */
285*f9fbec18Smcpowers mp_err
mp_bmul(const mp_int * a,const mp_int * b,mp_int * c)286*f9fbec18Smcpowers mp_bmul(const mp_int *a, const mp_int *b, mp_int *c)
287*f9fbec18Smcpowers {
288*f9fbec18Smcpowers mp_digit *pb, b_i;
289*f9fbec18Smcpowers mp_int tmp;
290*f9fbec18Smcpowers mp_size ib, a_used, b_used;
291*f9fbec18Smcpowers mp_err res = MP_OKAY;
292*f9fbec18Smcpowers
293*f9fbec18Smcpowers MP_DIGITS(&tmp) = 0;
294*f9fbec18Smcpowers
295*f9fbec18Smcpowers ARGCHK(a != NULL && b != NULL && c != NULL, MP_BADARG);
296*f9fbec18Smcpowers
297*f9fbec18Smcpowers if (a == c) {
298*f9fbec18Smcpowers MP_CHECKOK( mp_init_copy(&tmp, a) );
299*f9fbec18Smcpowers if (a == b)
300*f9fbec18Smcpowers b = &tmp;
301*f9fbec18Smcpowers a = &tmp;
302*f9fbec18Smcpowers } else if (b == c) {
303*f9fbec18Smcpowers MP_CHECKOK( mp_init_copy(&tmp, b) );
304*f9fbec18Smcpowers b = &tmp;
305*f9fbec18Smcpowers }
306*f9fbec18Smcpowers
307*f9fbec18Smcpowers if (MP_USED(a) < MP_USED(b)) {
308*f9fbec18Smcpowers const mp_int *xch = b; /* switch a and b if b longer */
309*f9fbec18Smcpowers b = a;
310*f9fbec18Smcpowers a = xch;
311*f9fbec18Smcpowers }
312*f9fbec18Smcpowers
313*f9fbec18Smcpowers MP_USED(c) = 1; MP_DIGIT(c, 0) = 0;
314*f9fbec18Smcpowers MP_CHECKOK( s_mp_pad(c, USED(a) + USED(b)) );
315*f9fbec18Smcpowers
316*f9fbec18Smcpowers pb = MP_DIGITS(b);
317*f9fbec18Smcpowers s_bmul_d(MP_DIGITS(a), MP_USED(a), *pb++, MP_DIGITS(c));
318*f9fbec18Smcpowers
319*f9fbec18Smcpowers /* Outer loop: Digits of b */
320*f9fbec18Smcpowers a_used = MP_USED(a);
321*f9fbec18Smcpowers b_used = MP_USED(b);
322*f9fbec18Smcpowers MP_USED(c) = a_used + b_used;
323*f9fbec18Smcpowers for (ib = 1; ib < b_used; ib++) {
324*f9fbec18Smcpowers b_i = *pb++;
325*f9fbec18Smcpowers
326*f9fbec18Smcpowers /* Inner product: Digits of a */
327*f9fbec18Smcpowers if (b_i)
328*f9fbec18Smcpowers s_bmul_d_add(MP_DIGITS(a), a_used, b_i, MP_DIGITS(c) + ib);
329*f9fbec18Smcpowers else
330*f9fbec18Smcpowers MP_DIGIT(c, ib + a_used) = b_i;
331*f9fbec18Smcpowers }
332*f9fbec18Smcpowers
333*f9fbec18Smcpowers s_mp_clamp(c);
334*f9fbec18Smcpowers
335*f9fbec18Smcpowers SIGN(c) = ZPOS;
336*f9fbec18Smcpowers
337*f9fbec18Smcpowers CLEANUP:
338*f9fbec18Smcpowers mp_clear(&tmp);
339*f9fbec18Smcpowers return res;
340*f9fbec18Smcpowers }
341*f9fbec18Smcpowers
342*f9fbec18Smcpowers
343*f9fbec18Smcpowers /* Compute modular reduction of a and store result in r.
344*f9fbec18Smcpowers * r could be a.
345*f9fbec18Smcpowers * For modular arithmetic, the irreducible polynomial f(t) is represented
346*f9fbec18Smcpowers * as an array of int[], where f(t) is of the form:
347*f9fbec18Smcpowers * f(t) = t^p[0] + t^p[1] + ... + t^p[k]
348*f9fbec18Smcpowers * where m = p[0] > p[1] > ... > p[k] = 0.
349*f9fbec18Smcpowers */
350*f9fbec18Smcpowers mp_err
mp_bmod(const mp_int * a,const unsigned int p[],mp_int * r)351*f9fbec18Smcpowers mp_bmod(const mp_int *a, const unsigned int p[], mp_int *r)
352*f9fbec18Smcpowers {
353*f9fbec18Smcpowers int j, k;
354*f9fbec18Smcpowers int n, dN, d0, d1;
355*f9fbec18Smcpowers mp_digit zz, *z, tmp;
356*f9fbec18Smcpowers mp_size used;
357*f9fbec18Smcpowers mp_err res = MP_OKAY;
358*f9fbec18Smcpowers
359*f9fbec18Smcpowers /* The algorithm does the reduction in place in r,
360*f9fbec18Smcpowers * if a != r, copy a into r first so reduction can be done in r
361*f9fbec18Smcpowers */
362*f9fbec18Smcpowers if (a != r) {
363*f9fbec18Smcpowers MP_CHECKOK( mp_copy(a, r) );
364*f9fbec18Smcpowers }
365*f9fbec18Smcpowers z = MP_DIGITS(r);
366*f9fbec18Smcpowers
367*f9fbec18Smcpowers /* start reduction */
368*f9fbec18Smcpowers dN = p[0] / MP_DIGIT_BITS;
369*f9fbec18Smcpowers used = MP_USED(r);
370*f9fbec18Smcpowers
371*f9fbec18Smcpowers for (j = used - 1; j > dN;) {
372*f9fbec18Smcpowers
373*f9fbec18Smcpowers zz = z[j];
374*f9fbec18Smcpowers if (zz == 0) {
375*f9fbec18Smcpowers j--; continue;
376*f9fbec18Smcpowers }
377*f9fbec18Smcpowers z[j] = 0;
378*f9fbec18Smcpowers
379*f9fbec18Smcpowers for (k = 1; p[k] > 0; k++) {
380*f9fbec18Smcpowers /* reducing component t^p[k] */
381*f9fbec18Smcpowers n = p[0] - p[k];
382*f9fbec18Smcpowers d0 = n % MP_DIGIT_BITS;
383*f9fbec18Smcpowers d1 = MP_DIGIT_BITS - d0;
384*f9fbec18Smcpowers n /= MP_DIGIT_BITS;
385*f9fbec18Smcpowers z[j-n] ^= (zz>>d0);
386*f9fbec18Smcpowers if (d0)
387*f9fbec18Smcpowers z[j-n-1] ^= (zz<<d1);
388*f9fbec18Smcpowers }
389*f9fbec18Smcpowers
390*f9fbec18Smcpowers /* reducing component t^0 */
391*f9fbec18Smcpowers n = dN;
392*f9fbec18Smcpowers d0 = p[0] % MP_DIGIT_BITS;
393*f9fbec18Smcpowers d1 = MP_DIGIT_BITS - d0;
394*f9fbec18Smcpowers z[j-n] ^= (zz >> d0);
395*f9fbec18Smcpowers if (d0)
396*f9fbec18Smcpowers z[j-n-1] ^= (zz << d1);
397*f9fbec18Smcpowers
398*f9fbec18Smcpowers }
399*f9fbec18Smcpowers
400*f9fbec18Smcpowers /* final round of reduction */
401*f9fbec18Smcpowers while (j == dN) {
402*f9fbec18Smcpowers
403*f9fbec18Smcpowers d0 = p[0] % MP_DIGIT_BITS;
404*f9fbec18Smcpowers zz = z[dN] >> d0;
405*f9fbec18Smcpowers if (zz == 0) break;
406*f9fbec18Smcpowers d1 = MP_DIGIT_BITS - d0;
407*f9fbec18Smcpowers
408*f9fbec18Smcpowers /* clear up the top d1 bits */
409*f9fbec18Smcpowers if (d0) z[dN] = (z[dN] << d1) >> d1;
410*f9fbec18Smcpowers *z ^= zz; /* reduction t^0 component */
411*f9fbec18Smcpowers
412*f9fbec18Smcpowers for (k = 1; p[k] > 0; k++) {
413*f9fbec18Smcpowers /* reducing component t^p[k]*/
414*f9fbec18Smcpowers n = p[k] / MP_DIGIT_BITS;
415*f9fbec18Smcpowers d0 = p[k] % MP_DIGIT_BITS;
416*f9fbec18Smcpowers d1 = MP_DIGIT_BITS - d0;
417*f9fbec18Smcpowers z[n] ^= (zz << d0);
418*f9fbec18Smcpowers tmp = zz >> d1;
419*f9fbec18Smcpowers if (d0 && tmp)
420*f9fbec18Smcpowers z[n+1] ^= tmp;
421*f9fbec18Smcpowers }
422*f9fbec18Smcpowers }
423*f9fbec18Smcpowers
424*f9fbec18Smcpowers s_mp_clamp(r);
425*f9fbec18Smcpowers CLEANUP:
426*f9fbec18Smcpowers return res;
427*f9fbec18Smcpowers }
428*f9fbec18Smcpowers
429*f9fbec18Smcpowers /* Compute the product of two polynomials a and b, reduce modulo p,
430*f9fbec18Smcpowers * Store the result in r. r could be a or b; a could be b.
431*f9fbec18Smcpowers */
432*f9fbec18Smcpowers mp_err
mp_bmulmod(const mp_int * a,const mp_int * b,const unsigned int p[],mp_int * r)433*f9fbec18Smcpowers mp_bmulmod(const mp_int *a, const mp_int *b, const unsigned int p[], mp_int *r)
434*f9fbec18Smcpowers {
435*f9fbec18Smcpowers mp_err res;
436*f9fbec18Smcpowers
437*f9fbec18Smcpowers if (a == b) return mp_bsqrmod(a, p, r);
438*f9fbec18Smcpowers if ((res = mp_bmul(a, b, r) ) != MP_OKAY)
439*f9fbec18Smcpowers return res;
440*f9fbec18Smcpowers return mp_bmod(r, p, r);
441*f9fbec18Smcpowers }
442*f9fbec18Smcpowers
443*f9fbec18Smcpowers /* Compute binary polynomial squaring c = a*a mod p .
444*f9fbec18Smcpowers * Parameter r and a can be identical.
445*f9fbec18Smcpowers */
446*f9fbec18Smcpowers
447*f9fbec18Smcpowers mp_err
mp_bsqrmod(const mp_int * a,const unsigned int p[],mp_int * r)448*f9fbec18Smcpowers mp_bsqrmod(const mp_int *a, const unsigned int p[], mp_int *r)
449*f9fbec18Smcpowers {
450*f9fbec18Smcpowers mp_digit *pa, *pr, a_i;
451*f9fbec18Smcpowers mp_int tmp;
452*f9fbec18Smcpowers mp_size ia, a_used;
453*f9fbec18Smcpowers mp_err res;
454*f9fbec18Smcpowers
455*f9fbec18Smcpowers ARGCHK(a != NULL && r != NULL, MP_BADARG);
456*f9fbec18Smcpowers MP_DIGITS(&tmp) = 0;
457*f9fbec18Smcpowers
458*f9fbec18Smcpowers if (a == r) {
459*f9fbec18Smcpowers MP_CHECKOK( mp_init_copy(&tmp, a) );
460*f9fbec18Smcpowers a = &tmp;
461*f9fbec18Smcpowers }
462*f9fbec18Smcpowers
463*f9fbec18Smcpowers MP_USED(r) = 1; MP_DIGIT(r, 0) = 0;
464*f9fbec18Smcpowers MP_CHECKOK( s_mp_pad(r, 2*USED(a)) );
465*f9fbec18Smcpowers
466*f9fbec18Smcpowers pa = MP_DIGITS(a);
467*f9fbec18Smcpowers pr = MP_DIGITS(r);
468*f9fbec18Smcpowers a_used = MP_USED(a);
469*f9fbec18Smcpowers MP_USED(r) = 2 * a_used;
470*f9fbec18Smcpowers
471*f9fbec18Smcpowers for (ia = 0; ia < a_used; ia++) {
472*f9fbec18Smcpowers a_i = *pa++;
473*f9fbec18Smcpowers *pr++ = gf2m_SQR0(a_i);
474*f9fbec18Smcpowers *pr++ = gf2m_SQR1(a_i);
475*f9fbec18Smcpowers }
476*f9fbec18Smcpowers
477*f9fbec18Smcpowers MP_CHECKOK( mp_bmod(r, p, r) );
478*f9fbec18Smcpowers s_mp_clamp(r);
479*f9fbec18Smcpowers SIGN(r) = ZPOS;
480*f9fbec18Smcpowers
481*f9fbec18Smcpowers CLEANUP:
482*f9fbec18Smcpowers mp_clear(&tmp);
483*f9fbec18Smcpowers return res;
484*f9fbec18Smcpowers }
485*f9fbec18Smcpowers
486*f9fbec18Smcpowers /* Compute binary polynomial y/x mod p, y divided by x, reduce modulo p.
487*f9fbec18Smcpowers * Store the result in r. r could be x or y, and x could equal y.
488*f9fbec18Smcpowers * Uses algorithm Modular_Division_GF(2^m) from
489*f9fbec18Smcpowers * Chang-Shantz, S. "From Euclid's GCD to Montgomery Multiplication to
490*f9fbec18Smcpowers * the Great Divide".
491*f9fbec18Smcpowers */
492*f9fbec18Smcpowers int
mp_bdivmod(const mp_int * y,const mp_int * x,const mp_int * pp,const unsigned int p[],mp_int * r)493*f9fbec18Smcpowers mp_bdivmod(const mp_int *y, const mp_int *x, const mp_int *pp,
494*f9fbec18Smcpowers const unsigned int p[], mp_int *r)
495*f9fbec18Smcpowers {
496*f9fbec18Smcpowers mp_int aa, bb, uu;
497*f9fbec18Smcpowers mp_int *a, *b, *u, *v;
498*f9fbec18Smcpowers mp_err res = MP_OKAY;
499*f9fbec18Smcpowers
500*f9fbec18Smcpowers MP_DIGITS(&aa) = 0;
501*f9fbec18Smcpowers MP_DIGITS(&bb) = 0;
502*f9fbec18Smcpowers MP_DIGITS(&uu) = 0;
503*f9fbec18Smcpowers
504*f9fbec18Smcpowers MP_CHECKOK( mp_init_copy(&aa, x) );
505*f9fbec18Smcpowers MP_CHECKOK( mp_init_copy(&uu, y) );
506*f9fbec18Smcpowers MP_CHECKOK( mp_init_copy(&bb, pp) );
507*f9fbec18Smcpowers MP_CHECKOK( s_mp_pad(r, USED(pp)) );
508*f9fbec18Smcpowers MP_USED(r) = 1; MP_DIGIT(r, 0) = 0;
509*f9fbec18Smcpowers
510*f9fbec18Smcpowers a = &aa; b= &bb; u=&uu; v=r;
511*f9fbec18Smcpowers /* reduce x and y mod p */
512*f9fbec18Smcpowers MP_CHECKOK( mp_bmod(a, p, a) );
513*f9fbec18Smcpowers MP_CHECKOK( mp_bmod(u, p, u) );
514*f9fbec18Smcpowers
515*f9fbec18Smcpowers while (!mp_isodd(a)) {
516*f9fbec18Smcpowers s_mp_div2(a);
517*f9fbec18Smcpowers if (mp_isodd(u)) {
518*f9fbec18Smcpowers MP_CHECKOK( mp_badd(u, pp, u) );
519*f9fbec18Smcpowers }
520*f9fbec18Smcpowers s_mp_div2(u);
521*f9fbec18Smcpowers }
522*f9fbec18Smcpowers
523*f9fbec18Smcpowers do {
524*f9fbec18Smcpowers if (mp_cmp_mag(b, a) > 0) {
525*f9fbec18Smcpowers MP_CHECKOK( mp_badd(b, a, b) );
526*f9fbec18Smcpowers MP_CHECKOK( mp_badd(v, u, v) );
527*f9fbec18Smcpowers do {
528*f9fbec18Smcpowers s_mp_div2(b);
529*f9fbec18Smcpowers if (mp_isodd(v)) {
530*f9fbec18Smcpowers MP_CHECKOK( mp_badd(v, pp, v) );
531*f9fbec18Smcpowers }
532*f9fbec18Smcpowers s_mp_div2(v);
533*f9fbec18Smcpowers } while (!mp_isodd(b));
534*f9fbec18Smcpowers }
535*f9fbec18Smcpowers else if ((MP_DIGIT(a,0) == 1) && (MP_USED(a) == 1))
536*f9fbec18Smcpowers break;
537*f9fbec18Smcpowers else {
538*f9fbec18Smcpowers MP_CHECKOK( mp_badd(a, b, a) );
539*f9fbec18Smcpowers MP_CHECKOK( mp_badd(u, v, u) );
540*f9fbec18Smcpowers do {
541*f9fbec18Smcpowers s_mp_div2(a);
542*f9fbec18Smcpowers if (mp_isodd(u)) {
543*f9fbec18Smcpowers MP_CHECKOK( mp_badd(u, pp, u) );
544*f9fbec18Smcpowers }
545*f9fbec18Smcpowers s_mp_div2(u);
546*f9fbec18Smcpowers } while (!mp_isodd(a));
547*f9fbec18Smcpowers }
548*f9fbec18Smcpowers } while (1);
549*f9fbec18Smcpowers
550*f9fbec18Smcpowers MP_CHECKOK( mp_copy(u, r) );
551*f9fbec18Smcpowers
552*f9fbec18Smcpowers CLEANUP:
553*f9fbec18Smcpowers /* XXX this appears to be a memory leak in the NSS code */
554*f9fbec18Smcpowers mp_clear(&aa);
555*f9fbec18Smcpowers mp_clear(&bb);
556*f9fbec18Smcpowers mp_clear(&uu);
557*f9fbec18Smcpowers return res;
558*f9fbec18Smcpowers
559*f9fbec18Smcpowers }
560*f9fbec18Smcpowers
561*f9fbec18Smcpowers /* Convert the bit-string representation of a polynomial a into an array
562*f9fbec18Smcpowers * of integers corresponding to the bits with non-zero coefficient.
563*f9fbec18Smcpowers * Up to max elements of the array will be filled. Return value is total
564*f9fbec18Smcpowers * number of coefficients that would be extracted if array was large enough.
565*f9fbec18Smcpowers */
566*f9fbec18Smcpowers int
mp_bpoly2arr(const mp_int * a,unsigned int p[],int max)567*f9fbec18Smcpowers mp_bpoly2arr(const mp_int *a, unsigned int p[], int max)
568*f9fbec18Smcpowers {
569*f9fbec18Smcpowers int i, j, k;
570*f9fbec18Smcpowers mp_digit top_bit, mask;
571*f9fbec18Smcpowers
572*f9fbec18Smcpowers top_bit = 1;
573*f9fbec18Smcpowers top_bit <<= MP_DIGIT_BIT - 1;
574*f9fbec18Smcpowers
575*f9fbec18Smcpowers for (k = 0; k < max; k++) p[k] = 0;
576*f9fbec18Smcpowers k = 0;
577*f9fbec18Smcpowers
578*f9fbec18Smcpowers for (i = MP_USED(a) - 1; i >= 0; i--) {
579*f9fbec18Smcpowers mask = top_bit;
580*f9fbec18Smcpowers for (j = MP_DIGIT_BIT - 1; j >= 0; j--) {
581*f9fbec18Smcpowers if (MP_DIGITS(a)[i] & mask) {
582*f9fbec18Smcpowers if (k < max) p[k] = MP_DIGIT_BIT * i + j;
583*f9fbec18Smcpowers k++;
584*f9fbec18Smcpowers }
585*f9fbec18Smcpowers mask >>= 1;
586*f9fbec18Smcpowers }
587*f9fbec18Smcpowers }
588*f9fbec18Smcpowers
589*f9fbec18Smcpowers return k;
590*f9fbec18Smcpowers }
591*f9fbec18Smcpowers
592*f9fbec18Smcpowers /* Convert the coefficient array representation of a polynomial to a
593*f9fbec18Smcpowers * bit-string. The array must be terminated by 0.
594*f9fbec18Smcpowers */
595*f9fbec18Smcpowers mp_err
mp_barr2poly(const unsigned int p[],mp_int * a)596*f9fbec18Smcpowers mp_barr2poly(const unsigned int p[], mp_int *a)
597*f9fbec18Smcpowers {
598*f9fbec18Smcpowers
599*f9fbec18Smcpowers mp_err res = MP_OKAY;
600*f9fbec18Smcpowers int i;
601*f9fbec18Smcpowers
602*f9fbec18Smcpowers mp_zero(a);
603*f9fbec18Smcpowers for (i = 0; p[i] > 0; i++) {
604*f9fbec18Smcpowers MP_CHECKOK( mpl_set_bit(a, p[i], 1) );
605*f9fbec18Smcpowers }
606*f9fbec18Smcpowers MP_CHECKOK( mpl_set_bit(a, 0, 1) );
607*f9fbec18Smcpowers
608*f9fbec18Smcpowers CLEANUP:
609*f9fbec18Smcpowers return res;
610*f9fbec18Smcpowers }
611