pd.if.org Git - zpackage/blob - tomsfastmath/src/mont/fp_montgomery_reduce.c

   1 /* TomsFastMath, a fast ISO C bignum library.
   2  *
   3  * This project is meant to fill in where LibTomMath
   4  * falls short.  That is speed ;-)
   5  *
   6  * This project is public domain and free for all purposes.
   7  *
   8  * Tom St Denis, tomstdenis@gmail.com
   9  */
  10 #include <tfm_private.h>
  11
  12 /******************************************************************/
  13 #if defined(TFM_X86) && !defined(TFM_SSE2)
  14 /* x86-32 code */
  15
  16 #define MONT_START
  17 #define MONT_FINI
  18 #define LOOP_END
  19 #define LOOP_START \
  20    mu = c[x] * mp
  21
  22 #define INNERMUL                                          \
  23 asm(                                                      \
  24    "movl %5,%%eax \n\t"                                   \
  25    "mull %4       \n\t"                                   \
  26    "addl %1,%%eax \n\t"                                   \
  27    "adcl $0,%%edx \n\t"                                   \
  28    "addl %%eax,%0 \n\t"                                   \
  29    "adcl $0,%%edx \n\t"                                   \
  30    "movl %%edx,%1 \n\t"                                   \
  31 :"=g"(_c[LO]), "=r"(cy)                                   \
  32 :"0"(_c[LO]), "1"(cy), "r"(mu), "r"(*tmpm++)              \
  33 : "%eax", "%edx", "cc")
  34
  35 #define PROPCARRY                           \
  36 asm(                                        \
  37    "addl   %1,%0    \n\t"                   \
  38    "setb   %%al     \n\t"                   \
  39    "movzbl %%al,%1 \n\t"                    \
  40 :"=g"(_c[LO]), "=r"(cy)                     \
  41 :"0"(_c[LO]), "1"(cy)                       \
  42 : "%eax", "cc")
  43
  44 /******************************************************************/
  45 #elif defined(TFM_X86_64)
  46 /* x86-64 code */
  47
  48 #define MONT_START
  49 #define MONT_FINI
  50 #define LOOP_END
  51 #define LOOP_START \
  52    mu = c[x] * mp
  53
  54 #define INNERMUL                                          \
  55 asm(                                                      \
  56    "movq %5,%%rax \n\t"                                   \
  57    "mulq %4       \n\t"                                   \
  58    "addq %1,%%rax \n\t"                                   \
  59    "adcq $0,%%rdx \n\t"                                   \
  60    "addq %%rax,%0 \n\t"                                   \
  61    "adcq $0,%%rdx \n\t"                                   \
  62    "movq %%rdx,%1 \n\t"                                   \
  63 :"=g"(_c[LO]), "=r"(cy)                                   \
  64 :"0"(_c[LO]), "1"(cy), "r"(mu), "r"(*tmpm++)              \
  65 : "%rax", "%rdx", "cc")
  66
  67 #define INNERMUL8 \
  68  asm(                  \
  69  "movq 0(%5),%%rax    \n\t"  \
  70  "movq 0(%2),%%r10    \n\t"  \
  71  "movq 0x8(%5),%%r11  \n\t"  \
  72  "mulq %4             \n\t"  \
  73  "addq %%r10,%%rax    \n\t"  \
  74  "adcq $0,%%rdx       \n\t"  \
  75  "movq 0x8(%2),%%r10  \n\t"  \
  76  "addq %3,%%rax       \n\t"  \
  77  "adcq $0,%%rdx       \n\t"  \
  78  "movq %%rax,0(%0)    \n\t"  \
  79  "movq %%rdx,%1       \n\t"  \
  80  \
  81  "movq %%r11,%%rax    \n\t"  \
  82  "movq 0x10(%5),%%r11 \n\t"  \
  83  "mulq %4             \n\t"  \
  84  "addq %%r10,%%rax    \n\t"  \
  85  "adcq $0,%%rdx       \n\t"  \
  86  "movq 0x10(%2),%%r10 \n\t"  \
  87  "addq %3,%%rax       \n\t"  \
  88  "adcq $0,%%rdx       \n\t"  \
  89  "movq %%rax,0x8(%0)  \n\t"  \
  90  "movq %%rdx,%1       \n\t"  \
  91  \
  92  "movq %%r11,%%rax    \n\t"  \
  93  "movq 0x18(%5),%%r11 \n\t"  \
  94  "mulq %4             \n\t"  \
  95  "addq %%r10,%%rax    \n\t"  \
  96  "adcq $0,%%rdx       \n\t"  \
  97  "movq 0x18(%2),%%r10 \n\t"  \
  98  "addq %3,%%rax       \n\t"  \
  99  "adcq $0,%%rdx       \n\t"  \
 100  "movq %%rax,0x10(%0) \n\t"  \
 101  "movq %%rdx,%1       \n\t"  \
 102  \
 103  "movq %%r11,%%rax    \n\t"  \
 104  "movq 0x20(%5),%%r11 \n\t"  \
 105  "mulq %4             \n\t"  \
 106  "addq %%r10,%%rax    \n\t"  \
 107  "adcq $0,%%rdx       \n\t"  \
 108  "movq 0x20(%2),%%r10 \n\t"  \
 109  "addq %3,%%rax       \n\t"  \
 110  "adcq $0,%%rdx       \n\t"  \
 111  "movq %%rax,0x18(%0) \n\t"  \
 112  "movq %%rdx,%1       \n\t"  \
 113  \
 114  "movq %%r11,%%rax    \n\t"  \
 115  "movq 0x28(%5),%%r11 \n\t"  \
 116  "mulq %4             \n\t"  \
 117  "addq %%r10,%%rax    \n\t"  \
 118  "adcq $0,%%rdx       \n\t"  \
 119  "movq 0x28(%2),%%r10 \n\t"  \
 120  "addq %3,%%rax       \n\t"  \
 121  "adcq $0,%%rdx       \n\t"  \
 122  "movq %%rax,0x20(%0) \n\t"  \
 123  "movq %%rdx,%1       \n\t"  \
 124  \
 125  "movq %%r11,%%rax    \n\t"  \
 126  "movq 0x30(%5),%%r11 \n\t"  \
 127  "mulq %4             \n\t"  \
 128  "addq %%r10,%%rax    \n\t"  \
 129  "adcq $0,%%rdx       \n\t"  \
 130  "movq 0x30(%2),%%r10 \n\t"  \
 131  "addq %3,%%rax       \n\t"  \
 132  "adcq $0,%%rdx       \n\t"  \
 133  "movq %%rax,0x28(%0) \n\t"  \
 134  "movq %%rdx,%1       \n\t"  \
 135  \
 136  "movq %%r11,%%rax    \n\t"  \
 137  "movq 0x38(%5),%%r11 \n\t"  \
 138  "mulq %4             \n\t"  \
 139  "addq %%r10,%%rax    \n\t"  \
 140  "adcq $0,%%rdx       \n\t"  \
 141  "movq 0x38(%2),%%r10 \n\t"  \
 142  "addq %3,%%rax       \n\t"  \
 143  "adcq $0,%%rdx       \n\t"  \
 144  "movq %%rax,0x30(%0) \n\t"  \
 145  "movq %%rdx,%1       \n\t"  \
 146  \
 147  "movq %%r11,%%rax    \n\t"  \
 148  "mulq %4             \n\t"  \
 149  "addq %%r10,%%rax    \n\t"  \
 150  "adcq $0,%%rdx       \n\t"  \
 151  "addq %3,%%rax       \n\t"  \
 152  "adcq $0,%%rdx       \n\t"  \
 153  "movq %%rax,0x38(%0) \n\t"  \
 154  "movq %%rdx,%1       \n\t"  \
 155  \
 156 :"=r"(_c), "=r"(cy)                    \
 157 : "0"(_c),  "1"(cy), "g"(mu), "r"(tmpm)\
 158 : "%rax", "%rdx", "%r10", "%r11", "cc")
 159
 160
 161 #define PROPCARRY                           \
 162 asm(                                        \
 163    "addq   %1,%0    \n\t"                   \
 164    "setb   %%al     \n\t"                   \
 165    "movzbq %%al,%1 \n\t"                    \
 166 :"=g"(_c[LO]), "=r"(cy)                     \
 167 :"0"(_c[LO]), "1"(cy)                       \
 168 : "%rax", "cc")
 169
 170 /******************************************************************/
 171 #elif defined(TFM_SSE2)
 172 /* SSE2 code (assumes 32-bit fp_digits) */
 173 /* XMM register assignments:
 174  * xmm0  *tmpm++, then Mu * (*tmpm++)
 175  * xmm1  c[x], then Mu
 176  * xmm2  mp
 177  * xmm3  cy
 178  * xmm4  _c[LO]
 179  */
 180
 181 #define MONT_START \
 182    asm("movd %0,%%mm2"::"g"(mp))
 183
 184 #define MONT_FINI \
 185    asm("emms")
 186
 187 #define LOOP_START          \
 188 asm(                        \
 189 "movd %0,%%mm1        \n\t" \
 190 "pxor %%mm3,%%mm3     \n\t" \
 191 "pmuludq %%mm2,%%mm1  \n\t" \
 192 :: "g"(c[x]))
 193
 194 /* pmuludq on mmx registers does a 32x32->64 multiply. */
 195 #define INNERMUL               \
 196 asm(                           \
 197    "movd %1,%%mm4        \n\t" \
 198    "movd %2,%%mm0        \n\t" \
 199    "paddq %%mm4,%%mm3    \n\t" \
 200    "pmuludq %%mm1,%%mm0  \n\t" \
 201    "paddq %%mm0,%%mm3    \n\t" \
 202    "movd %%mm3,%0        \n\t" \
 203    "psrlq $32, %%mm3     \n\t" \
 204 :"=g"(_c[LO]) : "0"(_c[LO]), "g"(*tmpm++) );
 205
 206 #define INNERMUL8 \
 207 asm(                           \
 208    "movd 0(%1),%%mm4     \n\t" \
 209    "movd 0(%2),%%mm0     \n\t" \
 210    "paddq %%mm4,%%mm3    \n\t" \
 211    "pmuludq %%mm1,%%mm0  \n\t" \
 212    "movd 4(%2),%%mm5     \n\t" \
 213    "paddq %%mm0,%%mm3    \n\t" \
 214    "movd 4(%1),%%mm6     \n\t" \
 215    "movd %%mm3,0(%0)     \n\t" \
 216    "psrlq $32, %%mm3     \n\t" \
 217 \
 218    "paddq %%mm6,%%mm3    \n\t" \
 219    "pmuludq %%mm1,%%mm5  \n\t" \
 220    "movd 8(%2),%%mm6     \n\t" \
 221    "paddq %%mm5,%%mm3    \n\t" \
 222    "movd 8(%1),%%mm7     \n\t" \
 223    "movd %%mm3,4(%0)     \n\t" \
 224    "psrlq $32, %%mm3     \n\t" \
 225 \
 226    "paddq %%mm7,%%mm3    \n\t" \
 227    "pmuludq %%mm1,%%mm6  \n\t" \
 228    "movd 12(%2),%%mm7    \n\t" \
 229    "paddq %%mm6,%%mm3    \n\t" \
 230    "movd 12(%1),%%mm5     \n\t" \
 231    "movd %%mm3,8(%0)     \n\t" \
 232    "psrlq $32, %%mm3     \n\t" \
 233 \
 234    "paddq %%mm5,%%mm3    \n\t" \
 235    "pmuludq %%mm1,%%mm7  \n\t" \
 236    "movd 16(%2),%%mm5    \n\t" \
 237    "paddq %%mm7,%%mm3    \n\t" \
 238    "movd 16(%1),%%mm6    \n\t" \
 239    "movd %%mm3,12(%0)    \n\t" \
 240    "psrlq $32, %%mm3     \n\t" \
 241 \
 242    "paddq %%mm6,%%mm3    \n\t" \
 243    "pmuludq %%mm1,%%mm5  \n\t" \
 244    "movd 20(%2),%%mm6    \n\t" \
 245    "paddq %%mm5,%%mm3    \n\t" \
 246    "movd 20(%1),%%mm7    \n\t" \
 247    "movd %%mm3,16(%0)    \n\t" \
 248    "psrlq $32, %%mm3     \n\t" \
 249 \
 250    "paddq %%mm7,%%mm3    \n\t" \
 251    "pmuludq %%mm1,%%mm6  \n\t" \
 252    "movd 24(%2),%%mm7    \n\t" \
 253    "paddq %%mm6,%%mm3    \n\t" \
 254    "movd 24(%1),%%mm5     \n\t" \
 255    "movd %%mm3,20(%0)    \n\t" \
 256    "psrlq $32, %%mm3     \n\t" \
 257 \
 258    "paddq %%mm5,%%mm3    \n\t" \
 259    "pmuludq %%mm1,%%mm7  \n\t" \
 260    "movd 28(%2),%%mm5    \n\t" \
 261    "paddq %%mm7,%%mm3    \n\t" \
 262    "movd 28(%1),%%mm6    \n\t" \
 263    "movd %%mm3,24(%0)    \n\t" \
 264    "psrlq $32, %%mm3     \n\t" \
 265 \
 266    "paddq %%mm6,%%mm3    \n\t" \
 267    "pmuludq %%mm1,%%mm5  \n\t" \
 268    "paddq %%mm5,%%mm3    \n\t" \
 269    "movd %%mm3,28(%0)    \n\t" \
 270    "psrlq $32, %%mm3     \n\t" \
 271 :"=r"(_c) : "0"(_c), "g"(tmpm) );
 272
 273 #define LOOP_END \
 274 asm( "movd %%mm3,%0  \n" :"=r"(cy))
 275
 276 #define PROPCARRY                           \
 277 asm(                                        \
 278    "addl   %1,%0    \n\t"                   \
 279    "setb   %%al     \n\t"                   \
 280    "movzbl %%al,%1 \n\t"                    \
 281 :"=g"(_c[LO]), "=r"(cy)                     \
 282 :"0"(_c[LO]), "1"(cy)                       \
 283 : "%eax", "cc")
 284
 285 /******************************************************************/
 286 #elif defined(TFM_ARM)
 287    /* ARMv4 code */
 288
 289 #define MONT_START
 290 #define MONT_FINI
 291 #define LOOP_END
 292 #define LOOP_START \
 293    mu = c[x] * mp
 294
 295 #define INNERMUL                    \
 296 asm(                                \
 297     " LDR    r0,%1            \n\t" \
 298     " ADDS   r0,r0,%0         \n\t" \
 299     " MOVCS  %0,#1            \n\t" \
 300     " MOVCC  %0,#0            \n\t" \
 301     " UMLAL  r0,%0,%3,%4      \n\t" \
 302     " STR    r0,%1            \n\t" \
 303 :"=r"(cy),"=m"(_c[0]):"0"(cy),"r"(mu),"r"(*tmpm++),"1"(_c[0]):"r0","cc");
 304
 305 #define PROPCARRY                  \
 306 asm(                               \
 307     " LDR   r0,%1            \n\t" \
 308     " ADDS  r0,r0,%0         \n\t" \
 309     " STR   r0,%1            \n\t" \
 310     " MOVCS %0,#1            \n\t" \
 311     " MOVCC %0,#0            \n\t" \
 312 :"=r"(cy),"=m"(_c[0]):"0"(cy),"1"(_c[0]):"r0","cc");
 313
 314 /******************************************************************/
 315 #elif defined(TFM_PPC32)
 316
 317 /* PPC32 */
 318 #define MONT_START
 319 #define MONT_FINI
 320 #define LOOP_END
 321 #define LOOP_START \
 322    mu = c[x] * mp
 323
 324 #define INNERMUL                     \
 325 asm(                                 \
 326    " mullw    16,%3,%4       \n\t"   \
 327    " mulhwu   17,%3,%4       \n\t"   \
 328    " addc     16,16,%2       \n\t"   \
 329    " addze    17,17          \n\t"   \
 330    " addc     %1,16,%5       \n\t"   \
 331    " addze    %0,17          \n\t"   \
 332 :"=r"(cy),"=r"(_c[0]):"0"(cy),"r"(mu),"r"(tmpm[0]),"1"(_c[0]):"16", "17", "cc"); ++tmpm;
 333
 334 #define PROPCARRY                    \
 335 asm(                                 \
 336    " addc     %1,%3,%2      \n\t"    \
 337    " xor      %0,%2,%2      \n\t"    \
 338    " addze    %0,%2         \n\t"    \
 339 :"=r"(cy),"=r"(_c[0]):"0"(cy),"1"(_c[0]):"cc");
 340
 341 /******************************************************************/
 342 #elif defined(TFM_PPC64)
 343
 344 /* PPC64 */
 345 #define MONT_START
 346 #define MONT_FINI
 347 #define LOOP_END
 348 #define LOOP_START \
 349    mu = c[x] * mp
 350
 351 #define INNERMUL                     \
 352 asm(                                 \
 353    " mulld    r16,%3,%4       \n\t"   \
 354    " mulhdu   r17,%3,%4       \n\t"   \
 355    " addc     r16,16,%0       \n\t"   \
 356    " addze    r17,r17          \n\t"   \
 357    " ldx      r18,0,%1        \n\t"   \
 358    " addc     r16,r16,r18       \n\t"   \
 359    " addze    %0,r17          \n\t"   \
 360    " sdx      r16,0,%1        \n\t"   \
 361 :"=r"(cy),"=m"(_c[0]):"0"(cy),"r"(mu),"r"(tmpm[0]),"1"(_c[0]):"r16", "r17", "r18","cc"); ++tmpm;
 362
 363 #define PROPCARRY                    \
 364 asm(                                 \
 365    " ldx      r16,0,%1       \n\t"    \
 366    " addc     r16,r16,%0      \n\t"    \
 367    " sdx      r16,0,%1       \n\t"    \
 368    " xor      %0,%0,%0      \n\t"    \
 369    " addze    %0,%0         \n\t"    \
 370 :"=r"(cy),"=m"(_c[0]):"0"(cy),"1"(_c[0]):"r16","cc");
 371
 372 /******************************************************************/
 373 #elif defined(TFM_AVR32)
 374
 375 /* AVR32 */
 376 #define MONT_START
 377 #define MONT_FINI
 378 #define LOOP_END
 379 #define LOOP_START \
 380    mu = c[x] * mp
 381
 382 #define INNERMUL                    \
 383 asm(                                \
 384     " ld.w   r2,%1            \n\t" \
 385     " add    r2,%0            \n\t" \
 386     " eor    r3,r3            \n\t" \
 387     " acr    r3               \n\t" \
 388     " macu.d r2,%3,%4         \n\t" \
 389     " st.w   %1,r2            \n\t" \
 390     " mov    %0,r3            \n\t" \
 391 :"=r"(cy),"=r"(_c):"0"(cy),"r"(mu),"r"(*tmpm++),"1"(_c):"r2","r3");
 392
 393 #define PROPCARRY                    \
 394 asm(                                 \
 395    " ld.w     r2,%1         \n\t"    \
 396    " add      r2,%0         \n\t"    \
 397    " st.w     %1,r2         \n\t"    \
 398    " eor      %0,%0         \n\t"    \
 399    " acr      %0            \n\t"    \
 400 :"=r"(cy),"=r"(&_c[0]):"0"(cy),"1"(&_c[0]):"r2","cc");
 401
 402 /******************************************************************/
 403 #elif defined(TFM_MIPS)
 404
 405 /* MIPS */
 406 #define MONT_START
 407 #define MONT_FINI
 408 #define LOOP_END
 409 #define LOOP_START \
 410    mu = c[x] * mp
 411
 412 #define INNERMUL                     \
 413 asm(                                 \
 414    " multu    %3,%4          \n\t"   \
 415    " mflo     $12            \n\t"   \
 416    " mfhi     $13            \n\t"   \
 417    " addu     $12,$12,%0     \n\t"   \
 418    " sltu     $10,$12,%0     \n\t"   \
 419    " addu     $13,$13,$10    \n\t"   \
 420    " lw       $10,%1         \n\t"   \
 421    " addu     $12,$12,$10    \n\t"   \
 422    " sltu     $10,$12,$10    \n\t"   \
 423    " addu     %0,$13,$10     \n\t"   \
 424    " sw       $12,%1         \n\t"   \
 425 :"=r"(cy),"=m"(_c[0]):"0"(cy),"r"(mu),"r"(tmpm[0]),"1"(_c[0]):"$10","$12","$13"); ++tmpm;
 426
 427 #define PROPCARRY                    \
 428 asm(                                 \
 429    " lw       $10,%1        \n\t"    \
 430    " addu     $10,$10,%0    \n\t"    \
 431    " sw       $10,%1        \n\t"    \
 432    " sltu     %0,$10,%0     \n\t"    \
 433 :"=r"(cy),"=m"(_c[0]):"0"(cy),"1"(_c[0]):"$10");
 434
 435 /******************************************************************/
 436 #else
 437
 438 /* ISO C code */
 439 #define MONT_START
 440 #define MONT_FINI
 441 #define LOOP_END
 442 #define LOOP_START \
 443    mu = c[x] * mp
 444
 445 #define INNERMUL                                      \
 446    do { fp_word t;                                    \
 447    _c[0] = t  = ((fp_word)_c[0] + (fp_word)cy) +      \
 448                 (((fp_word)mu) * ((fp_word)*tmpm++)); \
 449    cy = (t >> DIGIT_BIT);                             \
 450    } while (0)
 451
 452 #define PROPCARRY \
 453    do { fp_digit t = _c[0] += cy; cy = (t < cy); } while (0)
 454
 455 #endif
 456 /******************************************************************/
 457
 458
 459 #define LO  0
 460
 461 #ifdef TFM_SMALL_MONT_SET
 462 #include "fp_mont_small.i"
 463 #endif
 464
 465 /* computes x/R == x (mod N) via Montgomery Reduction */
 466 void fp_montgomery_reduce(fp_int *a, fp_int *m, fp_digit mp)
 467 {
 468    fp_digit c[FP_SIZE], *_c, *tmpm, mu;
 469    int      oldused, x, y, pa;
 470
 471    /* bail if too large */
 472    if (m->used > (FP_SIZE/2)) {
 473       return;
 474    }
 475
 476 #ifdef TFM_SMALL_MONT_SET
 477    if (m->used <= 16) {
 478       fp_montgomery_reduce_small(a, m, mp);
 479       return;
 480    }
 481 #endif
 482
 483 #if defined(USE_MEMSET)
 484    /* now zero the buff */
 485    memset(c, 0, sizeof c);
 486 #endif
 487    pa = m->used;
 488
 489    /* copy the input */
 490    oldused = a->used;
 491    for (x = 0; x < oldused; x++) {
 492        c[x] = a->dp[x];
 493    }
 494 #if !defined(USE_MEMSET)
 495    for (; x < 2*pa+1; x++) {
 496        c[x] = 0;
 497    }
 498 #endif
 499    MONT_START;
 500
 501    for (x = 0; x < pa; x++) {
 502        fp_digit cy = 0;
 503        /* get Mu for this round */
 504        LOOP_START;
 505        _c   = c + x;
 506        tmpm = m->dp;
 507        y = 0;
 508        #if defined(INNERMUL8)
 509         for (; y < (pa & ~7); y += 8) {
 510               INNERMUL8;
 511               _c   += 8;
 512               tmpm += 8;
 513            }
 514        #endif
 515
 516        for (; y < pa; y++) {
 517           INNERMUL;
 518           ++_c;
 519        }
 520        LOOP_END;
 521        while (cy) {
 522            PROPCARRY;
 523            ++_c;
 524        }
 525   }
 526
 527   /* now copy out */
 528   _c   = c + pa;
 529   tmpm = a->dp;
 530   for (x = 0; x < pa+1; x++) {
 531      *tmpm++ = *_c++;
 532   }
 533
 534   for (; x < oldused; x++)   {
 535      *tmpm++ = 0;
 536   }
 537
 538   MONT_FINI;
 539
 540   a->used = pa+1;
 541   fp_clamp(a);
 542
 543   /* if A >= m then A = A - m */
 544   if (fp_cmp_mag (a, m) != FP_LT) {
 545     s_fp_sub (a, m, a);
 546   }
 547 }
 548
 549
 550 /* $Source$ */
 551 /* $Revision$ */
 552 /* $Date$ */