Ashley Mills / cyassl-lib

This is a port of cyaSSL 2.7.0.

Dependents:   CyaSSL_DTLS_Cellular CyaSSL_DTLS_Ethernet

ctaocrypt/src/tfm.c@1:c0ce1562443a, 2013-09-05 (annotated)

Committer:
ashleymills
Date:
Thu Sep 05 15:55:50 2013 +0000
Revision:
1:c0ce1562443a
Parent:
0:714293de3836 
Nothing;

Who changed what in which revision?

UserRevisionLine numberNew contents of line
ashleymills 0:714293de3836 1 /* tfm.c
ashleymills 0:714293de3836 2 *
ashleymills 0:714293de3836 3 * Copyright (C) 2006-2013 wolfSSL Inc.
ashleymills 0:714293de3836 4 *
ashleymills 0:714293de3836 5 * This file is part of CyaSSL.
ashleymills 0:714293de3836 6 *
ashleymills 0:714293de3836 7 * CyaSSL is free software; you can redistribute it and/or modify
ashleymills 0:714293de3836 8 * it under the terms of the GNU General Public License as published by
ashleymills 0:714293de3836 9 * the Free Software Foundation; either version 2 of the License, or
ashleymills 0:714293de3836 10 * (at your option) any later version.
ashleymills 0:714293de3836 11 *
ashleymills 0:714293de3836 12 * CyaSSL is distributed in the hope that it will be useful,
ashleymills 0:714293de3836 13 * but WITHOUT ANY WARRANTY; without even the implied warranty of
ashleymills 0:714293de3836 14 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
ashleymills 0:714293de3836 15 * GNU General Public License for more details.
ashleymills 0:714293de3836 16 *
ashleymills 0:714293de3836 17 * You should have received a copy of the GNU General Public License
ashleymills 0:714293de3836 18 * along with this program; if not, write to the Free Software
ashleymills 0:714293de3836 19 * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA
ashleymills 0:714293de3836 20 */
ashleymills 0:714293de3836 21
ashleymills 0:714293de3836 22
ashleymills 0:714293de3836 23 /*
ashleymills 0:714293de3836 24 * Based on public domain TomsFastMath 0.10 by Tom St Denis, tomstdenis@iahu.ca,
ashleymills 0:714293de3836 25 * http://math.libtomcrypt.com
ashleymills 0:714293de3836 26 */
ashleymills 0:714293de3836 27
ashleymills 0:714293de3836 28 /**
ashleymills 0:714293de3836 29 * Edited by Moisés Guimarães (moisesguimaraesm@gmail.com)
ashleymills 0:714293de3836 30 * to fit CyaSSL's needs.
ashleymills 0:714293de3836 31 */
ashleymills 0:714293de3836 32
ashleymills 0:714293de3836 33 #ifdef HAVE_CONFIG_H
ashleymills 0:714293de3836 34 #include <config.h>
ashleymills 0:714293de3836 35 #endif
ashleymills 0:714293de3836 36
ashleymills 0:714293de3836 37 /* in case user set USE_FAST_MATH there */
ashleymills 0:714293de3836 38 #include <cyassl/ctaocrypt/settings.h>
ashleymills 0:714293de3836 39
ashleymills 0:714293de3836 40 #ifdef USE_FAST_MATH
ashleymills 0:714293de3836 41
ashleymills 0:714293de3836 42 #include <cyassl/ctaocrypt/tfm.h>
ashleymills 0:714293de3836 43 #include <ctaocrypt/src/asm.c> /* will define asm MACROS or C ones */
ashleymills 0:714293de3836 44
ashleymills 0:714293de3836 45
ashleymills 0:714293de3836 46 /* math settings check */
ashleymills 0:714293de3836 47 word32 CheckRunTimeSettings(void)
ashleymills 0:714293de3836 48 {
ashleymills 0:714293de3836 49 return CTC_SETTINGS;
ashleymills 0:714293de3836 50 }
ashleymills 0:714293de3836 51
ashleymills 0:714293de3836 52
ashleymills 0:714293de3836 53 /* math settings size check */
ashleymills 0:714293de3836 54 word32 CheckRunTimeFastMath(void)
ashleymills 0:714293de3836 55 {
ashleymills 0:714293de3836 56 return FP_SIZE;
ashleymills 0:714293de3836 57 }
ashleymills 0:714293de3836 58
ashleymills 0:714293de3836 59
ashleymills 0:714293de3836 60 /* Functions */
ashleymills 0:714293de3836 61
ashleymills 0:714293de3836 62 void fp_add(fp_int *a, fp_int *b, fp_int *c)
ashleymills 0:714293de3836 63 {
ashleymills 0:714293de3836 64 int sa, sb;
ashleymills 0:714293de3836 65
ashleymills 0:714293de3836 66 /* get sign of both inputs */
ashleymills 0:714293de3836 67 sa = a->sign;
ashleymills 0:714293de3836 68 sb = b->sign;
ashleymills 0:714293de3836 69
ashleymills 0:714293de3836 70 /* handle two cases, not four */
ashleymills 0:714293de3836 71 if (sa == sb) {
ashleymills 0:714293de3836 72 /* both positive or both negative */
ashleymills 0:714293de3836 73 /* add their magnitudes, copy the sign */
ashleymills 0:714293de3836 74 c->sign = sa;
ashleymills 0:714293de3836 75 s_fp_add (a, b, c);
ashleymills 0:714293de3836 76 } else {
ashleymills 0:714293de3836 77 /* one positive, the other negative */
ashleymills 0:714293de3836 78 /* subtract the one with the greater magnitude from */
ashleymills 0:714293de3836 79 /* the one of the lesser magnitude. The result gets */
ashleymills 0:714293de3836 80 /* the sign of the one with the greater magnitude. */
ashleymills 0:714293de3836 81 if (fp_cmp_mag (a, b) == FP_LT) {
ashleymills 0:714293de3836 82 c->sign = sb;
ashleymills 0:714293de3836 83 s_fp_sub (b, a, c);
ashleymills 0:714293de3836 84 } else {
ashleymills 0:714293de3836 85 c->sign = sa;
ashleymills 0:714293de3836 86 s_fp_sub (a, b, c);
ashleymills 0:714293de3836 87 }
ashleymills 0:714293de3836 88 }
ashleymills 0:714293de3836 89 }
ashleymills 0:714293de3836 90
ashleymills 0:714293de3836 91 /* unsigned addition */
ashleymills 0:714293de3836 92 void s_fp_add(fp_int *a, fp_int *b, fp_int *c)
ashleymills 0:714293de3836 93 {
ashleymills 0:714293de3836 94 int x, y, oldused;
ashleymills 0:714293de3836 95 register fp_word t;
ashleymills 0:714293de3836 96
ashleymills 0:714293de3836 97 y = MAX(a->used, b->used);
ashleymills 0:714293de3836 98 oldused = c->used;
ashleymills 0:714293de3836 99 c->used = y;
ashleymills 0:714293de3836 100
ashleymills 0:714293de3836 101 t = 0;
ashleymills 0:714293de3836 102 for (x = 0; x < y; x++) {
ashleymills 0:714293de3836 103 t += ((fp_word)a->dp[x]) + ((fp_word)b->dp[x]);
ashleymills 0:714293de3836 104 c->dp[x] = (fp_digit)t;
ashleymills 0:714293de3836 105 t >>= DIGIT_BIT;
ashleymills 0:714293de3836 106 }
ashleymills 0:714293de3836 107 if (t != 0 && x < FP_SIZE) {
ashleymills 0:714293de3836 108 c->dp[c->used++] = (fp_digit)t;
ashleymills 0:714293de3836 109 ++x;
ashleymills 0:714293de3836 110 }
ashleymills 0:714293de3836 111
ashleymills 0:714293de3836 112 c->used = x;
ashleymills 0:714293de3836 113 for (; x < oldused; x++) {
ashleymills 0:714293de3836 114 c->dp[x] = 0;
ashleymills 0:714293de3836 115 }
ashleymills 0:714293de3836 116 fp_clamp(c);
ashleymills 0:714293de3836 117 }
ashleymills 0:714293de3836 118
ashleymills 0:714293de3836 119 /* c = a - b */
ashleymills 0:714293de3836 120 void fp_sub(fp_int *a, fp_int *b, fp_int *c)
ashleymills 0:714293de3836 121 {
ashleymills 0:714293de3836 122 int sa, sb;
ashleymills 0:714293de3836 123
ashleymills 0:714293de3836 124 sa = a->sign;
ashleymills 0:714293de3836 125 sb = b->sign;
ashleymills 0:714293de3836 126
ashleymills 0:714293de3836 127 if (sa != sb) {
ashleymills 0:714293de3836 128 /* subtract a negative from a positive, OR */
ashleymills 0:714293de3836 129 /* subtract a positive from a negative. */
ashleymills 0:714293de3836 130 /* In either case, ADD their magnitudes, */
ashleymills 0:714293de3836 131 /* and use the sign of the first number. */
ashleymills 0:714293de3836 132 c->sign = sa;
ashleymills 0:714293de3836 133 s_fp_add (a, b, c);
ashleymills 0:714293de3836 134 } else {
ashleymills 0:714293de3836 135 /* subtract a positive from a positive, OR */
ashleymills 0:714293de3836 136 /* subtract a negative from a negative. */
ashleymills 0:714293de3836 137 /* First, take the difference between their */
ashleymills 0:714293de3836 138 /* magnitudes, then... */
ashleymills 0:714293de3836 139 if (fp_cmp_mag (a, b) != FP_LT) {
ashleymills 0:714293de3836 140 /* Copy the sign from the first */
ashleymills 0:714293de3836 141 c->sign = sa;
ashleymills 0:714293de3836 142 /* The first has a larger or equal magnitude */
ashleymills 0:714293de3836 143 s_fp_sub (a, b, c);
ashleymills 0:714293de3836 144 } else {
ashleymills 0:714293de3836 145 /* The result has the *opposite* sign from */
ashleymills 0:714293de3836 146 /* the first number. */
ashleymills 0:714293de3836 147 c->sign = (sa == FP_ZPOS) ? FP_NEG : FP_ZPOS;
ashleymills 0:714293de3836 148 /* The second has a larger magnitude */
ashleymills 0:714293de3836 149 s_fp_sub (b, a, c);
ashleymills 0:714293de3836 150 }
ashleymills 0:714293de3836 151 }
ashleymills 0:714293de3836 152 }
ashleymills 0:714293de3836 153
ashleymills 0:714293de3836 154 /* unsigned subtraction ||a|| >= ||b|| ALWAYS! */
ashleymills 0:714293de3836 155 void s_fp_sub(fp_int *a, fp_int *b, fp_int *c)
ashleymills 0:714293de3836 156 {
ashleymills 0:714293de3836 157 int x, oldbused, oldused;
ashleymills 0:714293de3836 158 fp_word t;
ashleymills 0:714293de3836 159
ashleymills 0:714293de3836 160 oldused = c->used;
ashleymills 0:714293de3836 161 oldbused = b->used;
ashleymills 0:714293de3836 162 c->used = a->used;
ashleymills 0:714293de3836 163 t = 0;
ashleymills 0:714293de3836 164 for (x = 0; x < oldbused; x++) {
ashleymills 0:714293de3836 165 t = ((fp_word)a->dp[x]) - (((fp_word)b->dp[x]) + t);
ashleymills 0:714293de3836 166 c->dp[x] = (fp_digit)t;
ashleymills 0:714293de3836 167 t = (t >> DIGIT_BIT)&1;
ashleymills 0:714293de3836 168 }
ashleymills 0:714293de3836 169 for (; x < a->used; x++) {
ashleymills 0:714293de3836 170 t = ((fp_word)a->dp[x]) - t;
ashleymills 0:714293de3836 171 c->dp[x] = (fp_digit)t;
ashleymills 0:714293de3836 172 t = (t >> DIGIT_BIT)&1;
ashleymills 0:714293de3836 173 }
ashleymills 0:714293de3836 174 for (; x < oldused; x++) {
ashleymills 0:714293de3836 175 c->dp[x] = 0;
ashleymills 0:714293de3836 176 }
ashleymills 0:714293de3836 177 fp_clamp(c);
ashleymills 0:714293de3836 178 }
ashleymills 0:714293de3836 179
ashleymills 0:714293de3836 180 /* c = a * b */
ashleymills 0:714293de3836 181 void fp_mul(fp_int *A, fp_int *B, fp_int *C)
ashleymills 0:714293de3836 182 {
ashleymills 0:714293de3836 183 int y, yy;
ashleymills 0:714293de3836 184
ashleymills 0:714293de3836 185 y = MAX(A->used, B->used);
ashleymills 0:714293de3836 186 yy = MIN(A->used, B->used);
ashleymills 0:714293de3836 187
ashleymills 0:714293de3836 188 /* call generic if we're out of range */
ashleymills 0:714293de3836 189 if (y + yy > FP_SIZE) {
ashleymills 0:714293de3836 190 fp_mul_comba(A, B, C);
ashleymills 0:714293de3836 191 return ;
ashleymills 0:714293de3836 192 }
ashleymills 0:714293de3836 193
ashleymills 0:714293de3836 194 /* pick a comba (unrolled 4/8/16/32 x or rolled) based on the size
ashleymills 0:714293de3836 195 of the largest input. We also want to avoid doing excess mults if the
ashleymills 0:714293de3836 196 inputs are not close to the next power of two. That is, for example,
ashleymills 0:714293de3836 197 if say y=17 then we would do (32-17)^2 = 225 unneeded multiplications
ashleymills 0:714293de3836 198 */
ashleymills 0:714293de3836 199
ashleymills 0:714293de3836 200 #ifdef TFM_MUL3
ashleymills 0:714293de3836 201 if (y <= 3) {
ashleymills 0:714293de3836 202 fp_mul_comba3(A,B,C);
ashleymills 0:714293de3836 203 return;
ashleymills 0:714293de3836 204 }
ashleymills 0:714293de3836 205 #endif
ashleymills 0:714293de3836 206 #ifdef TFM_MUL4
ashleymills 0:714293de3836 207 if (y == 4) {
ashleymills 0:714293de3836 208 fp_mul_comba4(A,B,C);
ashleymills 0:714293de3836 209 return;
ashleymills 0:714293de3836 210 }
ashleymills 0:714293de3836 211 #endif
ashleymills 0:714293de3836 212 #ifdef TFM_MUL6
ashleymills 0:714293de3836 213 if (y <= 6) {
ashleymills 0:714293de3836 214 fp_mul_comba6(A,B,C);
ashleymills 0:714293de3836 215 return;
ashleymills 0:714293de3836 216 }
ashleymills 0:714293de3836 217 #endif
ashleymills 0:714293de3836 218 #ifdef TFM_MUL7
ashleymills 0:714293de3836 219 if (y == 7) {
ashleymills 0:714293de3836 220 fp_mul_comba7(A,B,C);
ashleymills 0:714293de3836 221 return;
ashleymills 0:714293de3836 222 }
ashleymills 0:714293de3836 223 #endif
ashleymills 0:714293de3836 224 #ifdef TFM_MUL8
ashleymills 0:714293de3836 225 if (y == 8) {
ashleymills 0:714293de3836 226 fp_mul_comba8(A,B,C);
ashleymills 0:714293de3836 227 return;
ashleymills 0:714293de3836 228 }
ashleymills 0:714293de3836 229 #endif
ashleymills 0:714293de3836 230 #ifdef TFM_MUL9
ashleymills 0:714293de3836 231 if (y == 9) {
ashleymills 0:714293de3836 232 fp_mul_comba9(A,B,C);
ashleymills 0:714293de3836 233 return;
ashleymills 0:714293de3836 234 }
ashleymills 0:714293de3836 235 #endif
ashleymills 0:714293de3836 236 #ifdef TFM_MUL12
ashleymills 0:714293de3836 237 if (y <= 12) {
ashleymills 0:714293de3836 238 fp_mul_comba12(A,B,C);
ashleymills 0:714293de3836 239 return;
ashleymills 0:714293de3836 240 }
ashleymills 0:714293de3836 241 #endif
ashleymills 0:714293de3836 242 #ifdef TFM_MUL17
ashleymills 0:714293de3836 243 if (y <= 17) {
ashleymills 0:714293de3836 244 fp_mul_comba17(A,B,C);
ashleymills 0:714293de3836 245 return;
ashleymills 0:714293de3836 246 }
ashleymills 0:714293de3836 247 #endif
ashleymills 0:714293de3836 248
ashleymills 0:714293de3836 249 #ifdef TFM_SMALL_SET
ashleymills 0:714293de3836 250 if (y <= 16) {
ashleymills 0:714293de3836 251 fp_mul_comba_small(A,B,C);
ashleymills 0:714293de3836 252 return;
ashleymills 0:714293de3836 253 }
ashleymills 0:714293de3836 254 #endif
ashleymills 0:714293de3836 255 #if defined(TFM_MUL20)
ashleymills 0:714293de3836 256 if (y <= 20) {
ashleymills 0:714293de3836 257 fp_mul_comba20(A,B,C);
ashleymills 0:714293de3836 258 return;
ashleymills 0:714293de3836 259 }
ashleymills 0:714293de3836 260 #endif
ashleymills 0:714293de3836 261 #if defined(TFM_MUL24)
ashleymills 0:714293de3836 262 if (yy >= 16 && y <= 24) {
ashleymills 0:714293de3836 263 fp_mul_comba24(A,B,C);
ashleymills 0:714293de3836 264 return;
ashleymills 0:714293de3836 265 }
ashleymills 0:714293de3836 266 #endif
ashleymills 0:714293de3836 267 #if defined(TFM_MUL28)
ashleymills 0:714293de3836 268 if (yy >= 20 && y <= 28) {
ashleymills 0:714293de3836 269 fp_mul_comba28(A,B,C);
ashleymills 0:714293de3836 270 return;
ashleymills 0:714293de3836 271 }
ashleymills 0:714293de3836 272 #endif
ashleymills 0:714293de3836 273 #if defined(TFM_MUL32)
ashleymills 0:714293de3836 274 if (yy >= 24 && y <= 32) {
ashleymills 0:714293de3836 275 fp_mul_comba32(A,B,C);
ashleymills 0:714293de3836 276 return;
ashleymills 0:714293de3836 277 }
ashleymills 0:714293de3836 278 #endif
ashleymills 0:714293de3836 279 #if defined(TFM_MUL48)
ashleymills 0:714293de3836 280 if (yy >= 40 && y <= 48) {
ashleymills 0:714293de3836 281 fp_mul_comba48(A,B,C);
ashleymills 0:714293de3836 282 return;
ashleymills 0:714293de3836 283 }
ashleymills 0:714293de3836 284 #endif
ashleymills 0:714293de3836 285 #if defined(TFM_MUL64)
ashleymills 0:714293de3836 286 if (yy >= 56 && y <= 64) {
ashleymills 0:714293de3836 287 fp_mul_comba64(A,B,C);
ashleymills 0:714293de3836 288 return;
ashleymills 0:714293de3836 289 }
ashleymills 0:714293de3836 290 #endif
ashleymills 0:714293de3836 291 fp_mul_comba(A,B,C);
ashleymills 0:714293de3836 292 }
ashleymills 0:714293de3836 293
ashleymills 0:714293de3836 294 void fp_mul_2(fp_int * a, fp_int * b)
ashleymills 0:714293de3836 295 {
ashleymills 0:714293de3836 296 int x, oldused;
ashleymills 0:714293de3836 297
ashleymills 0:714293de3836 298 oldused = b->used;
ashleymills 0:714293de3836 299 b->used = a->used;
ashleymills 0:714293de3836 300
ashleymills 0:714293de3836 301 {
ashleymills 0:714293de3836 302 register fp_digit r, rr, *tmpa, *tmpb;
ashleymills 0:714293de3836 303
ashleymills 0:714293de3836 304 /* alias for source */
ashleymills 0:714293de3836 305 tmpa = a->dp;
ashleymills 0:714293de3836 306
ashleymills 0:714293de3836 307 /* alias for dest */
ashleymills 0:714293de3836 308 tmpb = b->dp;
ashleymills 0:714293de3836 309
ashleymills 0:714293de3836 310 /* carry */
ashleymills 0:714293de3836 311 r = 0;
ashleymills 0:714293de3836 312 for (x = 0; x < a->used; x++) {
ashleymills 0:714293de3836 313
ashleymills 0:714293de3836 314 /* get what will be the *next* carry bit from the
ashleymills 0:714293de3836 315 * MSB of the current digit
ashleymills 0:714293de3836 316 */
ashleymills 0:714293de3836 317 rr = *tmpa >> ((fp_digit)(DIGIT_BIT - 1));
ashleymills 0:714293de3836 318
ashleymills 0:714293de3836 319 /* now shift up this digit, add in the carry [from the previous] */
ashleymills 0:714293de3836 320 *tmpb++ = ((*tmpa++ << ((fp_digit)1)) | r);
ashleymills 0:714293de3836 321
ashleymills 0:714293de3836 322 /* copy the carry that would be from the source
ashleymills 0:714293de3836 323 * digit into the next iteration
ashleymills 0:714293de3836 324 */
ashleymills 0:714293de3836 325 r = rr;
ashleymills 0:714293de3836 326 }
ashleymills 0:714293de3836 327
ashleymills 0:714293de3836 328 /* new leading digit? */
ashleymills 0:714293de3836 329 if (r != 0 && b->used != (FP_SIZE-1)) {
ashleymills 0:714293de3836 330 /* add a MSB which is always 1 at this point */
ashleymills 0:714293de3836 331 *tmpb = 1;
ashleymills 0:714293de3836 332 ++(b->used);
ashleymills 0:714293de3836 333 }
ashleymills 0:714293de3836 334
ashleymills 0:714293de3836 335 /* now zero any excess digits on the destination
ashleymills 0:714293de3836 336 * that we didn't write to
ashleymills 0:714293de3836 337 */
ashleymills 0:714293de3836 338 tmpb = b->dp + b->used;
ashleymills 0:714293de3836 339 for (x = b->used; x < oldused; x++) {
ashleymills 0:714293de3836 340 *tmpb++ = 0;
ashleymills 0:714293de3836 341 }
ashleymills 0:714293de3836 342 }
ashleymills 0:714293de3836 343 b->sign = a->sign;
ashleymills 0:714293de3836 344 }
ashleymills 0:714293de3836 345
ashleymills 0:714293de3836 346 /* c = a * b */
ashleymills 0:714293de3836 347 void fp_mul_d(fp_int *a, fp_digit b, fp_int *c)
ashleymills 0:714293de3836 348 {
ashleymills 0:714293de3836 349 fp_word w;
ashleymills 0:714293de3836 350 int x, oldused;
ashleymills 0:714293de3836 351
ashleymills 0:714293de3836 352 oldused = c->used;
ashleymills 0:714293de3836 353 c->used = a->used;
ashleymills 0:714293de3836 354 c->sign = a->sign;
ashleymills 0:714293de3836 355 w = 0;
ashleymills 0:714293de3836 356 for (x = 0; x < a->used; x++) {
ashleymills 0:714293de3836 357 w = ((fp_word)a->dp[x]) * ((fp_word)b) + w;
ashleymills 0:714293de3836 358 c->dp[x] = (fp_digit)w;
ashleymills 0:714293de3836 359 w = w >> DIGIT_BIT;
ashleymills 0:714293de3836 360 }
ashleymills 0:714293de3836 361 if (w != 0 && (a->used != FP_SIZE)) {
ashleymills 0:714293de3836 362 c->dp[c->used++] = (fp_digit) w;
ashleymills 0:714293de3836 363 ++x;
ashleymills 0:714293de3836 364 }
ashleymills 0:714293de3836 365 for (; x < oldused; x++) {
ashleymills 0:714293de3836 366 c->dp[x] = 0;
ashleymills 0:714293de3836 367 }
ashleymills 0:714293de3836 368 fp_clamp(c);
ashleymills 0:714293de3836 369 }
ashleymills 0:714293de3836 370
ashleymills 0:714293de3836 371 /* c = a * 2**d */
ashleymills 0:714293de3836 372 void fp_mul_2d(fp_int *a, int b, fp_int *c)
ashleymills 0:714293de3836 373 {
ashleymills 0:714293de3836 374 fp_digit carry, carrytmp, shift;
ashleymills 0:714293de3836 375 int x;
ashleymills 0:714293de3836 376
ashleymills 0:714293de3836 377 /* copy it */
ashleymills 0:714293de3836 378 fp_copy(a, c);
ashleymills 0:714293de3836 379
ashleymills 0:714293de3836 380 /* handle whole digits */
ashleymills 0:714293de3836 381 if (b >= DIGIT_BIT) {
ashleymills 0:714293de3836 382 fp_lshd(c, b/DIGIT_BIT);
ashleymills 0:714293de3836 383 }
ashleymills 0:714293de3836 384 b %= DIGIT_BIT;
ashleymills 0:714293de3836 385
ashleymills 0:714293de3836 386 /* shift the digits */
ashleymills 0:714293de3836 387 if (b != 0) {
ashleymills 0:714293de3836 388 carry = 0;
ashleymills 0:714293de3836 389 shift = DIGIT_BIT - b;
ashleymills 0:714293de3836 390 for (x = 0; x < c->used; x++) {
ashleymills 0:714293de3836 391 carrytmp = c->dp[x] >> shift;
ashleymills 0:714293de3836 392 c->dp[x] = (c->dp[x] << b) + carry;
ashleymills 0:714293de3836 393 carry = carrytmp;
ashleymills 0:714293de3836 394 }
ashleymills 0:714293de3836 395 /* store last carry if room */
ashleymills 0:714293de3836 396 if (carry && x < FP_SIZE) {
ashleymills 0:714293de3836 397 c->dp[c->used++] = carry;
ashleymills 0:714293de3836 398 }
ashleymills 0:714293de3836 399 }
ashleymills 0:714293de3836 400 fp_clamp(c);
ashleymills 0:714293de3836 401 }
ashleymills 0:714293de3836 402
ashleymills 0:714293de3836 403 /* generic PxQ multiplier */
ashleymills 0:714293de3836 404 void fp_mul_comba(fp_int *A, fp_int *B, fp_int *C)
ashleymills 0:714293de3836 405 {
ashleymills 0:714293de3836 406 int ix, iy, iz, tx, ty, pa;
ashleymills 0:714293de3836 407 fp_digit c0, c1, c2, *tmpx, *tmpy;
ashleymills 0:714293de3836 408 fp_int tmp, *dst;
ashleymills 0:714293de3836 409
ashleymills 0:714293de3836 410 COMBA_START;
ashleymills 0:714293de3836 411 COMBA_CLEAR;
ashleymills 0:714293de3836 412
ashleymills 0:714293de3836 413 /* get size of output and trim */
ashleymills 0:714293de3836 414 pa = A->used + B->used;
ashleymills 0:714293de3836 415 if (pa >= FP_SIZE) {
ashleymills 0:714293de3836 416 pa = FP_SIZE-1;
ashleymills 0:714293de3836 417 }
ashleymills 0:714293de3836 418
ashleymills 0:714293de3836 419 if (A == C || B == C) {
ashleymills 0:714293de3836 420 fp_zero(&tmp);
ashleymills 0:714293de3836 421 dst = &tmp;
ashleymills 0:714293de3836 422 } else {
ashleymills 0:714293de3836 423 fp_zero(C);
ashleymills 0:714293de3836 424 dst = C;
ashleymills 0:714293de3836 425 }
ashleymills 0:714293de3836 426
ashleymills 0:714293de3836 427 for (ix = 0; ix < pa; ix++) {
ashleymills 0:714293de3836 428 /* get offsets into the two bignums */
ashleymills 0:714293de3836 429 ty = MIN(ix, B->used-1);
ashleymills 0:714293de3836 430 tx = ix - ty;
ashleymills 0:714293de3836 431
ashleymills 0:714293de3836 432 /* setup temp aliases */
ashleymills 0:714293de3836 433 tmpx = A->dp + tx;
ashleymills 0:714293de3836 434 tmpy = B->dp + ty;
ashleymills 0:714293de3836 435
ashleymills 0:714293de3836 436 /* this is the number of times the loop will iterrate, essentially its
ashleymills 0:714293de3836 437 while (tx++ < a->used && ty-- >= 0) { ... }
ashleymills 0:714293de3836 438 */
ashleymills 0:714293de3836 439 iy = MIN(A->used-tx, ty+1);
ashleymills 0:714293de3836 440
ashleymills 0:714293de3836 441 /* execute loop */
ashleymills 0:714293de3836 442 COMBA_FORWARD;
ashleymills 0:714293de3836 443 for (iz = 0; iz < iy; ++iz) {
ashleymills 0:714293de3836 444 /* TAO change COMBA_ADD back to MULADD */
ashleymills 0:714293de3836 445 MULADD(*tmpx++, *tmpy--);
ashleymills 0:714293de3836 446 }
ashleymills 0:714293de3836 447
ashleymills 0:714293de3836 448 /* store term */
ashleymills 0:714293de3836 449 COMBA_STORE(dst->dp[ix]);
ashleymills 0:714293de3836 450 }
ashleymills 0:714293de3836 451 COMBA_FINI;
ashleymills 0:714293de3836 452
ashleymills 0:714293de3836 453 dst->used = pa;
ashleymills 0:714293de3836 454 dst->sign = A->sign ^ B->sign;
ashleymills 0:714293de3836 455 fp_clamp(dst);
ashleymills 0:714293de3836 456 fp_copy(dst, C);
ashleymills 0:714293de3836 457 }
ashleymills 0:714293de3836 458
ashleymills 0:714293de3836 459 /* a/b => cb + d == a */
ashleymills 0:714293de3836 460 int fp_div(fp_int *a, fp_int *b, fp_int *c, fp_int *d)
ashleymills 0:714293de3836 461 {
ashleymills 0:714293de3836 462 fp_int q, x, y, t1, t2;
ashleymills 0:714293de3836 463 int n, t, i, norm, neg;
ashleymills 0:714293de3836 464
ashleymills 0:714293de3836 465 /* is divisor zero ? */
ashleymills 0:714293de3836 466 if (fp_iszero (b) == 1) {
ashleymills 0:714293de3836 467 return FP_VAL;
ashleymills 0:714293de3836 468 }
ashleymills 0:714293de3836 469
ashleymills 0:714293de3836 470 /* if a < b then q=0, r = a */
ashleymills 0:714293de3836 471 if (fp_cmp_mag (a, b) == FP_LT) {
ashleymills 0:714293de3836 472 if (d != NULL) {
ashleymills 0:714293de3836 473 fp_copy (a, d);
ashleymills 0:714293de3836 474 }
ashleymills 0:714293de3836 475 if (c != NULL) {
ashleymills 0:714293de3836 476 fp_zero (c);
ashleymills 0:714293de3836 477 }
ashleymills 0:714293de3836 478 return FP_OKAY;
ashleymills 0:714293de3836 479 }
ashleymills 0:714293de3836 480
ashleymills 0:714293de3836 481 fp_init(&q);
ashleymills 0:714293de3836 482 q.used = a->used + 2;
ashleymills 0:714293de3836 483
ashleymills 0:714293de3836 484 fp_init(&t1);
ashleymills 0:714293de3836 485 fp_init(&t2);
ashleymills 0:714293de3836 486 fp_init_copy(&x, a);
ashleymills 0:714293de3836 487 fp_init_copy(&y, b);
ashleymills 0:714293de3836 488
ashleymills 0:714293de3836 489 /* fix the sign */
ashleymills 0:714293de3836 490 neg = (a->sign == b->sign) ? FP_ZPOS : FP_NEG;
ashleymills 0:714293de3836 491 x.sign = y.sign = FP_ZPOS;
ashleymills 0:714293de3836 492
ashleymills 0:714293de3836 493 /* normalize both x and y, ensure that y >= b/2, [b == 2**DIGIT_BIT] */
ashleymills 0:714293de3836 494 norm = fp_count_bits(&y) % DIGIT_BIT;
ashleymills 0:714293de3836 495 if (norm < (int)(DIGIT_BIT-1)) {
ashleymills 0:714293de3836 496 norm = (DIGIT_BIT-1) - norm;
ashleymills 0:714293de3836 497 fp_mul_2d (&x, norm, &x);
ashleymills 0:714293de3836 498 fp_mul_2d (&y, norm, &y);
ashleymills 0:714293de3836 499 } else {
ashleymills 0:714293de3836 500 norm = 0;
ashleymills 0:714293de3836 501 }
ashleymills 0:714293de3836 502
ashleymills 0:714293de3836 503 /* note hac does 0 based, so if used==5 then its 0,1,2,3,4, e.g. use 4 */
ashleymills 0:714293de3836 504 n = x.used - 1;
ashleymills 0:714293de3836 505 t = y.used - 1;
ashleymills 0:714293de3836 506
ashleymills 0:714293de3836 507 /* while (x >= y*b**n-t) do { q[n-t] += 1; x -= y*b**{n-t} } */
ashleymills 0:714293de3836 508 fp_lshd (&y, n - t); /* y = y*b**{n-t} */
ashleymills 0:714293de3836 509
ashleymills 0:714293de3836 510 while (fp_cmp (&x, &y) != FP_LT) {
ashleymills 0:714293de3836 511 ++(q.dp[n - t]);
ashleymills 0:714293de3836 512 fp_sub (&x, &y, &x);
ashleymills 0:714293de3836 513 }
ashleymills 0:714293de3836 514
ashleymills 0:714293de3836 515 /* reset y by shifting it back down */
ashleymills 0:714293de3836 516 fp_rshd (&y, n - t);
ashleymills 0:714293de3836 517
ashleymills 0:714293de3836 518 /* step 3. for i from n down to (t + 1) */
ashleymills 0:714293de3836 519 for (i = n; i >= (t + 1); i--) {
ashleymills 0:714293de3836 520 if (i > x.used) {
ashleymills 0:714293de3836 521 continue;
ashleymills 0:714293de3836 522 }
ashleymills 0:714293de3836 523
ashleymills 0:714293de3836 524 /* step 3.1 if xi == yt then set q{i-t-1} to b-1,
ashleymills 0:714293de3836 525 * otherwise set q{i-t-1} to (xi*b + x{i-1})/yt */
ashleymills 0:714293de3836 526 if (x.dp[i] == y.dp[t]) {
ashleymills 0:714293de3836 527 q.dp[i - t - 1] = (fp_digit) ((((fp_word)1) << DIGIT_BIT) - 1);
ashleymills 0:714293de3836 528 } else {
ashleymills 0:714293de3836 529 fp_word tmp;
ashleymills 0:714293de3836 530 tmp = ((fp_word) x.dp[i]) << ((fp_word) DIGIT_BIT);
ashleymills 0:714293de3836 531 tmp |= ((fp_word) x.dp[i - 1]);
ashleymills 0:714293de3836 532 tmp /= ((fp_word)y.dp[t]);
ashleymills 0:714293de3836 533 q.dp[i - t - 1] = (fp_digit) (tmp);
ashleymills 0:714293de3836 534 }
ashleymills 0:714293de3836 535
ashleymills 0:714293de3836 536 /* while (q{i-t-1} * (yt * b + y{t-1})) >
ashleymills 0:714293de3836 537 xi * b**2 + xi-1 * b + xi-2
ashleymills 0:714293de3836 538
ashleymills 0:714293de3836 539 do q{i-t-1} -= 1;
ashleymills 0:714293de3836 540 */
ashleymills 0:714293de3836 541 q.dp[i - t - 1] = (q.dp[i - t - 1] + 1);
ashleymills 0:714293de3836 542 do {
ashleymills 0:714293de3836 543 q.dp[i - t - 1] = (q.dp[i - t - 1] - 1);
ashleymills 0:714293de3836 544
ashleymills 0:714293de3836 545 /* find left hand */
ashleymills 0:714293de3836 546 fp_zero (&t1);
ashleymills 0:714293de3836 547 t1.dp[0] = (t - 1 < 0) ? 0 : y.dp[t - 1];
ashleymills 0:714293de3836 548 t1.dp[1] = y.dp[t];
ashleymills 0:714293de3836 549 t1.used = 2;
ashleymills 0:714293de3836 550 fp_mul_d (&t1, q.dp[i - t - 1], &t1);
ashleymills 0:714293de3836 551
ashleymills 0:714293de3836 552 /* find right hand */
ashleymills 0:714293de3836 553 t2.dp[0] = (i - 2 < 0) ? 0 : x.dp[i - 2];
ashleymills 0:714293de3836 554 t2.dp[1] = (i - 1 < 0) ? 0 : x.dp[i - 1];
ashleymills 0:714293de3836 555 t2.dp[2] = x.dp[i];
ashleymills 0:714293de3836 556 t2.used = 3;
ashleymills 0:714293de3836 557 } while (fp_cmp_mag(&t1, &t2) == FP_GT);
ashleymills 0:714293de3836 558
ashleymills 0:714293de3836 559 /* step 3.3 x = x - q{i-t-1} * y * b**{i-t-1} */
ashleymills 0:714293de3836 560 fp_mul_d (&y, q.dp[i - t - 1], &t1);
ashleymills 0:714293de3836 561 fp_lshd (&t1, i - t - 1);
ashleymills 0:714293de3836 562 fp_sub (&x, &t1, &x);
ashleymills 0:714293de3836 563
ashleymills 0:714293de3836 564 /* if x < 0 then { x = x + y*b**{i-t-1}; q{i-t-1} -= 1; } */
ashleymills 0:714293de3836 565 if (x.sign == FP_NEG) {
ashleymills 0:714293de3836 566 fp_copy (&y, &t1);
ashleymills 0:714293de3836 567 fp_lshd (&t1, i - t - 1);
ashleymills 0:714293de3836 568 fp_add (&x, &t1, &x);
ashleymills 0:714293de3836 569 q.dp[i - t - 1] = q.dp[i - t - 1] - 1;
ashleymills 0:714293de3836 570 }
ashleymills 0:714293de3836 571 }
ashleymills 0:714293de3836 572
ashleymills 0:714293de3836 573 /* now q is the quotient and x is the remainder
ashleymills 0:714293de3836 574 * [which we have to normalize]
ashleymills 0:714293de3836 575 */
ashleymills 0:714293de3836 576
ashleymills 0:714293de3836 577 /* get sign before writing to c */
ashleymills 0:714293de3836 578 x.sign = x.used == 0 ? FP_ZPOS : a->sign;
ashleymills 0:714293de3836 579
ashleymills 0:714293de3836 580 if (c != NULL) {
ashleymills 0:714293de3836 581 fp_clamp (&q);
ashleymills 0:714293de3836 582 fp_copy (&q, c);
ashleymills 0:714293de3836 583 c->sign = neg;
ashleymills 0:714293de3836 584 }
ashleymills 0:714293de3836 585
ashleymills 0:714293de3836 586 if (d != NULL) {
ashleymills 0:714293de3836 587 fp_div_2d (&x, norm, &x, NULL);
ashleymills 0:714293de3836 588
ashleymills 0:714293de3836 589 /* the following is a kludge, essentially we were seeing the right remainder but
ashleymills 0:714293de3836 590 with excess digits that should have been zero
ashleymills 0:714293de3836 591 */
ashleymills 0:714293de3836 592 for (i = b->used; i < x.used; i++) {
ashleymills 0:714293de3836 593 x.dp[i] = 0;
ashleymills 0:714293de3836 594 }
ashleymills 0:714293de3836 595 fp_clamp(&x);
ashleymills 0:714293de3836 596 fp_copy (&x, d);
ashleymills 0:714293de3836 597 }
ashleymills 0:714293de3836 598
ashleymills 0:714293de3836 599 return FP_OKAY;
ashleymills 0:714293de3836 600 }
ashleymills 0:714293de3836 601
ashleymills 0:714293de3836 602 /* b = a/2 */
ashleymills 0:714293de3836 603 void fp_div_2(fp_int * a, fp_int * b)
ashleymills 0:714293de3836 604 {
ashleymills 0:714293de3836 605 int x, oldused;
ashleymills 0:714293de3836 606
ashleymills 0:714293de3836 607 oldused = b->used;
ashleymills 0:714293de3836 608 b->used = a->used;
ashleymills 0:714293de3836 609 {
ashleymills 0:714293de3836 610 register fp_digit r, rr, *tmpa, *tmpb;
ashleymills 0:714293de3836 611
ashleymills 0:714293de3836 612 /* source alias */
ashleymills 0:714293de3836 613 tmpa = a->dp + b->used - 1;
ashleymills 0:714293de3836 614
ashleymills 0:714293de3836 615 /* dest alias */
ashleymills 0:714293de3836 616 tmpb = b->dp + b->used - 1;
ashleymills 0:714293de3836 617
ashleymills 0:714293de3836 618 /* carry */
ashleymills 0:714293de3836 619 r = 0;
ashleymills 0:714293de3836 620 for (x = b->used - 1; x >= 0; x--) {
ashleymills 0:714293de3836 621 /* get the carry for the next iteration */
ashleymills 0:714293de3836 622 rr = *tmpa & 1;
ashleymills 0:714293de3836 623
ashleymills 0:714293de3836 624 /* shift the current digit, add in carry and store */
ashleymills 0:714293de3836 625 *tmpb-- = (*tmpa-- >> 1) | (r << (DIGIT_BIT - 1));
ashleymills 0:714293de3836 626
ashleymills 0:714293de3836 627 /* forward carry to next iteration */
ashleymills 0:714293de3836 628 r = rr;
ashleymills 0:714293de3836 629 }
ashleymills 0:714293de3836 630
ashleymills 0:714293de3836 631 /* zero excess digits */
ashleymills 0:714293de3836 632 tmpb = b->dp + b->used;
ashleymills 0:714293de3836 633 for (x = b->used; x < oldused; x++) {
ashleymills 0:714293de3836 634 *tmpb++ = 0;
ashleymills 0:714293de3836 635 }
ashleymills 0:714293de3836 636 }
ashleymills 0:714293de3836 637 b->sign = a->sign;
ashleymills 0:714293de3836 638 fp_clamp (b);
ashleymills 0:714293de3836 639 }
ashleymills 0:714293de3836 640
ashleymills 0:714293de3836 641 /* c = a / 2**b */
ashleymills 0:714293de3836 642 void fp_div_2d(fp_int *a, int b, fp_int *c, fp_int *d)
ashleymills 0:714293de3836 643 {
ashleymills 0:714293de3836 644 fp_digit D, r, rr;
ashleymills 0:714293de3836 645 int x;
ashleymills 0:714293de3836 646 fp_int t;
ashleymills 0:714293de3836 647
ashleymills 0:714293de3836 648 /* if the shift count is <= 0 then we do no work */
ashleymills 0:714293de3836 649 if (b <= 0) {
ashleymills 0:714293de3836 650 fp_copy (a, c);
ashleymills 0:714293de3836 651 if (d != NULL) {
ashleymills 0:714293de3836 652 fp_zero (d);
ashleymills 0:714293de3836 653 }
ashleymills 0:714293de3836 654 return;
ashleymills 0:714293de3836 655 }
ashleymills 0:714293de3836 656
ashleymills 0:714293de3836 657 fp_init(&t);
ashleymills 0:714293de3836 658
ashleymills 0:714293de3836 659 /* get the remainder */
ashleymills 0:714293de3836 660 if (d != NULL) {
ashleymills 0:714293de3836 661 fp_mod_2d (a, b, &t);
ashleymills 0:714293de3836 662 }
ashleymills 0:714293de3836 663
ashleymills 0:714293de3836 664 /* copy */
ashleymills 0:714293de3836 665 fp_copy(a, c);
ashleymills 0:714293de3836 666
ashleymills 0:714293de3836 667 /* shift by as many digits in the bit count */
ashleymills 0:714293de3836 668 if (b >= (int)DIGIT_BIT) {
ashleymills 0:714293de3836 669 fp_rshd (c, b / DIGIT_BIT);
ashleymills 0:714293de3836 670 }
ashleymills 0:714293de3836 671
ashleymills 0:714293de3836 672 /* shift any bit count < DIGIT_BIT */
ashleymills 0:714293de3836 673 D = (fp_digit) (b % DIGIT_BIT);
ashleymills 0:714293de3836 674 if (D != 0) {
ashleymills 0:714293de3836 675 register fp_digit *tmpc, mask, shift;
ashleymills 0:714293de3836 676
ashleymills 0:714293de3836 677 /* mask */
ashleymills 0:714293de3836 678 mask = (((fp_digit)1) << D) - 1;
ashleymills 0:714293de3836 679
ashleymills 0:714293de3836 680 /* shift for lsb */
ashleymills 0:714293de3836 681 shift = DIGIT_BIT - D;
ashleymills 0:714293de3836 682
ashleymills 0:714293de3836 683 /* alias */
ashleymills 0:714293de3836 684 tmpc = c->dp + (c->used - 1);
ashleymills 0:714293de3836 685
ashleymills 0:714293de3836 686 /* carry */
ashleymills 0:714293de3836 687 r = 0;
ashleymills 0:714293de3836 688 for (x = c->used - 1; x >= 0; x--) {
ashleymills 0:714293de3836 689 /* get the lower bits of this word in a temp */
ashleymills 0:714293de3836 690 rr = *tmpc & mask;
ashleymills 0:714293de3836 691
ashleymills 0:714293de3836 692 /* shift the current word and mix in the carry bits from the previous word */
ashleymills 0:714293de3836 693 *tmpc = (*tmpc >> D) | (r << shift);
ashleymills 0:714293de3836 694 --tmpc;
ashleymills 0:714293de3836 695
ashleymills 0:714293de3836 696 /* set the carry to the carry bits of the current word found above */
ashleymills 0:714293de3836 697 r = rr;
ashleymills 0:714293de3836 698 }
ashleymills 0:714293de3836 699 }
ashleymills 0:714293de3836 700 fp_clamp (c);
ashleymills 0:714293de3836 701 if (d != NULL) {
ashleymills 0:714293de3836 702 fp_copy (&t, d);
ashleymills 0:714293de3836 703 }
ashleymills 0:714293de3836 704 }
ashleymills 0:714293de3836 705
ashleymills 0:714293de3836 706 /* c = a mod b, 0 <= c < b */
ashleymills 0:714293de3836 707 int fp_mod(fp_int *a, fp_int *b, fp_int *c)
ashleymills 0:714293de3836 708 {
ashleymills 0:714293de3836 709 fp_int t;
ashleymills 0:714293de3836 710 int err;
ashleymills 0:714293de3836 711
ashleymills 0:714293de3836 712 fp_zero(&t);
ashleymills 0:714293de3836 713 if ((err = fp_div(a, b, NULL, &t)) != FP_OKAY) {
ashleymills 0:714293de3836 714 return err;
ashleymills 0:714293de3836 715 }
ashleymills 0:714293de3836 716 if (t.sign != b->sign) {
ashleymills 0:714293de3836 717 fp_add(&t, b, c);
ashleymills 0:714293de3836 718 } else {
ashleymills 0:714293de3836 719 fp_copy(&t, c);
ashleymills 0:714293de3836 720 }
ashleymills 0:714293de3836 721 return FP_OKAY;
ashleymills 0:714293de3836 722 }
ashleymills 0:714293de3836 723
ashleymills 0:714293de3836 724 /* c = a mod 2**d */
ashleymills 0:714293de3836 725 void fp_mod_2d(fp_int *a, int b, fp_int *c)
ashleymills 0:714293de3836 726 {
ashleymills 0:714293de3836 727 int x;
ashleymills 0:714293de3836 728
ashleymills 0:714293de3836 729 /* zero if count less than or equal to zero */
ashleymills 0:714293de3836 730 if (b <= 0) {
ashleymills 0:714293de3836 731 fp_zero(c);
ashleymills 0:714293de3836 732 return;
ashleymills 0:714293de3836 733 }
ashleymills 0:714293de3836 734
ashleymills 0:714293de3836 735 /* get copy of input */
ashleymills 0:714293de3836 736 fp_copy(a, c);
ashleymills 0:714293de3836 737
ashleymills 0:714293de3836 738 /* if 2**d is larger than we just return */
ashleymills 0:714293de3836 739 if (b >= (DIGIT_BIT * a->used)) {
ashleymills 0:714293de3836 740 return;
ashleymills 0:714293de3836 741 }
ashleymills 0:714293de3836 742
ashleymills 0:714293de3836 743 /* zero digits above the last digit of the modulus */
ashleymills 0:714293de3836 744 for (x = (b / DIGIT_BIT) + ((b % DIGIT_BIT) == 0 ? 0 : 1); x < c->used; x++) {
ashleymills 0:714293de3836 745 c->dp[x] = 0;
ashleymills 0:714293de3836 746 }
ashleymills 0:714293de3836 747 /* clear the digit that is not completely outside/inside the modulus */
ashleymills 0:714293de3836 748 c->dp[b / DIGIT_BIT] &= ~((fp_digit)0) >> (DIGIT_BIT - b);
ashleymills 0:714293de3836 749 fp_clamp (c);
ashleymills 0:714293de3836 750 }
ashleymills 0:714293de3836 751
ashleymills 0:714293de3836 752 static int fp_invmod_slow (fp_int * a, fp_int * b, fp_int * c)
ashleymills 0:714293de3836 753 {
ashleymills 0:714293de3836 754 fp_int x, y, u, v, A, B, C, D;
ashleymills 0:714293de3836 755 int res;
ashleymills 0:714293de3836 756
ashleymills 0:714293de3836 757 /* b cannot be negative */
ashleymills 0:714293de3836 758 if (b->sign == FP_NEG || fp_iszero(b) == 1) {
ashleymills 0:714293de3836 759 return FP_VAL;
ashleymills 0:714293de3836 760 }
ashleymills 0:714293de3836 761
ashleymills 0:714293de3836 762 /* init temps */
ashleymills 0:714293de3836 763 fp_init(&x); fp_init(&y);
ashleymills 0:714293de3836 764 fp_init(&u); fp_init(&v);
ashleymills 0:714293de3836 765 fp_init(&A); fp_init(&B);
ashleymills 0:714293de3836 766 fp_init(&C); fp_init(&D);
ashleymills 0:714293de3836 767
ashleymills 0:714293de3836 768 /* x = a, y = b */
ashleymills 0:714293de3836 769 if ((res = fp_mod(a, b, &x)) != FP_OKAY) {
ashleymills 0:714293de3836 770 return res;
ashleymills 0:714293de3836 771 }
ashleymills 0:714293de3836 772 fp_copy(b, &y);
ashleymills 0:714293de3836 773
ashleymills 0:714293de3836 774 /* 2. [modified] if x,y are both even then return an error! */
ashleymills 0:714293de3836 775 if (fp_iseven (&x) == 1 && fp_iseven (&y) == 1) {
ashleymills 0:714293de3836 776 return FP_VAL;
ashleymills 0:714293de3836 777 }
ashleymills 0:714293de3836 778
ashleymills 0:714293de3836 779 /* 3. u=x, v=y, A=1, B=0, C=0,D=1 */
ashleymills 0:714293de3836 780 fp_copy (&x, &u);
ashleymills 0:714293de3836 781 fp_copy (&y, &v);
ashleymills 0:714293de3836 782 fp_set (&A, 1);
ashleymills 0:714293de3836 783 fp_set (&D, 1);
ashleymills 0:714293de3836 784
ashleymills 0:714293de3836 785 top:
ashleymills 0:714293de3836 786 /* 4. while u is even do */
ashleymills 0:714293de3836 787 while (fp_iseven (&u) == 1) {
ashleymills 0:714293de3836 788 /* 4.1 u = u/2 */
ashleymills 0:714293de3836 789 fp_div_2 (&u, &u);
ashleymills 0:714293de3836 790
ashleymills 0:714293de3836 791 /* 4.2 if A or B is odd then */
ashleymills 0:714293de3836 792 if (fp_isodd (&A) == 1 || fp_isodd (&B) == 1) {
ashleymills 0:714293de3836 793 /* A = (A+y)/2, B = (B-x)/2 */
ashleymills 0:714293de3836 794 fp_add (&A, &y, &A);
ashleymills 0:714293de3836 795 fp_sub (&B, &x, &B);
ashleymills 0:714293de3836 796 }
ashleymills 0:714293de3836 797 /* A = A/2, B = B/2 */
ashleymills 0:714293de3836 798 fp_div_2 (&A, &A);
ashleymills 0:714293de3836 799 fp_div_2 (&B, &B);
ashleymills 0:714293de3836 800 }
ashleymills 0:714293de3836 801
ashleymills 0:714293de3836 802 /* 5. while v is even do */
ashleymills 0:714293de3836 803 while (fp_iseven (&v) == 1) {
ashleymills 0:714293de3836 804 /* 5.1 v = v/2 */
ashleymills 0:714293de3836 805 fp_div_2 (&v, &v);
ashleymills 0:714293de3836 806
ashleymills 0:714293de3836 807 /* 5.2 if C or D is odd then */
ashleymills 0:714293de3836 808 if (fp_isodd (&C) == 1 || fp_isodd (&D) == 1) {
ashleymills 0:714293de3836 809 /* C = (C+y)/2, D = (D-x)/2 */
ashleymills 0:714293de3836 810 fp_add (&C, &y, &C);
ashleymills 0:714293de3836 811 fp_sub (&D, &x, &D);
ashleymills 0:714293de3836 812 }
ashleymills 0:714293de3836 813 /* C = C/2, D = D/2 */
ashleymills 0:714293de3836 814 fp_div_2 (&C, &C);
ashleymills 0:714293de3836 815 fp_div_2 (&D, &D);
ashleymills 0:714293de3836 816 }
ashleymills 0:714293de3836 817
ashleymills 0:714293de3836 818 /* 6. if u >= v then */
ashleymills 0:714293de3836 819 if (fp_cmp (&u, &v) != FP_LT) {
ashleymills 0:714293de3836 820 /* u = u - v, A = A - C, B = B - D */
ashleymills 0:714293de3836 821 fp_sub (&u, &v, &u);
ashleymills 0:714293de3836 822 fp_sub (&A, &C, &A);
ashleymills 0:714293de3836 823 fp_sub (&B, &D, &B);
ashleymills 0:714293de3836 824 } else {
ashleymills 0:714293de3836 825 /* v - v - u, C = C - A, D = D - B */
ashleymills 0:714293de3836 826 fp_sub (&v, &u, &v);
ashleymills 0:714293de3836 827 fp_sub (&C, &A, &C);
ashleymills 0:714293de3836 828 fp_sub (&D, &B, &D);
ashleymills 0:714293de3836 829 }
ashleymills 0:714293de3836 830
ashleymills 0:714293de3836 831 /* if not zero goto step 4 */
ashleymills 0:714293de3836 832 if (fp_iszero (&u) == 0)
ashleymills 0:714293de3836 833 goto top;
ashleymills 0:714293de3836 834
ashleymills 0:714293de3836 835 /* now a = C, b = D, gcd == g*v */
ashleymills 0:714293de3836 836
ashleymills 0:714293de3836 837 /* if v != 1 then there is no inverse */
ashleymills 0:714293de3836 838 if (fp_cmp_d (&v, 1) != FP_EQ) {
ashleymills 0:714293de3836 839 return FP_VAL;
ashleymills 0:714293de3836 840 }
ashleymills 0:714293de3836 841
ashleymills 0:714293de3836 842 /* if its too low */
ashleymills 0:714293de3836 843 while (fp_cmp_d(&C, 0) == FP_LT) {
ashleymills 0:714293de3836 844 fp_add(&C, b, &C);
ashleymills 0:714293de3836 845 }
ashleymills 0:714293de3836 846
ashleymills 0:714293de3836 847 /* too big */
ashleymills 0:714293de3836 848 while (fp_cmp_mag(&C, b) != FP_LT) {
ashleymills 0:714293de3836 849 fp_sub(&C, b, &C);
ashleymills 0:714293de3836 850 }
ashleymills 0:714293de3836 851
ashleymills 0:714293de3836 852 /* C is now the inverse */
ashleymills 0:714293de3836 853 fp_copy(&C, c);
ashleymills 0:714293de3836 854 return FP_OKAY;
ashleymills 0:714293de3836 855 }
ashleymills 0:714293de3836 856
ashleymills 0:714293de3836 857 /* c = 1/a (mod b) for odd b only */
ashleymills 0:714293de3836 858 int fp_invmod(fp_int *a, fp_int *b, fp_int *c)
ashleymills 0:714293de3836 859 {
ashleymills 0:714293de3836 860 fp_int x, y, u, v, B, D;
ashleymills 0:714293de3836 861 int neg;
ashleymills 0:714293de3836 862
ashleymills 0:714293de3836 863 /* 2. [modified] b must be odd */
ashleymills 0:714293de3836 864 if (fp_iseven (b) == FP_YES) {
ashleymills 0:714293de3836 865 return fp_invmod_slow(a,b,c);
ashleymills 0:714293de3836 866 }
ashleymills 0:714293de3836 867
ashleymills 0:714293de3836 868 /* init all our temps */
ashleymills 0:714293de3836 869 fp_init(&x); fp_init(&y);
ashleymills 0:714293de3836 870 fp_init(&u); fp_init(&v);
ashleymills 0:714293de3836 871 fp_init(&B); fp_init(&D);
ashleymills 0:714293de3836 872
ashleymills 0:714293de3836 873 /* x == modulus, y == value to invert */
ashleymills 0:714293de3836 874 fp_copy(b, &x);
ashleymills 0:714293de3836 875
ashleymills 0:714293de3836 876 /* we need y = |a| */
ashleymills 0:714293de3836 877 fp_abs(a, &y);
ashleymills 0:714293de3836 878
ashleymills 0:714293de3836 879 /* 3. u=x, v=y, A=1, B=0, C=0,D=1 */
ashleymills 0:714293de3836 880 fp_copy(&x, &u);
ashleymills 0:714293de3836 881 fp_copy(&y, &v);
ashleymills 0:714293de3836 882 fp_set (&D, 1);
ashleymills 0:714293de3836 883
ashleymills 0:714293de3836 884 top:
ashleymills 0:714293de3836 885 /* 4. while u is even do */
ashleymills 0:714293de3836 886 while (fp_iseven (&u) == FP_YES) {
ashleymills 0:714293de3836 887 /* 4.1 u = u/2 */
ashleymills 0:714293de3836 888 fp_div_2 (&u, &u);
ashleymills 0:714293de3836 889
ashleymills 0:714293de3836 890 /* 4.2 if B is odd then */
ashleymills 0:714293de3836 891 if (fp_isodd (&B) == FP_YES) {
ashleymills 0:714293de3836 892 fp_sub (&B, &x, &B);
ashleymills 0:714293de3836 893 }
ashleymills 0:714293de3836 894 /* B = B/2 */
ashleymills 0:714293de3836 895 fp_div_2 (&B, &B);
ashleymills 0:714293de3836 896 }
ashleymills 0:714293de3836 897
ashleymills 0:714293de3836 898 /* 5. while v is even do */
ashleymills 0:714293de3836 899 while (fp_iseven (&v) == FP_YES) {
ashleymills 0:714293de3836 900 /* 5.1 v = v/2 */
ashleymills 0:714293de3836 901 fp_div_2 (&v, &v);
ashleymills 0:714293de3836 902
ashleymills 0:714293de3836 903 /* 5.2 if D is odd then */
ashleymills 0:714293de3836 904 if (fp_isodd (&D) == FP_YES) {
ashleymills 0:714293de3836 905 /* D = (D-x)/2 */
ashleymills 0:714293de3836 906 fp_sub (&D, &x, &D);
ashleymills 0:714293de3836 907 }
ashleymills 0:714293de3836 908 /* D = D/2 */
ashleymills 0:714293de3836 909 fp_div_2 (&D, &D);
ashleymills 0:714293de3836 910 }
ashleymills 0:714293de3836 911
ashleymills 0:714293de3836 912 /* 6. if u >= v then */
ashleymills 0:714293de3836 913 if (fp_cmp (&u, &v) != FP_LT) {
ashleymills 0:714293de3836 914 /* u = u - v, B = B - D */
ashleymills 0:714293de3836 915 fp_sub (&u, &v, &u);
ashleymills 0:714293de3836 916 fp_sub (&B, &D, &B);
ashleymills 0:714293de3836 917 } else {
ashleymills 0:714293de3836 918 /* v - v - u, D = D - B */
ashleymills 0:714293de3836 919 fp_sub (&v, &u, &v);
ashleymills 0:714293de3836 920 fp_sub (&D, &B, &D);
ashleymills 0:714293de3836 921 }
ashleymills 0:714293de3836 922
ashleymills 0:714293de3836 923 /* if not zero goto step 4 */
ashleymills 0:714293de3836 924 if (fp_iszero (&u) == FP_NO) {
ashleymills 0:714293de3836 925 goto top;
ashleymills 0:714293de3836 926 }
ashleymills 0:714293de3836 927
ashleymills 0:714293de3836 928 /* now a = C, b = D, gcd == g*v */
ashleymills 0:714293de3836 929
ashleymills 0:714293de3836 930 /* if v != 1 then there is no inverse */
ashleymills 0:714293de3836 931 if (fp_cmp_d (&v, 1) != FP_EQ) {
ashleymills 0:714293de3836 932 return FP_VAL;
ashleymills 0:714293de3836 933 }
ashleymills 0:714293de3836 934
ashleymills 0:714293de3836 935 /* b is now the inverse */
ashleymills 0:714293de3836 936 neg = a->sign;
ashleymills 0:714293de3836 937 while (D.sign == FP_NEG) {
ashleymills 0:714293de3836 938 fp_add (&D, b, &D);
ashleymills 0:714293de3836 939 }
ashleymills 0:714293de3836 940 fp_copy (&D, c);
ashleymills 0:714293de3836 941 c->sign = neg;
ashleymills 0:714293de3836 942 return FP_OKAY;
ashleymills 0:714293de3836 943 }
ashleymills 0:714293de3836 944
ashleymills 0:714293de3836 945 /* d = a * b (mod c) */
ashleymills 0:714293de3836 946 int fp_mulmod(fp_int *a, fp_int *b, fp_int *c, fp_int *d)
ashleymills 0:714293de3836 947 {
ashleymills 0:714293de3836 948 fp_int tmp;
ashleymills 0:714293de3836 949 fp_zero(&tmp);
ashleymills 0:714293de3836 950 fp_mul(a, b, &tmp);
ashleymills 0:714293de3836 951 return fp_mod(&tmp, c, d);
ashleymills 0:714293de3836 952 }
ashleymills 0:714293de3836 953
ashleymills 0:714293de3836 954 #ifdef TFM_TIMING_RESISTANT
ashleymills 0:714293de3836 955
ashleymills 0:714293de3836 956 /* timing resistant montgomery ladder based exptmod
ashleymills 0:714293de3836 957
ashleymills 0:714293de3836 958 Based on work by Marc Joye, Sung-Ming Yen, "The Montgomery Powering Ladder", Cryptographic Hardware and Embedded Systems, CHES 2002
ashleymills 0:714293de3836 959 */
ashleymills 0:714293de3836 960 static int _fp_exptmod(fp_int * G, fp_int * X, fp_int * P, fp_int * Y)
ashleymills 0:714293de3836 961 {
ashleymills 0:714293de3836 962 fp_int R[2];
ashleymills 0:714293de3836 963 fp_digit buf, mp;
ashleymills 0:714293de3836 964 int err, bitcnt, digidx, y;
ashleymills 0:714293de3836 965
ashleymills 0:714293de3836 966 /* now setup montgomery */
ashleymills 0:714293de3836 967 if ((err = fp_montgomery_setup (P, &mp)) != FP_OKAY) {
ashleymills 0:714293de3836 968 return err;
ashleymills 0:714293de3836 969 }
ashleymills 0:714293de3836 970
ashleymills 0:714293de3836 971 fp_init(&R[0]);
ashleymills 0:714293de3836 972 fp_init(&R[1]);
ashleymills 0:714293de3836 973
ashleymills 0:714293de3836 974 /* now we need R mod m */
ashleymills 0:714293de3836 975 fp_montgomery_calc_normalization (&R[0], P);
ashleymills 0:714293de3836 976
ashleymills 0:714293de3836 977 /* now set R[0][1] to G * R mod m */
ashleymills 0:714293de3836 978 if (fp_cmp_mag(P, G) != FP_GT) {
ashleymills 0:714293de3836 979 /* G > P so we reduce it first */
ashleymills 0:714293de3836 980 fp_mod(G, P, &R[1]);
ashleymills 0:714293de3836 981 } else {
ashleymills 0:714293de3836 982 fp_copy(G, &R[1]);
ashleymills 0:714293de3836 983 }
ashleymills 0:714293de3836 984 fp_mulmod (&R[1], &R[0], P, &R[1]);
ashleymills 0:714293de3836 985
ashleymills 0:714293de3836 986 /* for j = t-1 downto 0 do
ashleymills 0:714293de3836 987 r_!k = R0*R1; r_k = r_k^2
ashleymills 0:714293de3836 988 */
ashleymills 0:714293de3836 989
ashleymills 0:714293de3836 990 /* set initial mode and bit cnt */
ashleymills 0:714293de3836 991 bitcnt = 1;
ashleymills 0:714293de3836 992 buf = 0;
ashleymills 0:714293de3836 993 digidx = X->used - 1;
ashleymills 0:714293de3836 994
ashleymills 0:714293de3836 995 for (;;) {
ashleymills 0:714293de3836 996 /* grab next digit as required */
ashleymills 0:714293de3836 997 if (--bitcnt == 0) {
ashleymills 0:714293de3836 998 /* if digidx == -1 we are out of digits so break */
ashleymills 0:714293de3836 999 if (digidx == -1) {
ashleymills 0:714293de3836 1000 break;
ashleymills 0:714293de3836 1001 }
ashleymills 0:714293de3836 1002 /* read next digit and reset bitcnt */
ashleymills 0:714293de3836 1003 buf = X->dp[digidx--];
ashleymills 0:714293de3836 1004 bitcnt = (int)DIGIT_BIT;
ashleymills 0:714293de3836 1005 }
ashleymills 0:714293de3836 1006
ashleymills 0:714293de3836 1007 /* grab the next msb from the exponent */
ashleymills 0:714293de3836 1008 y = (int)(buf >> (DIGIT_BIT - 1)) & 1;
ashleymills 0:714293de3836 1009 buf <<= (fp_digit)1;
ashleymills 0:714293de3836 1010
ashleymills 0:714293de3836 1011 /* do ops */
ashleymills 0:714293de3836 1012 fp_mul(&R[0], &R[1], &R[y^1]); fp_montgomery_reduce(&R[y^1], P, mp);
ashleymills 0:714293de3836 1013 fp_sqr(&R[y], &R[y]); fp_montgomery_reduce(&R[y], P, mp);
ashleymills 0:714293de3836 1014 }
ashleymills 0:714293de3836 1015
ashleymills 0:714293de3836 1016 fp_montgomery_reduce(&R[0], P, mp);
ashleymills 0:714293de3836 1017 fp_copy(&R[0], Y);
ashleymills 0:714293de3836 1018 return FP_OKAY;
ashleymills 0:714293de3836 1019 }
ashleymills 0:714293de3836 1020
ashleymills 0:714293de3836 1021 #else
ashleymills 0:714293de3836 1022
ashleymills 0:714293de3836 1023 /* y = g**x (mod b)
ashleymills 0:714293de3836 1024 * Some restrictions... x must be positive and < b
ashleymills 0:714293de3836 1025 */
ashleymills 0:714293de3836 1026 static int _fp_exptmod(fp_int * G, fp_int * X, fp_int * P, fp_int * Y)
ashleymills 0:714293de3836 1027 {
ashleymills 0:714293de3836 1028 fp_int M[64], res;
ashleymills 0:714293de3836 1029 fp_digit buf, mp;
ashleymills 0:714293de3836 1030 int err, bitbuf, bitcpy, bitcnt, mode, digidx, x, y, winsize;
ashleymills 0:714293de3836 1031
ashleymills 0:714293de3836 1032 /* find window size */
ashleymills 0:714293de3836 1033 x = fp_count_bits (X);
ashleymills 0:714293de3836 1034 if (x <= 21) {
ashleymills 0:714293de3836 1035 winsize = 1;
ashleymills 0:714293de3836 1036 } else if (x <= 36) {
ashleymills 0:714293de3836 1037 winsize = 3;
ashleymills 0:714293de3836 1038 } else if (x <= 140) {
ashleymills 0:714293de3836 1039 winsize = 4;
ashleymills 0:714293de3836 1040 } else if (x <= 450) {
ashleymills 0:714293de3836 1041 winsize = 5;
ashleymills 0:714293de3836 1042 } else {
ashleymills 0:714293de3836 1043 winsize = 6;
ashleymills 0:714293de3836 1044 }
ashleymills 0:714293de3836 1045
ashleymills 0:714293de3836 1046 /* init M array */
ashleymills 0:714293de3836 1047 XMEMSET(M, 0, sizeof(M));
ashleymills 0:714293de3836 1048
ashleymills 0:714293de3836 1049 /* now setup montgomery */
ashleymills 0:714293de3836 1050 if ((err = fp_montgomery_setup (P, &mp)) != FP_OKAY) {
ashleymills 0:714293de3836 1051 return err;
ashleymills 0:714293de3836 1052 }
ashleymills 0:714293de3836 1053
ashleymills 0:714293de3836 1054 /* setup result */
ashleymills 0:714293de3836 1055 fp_init(&res);
ashleymills 0:714293de3836 1056
ashleymills 0:714293de3836 1057 /* create M table
ashleymills 0:714293de3836 1058 *
ashleymills 0:714293de3836 1059 * The M table contains powers of the input base, e.g. M[x] = G^x mod P
ashleymills 0:714293de3836 1060 *
ashleymills 0:714293de3836 1061 * The first half of the table is not computed though accept for M[0] and M[1]
ashleymills 0:714293de3836 1062 */
ashleymills 0:714293de3836 1063
ashleymills 0:714293de3836 1064 /* now we need R mod m */
ashleymills 0:714293de3836 1065 fp_montgomery_calc_normalization (&res, P);
ashleymills 0:714293de3836 1066
ashleymills 0:714293de3836 1067 /* now set M[1] to G * R mod m */
ashleymills 0:714293de3836 1068 if (fp_cmp_mag(P, G) != FP_GT) {
ashleymills 0:714293de3836 1069 /* G > P so we reduce it first */
ashleymills 0:714293de3836 1070 fp_mod(G, P, &M[1]);
ashleymills 0:714293de3836 1071 } else {
ashleymills 0:714293de3836 1072 fp_copy(G, &M[1]);
ashleymills 0:714293de3836 1073 }
ashleymills 0:714293de3836 1074 fp_mulmod (&M[1], &res, P, &M[1]);
ashleymills 0:714293de3836 1075
ashleymills 0:714293de3836 1076 /* compute the value at M[1<<(winsize-1)] by squaring M[1] (winsize-1) times */
ashleymills 0:714293de3836 1077 fp_copy (&M[1], &M[1 << (winsize - 1)]);
ashleymills 0:714293de3836 1078 for (x = 0; x < (winsize - 1); x++) {
ashleymills 0:714293de3836 1079 fp_sqr (&M[1 << (winsize - 1)], &M[1 << (winsize - 1)]);
ashleymills 0:714293de3836 1080 fp_montgomery_reduce (&M[1 << (winsize - 1)], P, mp);
ashleymills 0:714293de3836 1081 }
ashleymills 0:714293de3836 1082
ashleymills 0:714293de3836 1083 /* create upper table */
ashleymills 0:714293de3836 1084 for (x = (1 << (winsize - 1)) + 1; x < (1 << winsize); x++) {
ashleymills 0:714293de3836 1085 fp_mul(&M[x - 1], &M[1], &M[x]);
ashleymills 0:714293de3836 1086 fp_montgomery_reduce(&M[x], P, mp);
ashleymills 0:714293de3836 1087 }
ashleymills 0:714293de3836 1088
ashleymills 0:714293de3836 1089 /* set initial mode and bit cnt */
ashleymills 0:714293de3836 1090 mode = 0;
ashleymills 0:714293de3836 1091 bitcnt = 1;
ashleymills 0:714293de3836 1092 buf = 0;
ashleymills 0:714293de3836 1093 digidx = X->used - 1;
ashleymills 0:714293de3836 1094 bitcpy = 0;
ashleymills 0:714293de3836 1095 bitbuf = 0;
ashleymills 0:714293de3836 1096
ashleymills 0:714293de3836 1097 for (;;) {
ashleymills 0:714293de3836 1098 /* grab next digit as required */
ashleymills 0:714293de3836 1099 if (--bitcnt == 0) {
ashleymills 0:714293de3836 1100 /* if digidx == -1 we are out of digits so break */
ashleymills 0:714293de3836 1101 if (digidx == -1) {
ashleymills 0:714293de3836 1102 break;
ashleymills 0:714293de3836 1103 }
ashleymills 0:714293de3836 1104 /* read next digit and reset bitcnt */
ashleymills 0:714293de3836 1105 buf = X->dp[digidx--];
ashleymills 0:714293de3836 1106 bitcnt = (int)DIGIT_BIT;
ashleymills 0:714293de3836 1107 }
ashleymills 0:714293de3836 1108
ashleymills 0:714293de3836 1109 /* grab the next msb from the exponent */
ashleymills 0:714293de3836 1110 y = (int)(buf >> (DIGIT_BIT - 1)) & 1;
ashleymills 0:714293de3836 1111 buf <<= (fp_digit)1;
ashleymills 0:714293de3836 1112
ashleymills 0:714293de3836 1113 /* if the bit is zero and mode == 0 then we ignore it
ashleymills 0:714293de3836 1114 * These represent the leading zero bits before the first 1 bit
ashleymills 0:714293de3836 1115 * in the exponent. Technically this opt is not required but it
ashleymills 0:714293de3836 1116 * does lower the # of trivial squaring/reductions used
ashleymills 0:714293de3836 1117 */
ashleymills 0:714293de3836 1118 if (mode == 0 && y == 0) {
ashleymills 0:714293de3836 1119 continue;
ashleymills 0:714293de3836 1120 }
ashleymills 0:714293de3836 1121
ashleymills 0:714293de3836 1122 /* if the bit is zero and mode == 1 then we square */
ashleymills 0:714293de3836 1123 if (mode == 1 && y == 0) {
ashleymills 0:714293de3836 1124 fp_sqr(&res, &res);
ashleymills 0:714293de3836 1125 fp_montgomery_reduce(&res, P, mp);
ashleymills 0:714293de3836 1126 continue;
ashleymills 0:714293de3836 1127 }
ashleymills 0:714293de3836 1128
ashleymills 0:714293de3836 1129 /* else we add it to the window */
ashleymills 0:714293de3836 1130 bitbuf |= (y << (winsize - ++bitcpy));
ashleymills 0:714293de3836 1131 mode = 2;
ashleymills 0:714293de3836 1132
ashleymills 0:714293de3836 1133 if (bitcpy == winsize) {
ashleymills 0:714293de3836 1134 /* ok window is filled so square as required and multiply */
ashleymills 0:714293de3836 1135 /* square first */
ashleymills 0:714293de3836 1136 for (x = 0; x < winsize; x++) {
ashleymills 0:714293de3836 1137 fp_sqr(&res, &res);
ashleymills 0:714293de3836 1138 fp_montgomery_reduce(&res, P, mp);
ashleymills 0:714293de3836 1139 }
ashleymills 0:714293de3836 1140
ashleymills 0:714293de3836 1141 /* then multiply */
ashleymills 0:714293de3836 1142 fp_mul(&res, &M[bitbuf], &res);
ashleymills 0:714293de3836 1143 fp_montgomery_reduce(&res, P, mp);
ashleymills 0:714293de3836 1144
ashleymills 0:714293de3836 1145 /* empty window and reset */
ashleymills 0:714293de3836 1146 bitcpy = 0;
ashleymills 0:714293de3836 1147 bitbuf = 0;
ashleymills 0:714293de3836 1148 mode = 1;
ashleymills 0:714293de3836 1149 }
ashleymills 0:714293de3836 1150 }
ashleymills 0:714293de3836 1151
ashleymills 0:714293de3836 1152 /* if bits remain then square/multiply */
ashleymills 0:714293de3836 1153 if (mode == 2 && bitcpy > 0) {
ashleymills 0:714293de3836 1154 /* square then multiply if the bit is set */
ashleymills 0:714293de3836 1155 for (x = 0; x < bitcpy; x++) {
ashleymills 0:714293de3836 1156 fp_sqr(&res, &res);
ashleymills 0:714293de3836 1157 fp_montgomery_reduce(&res, P, mp);
ashleymills 0:714293de3836 1158
ashleymills 0:714293de3836 1159 /* get next bit of the window */
ashleymills 0:714293de3836 1160 bitbuf <<= 1;
ashleymills 0:714293de3836 1161 if ((bitbuf & (1 << winsize)) != 0) {
ashleymills 0:714293de3836 1162 /* then multiply */
ashleymills 0:714293de3836 1163 fp_mul(&res, &M[1], &res);
ashleymills 0:714293de3836 1164 fp_montgomery_reduce(&res, P, mp);
ashleymills 0:714293de3836 1165 }
ashleymills 0:714293de3836 1166 }
ashleymills 0:714293de3836 1167 }
ashleymills 0:714293de3836 1168
ashleymills 0:714293de3836 1169 /* fixup result if Montgomery reduction is used
ashleymills 0:714293de3836 1170 * recall that any value in a Montgomery system is
ashleymills 0:714293de3836 1171 * actually multiplied by R mod n. So we have
ashleymills 0:714293de3836 1172 * to reduce one more time to cancel out the factor
ashleymills 0:714293de3836 1173 * of R.
ashleymills 0:714293de3836 1174 */
ashleymills 0:714293de3836 1175 fp_montgomery_reduce(&res, P, mp);
ashleymills 0:714293de3836 1176
ashleymills 0:714293de3836 1177 /* swap res with Y */
ashleymills 0:714293de3836 1178 fp_copy (&res, Y);
ashleymills 0:714293de3836 1179 return FP_OKAY;
ashleymills 0:714293de3836 1180 }
ashleymills 0:714293de3836 1181
ashleymills 0:714293de3836 1182 #endif
ashleymills 0:714293de3836 1183
ashleymills 0:714293de3836 1184 int fp_exptmod(fp_int * G, fp_int * X, fp_int * P, fp_int * Y)
ashleymills 0:714293de3836 1185 {
ashleymills 0:714293de3836 1186 /* prevent overflows */
ashleymills 0:714293de3836 1187 if (P->used > (FP_SIZE/2)) {
ashleymills 0:714293de3836 1188 return FP_VAL;
ashleymills 0:714293de3836 1189 }
ashleymills 0:714293de3836 1190
ashleymills 0:714293de3836 1191 if (X->sign == FP_NEG) {
ashleymills 0:714293de3836 1192 #ifndef POSITIVE_EXP_ONLY /* reduce stack if assume no negatives */
ashleymills 0:714293de3836 1193 int err;
ashleymills 0:714293de3836 1194 fp_int tmp;
ashleymills 0:714293de3836 1195
ashleymills 0:714293de3836 1196 /* yes, copy G and invmod it */
ashleymills 0:714293de3836 1197 fp_copy(G, &tmp);
ashleymills 0:714293de3836 1198 if ((err = fp_invmod(&tmp, P, &tmp)) != FP_OKAY) {
ashleymills 0:714293de3836 1199 return err;
ashleymills 0:714293de3836 1200 }
ashleymills 0:714293de3836 1201 X->sign = FP_ZPOS;
ashleymills 0:714293de3836 1202 err = _fp_exptmod(&tmp, X, P, Y);
ashleymills 0:714293de3836 1203 if (X != Y) {
ashleymills 0:714293de3836 1204 X->sign = FP_NEG;
ashleymills 0:714293de3836 1205 }
ashleymills 0:714293de3836 1206 return err;
ashleymills 0:714293de3836 1207 #else
ashleymills 0:714293de3836 1208 return FP_VAL;
ashleymills 0:714293de3836 1209 #endif
ashleymills 0:714293de3836 1210 }
ashleymills 0:714293de3836 1211 else {
ashleymills 0:714293de3836 1212 /* Positive exponent so just exptmod */
ashleymills 0:714293de3836 1213 return _fp_exptmod(G, X, P, Y);
ashleymills 0:714293de3836 1214 }
ashleymills 0:714293de3836 1215 }
ashleymills 0:714293de3836 1216
ashleymills 0:714293de3836 1217 /* computes a = 2**b */
ashleymills 0:714293de3836 1218 void fp_2expt(fp_int *a, int b)
ashleymills 0:714293de3836 1219 {
ashleymills 0:714293de3836 1220 int z;
ashleymills 0:714293de3836 1221
ashleymills 0:714293de3836 1222 /* zero a as per default */
ashleymills 0:714293de3836 1223 fp_zero (a);
ashleymills 0:714293de3836 1224
ashleymills 0:714293de3836 1225 if (b < 0) {
ashleymills 0:714293de3836 1226 return;
ashleymills 0:714293de3836 1227 }
ashleymills 0:714293de3836 1228
ashleymills 0:714293de3836 1229 z = b / DIGIT_BIT;
ashleymills 0:714293de3836 1230 if (z >= FP_SIZE) {
ashleymills 0:714293de3836 1231 return;
ashleymills 0:714293de3836 1232 }
ashleymills 0:714293de3836 1233
ashleymills 0:714293de3836 1234 /* set the used count of where the bit will go */
ashleymills 0:714293de3836 1235 a->used = z + 1;
ashleymills 0:714293de3836 1236
ashleymills 0:714293de3836 1237 /* put the single bit in its place */
ashleymills 0:714293de3836 1238 a->dp[z] = ((fp_digit)1) << (b % DIGIT_BIT);
ashleymills 0:714293de3836 1239 }
ashleymills 0:714293de3836 1240
ashleymills 0:714293de3836 1241 /* b = a*a */
ashleymills 0:714293de3836 1242 void fp_sqr(fp_int *A, fp_int *B)
ashleymills 0:714293de3836 1243 {
ashleymills 0:714293de3836 1244 int y = A->used;
ashleymills 0:714293de3836 1245
ashleymills 0:714293de3836 1246 /* call generic if we're out of range */
ashleymills 0:714293de3836 1247 if (y + y > FP_SIZE) {
ashleymills 0:714293de3836 1248 fp_sqr_comba(A, B);
ashleymills 0:714293de3836 1249 return ;
ashleymills 0:714293de3836 1250 }
ashleymills 0:714293de3836 1251
ashleymills 0:714293de3836 1252 #if defined(TFM_SQR3)
ashleymills 0:714293de3836 1253 if (y <= 3) {
ashleymills 0:714293de3836 1254 fp_sqr_comba3(A,B);
ashleymills 0:714293de3836 1255 return;
ashleymills 0:714293de3836 1256 }
ashleymills 0:714293de3836 1257 #endif
ashleymills 0:714293de3836 1258 #if defined(TFM_SQR4)
ashleymills 0:714293de3836 1259 if (y == 4) {
ashleymills 0:714293de3836 1260 fp_sqr_comba4(A,B);
ashleymills 0:714293de3836 1261 return;
ashleymills 0:714293de3836 1262 }
ashleymills 0:714293de3836 1263 #endif
ashleymills 0:714293de3836 1264 #if defined(TFM_SQR6)
ashleymills 0:714293de3836 1265 if (y <= 6) {
ashleymills 0:714293de3836 1266 fp_sqr_comba6(A,B);
ashleymills 0:714293de3836 1267 return;
ashleymills 0:714293de3836 1268 }
ashleymills 0:714293de3836 1269 #endif
ashleymills 0:714293de3836 1270 #if defined(TFM_SQR7)
ashleymills 0:714293de3836 1271 if (y == 7) {
ashleymills 0:714293de3836 1272 fp_sqr_comba7(A,B);
ashleymills 0:714293de3836 1273 return;
ashleymills 0:714293de3836 1274 }
ashleymills 0:714293de3836 1275 #endif
ashleymills 0:714293de3836 1276 #if defined(TFM_SQR8)
ashleymills 0:714293de3836 1277 if (y == 8) {
ashleymills 0:714293de3836 1278 fp_sqr_comba8(A,B);
ashleymills 0:714293de3836 1279 return;
ashleymills 0:714293de3836 1280 }
ashleymills 0:714293de3836 1281 #endif
ashleymills 0:714293de3836 1282 #if defined(TFM_SQR9)
ashleymills 0:714293de3836 1283 if (y == 9) {
ashleymills 0:714293de3836 1284 fp_sqr_comba9(A,B);
ashleymills 0:714293de3836 1285 return;
ashleymills 0:714293de3836 1286 }
ashleymills 0:714293de3836 1287 #endif
ashleymills 0:714293de3836 1288 #if defined(TFM_SQR12)
ashleymills 0:714293de3836 1289 if (y <= 12) {
ashleymills 0:714293de3836 1290 fp_sqr_comba12(A,B);
ashleymills 0:714293de3836 1291 return;
ashleymills 0:714293de3836 1292 }
ashleymills 0:714293de3836 1293 #endif
ashleymills 0:714293de3836 1294 #if defined(TFM_SQR17)
ashleymills 0:714293de3836 1295 if (y <= 17) {
ashleymills 0:714293de3836 1296 fp_sqr_comba17(A,B);
ashleymills 0:714293de3836 1297 return;
ashleymills 0:714293de3836 1298 }
ashleymills 0:714293de3836 1299 #endif
ashleymills 0:714293de3836 1300 #if defined(TFM_SMALL_SET)
ashleymills 0:714293de3836 1301 if (y <= 16) {
ashleymills 0:714293de3836 1302 fp_sqr_comba_small(A,B);
ashleymills 0:714293de3836 1303 return;
ashleymills 0:714293de3836 1304 }
ashleymills 0:714293de3836 1305 #endif
ashleymills 0:714293de3836 1306 #if defined(TFM_SQR20)
ashleymills 0:714293de3836 1307 if (y <= 20) {
ashleymills 0:714293de3836 1308 fp_sqr_comba20(A,B);
ashleymills 0:714293de3836 1309 return;
ashleymills 0:714293de3836 1310 }
ashleymills 0:714293de3836 1311 #endif
ashleymills 0:714293de3836 1312 #if defined(TFM_SQR24)
ashleymills 0:714293de3836 1313 if (y <= 24) {
ashleymills 0:714293de3836 1314 fp_sqr_comba24(A,B);
ashleymills 0:714293de3836 1315 return;
ashleymills 0:714293de3836 1316 }
ashleymills 0:714293de3836 1317 #endif
ashleymills 0:714293de3836 1318 #if defined(TFM_SQR28)
ashleymills 0:714293de3836 1319 if (y <= 28) {
ashleymills 0:714293de3836 1320 fp_sqr_comba28(A,B);
ashleymills 0:714293de3836 1321 return;
ashleymills 0:714293de3836 1322 }
ashleymills 0:714293de3836 1323 #endif
ashleymills 0:714293de3836 1324 #if defined(TFM_SQR32)
ashleymills 0:714293de3836 1325 if (y <= 32) {
ashleymills 0:714293de3836 1326 fp_sqr_comba32(A,B);
ashleymills 0:714293de3836 1327 return;
ashleymills 0:714293de3836 1328 }
ashleymills 0:714293de3836 1329 #endif
ashleymills 0:714293de3836 1330 #if defined(TFM_SQR48)
ashleymills 0:714293de3836 1331 if (y <= 48) {
ashleymills 0:714293de3836 1332 fp_sqr_comba48(A,B);
ashleymills 0:714293de3836 1333 return;
ashleymills 0:714293de3836 1334 }
ashleymills 0:714293de3836 1335 #endif
ashleymills 0:714293de3836 1336 #if defined(TFM_SQR64)
ashleymills 0:714293de3836 1337 if (y <= 64) {
ashleymills 0:714293de3836 1338 fp_sqr_comba64(A,B);
ashleymills 0:714293de3836 1339 return;
ashleymills 0:714293de3836 1340 }
ashleymills 0:714293de3836 1341 #endif
ashleymills 0:714293de3836 1342 fp_sqr_comba(A, B);
ashleymills 0:714293de3836 1343 }
ashleymills 0:714293de3836 1344
ashleymills 0:714293de3836 1345 /* generic comba squarer */
ashleymills 0:714293de3836 1346 void fp_sqr_comba(fp_int *A, fp_int *B)
ashleymills 0:714293de3836 1347 {
ashleymills 0:714293de3836 1348 int pa, ix, iz;
ashleymills 0:714293de3836 1349 fp_digit c0, c1, c2;
ashleymills 0:714293de3836 1350 fp_int tmp, *dst;
ashleymills 0:714293de3836 1351 #ifdef TFM_ISO
ashleymills 0:714293de3836 1352 fp_word tt;
ashleymills 0:714293de3836 1353 #endif
ashleymills 0:714293de3836 1354
ashleymills 0:714293de3836 1355 /* get size of output and trim */
ashleymills 0:714293de3836 1356 pa = A->used + A->used;
ashleymills 0:714293de3836 1357 if (pa >= FP_SIZE) {
ashleymills 0:714293de3836 1358 pa = FP_SIZE-1;
ashleymills 0:714293de3836 1359 }
ashleymills 0:714293de3836 1360
ashleymills 0:714293de3836 1361 /* number of output digits to produce */
ashleymills 0:714293de3836 1362 COMBA_START;
ashleymills 0:714293de3836 1363 COMBA_CLEAR;
ashleymills 0:714293de3836 1364
ashleymills 0:714293de3836 1365 if (A == B) {
ashleymills 0:714293de3836 1366 fp_zero(&tmp);
ashleymills 0:714293de3836 1367 dst = &tmp;
ashleymills 0:714293de3836 1368 } else {
ashleymills 0:714293de3836 1369 fp_zero(B);
ashleymills 0:714293de3836 1370 dst = B;
ashleymills 0:714293de3836 1371 }
ashleymills 0:714293de3836 1372
ashleymills 0:714293de3836 1373 for (ix = 0; ix < pa; ix++) {
ashleymills 0:714293de3836 1374 int tx, ty, iy;
ashleymills 0:714293de3836 1375 fp_digit *tmpy, *tmpx;
ashleymills 0:714293de3836 1376
ashleymills 0:714293de3836 1377 /* get offsets into the two bignums */
ashleymills 0:714293de3836 1378 ty = MIN(A->used-1, ix);
ashleymills 0:714293de3836 1379 tx = ix - ty;
ashleymills 0:714293de3836 1380
ashleymills 0:714293de3836 1381 /* setup temp aliases */
ashleymills 0:714293de3836 1382 tmpx = A->dp + tx;
ashleymills 0:714293de3836 1383 tmpy = A->dp + ty;
ashleymills 0:714293de3836 1384
ashleymills 0:714293de3836 1385 /* this is the number of times the loop will iterrate,
ashleymills 0:714293de3836 1386 while (tx++ < a->used && ty-- >= 0) { ... }
ashleymills 0:714293de3836 1387 */
ashleymills 0:714293de3836 1388 iy = MIN(A->used-tx, ty+1);
ashleymills 0:714293de3836 1389
ashleymills 0:714293de3836 1390 /* now for squaring tx can never equal ty
ashleymills 0:714293de3836 1391 * we halve the distance since they approach
ashleymills 0:714293de3836 1392 * at a rate of 2x and we have to round because
ashleymills 0:714293de3836 1393 * odd cases need to be executed
ashleymills 0:714293de3836 1394 */
ashleymills 0:714293de3836 1395 iy = MIN(iy, (ty-tx+1)>>1);
ashleymills 0:714293de3836 1396
ashleymills 0:714293de3836 1397 /* forward carries */
ashleymills 0:714293de3836 1398 COMBA_FORWARD;
ashleymills 0:714293de3836 1399
ashleymills 0:714293de3836 1400 /* execute loop */
ashleymills 0:714293de3836 1401 for (iz = 0; iz < iy; iz++) {
ashleymills 0:714293de3836 1402 SQRADD2(*tmpx++, *tmpy--);
ashleymills 0:714293de3836 1403 }
ashleymills 0:714293de3836 1404
ashleymills 0:714293de3836 1405 /* even columns have the square term in them */
ashleymills 0:714293de3836 1406 if ((ix&1) == 0) {
ashleymills 0:714293de3836 1407 /* TAO change COMBA_ADD back to SQRADD */
ashleymills 0:714293de3836 1408 SQRADD(A->dp[ix>>1], A->dp[ix>>1]);
ashleymills 0:714293de3836 1409 }
ashleymills 0:714293de3836 1410
ashleymills 0:714293de3836 1411 /* store it */
ashleymills 0:714293de3836 1412 COMBA_STORE(dst->dp[ix]);
ashleymills 0:714293de3836 1413 }
ashleymills 0:714293de3836 1414
ashleymills 0:714293de3836 1415 COMBA_FINI;
ashleymills 0:714293de3836 1416
ashleymills 0:714293de3836 1417 /* setup dest */
ashleymills 0:714293de3836 1418 dst->used = pa;
ashleymills 0:714293de3836 1419 fp_clamp (dst);
ashleymills 0:714293de3836 1420 if (dst != B) {
ashleymills 0:714293de3836 1421 fp_copy(dst, B);
ashleymills 0:714293de3836 1422 }
ashleymills 0:714293de3836 1423 }
ashleymills 0:714293de3836 1424
ashleymills 0:714293de3836 1425 int fp_cmp(fp_int *a, fp_int *b)
ashleymills 0:714293de3836 1426 {
ashleymills 0:714293de3836 1427 if (a->sign == FP_NEG && b->sign == FP_ZPOS) {
ashleymills 0:714293de3836 1428 return FP_LT;
ashleymills 0:714293de3836 1429 } else if (a->sign == FP_ZPOS && b->sign == FP_NEG) {
ashleymills 0:714293de3836 1430 return FP_GT;
ashleymills 0:714293de3836 1431 } else {
ashleymills 0:714293de3836 1432 /* compare digits */
ashleymills 0:714293de3836 1433 if (a->sign == FP_NEG) {
ashleymills 0:714293de3836 1434 /* if negative compare opposite direction */
ashleymills 0:714293de3836 1435 return fp_cmp_mag(b, a);
ashleymills 0:714293de3836 1436 } else {
ashleymills 0:714293de3836 1437 return fp_cmp_mag(a, b);
ashleymills 0:714293de3836 1438 }
ashleymills 0:714293de3836 1439 }
ashleymills 0:714293de3836 1440 }
ashleymills 0:714293de3836 1441
ashleymills 0:714293de3836 1442 /* compare against a single digit */
ashleymills 0:714293de3836 1443 int fp_cmp_d(fp_int *a, fp_digit b)
ashleymills 0:714293de3836 1444 {
ashleymills 0:714293de3836 1445 /* compare based on sign */
ashleymills 0:714293de3836 1446 if ((b && a->used == 0) || a->sign == FP_NEG) {
ashleymills 0:714293de3836 1447 return FP_LT;
ashleymills 0:714293de3836 1448 }
ashleymills 0:714293de3836 1449
ashleymills 0:714293de3836 1450 /* compare based on magnitude */
ashleymills 0:714293de3836 1451 if (a->used > 1) {
ashleymills 0:714293de3836 1452 return FP_GT;
ashleymills 0:714293de3836 1453 }
ashleymills 0:714293de3836 1454
ashleymills 0:714293de3836 1455 /* compare the only digit of a to b */
ashleymills 0:714293de3836 1456 if (a->dp[0] > b) {
ashleymills 0:714293de3836 1457 return FP_GT;
ashleymills 0:714293de3836 1458 } else if (a->dp[0] < b) {
ashleymills 0:714293de3836 1459 return FP_LT;
ashleymills 0:714293de3836 1460 } else {
ashleymills 0:714293de3836 1461 return FP_EQ;
ashleymills 0:714293de3836 1462 }
ashleymills 0:714293de3836 1463
ashleymills 0:714293de3836 1464 }
ashleymills 0:714293de3836 1465
ashleymills 0:714293de3836 1466 int fp_cmp_mag(fp_int *a, fp_int *b)
ashleymills 0:714293de3836 1467 {
ashleymills 0:714293de3836 1468 int x;
ashleymills 0:714293de3836 1469
ashleymills 0:714293de3836 1470 if (a->used > b->used) {
ashleymills 0:714293de3836 1471 return FP_GT;
ashleymills 0:714293de3836 1472 } else if (a->used < b->used) {
ashleymills 0:714293de3836 1473 return FP_LT;
ashleymills 0:714293de3836 1474 } else {
ashleymills 0:714293de3836 1475 for (x = a->used - 1; x >= 0; x--) {
ashleymills 0:714293de3836 1476 if (a->dp[x] > b->dp[x]) {
ashleymills 0:714293de3836 1477 return FP_GT;
ashleymills 0:714293de3836 1478 } else if (a->dp[x] < b->dp[x]) {
ashleymills 0:714293de3836 1479 return FP_LT;
ashleymills 0:714293de3836 1480 }
ashleymills 0:714293de3836 1481 }
ashleymills 0:714293de3836 1482 }
ashleymills 0:714293de3836 1483 return FP_EQ;
ashleymills 0:714293de3836 1484 }
ashleymills 0:714293de3836 1485
ashleymills 0:714293de3836 1486 /* setups the montgomery reduction */
ashleymills 0:714293de3836 1487 int fp_montgomery_setup(fp_int *a, fp_digit *rho)
ashleymills 0:714293de3836 1488 {
ashleymills 0:714293de3836 1489 fp_digit x, b;
ashleymills 0:714293de3836 1490
ashleymills 0:714293de3836 1491 /* fast inversion mod 2**k
ashleymills 0:714293de3836 1492 *
ashleymills 0:714293de3836 1493 * Based on the fact that
ashleymills 0:714293de3836 1494 *
ashleymills 0:714293de3836 1495 * XA = 1 (mod 2**n) => (X(2-XA)) A = 1 (mod 2**2n)
ashleymills 0:714293de3836 1496 * => 2*X*A - X*X*A*A = 1
ashleymills 0:714293de3836 1497 * => 2*(1) - (1) = 1
ashleymills 0:714293de3836 1498 */
ashleymills 0:714293de3836 1499 b = a->dp[0];
ashleymills 0:714293de3836 1500
ashleymills 0:714293de3836 1501 if ((b & 1) == 0) {
ashleymills 0:714293de3836 1502 return FP_VAL;
ashleymills 0:714293de3836 1503 }
ashleymills 0:714293de3836 1504
ashleymills 0:714293de3836 1505 x = (((b + 2) & 4) << 1) + b; /* here x*a==1 mod 2**4 */
ashleymills 0:714293de3836 1506 x *= 2 - b * x; /* here x*a==1 mod 2**8 */
ashleymills 0:714293de3836 1507 x *= 2 - b * x; /* here x*a==1 mod 2**16 */
ashleymills 0:714293de3836 1508 x *= 2 - b * x; /* here x*a==1 mod 2**32 */
ashleymills 0:714293de3836 1509 #ifdef FP_64BIT
ashleymills 0:714293de3836 1510 x *= 2 - b * x; /* here x*a==1 mod 2**64 */
ashleymills 0:714293de3836 1511 #endif
ashleymills 0:714293de3836 1512
ashleymills 0:714293de3836 1513 /* rho = -1/m mod b */
ashleymills 0:714293de3836 1514 *rho = (fp_digit) (((fp_word) 1 << ((fp_word) DIGIT_BIT)) - ((fp_word)x));
ashleymills 0:714293de3836 1515
ashleymills 0:714293de3836 1516 return FP_OKAY;
ashleymills 0:714293de3836 1517 }
ashleymills 0:714293de3836 1518
ashleymills 0:714293de3836 1519 /* computes a = B**n mod b without division or multiplication useful for
ashleymills 0:714293de3836 1520 * normalizing numbers in a Montgomery system.
ashleymills 0:714293de3836 1521 */
ashleymills 0:714293de3836 1522 void fp_montgomery_calc_normalization(fp_int *a, fp_int *b)
ashleymills 0:714293de3836 1523 {
ashleymills 0:714293de3836 1524 int x, bits;
ashleymills 0:714293de3836 1525
ashleymills 0:714293de3836 1526 /* how many bits of last digit does b use */
ashleymills 0:714293de3836 1527 bits = fp_count_bits (b) % DIGIT_BIT;
ashleymills 0:714293de3836 1528 if (!bits) bits = DIGIT_BIT;
ashleymills 0:714293de3836 1529
ashleymills 0:714293de3836 1530 /* compute A = B^(n-1) * 2^(bits-1) */
ashleymills 0:714293de3836 1531 if (b->used > 1) {
ashleymills 0:714293de3836 1532 fp_2expt (a, (b->used - 1) * DIGIT_BIT + bits - 1);
ashleymills 0:714293de3836 1533 } else {
ashleymills 0:714293de3836 1534 fp_set(a, 1);
ashleymills 0:714293de3836 1535 bits = 1;
ashleymills 0:714293de3836 1536 }
ashleymills 0:714293de3836 1537
ashleymills 0:714293de3836 1538 /* now compute C = A * B mod b */
ashleymills 0:714293de3836 1539 for (x = bits - 1; x < (int)DIGIT_BIT; x++) {
ashleymills 0:714293de3836 1540 fp_mul_2 (a, a);
ashleymills 0:714293de3836 1541 if (fp_cmp_mag (a, b) != FP_LT) {
ashleymills 0:714293de3836 1542 s_fp_sub (a, b, a);
ashleymills 0:714293de3836 1543 }
ashleymills 0:714293de3836 1544 }
ashleymills 0:714293de3836 1545 }
ashleymills 0:714293de3836 1546
ashleymills 0:714293de3836 1547
ashleymills 0:714293de3836 1548 #ifdef TFM_SMALL_MONT_SET
ashleymills 0:714293de3836 1549 #include "fp_mont_small.i"
ashleymills 0:714293de3836 1550 #endif
ashleymills 0:714293de3836 1551
ashleymills 0:714293de3836 1552 /* computes x/R == x (mod N) via Montgomery Reduction */
ashleymills 0:714293de3836 1553 void fp_montgomery_reduce(fp_int *a, fp_int *m, fp_digit mp)
ashleymills 0:714293de3836 1554 {
ashleymills 0:714293de3836 1555 fp_digit c[FP_SIZE], *_c, *tmpm, mu = 0;
ashleymills 0:714293de3836 1556 int oldused, x, y, pa;
ashleymills 0:714293de3836 1557
ashleymills 0:714293de3836 1558 /* bail if too large */
ashleymills 0:714293de3836 1559 if (m->used > (FP_SIZE/2)) {
ashleymills 0:714293de3836 1560 (void)mu; /* shut up compiler */
ashleymills 0:714293de3836 1561 return;
ashleymills 0:714293de3836 1562 }
ashleymills 0:714293de3836 1563
ashleymills 0:714293de3836 1564 #ifdef TFM_SMALL_MONT_SET
ashleymills 0:714293de3836 1565 if (m->used <= 16) {
ashleymills 0:714293de3836 1566 fp_montgomery_reduce_small(a, m, mp);
ashleymills 0:714293de3836 1567 return;
ashleymills 0:714293de3836 1568 }
ashleymills 0:714293de3836 1569 #endif
ashleymills 0:714293de3836 1570
ashleymills 0:714293de3836 1571
ashleymills 0:714293de3836 1572 /* now zero the buff */
ashleymills 0:714293de3836 1573 XMEMSET(c, 0, sizeof c);
ashleymills 0:714293de3836 1574 pa = m->used;
ashleymills 0:714293de3836 1575
ashleymills 0:714293de3836 1576 /* copy the input */
ashleymills 0:714293de3836 1577 oldused = a->used;
ashleymills 0:714293de3836 1578 for (x = 0; x < oldused; x++) {
ashleymills 0:714293de3836 1579 c[x] = a->dp[x];
ashleymills 0:714293de3836 1580 }
ashleymills 0:714293de3836 1581 MONT_START;
ashleymills 0:714293de3836 1582
ashleymills 0:714293de3836 1583 for (x = 0; x < pa; x++) {
ashleymills 0:714293de3836 1584 fp_digit cy = 0;
ashleymills 0:714293de3836 1585 /* get Mu for this round */
ashleymills 0:714293de3836 1586 LOOP_START;
ashleymills 0:714293de3836 1587 _c = c + x;
ashleymills 0:714293de3836 1588 tmpm = m->dp;
ashleymills 0:714293de3836 1589 y = 0;
ashleymills 0:714293de3836 1590 #if (defined(TFM_SSE2) || defined(TFM_X86_64))
ashleymills 0:714293de3836 1591 for (; y < (pa & ~7); y += 8) {
ashleymills 0:714293de3836 1592 INNERMUL8;
ashleymills 0:714293de3836 1593 _c += 8;
ashleymills 0:714293de3836 1594 tmpm += 8;
ashleymills 0:714293de3836 1595 }
ashleymills 0:714293de3836 1596 #endif
ashleymills 0:714293de3836 1597
ashleymills 0:714293de3836 1598 for (; y < pa; y++) {
ashleymills 0:714293de3836 1599 INNERMUL;
ashleymills 0:714293de3836 1600 ++_c;
ashleymills 0:714293de3836 1601 }
ashleymills 0:714293de3836 1602 LOOP_END;
ashleymills 0:714293de3836 1603 while (cy) {
ashleymills 0:714293de3836 1604 PROPCARRY;
ashleymills 0:714293de3836 1605 ++_c;
ashleymills 0:714293de3836 1606 }
ashleymills 0:714293de3836 1607 }
ashleymills 0:714293de3836 1608
ashleymills 0:714293de3836 1609 /* now copy out */
ashleymills 0:714293de3836 1610 _c = c + pa;
ashleymills 0:714293de3836 1611 tmpm = a->dp;
ashleymills 0:714293de3836 1612 for (x = 0; x < pa+1; x++) {
ashleymills 0:714293de3836 1613 *tmpm++ = *_c++;
ashleymills 0:714293de3836 1614 }
ashleymills 0:714293de3836 1615
ashleymills 0:714293de3836 1616 for (; x < oldused; x++) {
ashleymills 0:714293de3836 1617 *tmpm++ = 0;
ashleymills 0:714293de3836 1618 }
ashleymills 0:714293de3836 1619
ashleymills 0:714293de3836 1620 MONT_FINI;
ashleymills 0:714293de3836 1621
ashleymills 0:714293de3836 1622 a->used = pa+1;
ashleymills 0:714293de3836 1623 fp_clamp(a);
ashleymills 0:714293de3836 1624
ashleymills 0:714293de3836 1625 /* if A >= m then A = A - m */
ashleymills 0:714293de3836 1626 if (fp_cmp_mag (a, m) != FP_LT) {
ashleymills 0:714293de3836 1627 s_fp_sub (a, m, a);
ashleymills 0:714293de3836 1628 }
ashleymills 0:714293de3836 1629 }
ashleymills 0:714293de3836 1630
ashleymills 0:714293de3836 1631 void fp_read_unsigned_bin(fp_int *a, unsigned char *b, int c)
ashleymills 0:714293de3836 1632 {
ashleymills 0:714293de3836 1633 /* zero the int */
ashleymills 0:714293de3836 1634 fp_zero (a);
ashleymills 0:714293de3836 1635
ashleymills 0:714293de3836 1636 /* If we know the endianness of this architecture, and we're using
ashleymills 0:714293de3836 1637 32-bit fp_digits, we can optimize this */
ashleymills 0:714293de3836 1638 #if (defined(LITTLE_ENDIAN_ORDER) || defined(BIG_ENDIAN_ORDER)) && !defined(FP_64BIT)
ashleymills 0:714293de3836 1639 /* But not for both simultaneously */
ashleymills 0:714293de3836 1640 #if defined(LITTLE_ENDIAN_ORDER) && defined(BIG_ENDIAN_ORDER)
ashleymills 0:714293de3836 1641 #error Both LITTLE_ENDIAN_ORDER and BIG_ENDIAN_ORDER defined.
ashleymills 0:714293de3836 1642 #endif
ashleymills 0:714293de3836 1643 {
ashleymills 0:714293de3836 1644 unsigned char *pd = (unsigned char *)a->dp;
ashleymills 0:714293de3836 1645
ashleymills 0:714293de3836 1646 if ((unsigned)c > (FP_SIZE * sizeof(fp_digit))) {
ashleymills 0:714293de3836 1647 int excess = c - (FP_SIZE * sizeof(fp_digit));
ashleymills 0:714293de3836 1648 c -= excess;
ashleymills 0:714293de3836 1649 b += excess;
ashleymills 0:714293de3836 1650 }
ashleymills 0:714293de3836 1651 a->used = (c + sizeof(fp_digit) - 1)/sizeof(fp_digit);
ashleymills 0:714293de3836 1652 /* read the bytes in */
ashleymills 0:714293de3836 1653 #ifdef BIG_ENDIAN_ORDER
ashleymills 0:714293de3836 1654 {
ashleymills 0:714293de3836 1655 /* Use Duff's device to unroll the loop. */
ashleymills 0:714293de3836 1656 int idx = (c - 1) & ~3;
ashleymills 0:714293de3836 1657 switch (c % 4) {
ashleymills 0:714293de3836 1658 case 0: do { pd[idx+0] = *b++;
ashleymills 0:714293de3836 1659 case 3: pd[idx+1] = *b++;
ashleymills 0:714293de3836 1660 case 2: pd[idx+2] = *b++;
ashleymills 0:714293de3836 1661 case 1: pd[idx+3] = *b++;
ashleymills 0:714293de3836 1662 idx -= 4;
ashleymills 0:714293de3836 1663 } while ((c -= 4) > 0);
ashleymills 0:714293de3836 1664 }
ashleymills 0:714293de3836 1665 }
ashleymills 0:714293de3836 1666 #else
ashleymills 0:714293de3836 1667 for (c -= 1; c >= 0; c -= 1) {
ashleymills 0:714293de3836 1668 pd[c] = *b++;
ashleymills 0:714293de3836 1669 }
ashleymills 0:714293de3836 1670 #endif
ashleymills 0:714293de3836 1671 }
ashleymills 0:714293de3836 1672 #else
ashleymills 0:714293de3836 1673 /* read the bytes in */
ashleymills 0:714293de3836 1674 for (; c > 0; c--) {
ashleymills 0:714293de3836 1675 fp_mul_2d (a, 8, a);
ashleymills 0:714293de3836 1676 a->dp[0] |= *b++;
ashleymills 0:714293de3836 1677 a->used += 1;
ashleymills 0:714293de3836 1678 }
ashleymills 0:714293de3836 1679 #endif
ashleymills 0:714293de3836 1680 fp_clamp (a);
ashleymills 0:714293de3836 1681 }
ashleymills 0:714293de3836 1682
ashleymills 0:714293de3836 1683 void fp_to_unsigned_bin(fp_int *a, unsigned char *b)
ashleymills 0:714293de3836 1684 {
ashleymills 0:714293de3836 1685 int x;
ashleymills 0:714293de3836 1686 fp_int t;
ashleymills 0:714293de3836 1687
ashleymills 0:714293de3836 1688 fp_init_copy(&t, a);
ashleymills 0:714293de3836 1689
ashleymills 0:714293de3836 1690 x = 0;
ashleymills 0:714293de3836 1691 while (fp_iszero (&t) == FP_NO) {
ashleymills 0:714293de3836 1692 b[x++] = (unsigned char) (t.dp[0] & 255);
ashleymills 0:714293de3836 1693 fp_div_2d (&t, 8, &t, NULL);
ashleymills 0:714293de3836 1694 }
ashleymills 0:714293de3836 1695 fp_reverse (b, x);
ashleymills 0:714293de3836 1696 }
ashleymills 0:714293de3836 1697
ashleymills 0:714293de3836 1698 int fp_unsigned_bin_size(fp_int *a)
ashleymills 0:714293de3836 1699 {
ashleymills 0:714293de3836 1700 int size = fp_count_bits (a);
ashleymills 0:714293de3836 1701 return (size / 8 + ((size & 7) != 0 ? 1 : 0));
ashleymills 0:714293de3836 1702 }
ashleymills 0:714293de3836 1703
ashleymills 0:714293de3836 1704 void fp_set(fp_int *a, fp_digit b)
ashleymills 0:714293de3836 1705 {
ashleymills 0:714293de3836 1706 fp_zero(a);
ashleymills 0:714293de3836 1707 a->dp[0] = b;
ashleymills 0:714293de3836 1708 a->used = a->dp[0] ? 1 : 0;
ashleymills 0:714293de3836 1709 }
ashleymills 0:714293de3836 1710
ashleymills 0:714293de3836 1711 int fp_count_bits (fp_int * a)
ashleymills 0:714293de3836 1712 {
ashleymills 0:714293de3836 1713 int r;
ashleymills 0:714293de3836 1714 fp_digit q;
ashleymills 0:714293de3836 1715
ashleymills 0:714293de3836 1716 /* shortcut */
ashleymills 0:714293de3836 1717 if (a->used == 0) {
ashleymills 0:714293de3836 1718 return 0;
ashleymills 0:714293de3836 1719 }
ashleymills 0:714293de3836 1720
ashleymills 0:714293de3836 1721 /* get number of digits and add that */
ashleymills 0:714293de3836 1722 r = (a->used - 1) * DIGIT_BIT;
ashleymills 0:714293de3836 1723
ashleymills 0:714293de3836 1724 /* take the last digit and count the bits in it */
ashleymills 0:714293de3836 1725 q = a->dp[a->used - 1];
ashleymills 0:714293de3836 1726 while (q > ((fp_digit) 0)) {
ashleymills 0:714293de3836 1727 ++r;
ashleymills 0:714293de3836 1728 q >>= ((fp_digit) 1);
ashleymills 0:714293de3836 1729 }
ashleymills 0:714293de3836 1730 return r;
ashleymills 0:714293de3836 1731 }
ashleymills 0:714293de3836 1732
ashleymills 0:714293de3836 1733 void fp_lshd(fp_int *a, int x)
ashleymills 0:714293de3836 1734 {
ashleymills 0:714293de3836 1735 int y;
ashleymills 0:714293de3836 1736
ashleymills 0:714293de3836 1737 /* move up and truncate as required */
ashleymills 0:714293de3836 1738 y = MIN(a->used + x - 1, (int)(FP_SIZE-1));
ashleymills 0:714293de3836 1739
ashleymills 0:714293de3836 1740 /* store new size */
ashleymills 0:714293de3836 1741 a->used = y + 1;
ashleymills 0:714293de3836 1742
ashleymills 0:714293de3836 1743 /* move digits */
ashleymills 0:714293de3836 1744 for (; y >= x; y--) {
ashleymills 0:714293de3836 1745 a->dp[y] = a->dp[y-x];
ashleymills 0:714293de3836 1746 }
ashleymills 0:714293de3836 1747
ashleymills 0:714293de3836 1748 /* zero lower digits */
ashleymills 0:714293de3836 1749 for (; y >= 0; y--) {
ashleymills 0:714293de3836 1750 a->dp[y] = 0;
ashleymills 0:714293de3836 1751 }
ashleymills 0:714293de3836 1752
ashleymills 0:714293de3836 1753 /* clamp digits */
ashleymills 0:714293de3836 1754 fp_clamp(a);
ashleymills 0:714293de3836 1755 }
ashleymills 0:714293de3836 1756
ashleymills 0:714293de3836 1757 void fp_rshd(fp_int *a, int x)
ashleymills 0:714293de3836 1758 {
ashleymills 0:714293de3836 1759 int y;
ashleymills 0:714293de3836 1760
ashleymills 0:714293de3836 1761 /* too many digits just zero and return */
ashleymills 0:714293de3836 1762 if (x >= a->used) {
ashleymills 0:714293de3836 1763 fp_zero(a);
ashleymills 0:714293de3836 1764 return;
ashleymills 0:714293de3836 1765 }
ashleymills 0:714293de3836 1766
ashleymills 0:714293de3836 1767 /* shift */
ashleymills 0:714293de3836 1768 for (y = 0; y < a->used - x; y++) {
ashleymills 0:714293de3836 1769 a->dp[y] = a->dp[y+x];
ashleymills 0:714293de3836 1770 }
ashleymills 0:714293de3836 1771
ashleymills 0:714293de3836 1772 /* zero rest */
ashleymills 0:714293de3836 1773 for (; y < a->used; y++) {
ashleymills 0:714293de3836 1774 a->dp[y] = 0;
ashleymills 0:714293de3836 1775 }
ashleymills 0:714293de3836 1776
ashleymills 0:714293de3836 1777 /* decrement count */
ashleymills 0:714293de3836 1778 a->used -= x;
ashleymills 0:714293de3836 1779 fp_clamp(a);
ashleymills 0:714293de3836 1780 }
ashleymills 0:714293de3836 1781
ashleymills 0:714293de3836 1782 /* reverse an array, used for radix code */
ashleymills 0:714293de3836 1783 void fp_reverse (unsigned char *s, int len)
ashleymills 0:714293de3836 1784 {
ashleymills 0:714293de3836 1785 int ix, iy;
ashleymills 0:714293de3836 1786 unsigned char t;
ashleymills 0:714293de3836 1787
ashleymills 0:714293de3836 1788 ix = 0;
ashleymills 0:714293de3836 1789 iy = len - 1;
ashleymills 0:714293de3836 1790 while (ix < iy) {
ashleymills 0:714293de3836 1791 t = s[ix];
ashleymills 0:714293de3836 1792 s[ix] = s[iy];
ashleymills 0:714293de3836 1793 s[iy] = t;
ashleymills 0:714293de3836 1794 ++ix;
ashleymills 0:714293de3836 1795 --iy;
ashleymills 0:714293de3836 1796 }
ashleymills 0:714293de3836 1797 }
ashleymills 0:714293de3836 1798
ashleymills 0:714293de3836 1799
ashleymills 0:714293de3836 1800 /* c = a - b */
ashleymills 0:714293de3836 1801 void fp_sub_d(fp_int *a, fp_digit b, fp_int *c)
ashleymills 0:714293de3836 1802 {
ashleymills 0:714293de3836 1803 fp_int tmp;
ashleymills 0:714293de3836 1804 fp_set(&tmp, b);
ashleymills 0:714293de3836 1805 fp_sub(a, &tmp, c);
ashleymills 0:714293de3836 1806 }
ashleymills 0:714293de3836 1807
ashleymills 0:714293de3836 1808
ashleymills 0:714293de3836 1809 /* CyaSSL callers from normal lib */
ashleymills 0:714293de3836 1810
ashleymills 0:714293de3836 1811 /* init a new mp_int */
ashleymills 0:714293de3836 1812 int mp_init (mp_int * a)
ashleymills 0:714293de3836 1813 {
ashleymills 0:714293de3836 1814 if (a)
ashleymills 0:714293de3836 1815 fp_init(a);
ashleymills 0:714293de3836 1816 return MP_OKAY;
ashleymills 0:714293de3836 1817 }
ashleymills 0:714293de3836 1818
ashleymills 0:714293de3836 1819 /* clear one (frees) */
ashleymills 0:714293de3836 1820 void mp_clear (mp_int * a)
ashleymills 0:714293de3836 1821 {
ashleymills 0:714293de3836 1822 fp_zero(a);
ashleymills 0:714293de3836 1823 }
ashleymills 0:714293de3836 1824
ashleymills 0:714293de3836 1825 /* handle up to 6 inits */
ashleymills 0:714293de3836 1826 int mp_init_multi(mp_int* a, mp_int* b, mp_int* c, mp_int* d, mp_int* e, mp_int* f)
ashleymills 0:714293de3836 1827 {
ashleymills 0:714293de3836 1828 if (a)
ashleymills 0:714293de3836 1829 fp_init(a);
ashleymills 0:714293de3836 1830 if (b)
ashleymills 0:714293de3836 1831 fp_init(b);
ashleymills 0:714293de3836 1832 if (c)
ashleymills 0:714293de3836 1833 fp_init(c);
ashleymills 0:714293de3836 1834 if (d)
ashleymills 0:714293de3836 1835 fp_init(d);
ashleymills 0:714293de3836 1836 if (e)
ashleymills 0:714293de3836 1837 fp_init(e);
ashleymills 0:714293de3836 1838 if (f)
ashleymills 0:714293de3836 1839 fp_init(f);
ashleymills 0:714293de3836 1840
ashleymills 0:714293de3836 1841 return MP_OKAY;
ashleymills 0:714293de3836 1842 }
ashleymills 0:714293de3836 1843
ashleymills 0:714293de3836 1844 /* high level addition (handles signs) */
ashleymills 0:714293de3836 1845 int mp_add (mp_int * a, mp_int * b, mp_int * c)
ashleymills 0:714293de3836 1846 {
ashleymills 0:714293de3836 1847 fp_add(a, b, c);
ashleymills 0:714293de3836 1848 return MP_OKAY;
ashleymills 0:714293de3836 1849 }
ashleymills 0:714293de3836 1850
ashleymills 0:714293de3836 1851 /* high level subtraction (handles signs) */
ashleymills 0:714293de3836 1852 int mp_sub (mp_int * a, mp_int * b, mp_int * c)
ashleymills 0:714293de3836 1853 {
ashleymills 0:714293de3836 1854 fp_sub(a, b, c);
ashleymills 0:714293de3836 1855 return MP_OKAY;
ashleymills 0:714293de3836 1856 }
ashleymills 0:714293de3836 1857
ashleymills 0:714293de3836 1858 /* high level multiplication (handles sign) */
ashleymills 0:714293de3836 1859 int mp_mul (mp_int * a, mp_int * b, mp_int * c)
ashleymills 0:714293de3836 1860 {
ashleymills 0:714293de3836 1861 fp_mul(a, b, c);
ashleymills 0:714293de3836 1862 return MP_OKAY;
ashleymills 0:714293de3836 1863 }
ashleymills 0:714293de3836 1864
ashleymills 0:714293de3836 1865 /* d = a * b (mod c) */
ashleymills 0:714293de3836 1866 int mp_mulmod (mp_int * a, mp_int * b, mp_int * c, mp_int * d)
ashleymills 0:714293de3836 1867 {
ashleymills 0:714293de3836 1868 return fp_mulmod(a, b, c, d);
ashleymills 0:714293de3836 1869 }
ashleymills 0:714293de3836 1870
ashleymills 0:714293de3836 1871 /* c = a mod b, 0 <= c < b */
ashleymills 0:714293de3836 1872 int mp_mod (mp_int * a, mp_int * b, mp_int * c)
ashleymills 0:714293de3836 1873 {
ashleymills 0:714293de3836 1874 return fp_mod (a, b, c);
ashleymills 0:714293de3836 1875 }
ashleymills 0:714293de3836 1876
ashleymills 0:714293de3836 1877 /* hac 14.61, pp608 */
ashleymills 0:714293de3836 1878 int mp_invmod (mp_int * a, mp_int * b, mp_int * c)
ashleymills 0:714293de3836 1879 {
ashleymills 0:714293de3836 1880 return fp_invmod(a, b, c);
ashleymills 0:714293de3836 1881 }
ashleymills 0:714293de3836 1882
ashleymills 0:714293de3836 1883 /* this is a shell function that calls either the normal or Montgomery
ashleymills 0:714293de3836 1884 * exptmod functions. Originally the call to the montgomery code was
ashleymills 0:714293de3836 1885 * embedded in the normal function but that wasted alot of stack space
ashleymills 0:714293de3836 1886 * for nothing (since 99% of the time the Montgomery code would be called)
ashleymills 0:714293de3836 1887 */
ashleymills 0:714293de3836 1888 int mp_exptmod (mp_int * G, mp_int * X, mp_int * P, mp_int * Y)
ashleymills 0:714293de3836 1889 {
ashleymills 0:714293de3836 1890 return fp_exptmod(G, X, P, Y);
ashleymills 0:714293de3836 1891 }
ashleymills 0:714293de3836 1892
ashleymills 0:714293de3836 1893 /* compare two ints (signed)*/
ashleymills 0:714293de3836 1894 int mp_cmp (mp_int * a, mp_int * b)
ashleymills 0:714293de3836 1895 {
ashleymills 0:714293de3836 1896 return fp_cmp(a, b);
ashleymills 0:714293de3836 1897 }
ashleymills 0:714293de3836 1898
ashleymills 0:714293de3836 1899 /* compare a digit */
ashleymills 0:714293de3836 1900 int mp_cmp_d(mp_int * a, mp_digit b)
ashleymills 0:714293de3836 1901 {
ashleymills 0:714293de3836 1902 return fp_cmp_d(a, b);
ashleymills 0:714293de3836 1903 }
ashleymills 0:714293de3836 1904
ashleymills 0:714293de3836 1905 /* get the size for an unsigned equivalent */
ashleymills 0:714293de3836 1906 int mp_unsigned_bin_size (mp_int * a)
ashleymills 0:714293de3836 1907 {
ashleymills 0:714293de3836 1908 return fp_unsigned_bin_size(a);
ashleymills 0:714293de3836 1909 }
ashleymills 0:714293de3836 1910
ashleymills 0:714293de3836 1911 /* store in unsigned [big endian] format */
ashleymills 0:714293de3836 1912 int mp_to_unsigned_bin (mp_int * a, unsigned char *b)
ashleymills 0:714293de3836 1913 {
ashleymills 0:714293de3836 1914 fp_to_unsigned_bin(a,b);
ashleymills 0:714293de3836 1915 return MP_OKAY;
ashleymills 0:714293de3836 1916 }
ashleymills 0:714293de3836 1917
ashleymills 0:714293de3836 1918 /* reads a unsigned char array, assumes the msb is stored first [big endian] */
ashleymills 0:714293de3836 1919 int mp_read_unsigned_bin (mp_int * a, const unsigned char *b, int c)
ashleymills 0:714293de3836 1920 {
ashleymills 0:714293de3836 1921 fp_read_unsigned_bin(a, (unsigned char *)b, c);
ashleymills 0:714293de3836 1922 return MP_OKAY;
ashleymills 0:714293de3836 1923 }
ashleymills 0:714293de3836 1924
ashleymills 0:714293de3836 1925
ashleymills 0:714293de3836 1926 int mp_sub_d(fp_int *a, fp_digit b, fp_int *c)
ashleymills 0:714293de3836 1927 {
ashleymills 0:714293de3836 1928 fp_sub_d(a, b, c);
ashleymills 0:714293de3836 1929 return MP_OKAY;
ashleymills 0:714293de3836 1930 }
ashleymills 0:714293de3836 1931
ashleymills 0:714293de3836 1932
ashleymills 0:714293de3836 1933 /* fast math conversion */
ashleymills 0:714293de3836 1934 int mp_copy(fp_int* a, fp_int* b)
ashleymills 0:714293de3836 1935 {
ashleymills 0:714293de3836 1936 fp_copy(a, b);
ashleymills 0:714293de3836 1937 return MP_OKAY;
ashleymills 0:714293de3836 1938 }
ashleymills 0:714293de3836 1939
ashleymills 0:714293de3836 1940
ashleymills 0:714293de3836 1941 /* fast math conversion */
ashleymills 0:714293de3836 1942 int mp_isodd(mp_int* a)
ashleymills 0:714293de3836 1943 {
ashleymills 0:714293de3836 1944 return fp_isodd(a);
ashleymills 0:714293de3836 1945 }
ashleymills 0:714293de3836 1946
ashleymills 0:714293de3836 1947
ashleymills 0:714293de3836 1948 /* fast math conversion */
ashleymills 0:714293de3836 1949 int mp_iszero(mp_int* a)
ashleymills 0:714293de3836 1950 {
ashleymills 0:714293de3836 1951 return fp_iszero(a);
ashleymills 0:714293de3836 1952 }
ashleymills 0:714293de3836 1953
ashleymills 0:714293de3836 1954
ashleymills 0:714293de3836 1955 /* fast math conversion */
ashleymills 0:714293de3836 1956 int mp_count_bits (mp_int* a)
ashleymills 0:714293de3836 1957 {
ashleymills 0:714293de3836 1958 return fp_count_bits(a);
ashleymills 0:714293de3836 1959 }
ashleymills 0:714293de3836 1960
ashleymills 0:714293de3836 1961
ashleymills 0:714293de3836 1962 /* fast math wrappers */
ashleymills 0:714293de3836 1963 int mp_set_int(fp_int *a, fp_digit b)
ashleymills 0:714293de3836 1964 {
ashleymills 0:714293de3836 1965 fp_set(a, b);
ashleymills 0:714293de3836 1966 return MP_OKAY;
ashleymills 0:714293de3836 1967 }
ashleymills 0:714293de3836 1968
ashleymills 0:714293de3836 1969
ashleymills 0:714293de3836 1970 #if defined(CYASSL_KEY_GEN) || defined (HAVE_ECC)
ashleymills 0:714293de3836 1971
ashleymills 0:714293de3836 1972 /* c = a * a (mod b) */
ashleymills 0:714293de3836 1973 int fp_sqrmod(fp_int *a, fp_int *b, fp_int *c)
ashleymills 0:714293de3836 1974 {
ashleymills 0:714293de3836 1975 fp_int tmp;
ashleymills 0:714293de3836 1976 fp_zero(&tmp);
ashleymills 0:714293de3836 1977 fp_sqr(a, &tmp);
ashleymills 0:714293de3836 1978 return fp_mod(&tmp, b, c);
ashleymills 0:714293de3836 1979 }
ashleymills 0:714293de3836 1980
ashleymills 0:714293de3836 1981 /* fast math conversion */
ashleymills 0:714293de3836 1982 int mp_sqrmod(mp_int *a, mp_int *b, mp_int *c)
ashleymills 0:714293de3836 1983 {
ashleymills 0:714293de3836 1984 return fp_sqrmod(a, b, c);
ashleymills 0:714293de3836 1985 }
ashleymills 0:714293de3836 1986
ashleymills 0:714293de3836 1987 /* fast math conversion */
ashleymills 0:714293de3836 1988 int mp_montgomery_calc_normalization(mp_int *a, mp_int *b)
ashleymills 0:714293de3836 1989 {
ashleymills 0:714293de3836 1990 fp_montgomery_calc_normalization(a, b);
ashleymills 0:714293de3836 1991 return MP_OKAY;
ashleymills 0:714293de3836 1992 }
ashleymills 0:714293de3836 1993
ashleymills 0:714293de3836 1994 #endif /* CYASSL_KEYGEN || HAVE_ECC */
ashleymills 0:714293de3836 1995
ashleymills 0:714293de3836 1996
ashleymills 0:714293de3836 1997 #ifdef CYASSL_KEY_GEN
ashleymills 0:714293de3836 1998
ashleymills 0:714293de3836 1999 void fp_gcd(fp_int *a, fp_int *b, fp_int *c);
ashleymills 0:714293de3836 2000 void fp_lcm(fp_int *a, fp_int *b, fp_int *c);
ashleymills 0:714293de3836 2001 int fp_isprime(fp_int *a);
ashleymills 0:714293de3836 2002 int fp_cnt_lsb(fp_int *a);
ashleymills 0:714293de3836 2003
ashleymills 0:714293de3836 2004 int mp_gcd(fp_int *a, fp_int *b, fp_int *c)
ashleymills 0:714293de3836 2005 {
ashleymills 0:714293de3836 2006 fp_gcd(a, b, c);
ashleymills 0:714293de3836 2007 return MP_OKAY;
ashleymills 0:714293de3836 2008 }
ashleymills 0:714293de3836 2009
ashleymills 0:714293de3836 2010
ashleymills 0:714293de3836 2011 int mp_lcm(fp_int *a, fp_int *b, fp_int *c)
ashleymills 0:714293de3836 2012 {
ashleymills 0:714293de3836 2013 fp_lcm(a, b, c);
ashleymills 0:714293de3836 2014 return MP_OKAY;
ashleymills 0:714293de3836 2015 }
ashleymills 0:714293de3836 2016
ashleymills 0:714293de3836 2017
ashleymills 0:714293de3836 2018 int mp_prime_is_prime(mp_int* a, int t, int* result)
ashleymills 0:714293de3836 2019 {
ashleymills 0:714293de3836 2020 (void)t;
ashleymills 0:714293de3836 2021 *result = fp_isprime(a);
ashleymills 0:714293de3836 2022 return MP_OKAY;
ashleymills 0:714293de3836 2023 }
ashleymills 0:714293de3836 2024
ashleymills 0:714293de3836 2025
ashleymills 0:714293de3836 2026
ashleymills 0:714293de3836 2027 static int s_is_power_of_two(fp_digit b, int *p)
ashleymills 0:714293de3836 2028 {
ashleymills 0:714293de3836 2029 int x;
ashleymills 0:714293de3836 2030
ashleymills 0:714293de3836 2031 /* fast return if no power of two */
ashleymills 0:714293de3836 2032 if ((b==0) || (b & (b-1))) {
ashleymills 0:714293de3836 2033 return 0;
ashleymills 0:714293de3836 2034 }
ashleymills 0:714293de3836 2035
ashleymills 0:714293de3836 2036 for (x = 0; x < DIGIT_BIT; x++) {
ashleymills 0:714293de3836 2037 if (b == (((fp_digit)1)<<x)) {
ashleymills 0:714293de3836 2038 *p = x;
ashleymills 0:714293de3836 2039 return 1;
ashleymills 0:714293de3836 2040 }
ashleymills 0:714293de3836 2041 }
ashleymills 0:714293de3836 2042 return 0;
ashleymills 0:714293de3836 2043 }
ashleymills 0:714293de3836 2044
ashleymills 0:714293de3836 2045 /* a/b => cb + d == a */
ashleymills 0:714293de3836 2046 static int fp_div_d(fp_int *a, fp_digit b, fp_int *c, fp_digit *d)
ashleymills 0:714293de3836 2047 {
ashleymills 0:714293de3836 2048 fp_int q;
ashleymills 0:714293de3836 2049 fp_word w;
ashleymills 0:714293de3836 2050 fp_digit t;
ashleymills 0:714293de3836 2051 int ix;
ashleymills 0:714293de3836 2052
ashleymills 0:714293de3836 2053 /* cannot divide by zero */
ashleymills 0:714293de3836 2054 if (b == 0) {
ashleymills 0:714293de3836 2055 return FP_VAL;
ashleymills 0:714293de3836 2056 }
ashleymills 0:714293de3836 2057
ashleymills 0:714293de3836 2058 /* quick outs */
ashleymills 0:714293de3836 2059 if (b == 1 || fp_iszero(a) == 1) {
ashleymills 0:714293de3836 2060 if (d != NULL) {
ashleymills 0:714293de3836 2061 *d = 0;
ashleymills 0:714293de3836 2062 }
ashleymills 0:714293de3836 2063 if (c != NULL) {
ashleymills 0:714293de3836 2064 fp_copy(a, c);
ashleymills 0:714293de3836 2065 }
ashleymills 0:714293de3836 2066 return FP_OKAY;
ashleymills 0:714293de3836 2067 }
ashleymills 0:714293de3836 2068
ashleymills 0:714293de3836 2069 /* power of two ? */
ashleymills 0:714293de3836 2070 if (s_is_power_of_two(b, &ix) == 1) {
ashleymills 0:714293de3836 2071 if (d != NULL) {
ashleymills 0:714293de3836 2072 *d = a->dp[0] & ((((fp_digit)1)<<ix) - 1);
ashleymills 0:714293de3836 2073 }
ashleymills 0:714293de3836 2074 if (c != NULL) {
ashleymills 0:714293de3836 2075 fp_div_2d(a, ix, c, NULL);
ashleymills 0:714293de3836 2076 }
ashleymills 0:714293de3836 2077 return FP_OKAY;
ashleymills 0:714293de3836 2078 }
ashleymills 0:714293de3836 2079
ashleymills 0:714293de3836 2080 /* no easy answer [c'est la vie]. Just division */
ashleymills 0:714293de3836 2081 fp_init(&q);
ashleymills 0:714293de3836 2082
ashleymills 0:714293de3836 2083 q.used = a->used;
ashleymills 0:714293de3836 2084 q.sign = a->sign;
ashleymills 0:714293de3836 2085 w = 0;
ashleymills 0:714293de3836 2086 for (ix = a->used - 1; ix >= 0; ix--) {
ashleymills 0:714293de3836 2087 w = (w << ((fp_word)DIGIT_BIT)) | ((fp_word)a->dp[ix]);
ashleymills 0:714293de3836 2088
ashleymills 0:714293de3836 2089 if (w >= b) {
ashleymills 0:714293de3836 2090 t = (fp_digit)(w / b);
ashleymills 0:714293de3836 2091 w -= ((fp_word)t) * ((fp_word)b);
ashleymills 0:714293de3836 2092 } else {
ashleymills 0:714293de3836 2093 t = 0;
ashleymills 0:714293de3836 2094 }
ashleymills 0:714293de3836 2095 q.dp[ix] = (fp_digit)t;
ashleymills 0:714293de3836 2096 }
ashleymills 0:714293de3836 2097
ashleymills 0:714293de3836 2098 if (d != NULL) {
ashleymills 0:714293de3836 2099 *d = (fp_digit)w;
ashleymills 0:714293de3836 2100 }
ashleymills 0:714293de3836 2101
ashleymills 0:714293de3836 2102 if (c != NULL) {
ashleymills 0:714293de3836 2103 fp_clamp(&q);
ashleymills 0:714293de3836 2104 fp_copy(&q, c);
ashleymills 0:714293de3836 2105 }
ashleymills 0:714293de3836 2106
ashleymills 0:714293de3836 2107 return FP_OKAY;
ashleymills 0:714293de3836 2108 }
ashleymills 0:714293de3836 2109
ashleymills 0:714293de3836 2110
ashleymills 0:714293de3836 2111 /* c = a mod b, 0 <= c < b */
ashleymills 0:714293de3836 2112 static int fp_mod_d(fp_int *a, fp_digit b, fp_digit *c)
ashleymills 0:714293de3836 2113 {
ashleymills 0:714293de3836 2114 return fp_div_d(a, b, NULL, c);
ashleymills 0:714293de3836 2115 }
ashleymills 0:714293de3836 2116
ashleymills 0:714293de3836 2117
ashleymills 0:714293de3836 2118 /* Miller-Rabin test of "a" to the base of "b" as described in
ashleymills 0:714293de3836 2119 * HAC pp. 139 Algorithm 4.24
ashleymills 0:714293de3836 2120 *
ashleymills 0:714293de3836 2121 * Sets result to 0 if definitely composite or 1 if probably prime.
ashleymills 0:714293de3836 2122 * Randomly the chance of error is no more than 1/4 and often
ashleymills 0:714293de3836 2123 * very much lower.
ashleymills 0:714293de3836 2124 */
ashleymills 0:714293de3836 2125 static void fp_prime_miller_rabin (fp_int * a, fp_int * b, int *result)
ashleymills 0:714293de3836 2126 {
ashleymills 0:714293de3836 2127 fp_int n1, y, r;
ashleymills 0:714293de3836 2128 int s, j;
ashleymills 0:714293de3836 2129
ashleymills 0:714293de3836 2130 /* default */
ashleymills 0:714293de3836 2131 *result = FP_NO;
ashleymills 0:714293de3836 2132
ashleymills 0:714293de3836 2133 /* ensure b > 1 */
ashleymills 0:714293de3836 2134 if (fp_cmp_d(b, 1) != FP_GT) {
ashleymills 0:714293de3836 2135 return;
ashleymills 0:714293de3836 2136 }
ashleymills 0:714293de3836 2137
ashleymills 0:714293de3836 2138 /* get n1 = a - 1 */
ashleymills 0:714293de3836 2139 fp_init_copy(&n1, a);
ashleymills 0:714293de3836 2140 fp_sub_d(&n1, 1, &n1);
ashleymills 0:714293de3836 2141
ashleymills 0:714293de3836 2142 /* set 2**s * r = n1 */
ashleymills 0:714293de3836 2143 fp_init_copy(&r, &n1);
ashleymills 0:714293de3836 2144
ashleymills 0:714293de3836 2145 /* count the number of least significant bits
ashleymills 0:714293de3836 2146 * which are zero
ashleymills 0:714293de3836 2147 */
ashleymills 0:714293de3836 2148 s = fp_cnt_lsb(&r);
ashleymills 0:714293de3836 2149
ashleymills 0:714293de3836 2150 /* now divide n - 1 by 2**s */
ashleymills 0:714293de3836 2151 fp_div_2d (&r, s, &r, NULL);
ashleymills 0:714293de3836 2152
ashleymills 0:714293de3836 2153 /* compute y = b**r mod a */
ashleymills 0:714293de3836 2154 fp_init(&y);
ashleymills 0:714293de3836 2155 fp_exptmod(b, &r, a, &y);
ashleymills 0:714293de3836 2156
ashleymills 0:714293de3836 2157 /* if y != 1 and y != n1 do */
ashleymills 0:714293de3836 2158 if (fp_cmp_d (&y, 1) != FP_EQ && fp_cmp (&y, &n1) != FP_EQ) {
ashleymills 0:714293de3836 2159 j = 1;
ashleymills 0:714293de3836 2160 /* while j <= s-1 and y != n1 */
ashleymills 0:714293de3836 2161 while ((j <= (s - 1)) && fp_cmp (&y, &n1) != FP_EQ) {
ashleymills 0:714293de3836 2162 fp_sqrmod (&y, a, &y);
ashleymills 0:714293de3836 2163
ashleymills 0:714293de3836 2164 /* if y == 1 then composite */
ashleymills 0:714293de3836 2165 if (fp_cmp_d (&y, 1) == FP_EQ) {
ashleymills 0:714293de3836 2166 return;
ashleymills 0:714293de3836 2167 }
ashleymills 0:714293de3836 2168 ++j;
ashleymills 0:714293de3836 2169 }
ashleymills 0:714293de3836 2170
ashleymills 0:714293de3836 2171 /* if y != n1 then composite */
ashleymills 0:714293de3836 2172 if (fp_cmp (&y, &n1) != FP_EQ) {
ashleymills 0:714293de3836 2173 return;
ashleymills 0:714293de3836 2174 }
ashleymills 0:714293de3836 2175 }
ashleymills 0:714293de3836 2176
ashleymills 0:714293de3836 2177 /* probably prime now */
ashleymills 0:714293de3836 2178 *result = FP_YES;
ashleymills 0:714293de3836 2179 }
ashleymills 0:714293de3836 2180
ashleymills 0:714293de3836 2181
ashleymills 0:714293de3836 2182 /* a few primes */
ashleymills 0:714293de3836 2183 static const fp_digit primes[256] = {
ashleymills 0:714293de3836 2184 0x0002, 0x0003, 0x0005, 0x0007, 0x000B, 0x000D, 0x0011, 0x0013,
ashleymills 0:714293de3836 2185 0x0017, 0x001D, 0x001F, 0x0025, 0x0029, 0x002B, 0x002F, 0x0035,
ashleymills 0:714293de3836 2186 0x003B, 0x003D, 0x0043, 0x0047, 0x0049, 0x004F, 0x0053, 0x0059,
ashleymills 0:714293de3836 2187 0x0061, 0x0065, 0x0067, 0x006B, 0x006D, 0x0071, 0x007F, 0x0083,
ashleymills 0:714293de3836 2188 0x0089, 0x008B, 0x0095, 0x0097, 0x009D, 0x00A3, 0x00A7, 0x00AD,
ashleymills 0:714293de3836 2189 0x00B3, 0x00B5, 0x00BF, 0x00C1, 0x00C5, 0x00C7, 0x00D3, 0x00DF,
ashleymills 0:714293de3836 2190 0x00E3, 0x00E5, 0x00E9, 0x00EF, 0x00F1, 0x00FB, 0x0101, 0x0107,
ashleymills 0:714293de3836 2191 0x010D, 0x010F, 0x0115, 0x0119, 0x011B, 0x0125, 0x0133, 0x0137,
ashleymills 0:714293de3836 2192
ashleymills 0:714293de3836 2193 0x0139, 0x013D, 0x014B, 0x0151, 0x015B, 0x015D, 0x0161, 0x0167,
ashleymills 0:714293de3836 2194 0x016F, 0x0175, 0x017B, 0x017F, 0x0185, 0x018D, 0x0191, 0x0199,
ashleymills 0:714293de3836 2195 0x01A3, 0x01A5, 0x01AF, 0x01B1, 0x01B7, 0x01BB, 0x01C1, 0x01C9,
ashleymills 0:714293de3836 2196 0x01CD, 0x01CF, 0x01D3, 0x01DF, 0x01E7, 0x01EB, 0x01F3, 0x01F7,
ashleymills 0:714293de3836 2197 0x01FD, 0x0209, 0x020B, 0x021D, 0x0223, 0x022D, 0x0233, 0x0239,
ashleymills 0:714293de3836 2198 0x023B, 0x0241, 0x024B, 0x0251, 0x0257, 0x0259, 0x025F, 0x0265,
ashleymills 0:714293de3836 2199 0x0269, 0x026B, 0x0277, 0x0281, 0x0283, 0x0287, 0x028D, 0x0293,
ashleymills 0:714293de3836 2200 0x0295, 0x02A1, 0x02A5, 0x02AB, 0x02B3, 0x02BD, 0x02C5, 0x02CF,
ashleymills 0:714293de3836 2201
ashleymills 0:714293de3836 2202 0x02D7, 0x02DD, 0x02E3, 0x02E7, 0x02EF, 0x02F5, 0x02F9, 0x0301,
ashleymills 0:714293de3836 2203 0x0305, 0x0313, 0x031D, 0x0329, 0x032B, 0x0335, 0x0337, 0x033B,
ashleymills 0:714293de3836 2204 0x033D, 0x0347, 0x0355, 0x0359, 0x035B, 0x035F, 0x036D, 0x0371,
ashleymills 0:714293de3836 2205 0x0373, 0x0377, 0x038B, 0x038F, 0x0397, 0x03A1, 0x03A9, 0x03AD,
ashleymills 0:714293de3836 2206 0x03B3, 0x03B9, 0x03C7, 0x03CB, 0x03D1, 0x03D7, 0x03DF, 0x03E5,
ashleymills 0:714293de3836 2207 0x03F1, 0x03F5, 0x03FB, 0x03FD, 0x0407, 0x0409, 0x040F, 0x0419,
ashleymills 0:714293de3836 2208 0x041B, 0x0425, 0x0427, 0x042D, 0x043F, 0x0443, 0x0445, 0x0449,
ashleymills 0:714293de3836 2209 0x044F, 0x0455, 0x045D, 0x0463, 0x0469, 0x047F, 0x0481, 0x048B,
ashleymills 0:714293de3836 2210
ashleymills 0:714293de3836 2211 0x0493, 0x049D, 0x04A3, 0x04A9, 0x04B1, 0x04BD, 0x04C1, 0x04C7,
ashleymills 0:714293de3836 2212 0x04CD, 0x04CF, 0x04D5, 0x04E1, 0x04EB, 0x04FD, 0x04FF, 0x0503,
ashleymills 0:714293de3836 2213 0x0509, 0x050B, 0x0511, 0x0515, 0x0517, 0x051B, 0x0527, 0x0529,
ashleymills 0:714293de3836 2214 0x052F, 0x0551, 0x0557, 0x055D, 0x0565, 0x0577, 0x0581, 0x058F,
ashleymills 0:714293de3836 2215 0x0593, 0x0595, 0x0599, 0x059F, 0x05A7, 0x05AB, 0x05AD, 0x05B3,
ashleymills 0:714293de3836 2216 0x05BF, 0x05C9, 0x05CB, 0x05CF, 0x05D1, 0x05D5, 0x05DB, 0x05E7,
ashleymills 0:714293de3836 2217 0x05F3, 0x05FB, 0x0607, 0x060D, 0x0611, 0x0617, 0x061F, 0x0623,
ashleymills 0:714293de3836 2218 0x062B, 0x062F, 0x063D, 0x0641, 0x0647, 0x0649, 0x064D, 0x0653
ashleymills 0:714293de3836 2219 };
ashleymills 0:714293de3836 2220
ashleymills 0:714293de3836 2221 int fp_isprime(fp_int *a)
ashleymills 0:714293de3836 2222 {
ashleymills 0:714293de3836 2223 fp_int b;
ashleymills 0:714293de3836 2224 fp_digit d = 0;
ashleymills 0:714293de3836 2225 int r, res;
ashleymills 0:714293de3836 2226
ashleymills 0:714293de3836 2227 /* do trial division */
ashleymills 0:714293de3836 2228 for (r = 0; r < 256; r++) {
ashleymills 0:714293de3836 2229 fp_mod_d(a, primes[r], &d);
ashleymills 0:714293de3836 2230 if (d == 0) {
ashleymills 0:714293de3836 2231 return FP_NO;
ashleymills 0:714293de3836 2232 }
ashleymills 0:714293de3836 2233 }
ashleymills 0:714293de3836 2234
ashleymills 0:714293de3836 2235 /* now do 8 miller rabins */
ashleymills 0:714293de3836 2236 fp_init(&b);
ashleymills 0:714293de3836 2237 for (r = 0; r < 8; r++) {
ashleymills 0:714293de3836 2238 fp_set(&b, primes[r]);
ashleymills 0:714293de3836 2239 fp_prime_miller_rabin(a, &b, &res);
ashleymills 0:714293de3836 2240 if (res == FP_NO) {
ashleymills 0:714293de3836 2241 return FP_NO;
ashleymills 0:714293de3836 2242 }
ashleymills 0:714293de3836 2243 }
ashleymills 0:714293de3836 2244 return FP_YES;
ashleymills 0:714293de3836 2245 }
ashleymills 0:714293de3836 2246
ashleymills 0:714293de3836 2247
ashleymills 0:714293de3836 2248 /* c = [a, b] */
ashleymills 0:714293de3836 2249 void fp_lcm(fp_int *a, fp_int *b, fp_int *c)
ashleymills 0:714293de3836 2250 {
ashleymills 0:714293de3836 2251 fp_int t1, t2;
ashleymills 0:714293de3836 2252
ashleymills 0:714293de3836 2253 fp_init(&t1);
ashleymills 0:714293de3836 2254 fp_init(&t2);
ashleymills 0:714293de3836 2255 fp_gcd(a, b, &t1);
ashleymills 0:714293de3836 2256 if (fp_cmp_mag(a, b) == FP_GT) {
ashleymills 0:714293de3836 2257 fp_div(a, &t1, &t2, NULL);
ashleymills 0:714293de3836 2258 fp_mul(b, &t2, c);
ashleymills 0:714293de3836 2259 } else {
ashleymills 0:714293de3836 2260 fp_div(b, &t1, &t2, NULL);
ashleymills 0:714293de3836 2261 fp_mul(a, &t2, c);
ashleymills 0:714293de3836 2262 }
ashleymills 0:714293de3836 2263 }
ashleymills 0:714293de3836 2264
ashleymills 0:714293de3836 2265
ashleymills 0:714293de3836 2266 static const int lnz[16] = {
ashleymills 0:714293de3836 2267 4, 0, 1, 0, 2, 0, 1, 0, 3, 0, 1, 0, 2, 0, 1, 0
ashleymills 0:714293de3836 2268 };
ashleymills 0:714293de3836 2269
ashleymills 0:714293de3836 2270 /* Counts the number of lsbs which are zero before the first zero bit */
ashleymills 0:714293de3836 2271 int fp_cnt_lsb(fp_int *a)
ashleymills 0:714293de3836 2272 {
ashleymills 0:714293de3836 2273 int x;
ashleymills 0:714293de3836 2274 fp_digit q, qq;
ashleymills 0:714293de3836 2275
ashleymills 0:714293de3836 2276 /* easy out */
ashleymills 0:714293de3836 2277 if (fp_iszero(a) == 1) {
ashleymills 0:714293de3836 2278 return 0;
ashleymills 0:714293de3836 2279 }
ashleymills 0:714293de3836 2280
ashleymills 0:714293de3836 2281 /* scan lower digits until non-zero */
ashleymills 0:714293de3836 2282 for (x = 0; x < a->used && a->dp[x] == 0; x++);
ashleymills 0:714293de3836 2283 q = a->dp[x];
ashleymills 0:714293de3836 2284 x *= DIGIT_BIT;
ashleymills 0:714293de3836 2285
ashleymills 0:714293de3836 2286 /* now scan this digit until a 1 is found */
ashleymills 0:714293de3836 2287 if ((q & 1) == 0) {
ashleymills 0:714293de3836 2288 do {
ashleymills 0:714293de3836 2289 qq = q & 15;
ashleymills 0:714293de3836 2290 x += lnz[qq];
ashleymills 0:714293de3836 2291 q >>= 4;
ashleymills 0:714293de3836 2292 } while (qq == 0);
ashleymills 0:714293de3836 2293 }
ashleymills 0:714293de3836 2294 return x;
ashleymills 0:714293de3836 2295 }
ashleymills 0:714293de3836 2296
ashleymills 0:714293de3836 2297
ashleymills 0:714293de3836 2298 /* c = (a, b) */
ashleymills 0:714293de3836 2299 void fp_gcd(fp_int *a, fp_int *b, fp_int *c)
ashleymills 0:714293de3836 2300 {
ashleymills 0:714293de3836 2301 fp_int u, v, r;
ashleymills 0:714293de3836 2302
ashleymills 0:714293de3836 2303 /* either zero than gcd is the largest */
ashleymills 0:714293de3836 2304 if (fp_iszero (a) == 1 && fp_iszero (b) == 0) {
ashleymills 0:714293de3836 2305 fp_abs (b, c);
ashleymills 0:714293de3836 2306 return;
ashleymills 0:714293de3836 2307 }
ashleymills 0:714293de3836 2308 if (fp_iszero (a) == 0 && fp_iszero (b) == 1) {
ashleymills 0:714293de3836 2309 fp_abs (a, c);
ashleymills 0:714293de3836 2310 return;
ashleymills 0:714293de3836 2311 }
ashleymills 0:714293de3836 2312
ashleymills 0:714293de3836 2313 /* optimized. At this point if a == 0 then
ashleymills 0:714293de3836 2314 * b must equal zero too
ashleymills 0:714293de3836 2315 */
ashleymills 0:714293de3836 2316 if (fp_iszero (a) == 1) {
ashleymills 0:714293de3836 2317 fp_zero(c);
ashleymills 0:714293de3836 2318 return;
ashleymills 0:714293de3836 2319 }
ashleymills 0:714293de3836 2320
ashleymills 0:714293de3836 2321 /* sort inputs */
ashleymills 0:714293de3836 2322 if (fp_cmp_mag(a, b) != FP_LT) {
ashleymills 0:714293de3836 2323 fp_init_copy(&u, a);
ashleymills 0:714293de3836 2324 fp_init_copy(&v, b);
ashleymills 0:714293de3836 2325 } else {
ashleymills 0:714293de3836 2326 fp_init_copy(&u, b);
ashleymills 0:714293de3836 2327 fp_init_copy(&v, a);
ashleymills 0:714293de3836 2328 }
ashleymills 0:714293de3836 2329
ashleymills 0:714293de3836 2330 fp_zero(&r);
ashleymills 0:714293de3836 2331 while (fp_iszero(&v) == FP_NO) {
ashleymills 0:714293de3836 2332 fp_mod(&u, &v, &r);
ashleymills 0:714293de3836 2333 fp_copy(&v, &u);
ashleymills 0:714293de3836 2334 fp_copy(&r, &v);
ashleymills 0:714293de3836 2335 }
ashleymills 0:714293de3836 2336 fp_copy(&u, c);
ashleymills 0:714293de3836 2337 }
ashleymills 0:714293de3836 2338
ashleymills 0:714293de3836 2339 #endif /* CYASSL_KEY_GEN */
ashleymills 0:714293de3836 2340
ashleymills 0:714293de3836 2341
ashleymills 0:714293de3836 2342 #if defined(HAVE_ECC) || !defined(NO_PWDBASED)
ashleymills 0:714293de3836 2343 /* c = a + b */
ashleymills 0:714293de3836 2344 void fp_add_d(fp_int *a, fp_digit b, fp_int *c)
ashleymills 0:714293de3836 2345 {
ashleymills 0:714293de3836 2346 fp_int tmp;
ashleymills 0:714293de3836 2347 fp_set(&tmp, b);
ashleymills 0:714293de3836 2348 fp_add(a,&tmp,c);
ashleymills 0:714293de3836 2349 }
ashleymills 0:714293de3836 2350
ashleymills 0:714293de3836 2351 /* external compatibility */
ashleymills 0:714293de3836 2352 int mp_add_d(fp_int *a, fp_digit b, fp_int *c)
ashleymills 0:714293de3836 2353 {
ashleymills 0:714293de3836 2354 fp_add_d(a, b, c);
ashleymills 0:714293de3836 2355 return MP_OKAY;
ashleymills 0:714293de3836 2356 }
ashleymills 0:714293de3836 2357
ashleymills 0:714293de3836 2358 #endif /* HAVE_ECC || !NO_PWDBASED */
ashleymills 0:714293de3836 2359
ashleymills 0:714293de3836 2360
ashleymills 0:714293de3836 2361 #ifdef HAVE_ECC
ashleymills 0:714293de3836 2362
ashleymills 0:714293de3836 2363 /* chars used in radix conversions */
ashleymills 0:714293de3836 2364 const char *fp_s_rmap = "0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZabcdefghijklmnopqrstuvwxyz+/";
ashleymills 0:714293de3836 2365
ashleymills 0:714293de3836 2366 static int fp_read_radix(fp_int *a, const char *str, int radix)
ashleymills 0:714293de3836 2367 {
ashleymills 0:714293de3836 2368 int y, neg;
ashleymills 0:714293de3836 2369 char ch;
ashleymills 0:714293de3836 2370
ashleymills 0:714293de3836 2371 /* make sure the radix is ok */
ashleymills 0:714293de3836 2372 if (radix < 2 || radix > 64) {
ashleymills 0:714293de3836 2373 return FP_VAL;
ashleymills 0:714293de3836 2374 }
ashleymills 0:714293de3836 2375
ashleymills 0:714293de3836 2376 /* if the leading digit is a
ashleymills 0:714293de3836 2377 * minus set the sign to negative.
ashleymills 0:714293de3836 2378 */
ashleymills 0:714293de3836 2379 if (*str == '-') {
ashleymills 0:714293de3836 2380 ++str;
ashleymills 0:714293de3836 2381 neg = FP_NEG;
ashleymills 0:714293de3836 2382 } else {
ashleymills 0:714293de3836 2383 neg = FP_ZPOS;
ashleymills 0:714293de3836 2384 }
ashleymills 0:714293de3836 2385
ashleymills 0:714293de3836 2386 /* set the integer to the default of zero */
ashleymills 0:714293de3836 2387 fp_zero (a);
ashleymills 0:714293de3836 2388
ashleymills 0:714293de3836 2389 /* process each digit of the string */
ashleymills 0:714293de3836 2390 while (*str) {
ashleymills 0:714293de3836 2391 /* if the radix < 36 the conversion is case insensitive
ashleymills 0:714293de3836 2392 * this allows numbers like 1AB and 1ab to represent the same value
ashleymills 0:714293de3836 2393 * [e.g. in hex]
ashleymills 0:714293de3836 2394 */
ashleymills 0:714293de3836 2395 ch = (char) ((radix < 36) ? XTOUPPER(*str) : *str);
ashleymills 0:714293de3836 2396 for (y = 0; y < 64; y++) {
ashleymills 0:714293de3836 2397 if (ch == fp_s_rmap[y]) {
ashleymills 0:714293de3836 2398 break;
ashleymills 0:714293de3836 2399 }
ashleymills 0:714293de3836 2400 }
ashleymills 0:714293de3836 2401
ashleymills 0:714293de3836 2402 /* if the char was found in the map
ashleymills 0:714293de3836 2403 * and is less than the given radix add it
ashleymills 0:714293de3836 2404 * to the number, otherwise exit the loop.
ashleymills 0:714293de3836 2405 */
ashleymills 0:714293de3836 2406 if (y < radix) {
ashleymills 0:714293de3836 2407 fp_mul_d (a, (fp_digit) radix, a);
ashleymills 0:714293de3836 2408 fp_add_d (a, (fp_digit) y, a);
ashleymills 0:714293de3836 2409 } else {
ashleymills 0:714293de3836 2410 break;
ashleymills 0:714293de3836 2411 }
ashleymills 0:714293de3836 2412 ++str;
ashleymills 0:714293de3836 2413 }
ashleymills 0:714293de3836 2414
ashleymills 0:714293de3836 2415 /* set the sign only if a != 0 */
ashleymills 0:714293de3836 2416 if (fp_iszero(a) != FP_YES) {
ashleymills 0:714293de3836 2417 a->sign = neg;
ashleymills 0:714293de3836 2418 }
ashleymills 0:714293de3836 2419 return FP_OKAY;
ashleymills 0:714293de3836 2420 }
ashleymills 0:714293de3836 2421
ashleymills 0:714293de3836 2422 /* fast math conversion */
ashleymills 0:714293de3836 2423 int mp_read_radix(mp_int *a, const char *str, int radix)
ashleymills 0:714293de3836 2424 {
ashleymills 0:714293de3836 2425 return fp_read_radix(a, str, radix);
ashleymills 0:714293de3836 2426 }
ashleymills 0:714293de3836 2427
ashleymills 0:714293de3836 2428 /* fast math conversion */
ashleymills 0:714293de3836 2429 int mp_set(fp_int *a, fp_digit b)
ashleymills 0:714293de3836 2430 {
ashleymills 0:714293de3836 2431 fp_set(a,b);
ashleymills 0:714293de3836 2432 return MP_OKAY;
ashleymills 0:714293de3836 2433 }
ashleymills 0:714293de3836 2434
ashleymills 0:714293de3836 2435 /* fast math conversion */
ashleymills 0:714293de3836 2436 int mp_sqr(fp_int *A, fp_int *B)
ashleymills 0:714293de3836 2437 {
ashleymills 0:714293de3836 2438 fp_sqr(A, B);
ashleymills 0:714293de3836 2439 return MP_OKAY;
ashleymills 0:714293de3836 2440 }
ashleymills 0:714293de3836 2441
ashleymills 0:714293de3836 2442 /* fast math conversion */
ashleymills 0:714293de3836 2443 int mp_montgomery_reduce(fp_int *a, fp_int *m, fp_digit mp)
ashleymills 0:714293de3836 2444 {
ashleymills 0:714293de3836 2445 fp_montgomery_reduce(a, m, mp);
ashleymills 0:714293de3836 2446 return MP_OKAY;
ashleymills 0:714293de3836 2447 }
ashleymills 0:714293de3836 2448
ashleymills 0:714293de3836 2449
ashleymills 0:714293de3836 2450 /* fast math conversion */
ashleymills 0:714293de3836 2451 int mp_montgomery_setup(fp_int *a, fp_digit *rho)
ashleymills 0:714293de3836 2452 {
ashleymills 0:714293de3836 2453 return fp_montgomery_setup(a, rho);
ashleymills 0:714293de3836 2454 }
ashleymills 0:714293de3836 2455
ashleymills 0:714293de3836 2456 int mp_div_2(fp_int * a, fp_int * b)
ashleymills 0:714293de3836 2457 {
ashleymills 0:714293de3836 2458 fp_div_2(a, b);
ashleymills 0:714293de3836 2459 return MP_OKAY;
ashleymills 0:714293de3836 2460 }
ashleymills 0:714293de3836 2461
ashleymills 0:714293de3836 2462
ashleymills 0:714293de3836 2463 int mp_init_copy(fp_int * a, fp_int * b)
ashleymills 0:714293de3836 2464 {
ashleymills 0:714293de3836 2465 fp_init_copy(a, b);
ashleymills 0:714293de3836 2466 return MP_OKAY;
ashleymills 0:714293de3836 2467 }
ashleymills 0:714293de3836 2468
ashleymills 0:714293de3836 2469
ashleymills 0:714293de3836 2470
ashleymills 0:714293de3836 2471 #endif /* HAVE_ECC */
ashleymills 0:714293de3836 2472
ashleymills 0:714293de3836 2473 #endif /* USE_FAST_MATH */