1/* Copyright (c) 1998,2011-2014 Apple Inc. All Rights Reserved. 2 * 3 * NOTICE: USE OF THE MATERIALS ACCOMPANYING THIS NOTICE IS SUBJECT 4 * TO THE TERMS OF THE SIGNED "FAST ELLIPTIC ENCRYPTION (FEE) REFERENCE 5 * SOURCE CODE EVALUATION AGREEMENT" BETWEEN APPLE, INC. AND THE 6 * ORIGINAL LICENSEE THAT OBTAINED THESE MATERIALS FROM APPLE, 7 * INC. ANY USE OF THESE MATERIALS NOT PERMITTED BY SUCH AGREEMENT WILL 8 * EXPOSE YOU TO LIABILITY. 9 *************************************************************************** 10 11 giantIntegers.c - library for large-integer arithmetic. 12 13 Revision History 14 ---------------- 15 Fixed a==b bug in addg(). 16 10/06/98 ap 17 Changed to compile with C++. 18 13 Apr 98 Fixed shiftright(1) bug in modg_via_recip. 19 09 Apr 98 at Apple 20 Major rewrite of core arithmetic routines to make this module 21 independent of size of giantDigit. 22 Removed idivg() and radixdiv(). 23 20 Jan 98 at Apple 24 Deleted FFT arithmetic; simplified mulg(). 25 09 Jan 98 at Apple 26 gshiftright() optimization. 27 08 Jan 98 at Apple 28 newGiant() returns NULL on malloc failure 29 24 Dec 97 at Apple 30 New grammarSquare(); optimized modg_via_recip() 31 11 Jun 97 at Apple 32 Added modg_via_recip(), divg_via_recip(), make_recip() 33 Added new multiple giant stack mechanism 34 Fixed potential packing/alignment bug in copyGiant() 35 Added profiling for borrowGiant(), returnGiant() 36 Deleted obsolete ifdef'd code 37 Deleted newgiant() 38 All calls to borrowGiant() now specify required size (no more 39 borrowGiant(0) calls) 40 08 May 97 at Apple 41 Changed size of giantstruct.n to 1 for Mac build 42 05 Feb 97 at Apple 43 newGiant() no longer modifies CurrentMaxShorts or giant stack 44 Added modg profiling 45 01 Feb 97 at NeXT 46 Added iszero() check in gcompg 47 17 Jan 97 at NeXT 48 Fixed negation bug in gmersennemod() 49 Fixed n[words-1] == 0 bug in extractbits() 50 Cleaned up lots of static declarations 51 19 Sep 96 at NeXT 52 Fixed --size underflow bug in normal_subg(). 53 4 Sep 96 at NeXT 54 Fixed (b<n), (sign<0) case in gmersennemod() to allow for arbitrary n. 55 9 Aug 96 at NeXT 56 Fixed sign-extend bug in data_to_giant(). 57 7 Aug 96 at NeXT 58 Changed precision in newtondivide(). 59 Removed #ifdef UNUSED code. 60 24 Jul 96 at NeXT 61 Added scompg(). 62 Created. 63 64*/ 65 66#include <stdio.h> 67#include <stdlib.h> 68#include <string.h> 69#include "platform.h" 70#include "giantIntegers.h" 71#include "feeDebug.h" 72#include "ckconfig.h" 73#include "ellipticMeasure.h" 74#include "falloc.h" 75#include "giantPortCommon.h" 76 77#ifdef FEE_DEBUG 78#if (GIANT_LOG2_BITS_PER_DIGIT == 4) 79#warning Compiling with two-byte giantDigits 80#endif 81#endif 82 83#if 0 84#if FEE_DEBUG 85char printbuf1[200]; 86char printbuf2[200]; 87char printbuf3[200]; 88void printGiantBuf(giant x) 89{ 90 int i; 91 92 sprintf(printbuf2, "sign=%d cap=%d n[]=", x->sign, x->capacity); 93 for(i=0; i<abs(x->sign); i++) { 94 sprintf(printbuf3 + 10*i, "%u:", x->n[i]); 95 } 96} 97 98char printbuf4[200]; 99char printbuf5[200]; 100void printGiantBuf2(giant x) 101{ 102 int i; 103 104 sprintf(printbuf4, "sign=%d cap=%d n[]=", x->sign, x->capacity); 105 for(i=0; i<abs(x->sign); i++) { 106 sprintf(printbuf5 + 10*i, "%u:", x->n[i]); 107 } 108} 109#endif /* FEE_DEBUG */ 110#endif /* 0 */ 111 112/******** debugging flags *********/ 113 114/* 115 * Flag use of unoptimized divg, modg, binvg 116 */ 117//#define WARN_UNOPTIMIZE FEE_DEBUG 118#define WARN_UNOPTIMIZE 0 119 120/* 121 * Log interesting giant stack events 122 */ 123#define LOG_GIANT_STACK 0 124 125/* 126 * Log allocation of giant larger than stack size 127 */ 128#define LOG_GIANT_STACK_OVERFLOW 1 129 130/* 131 * Flag newGiant(0) and borrowGiant(0) calls 132 */ 133#define WARN_ZERO_GIANT_SIZE FEE_DEBUG 134 135/* temp mac-only giant initialization debug */ 136#define GIANT_MAC_DEBUG 0 137#if GIANT_MAC_DEBUG 138 139#include <string.h> 140#include <TextUtils.h> 141 142/* this one needs a writable string */ 143static void logCom(unsigned char *str) { 144 c2pstr((char *)str); 145 DebugStr(str); 146} 147 148/* constant strings */ 149void dblog0(const char *str) { 150 Str255 outStr; 151 strcpy((char *)outStr, str); 152 logCom(outStr); 153} 154 155#else 156#define dblog0(s) 157 158#endif /* GIANT_MAC_DEBUG */ 159 160#ifndef min 161#define min(a,b) ((a)<(b)? (a) : (b)) 162#endif // min 163#ifndef max 164#define max(a,b) ((a)>(b)? (a) : (b)) 165#endif // max 166 167#ifndef TRUE 168#define TRUE 1 169#endif // TRUE 170#ifndef FALSE 171#define FALSE 0 172#endif // FALSE 173 174static void absg(giant g); /* g := |g|. */ 175 176/************** globals *******************/ 177 178 179/* ------ giant stack package ------ */ 180 181/* 182 * The giant stack package is a local cache which allows us to avoid calls 183 * to malloc() for borrowGiant(). On a 90 Mhz Pentium, enabling the 184 * giant stack package shows about a 1.35 speedup factor over an identical 185 * CryptKit without the giant stacks enabled. 186 */ 187 188#if GIANTS_VIA_STACK 189 190#if LOG_GIANT_STACK 191#define gstackDbg(x) printf x 192#else // LOG_GIANT_STACK 193#define gstackDbg(x) 194#endif // LOG_GIANT_STACK 195 196typedef struct { 197 unsigned numDigits; // capacity of giants in this stack 198 unsigned numFree; // number of free giants in stack 199 unsigned totalGiants; // total number in *stack 200 giant *stack; 201} gstack; 202 203static gstack *gstacks = NULL; // array of stacks 204static unsigned numGstacks = 0; // # of elements in gstacks 205static int gstackInitd = 0; // this module has been init'd 206 207#define INIT_NUM_GIANTS 16 /* initial # of giants / stack */ 208#define MIN_GIANT_SIZE 4 /* numDigits for gstack[0] */ 209#define GIANT_SIZE_INCR 2 /* in << bits */ 210 211/* 212 * Initialize giant stacks, with up to specified max giant size. 213 */ 214void initGiantStacks(unsigned maxDigits) 215{ 216 unsigned curSize = MIN_GIANT_SIZE; 217 unsigned sz; 218 unsigned i; 219 220 dblog0("initGiantStacks\n"); 221 222 if(gstackInitd) { 223 /* 224 * Shouldn't be called more than once... 225 */ 226 printf("multiple initGiantStacks calls\n"); 227 return; 228 } 229 gstackDbg(("initGiantStacks(%d)\n", maxDigits)); 230 231 /* 232 * How many stacks? 233 */ 234 numGstacks = 1; 235 while(curSize<=maxDigits) { 236 curSize <<= GIANT_SIZE_INCR; 237 numGstacks++; 238 } 239 240 sz = sizeof(gstack) * numGstacks; 241 gstacks = (gstack*) fmalloc(sz); 242 bzero(gstacks, sz); 243 244 curSize = MIN_GIANT_SIZE; 245 for(i=0; i<numGstacks; i++) { 246 gstacks[i].numDigits = curSize; 247 curSize <<= GIANT_SIZE_INCR; 248 } 249 250 gstackInitd = 1; 251} 252 253/* called at shut down - free resources */ 254void freeGiantStacks(void) 255{ 256 int i; 257 int j; 258 gstack *gs; 259 260 if(!gstackInitd) { 261 return; 262 } 263 for(i=0; i<numGstacks; i++) { 264 gs = &gstacks[i]; 265 for(j=0; j<gs->numFree; j++) { 266 freeGiant(gs->stack[j]); 267 gs->stack[j] = NULL; 268 } 269 /* and the stack itself - may be null if this was never used */ 270 if(gs->stack != NULL) { 271 ffree(gs->stack); 272 gs->stack = NULL; 273 } 274 } 275 ffree(gstacks); 276 gstacks = NULL; 277 gstackInitd = 0; 278} 279 280#endif // GIANTS_VIA_STACK 281 282giant borrowGiant(unsigned numDigits) 283{ 284 giant result; 285 286 #if GIANTS_VIA_STACK 287 288 unsigned stackNum; 289 gstack *gs = gstacks; 290 291 #if WARN_ZERO_GIANT_SIZE 292 if(numDigits == 0) { 293 printf("borrowGiant(0)\n"); 294 numDigits = gstacks[numGstacks-1].numDigits; 295 } 296 #endif // WARN_ZERO_GIANT_SIZE 297 298 /* 299 * Find appropriate stack 300 */ 301 if(numDigits <= MIN_GIANT_SIZE) 302 stackNum = 0; 303 else if (numDigits <= (MIN_GIANT_SIZE << GIANT_SIZE_INCR)) 304 stackNum = 1; 305 else if (numDigits <= (MIN_GIANT_SIZE << (2 * GIANT_SIZE_INCR))) 306 stackNum = 2; 307 else if (numDigits <= (MIN_GIANT_SIZE << (3 * GIANT_SIZE_INCR))) 308 stackNum = 3; 309 else if (numDigits <= (MIN_GIANT_SIZE << (4 * GIANT_SIZE_INCR))) 310 stackNum = 4; 311 else 312 stackNum = numGstacks; 313 314 if(stackNum >= numGstacks) { 315 /* 316 * out of bounds; just malloc 317 */ 318 #if LOG_GIANT_STACK_OVERFLOW 319 gstackDbg(("giantFromStack overflow; numDigits %d\n", 320 numDigits)); 321 #endif // LOG_GIANT_STACK_OVERFLOW 322 return newGiant(numDigits); 323 } 324 gs = &gstacks[stackNum]; 325 326 #if GIANT_MAC_DEBUG 327 if((gs->numFree != 0) && (gs->stack == NULL)) { 328 dblog0("borrowGiant: null stack!\n"); 329 } 330 #endif 331 332 if(gs->numFree != 0) { 333 result = gs->stack[--gs->numFree]; 334 } 335 else { 336 /* 337 * Stack empty; malloc 338 */ 339 result = newGiant(gs->numDigits); 340 } 341 342 #else /* GIANTS_VIA_STACK */ 343 344 result = newGiant(numDigits); 345 346 #endif /* GIANTS_VIA_STACK */ 347 348 PROF_INCR(numBorrows); 349 return result; 350} 351 352void returnGiant(giant g) 353{ 354 355 #if GIANTS_VIA_STACK 356 357 unsigned stackNum; 358 gstack *gs; 359 unsigned cap = g->capacity; 360 361 362 #if FEE_DEBUG 363 if(!gstackInitd) { 364 CKRaise("returnGiant before stacks initialized!"); 365 } 366 #endif // FEE_DEBUG 367 368 #if GIANT_MAC_DEBUG 369 if(g == NULL) { 370 dblog0("returnGiant: null g!\n"); 371 } 372 #endif 373 374 /* 375 * Find appropriate stack. Note we expect exact match of 376 * capacity and stack's giant size. 377 */ 378 /* 379 * Optimized unrolled loop. Just make sure there are enough cases 380 * to handle all of the stacks. Errors in this case will be flagged 381 * via LOG_GIANT_STACK_OVERFLOW. 382 */ 383 switch(cap) { 384 case MIN_GIANT_SIZE: 385 stackNum = 0; 386 break; 387 case MIN_GIANT_SIZE << GIANT_SIZE_INCR: 388 stackNum = 1; 389 break; 390 case MIN_GIANT_SIZE << (2 * GIANT_SIZE_INCR): 391 stackNum = 2; 392 break; 393 case MIN_GIANT_SIZE << (3 * GIANT_SIZE_INCR): 394 stackNum = 3; 395 break; 396 case MIN_GIANT_SIZE << (4 * GIANT_SIZE_INCR): 397 stackNum = 4; 398 break; 399 default: 400 stackNum = numGstacks; 401 break; 402 } 403 404 if(stackNum >= numGstacks) { 405 /* 406 * out of bounds; just free 407 */ 408 #if LOG_GIANT_STACK_OVERFLOW 409 gstackDbg(("giantToStack overflow; numDigits %d\n", cap)); 410 #endif // LOG_GIANT_STACK_OVERFLOW 411 freeGiant(g); 412 return; 413 } 414 gs = &gstacks[stackNum]; 415 if(gs->numFree == gs->totalGiants) { 416 if(gs->totalGiants == 0) { 417 gstackDbg(("Initial alloc of gstack(%d)\n", 418 gs->numDigits)); 419 gs->totalGiants = INIT_NUM_GIANTS; 420 } 421 else { 422 gs->totalGiants *= 2; 423 gstackDbg(("Bumping gstack(%d) to %d\n", 424 gs->numDigits, gs->totalGiants)); 425 } 426 gs->stack = (giantstruct**) frealloc(gs->stack, gs->totalGiants*sizeof(giant)); 427 } 428 g->sign = 0; // not sure this is important... 429 gs->stack[gs->numFree++] = g; 430 431 #if GIANT_MAC_DEBUG 432 if((gs->numFree != 0) && (gs->stack == NULL)) { 433 dblog0("borrowGiant: null stack!\n"); 434 } 435 #endif 436 437 #else /* GIANTS_VIA_STACK */ 438 439 freeGiant(g); 440 441 #endif /* GIANTS_VIA_STACK */ 442} 443 444void freeGiant(giant x) { 445 ffree(x); 446} 447 448giant newGiant(unsigned numDigits) { 449 // giant sufficient for 2^numbits+16 sized ops 450 int size; 451 giant result; 452 453 #if WARN_ZERO_GIANT_SIZE 454 if(numDigits == 0) { 455 printf("newGiant(0)\n"); 456 #if GIANTS_VIA_STACK 457 numDigits = gstacks[numGstacks-1].totalGiants; 458 #else 459 /* HACK */ 460 numDigits = 20; 461 #endif 462 } 463 #endif // WARN_ZERO_GIANT_SIZE 464 465 size = (numDigits-1) * GIANT_BYTES_PER_DIGIT + sizeof(giantstruct); 466 result = (giant)fmalloc(size); 467 if(result == NULL) { 468 return NULL; 469 } 470 result->sign = 0; 471 result->capacity = numDigits; 472 return result; 473} 474 475giant copyGiant(giant x) 476{ 477 int bytes; 478 479 giant result = newGiant(x->capacity); 480 481 /* 482 * 13 Jun 1997 483 * NO! this assumes packed alignment 484 */ 485 bytes = sizeof(giantstruct) + 486 ((x->capacity - 1) * GIANT_BYTES_PER_DIGIT); 487 bcopy(x, result, bytes); 488 return result; 489} 490 491/* ------ initialization and utility routines ------ */ 492 493 494unsigned bitlen(giant n) { 495 unsigned b = GIANT_BITS_PER_DIGIT; 496 giantDigit c = 1 << (GIANT_BITS_PER_DIGIT - 1); 497 giantDigit w; 498 499 if (isZero(n)) { 500 return(0); 501 } 502 w = n->n[abs(n->sign) - 1]; 503 if (!w) { 504 CKRaise("bitlen - no bit set!"); 505 } 506 while((w&c) == 0) { 507 b--; 508 c >>= 1; 509 } 510 return(GIANT_BITS_PER_DIGIT * (abs(n->sign)-1) + b); 511} 512 513int bitval(giant n, int pos) { 514 int i = abs(pos) >> GIANT_LOG2_BITS_PER_DIGIT; 515 giantDigit c = 1 << (pos & (GIANT_BITS_PER_DIGIT - 1)); 516 517 return((n->n[i]) & c); 518} 519 520int gsign(giant g) 521/* returns the sign of g */ 522{ 523 if (isZero(g)) return(0); 524 if (g->sign > 0) return(1); 525 return(-1); 526} 527 528/* 529 * Adjust sign for possible leading (m.s.) zero digits 530 */ 531void gtrimSign(giant g) 532{ 533 int numDigits = abs(g->sign); 534 int i; 535 536 for(i=numDigits-1; i>=0; i--) { 537 if(g->n[i] == 0) { 538 numDigits--; 539 } 540 else { 541 break; 542 } 543 } 544 if(g->sign < 0) { 545 g->sign = -numDigits; 546 } 547 else { 548 g->sign = numDigits; 549 } 550} 551 552 553int isone(giant g) { 554 return((g->sign==1)&&(g->n[0]==1)); 555} 556 557int isZero(giant thegiant) { 558/* Returns TRUE if thegiant == 0. */ 559 int count; 560 int length = abs(thegiant->sign); 561 giantDigit *numpointer; 562 563 if (length) { 564 numpointer = thegiant->n; 565 566 for(count = 0; count<length; ++count,++numpointer) 567 if (*numpointer != 0 ) return(FALSE); 568 } 569 return(TRUE); 570} 571 572int gcompg(giant a, giant b) 573/* returns -1,0,1 if a<b, a=b, a>b, respectively */ 574{ 575 int sa = a->sign; 576 int j; 577 int sb = b->sign; 578 giantDigit va; 579 giantDigit vb; 580 int sgn; 581 582 if(isZero(a) && isZero(b)) return 0; 583 if(sa > sb) return(1); 584 if(sa < sb) return(-1); 585 if(sa < 0) { 586 sa = -sa; /* Take absolute value of sa */ 587 sgn = -1; 588 } else sgn = 1; 589 for(j = sa-1; j >= 0; j--) { 590 va = a->n[j]; vb = b->n[j]; 591 if (va > vb) return(sgn); 592 if (va < vb) return(-sgn); 593 } 594 return(0); 595} 596 597/* destgiant becomes equal to srcgiant */ 598void gtog(giant srcgiant, giant destgiant) { 599 600 int numbytes; 601 602 CKASSERT(srcgiant != NULL); 603 numbytes = abs(srcgiant->sign) * GIANT_BYTES_PER_DIGIT; 604 if (destgiant->capacity < abs(srcgiant->sign)) 605 CKRaise("gtog overflow!!"); 606 memcpy((char *)destgiant->n, (char *)srcgiant->n, numbytes); 607 destgiant->sign = srcgiant->sign; 608} 609 610void int_to_giant(int i, giant g) { 611/* The giant g becomes set to the integer value i. */ 612 int isneg = (i<0); 613 unsigned int j = abs(i); 614 unsigned dex; 615 616 g->sign = 0; 617 if (i==0) { 618 g->n[0] = 0; 619 return; 620 } 621 622 if(GIANT_BYTES_PER_DIGIT == sizeof(int)) { 623 g->n[0] = j; 624 g->sign = 1; 625 } 626 else { 627 /* one loop per digit */ 628 unsigned scnt = GIANT_BITS_PER_DIGIT; // fool compiler 629 630 for(dex=0; dex<sizeof(int); dex++) { 631 g->n[dex] = j & GIANT_DIGIT_MASK; 632 j >>= scnt; 633 g->sign++; 634 if(j == 0) { 635 break; 636 } 637 } 638 } 639 if (isneg) { 640 g->sign = -(g->sign); 641 } 642} 643 644/*------------- Arithmetic --------------*/ 645 646void negg(giant g) { 647/* g becomes -g */ 648 g->sign = -g->sign; 649} 650 651void iaddg(int i, giant g) { /* positive g becomes g + (int)i */ 652 int j; 653 giantDigit carry; 654 int size = abs(g->sign); 655 656 if (isZero(g)) { 657 int_to_giant(i,g); 658 } 659 else { 660 carry = i; 661 for(j=0; ((j<size) && (carry != 0)); j++) { 662 g->n[j] = giantAddDigits(g->n[j], carry, &carry); 663 } 664 if(carry) { 665 ++g->sign; 666 // realloc 667 if (g->sign > (int)g->capacity) CKRaise("iaddg overflow!"); 668 g->n[size] = carry; 669 } 670 } 671} 672 673/* 674 * g *= (int n) 675 * 676 * FIXME - we can improve this... 677 */ 678void imulg(unsigned n, giant g) 679{ 680 giant tmp = borrowGiant(abs(g->sign) + sizeof(int)); 681 682 int_to_giant(n, tmp); 683 mulg(tmp, g); 684 returnGiant(tmp); 685} 686 687static void normal_addg(giant a, giant b) 688/* b := a + b, both a,b assumed non-negative. */ 689{ 690 giantDigit carry1 = 0; 691 giantDigit carry2 = 0; 692 int asize = a->sign, bsize = b->sign; 693 giantDigit *an = a->n; 694 giantDigit *bn = b->n; 695 giantDigit tmp; 696 int j; 697 int comSize; 698 int maxSize; 699 700 if(asize < bsize) { 701 comSize = asize; 702 maxSize = bsize; 703 } 704 else { 705 comSize = bsize; 706 maxSize = asize; 707 } 708 709 /* first handle the common digits */ 710 for(j=0; j<comSize; j++) { 711 /* 712 * first add the carry, then an[j] - either add could result 713 * in another carry 714 */ 715 if(carry1 || carry2) { 716 tmp = giantAddDigits(bn[j], (giantDigit)1, &carry1); 717 } 718 else { 719 carry1 = 0; 720 tmp = bn[j]; 721 } 722 bn[j] = giantAddDigits(tmp, an[j], &carry2); 723 } 724 725 if(asize < bsize) { 726 727 /* now propagate remaining carry beyond asize */ 728 if(carry2) { 729 carry1 = 1; 730 } 731 if(carry1) { 732 for(; j<bsize; j++) { 733 bn[j] = giantAddDigits(bn[j], (giantDigit)1, &carry1); 734 if(carry1 == 0) { 735 break; 736 } 737 } 738 } 739 } else { 740 /* now propagate remaining an[] and carry beyond bsize */ 741 if(carry2) { 742 carry1 = 1; 743 } 744 for(; j<asize; j++) { 745 if(carry1) { 746 bn[j] = giantAddDigits(an[j], (giantDigit)1, &carry1); 747 } 748 else { 749 bn[j] = an[j]; 750 carry1 = 0; 751 } 752 } 753 } 754 b->sign = maxSize; 755 if(carry1) { 756 // realloc? 757 bn[j] = 1; 758 b->sign++; 759 if (b->sign > (int)b->capacity) CKRaise("iaddg overflow!"); 760 } 761 762} 763 764static void normal_subg(giant a, giant b) 765/* b := b - a; requires b, a non-negative and b >= a. */ 766{ 767 int j; 768 int size = b->sign; 769 giantDigit tmp; 770 giantDigit borrow1 = 0; 771 giantDigit borrow2 = 0; 772 giantDigit *an = a->n; 773 giantDigit *bn = b->n; 774 775 if(a->sign == 0) { 776 return; 777 } 778 779 for (j=0; j<a->sign; ++j) { 780 if(borrow1 || borrow2) { 781 tmp = giantSubDigits(bn[j], (giantDigit)1, &borrow1); 782 } 783 else { 784 tmp = bn[j]; 785 borrow1 = 0; 786 } 787 bn[j] = giantSubDigits(tmp, an[j], &borrow2); 788 } 789 if(borrow1 || borrow2) { 790 /* propagate borrow thru remainder of bn[] */ 791 borrow1 = 1; 792 for (j=a->sign; j<size; ++j) { 793 if(borrow1) { 794 bn[j] = giantSubDigits(bn[j], (giantDigit)1, &borrow1); 795 } 796 else { 797 break; 798 } 799 } 800 } 801 802 /* adjust sign for leading zero digits */ 803 while((size-- > 0) && (b->n[size] == 0)) 804 ; 805 b->sign = (b->n[size] == 0)? 0 : size+1; 806} 807 808static void reverse_subg(giant a, giant b) 809/* b := a - b; requires b, a non-negative and a >= b. */ 810{ 811 int j; 812 int size = a->sign; 813 giantDigit tmp; 814 giantDigit borrow1 = 0; 815 giantDigit borrow2 = 0; 816 giantDigit *an = a->n; 817 giantDigit *bn = b->n; 818 819 if(b->sign == 0) { 820 gtog(a, b); 821 return; 822 } 823 for (j=0; j<b->sign; ++j) { 824 if(borrow1 || borrow2) { 825 tmp = giantSubDigits(an[j], (giantDigit)1, &borrow1); 826 } 827 else { 828 tmp = an[j]; 829 borrow1 = 0; 830 } 831 bn[j] = giantSubDigits(tmp, bn[j], &borrow2); 832 } 833 if(borrow1 || borrow2) { 834 /* propagate borrow thru remainder of bn[] */ 835 borrow1 = 1; 836 } 837 for (j=b->sign; j<size; ++j) { 838 if(borrow1) { 839 bn[j] = giantSubDigits(an[j], (giantDigit)1, &borrow1); 840 } 841 else { 842 bn[j] = an[j]; 843 borrow1 = 0; 844 } 845 } 846 847 b->sign = size; /* REC, 21 Apr 1996. */ 848 while(!b->n[--size]); 849 b->sign = size+1; 850} 851 852 853void addg(giant a, giant b) 854/* b := b + a, any signs any result. */ 855{ int asgn = a->sign, bsgn = b->sign; 856 if(asgn == 0) return; 857 if(bsgn == 0) { 858 gtog(a,b); 859 return; 860 } 861 if((asgn < 0) == (bsgn < 0)) { 862 if(bsgn > 0) { 863 normal_addg(a,b); 864 return; 865 } 866 negg(a); if(a != b) negg(b); normal_addg(a,b); /* Fix REC 1 Dec 98. */ 867 negg(a); if(a != b) negg(b); return; /* Fix REC 1 Dec 98. */ 868 } 869 if(bsgn > 0) { 870 negg(a); 871 if(gcompg(b,a) >= 0) { 872 normal_subg(a,b); 873 negg(a); 874 return; 875 } 876 reverse_subg(a,b); 877 negg(a); 878 negg(b); 879 return; 880 } 881 negg(b); 882 if(gcompg(b,a) < 0) { 883 reverse_subg(a,b); 884 return; 885 } 886 normal_subg(a,b); 887 negg(b); 888 return; 889} 890 891void subg(giant a, giant b) 892/* b := b - a, any signs, any result. */ 893{ 894 int asgn = a->sign, bsgn = b->sign; 895 if(asgn == 0) return; 896 if(bsgn == 0) { 897 gtog(a,b); 898 negg(b); 899 return; 900 } 901 if((asgn < 0) != (bsgn < 0)) { 902 if(bsgn > 0) { 903 negg(a); 904 normal_addg(a,b); 905 negg(a); 906 return; 907 } 908 negg(b); 909 normal_addg(a,b); 910 negg(b); 911 return; 912 } 913 if(bsgn > 0) { 914 if(gcompg(b,a) >= 0) { 915 normal_subg(a,b); 916 return; 917 } 918 reverse_subg(a,b); 919 negg(b); 920 return; 921 } 922 negg(a); negg(b); 923 if(gcompg(b,a) >= 0) { 924 normal_subg(a,b); 925 negg(a); 926 negg(b); 927 return; 928 } 929 reverse_subg(a,b); 930 negg(a); 931 return; 932} 933 934static void bdivg(giant v, giant u) 935/* u becomes greatest power of two not exceeding u/v. */ 936{ 937 int diff = bitlen(u) - bitlen(v); 938 giant scratch7; 939 940 if (diff<0) { 941 int_to_giant(0,u); 942 return; 943 } 944 scratch7 = borrowGiant(u->capacity); 945 gtog(v, scratch7); 946 gshiftleft(diff,scratch7); 947 if(gcompg(u,scratch7) < 0) diff--; 948 if(diff<0) { 949 int_to_giant(0,u); 950 returnGiant(scratch7); 951 return; 952 } 953 int_to_giant(1,u); 954 gshiftleft(diff,u); 955 returnGiant(scratch7); 956} 957 958int binvaux(giant p, giant x) 959/* Binary inverse method. 960 Returns zero if no inverse exists, in which case x becomes 961 GCD(x,p). */ 962{ 963 giant scratch7; 964 giant u0; 965 giant u1; 966 giant v0; 967 giant v1; 968 int result = 1; 969 int giantSize; 970 PROF_START; 971 972 if(isone(x)) return(result); 973 giantSize = 4 * abs(p->sign); 974 scratch7 = borrowGiant(giantSize); 975 u0 = borrowGiant(giantSize); 976 u1 = borrowGiant(giantSize); 977 v0 = borrowGiant(giantSize); 978 v1 = borrowGiant(giantSize); 979 int_to_giant(1, v0); gtog(x, v1); 980 int_to_giant(0,x); gtog(p, u1); 981 while(!isZero(v1)) { 982 gtog(u1, u0); bdivg(v1, u0); 983 gtog(x, scratch7); 984 gtog(v0, x); 985 mulg(u0, v0); 986 subg(v0,scratch7); 987 gtog(scratch7, v0); 988 989 gtog(u1, scratch7); 990 gtog(v1, u1); 991 mulg(u0, v1); 992 subg(v1,scratch7); 993 gtog(scratch7, v1); 994 } 995 if (!isone(u1)) { 996 gtog(u1,x); 997 if(x->sign<0) addg(p, x); 998 result = 0; 999 goto done; 1000 } 1001 if (x->sign<0) addg(p, x); 1002 done: 1003 returnGiant(scratch7); 1004 returnGiant(u0); 1005 returnGiant(u1); 1006 returnGiant(v0); 1007 returnGiant(v1); 1008 PROF_END(binvauxTime); 1009 return(result); 1010} 1011 1012/* 1013 * Superceded by binvg_cp() 1014 */ 1015#if 0 1016int binvg(giant p, giant x) 1017{ 1018 modg(p, x); 1019 return(binvaux(p,x)); 1020} 1021#endif 1022 1023static void absg(giant g) { 1024/* g becomes the absolute value of g */ 1025 if (g->sign < 0) g->sign = -g->sign; 1026} 1027 1028void gshiftleft(int bits, giant g) { 1029/* shift g left bits bits. Equivalent to g = g*2^bits */ 1030 int rem = bits & (GIANT_BITS_PER_DIGIT - 1); 1031 int crem = GIANT_BITS_PER_DIGIT - rem; 1032 int digits = 1 + (bits >> GIANT_LOG2_BITS_PER_DIGIT); 1033 int size = abs(g->sign); 1034 int j; 1035 int k; 1036 int sign = gsign(g); 1037 giantDigit carry; 1038 giantDigit dat; 1039 1040 #if FEE_DEBUG 1041 if(bits < 0) { 1042 CKRaise("gshiftleft(-bits)\n"); 1043 } 1044 #endif /* FEE_DEBUG */ 1045 1046 if(!bits) return; 1047 if(!size) return; 1048 if((size+digits) > (int)g->capacity) { 1049 CKRaise("gshiftleft overflow"); 1050 return; 1051 } 1052 k = size - 1 + digits; // (MSD of result + 1) 1053 carry = 0; 1054 1055 /* bug fix for 32-bit giantDigits; this is also an optimization for 1056 * other sizes. rem=0 means we're shifting strictly by digits, no 1057 * bit shifts. */ 1058 if(rem == 0) { 1059 g->n[k] = 0; // XXX hack - for sign fixup 1060 for(j=size-1; j>=0; j--) { 1061 g->n[--k] = g->n[j]; 1062 } 1063 do{ 1064 g->n[--k] = 0; 1065 } while(k>0); 1066 } 1067 else { 1068 /* 1069 * normal unaligned case 1070 * FIXME - this writes past g->n[size-1] the first time thru! 1071 */ 1072 for(j=size-1; j>=0; j--) { 1073 dat = g->n[j]; 1074 g->n[k--] = (dat >> crem) | carry; 1075 carry = (dat << rem); 1076 } 1077 do{ 1078 g->n[k--] = carry; 1079 carry = 0; 1080 } while(k>=0); 1081 } 1082 k = size - 1 + digits; 1083 if(g->n[k] == 0) --k; 1084 g->sign = sign * (k+1); 1085 if (abs(g->sign) > g->capacity) { 1086 CKRaise("gshiftleft overflow"); 1087 } 1088} 1089 1090void gshiftright(int bits, giant g) { 1091/* shift g right bits bits. Equivalent to g = g/2^bits */ 1092 int j; 1093 int size=abs(g->sign); 1094 giantDigit carry; 1095 int digits = bits >> GIANT_LOG2_BITS_PER_DIGIT; 1096 int remain = bits & (GIANT_BITS_PER_DIGIT - 1); 1097 int cremain = GIANT_BITS_PER_DIGIT - remain; 1098 1099 #if FEE_DEBUG 1100 if(bits < 0) { 1101 CKRaise("gshiftright(-bits)\n"); 1102 } 1103 #endif /* FEE_DEBUG */ 1104 if(bits==0) return; 1105 if(isZero(g)) return; 1106 if (digits >= size) { 1107 g->sign = 0; 1108 return; 1109 } 1110 1111 size -= digits; 1112 1113/* Begin OPT: 9 Jan 98 REC. */ 1114 if(remain == 0) { 1115 if(g->sign > 0) { 1116 g->sign = size; 1117 } 1118 else { 1119 g->sign = -size; 1120 } 1121 for(j=0; j < size; j++) { 1122 g->n[j] = g->n[j+digits]; 1123 } 1124 return; 1125 } 1126/* End OPT: 9 Jan 98 REC. */ 1127 1128 for(j=0;j<size;++j) { 1129 if (j==size-1) { 1130 carry = 0; 1131 } 1132 else { 1133 carry = (g->n[j+digits+1]) << cremain; 1134 } 1135 g->n[j] = ((g->n[j+digits]) >> remain ) | carry; 1136 } 1137 if (g->n[size-1] == 0) { 1138 --size; 1139 } 1140 if(g->sign > 0) { 1141 g->sign = size; 1142 } 1143 else { 1144 g->sign = -size; 1145 } 1146 if (abs(g->sign) > g->capacity) { 1147 CKRaise("gshiftright overflow"); 1148 } 1149} 1150 1151 1152void extractbits(unsigned n, giant src, giant dest) { 1153/* dest becomes lowermost n bits of src. Equivalent to dest = src % 2^n */ 1154 int digits = n >> GIANT_LOG2_BITS_PER_DIGIT; 1155 int numbytes = digits * GIANT_BYTES_PER_DIGIT; 1156 int bits = n & (GIANT_BITS_PER_DIGIT - 1); 1157 1158 if (n <= 0) { 1159 return; 1160 } 1161 if (dest->capacity * 8 * GIANT_BYTES_PER_DIGIT < n) { 1162 CKRaise("extractbits - not enough room"); 1163 } 1164 if (digits >= abs(src->sign)) { 1165 gtog(src,dest); 1166 } 1167 else { 1168 memcpy((char *)(dest->n), (char *)(src->n), numbytes); 1169 if (bits) { 1170 dest->n[digits] = src->n[digits] & ((1<<bits)-1); 1171 ++digits; 1172 } 1173 /* Next, fix by REC, 12 Jan 97. */ 1174 // while((dest->n[words-1] == 0) && (words > 0)) --words; 1175 while((digits > 0) && (dest->n[digits-1] == 0)) { 1176 --digits; 1177 } 1178 if(src->sign < 0) { 1179 dest->sign = -digits; 1180 } 1181 else { 1182 dest->sign = digits; 1183 } 1184 } 1185 if (abs(dest->sign) > dest->capacity) { 1186 CKRaise("extractbits overflow"); 1187 } 1188} 1189 1190#define NEW_MERSENNE 0 1191 1192/* 1193 * New gmersennemod, 24 Dec 1997. This runs significantly slower than the 1194 * original. 1195 */ 1196#if NEW_MERSENNE 1197 1198void 1199gmersennemod( 1200 int n, 1201 giant g 1202) 1203/* g := g (mod 2^n - 1) */ 1204{ 1205 int the_sign; 1206 giant scratch3 = borrowGiant(g->capacity); 1207 giant scratch4 = borrowGiant(1); 1208 1209 if ((the_sign = gsign(g)) < 0) absg(g); 1210 while (bitlen(g) > n) { 1211 gtog(g,scratch3); 1212 gshiftright(n,scratch3); 1213 addg(scratch3,g); 1214 gshiftleft(n,scratch3); 1215 subg(scratch3,g); 1216 } 1217 if(isZero(g)) goto out; 1218 int_to_giant(1,scratch3); 1219 gshiftleft(n,scratch3); 1220 int_to_giant(1,scratch4); 1221 subg(scratch4,scratch3); 1222 if(gcompg(g,scratch3) >= 0) subg(scratch3,g); 1223 if (the_sign < 0) { 1224 g->sign = -g->sign; 1225 addg(scratch3,g); 1226 } 1227out: 1228 returnGiant(scratch3); 1229 returnGiant(scratch4); 1230} 1231 1232#else /* NEW_MERSENNE */ 1233 1234void gmersennemod(int n, giant g) { 1235/* g becomes g mod ((2^n)-1) 1236 31 Jul 96 modified REC. 1237 17 Jan 97 modified REC. 1238*/ 1239 unsigned bits = n & (GIANT_BITS_PER_DIGIT - 1); 1240 unsigned digits = 1 + ((n-1) >> GIANT_LOG2_BITS_PER_DIGIT); 1241 int isPositive = (g->sign > 0); 1242 int j; 1243 int b; 1244 int size; 1245 int foundzero; 1246 giantDigit mask = (bits == 0) ? GIANT_DIGIT_MASK : (giantDigit)((1<<bits)-1); 1247 giant scratch1; 1248 1249 b = bitlen(g); 1250 if(b < n) { 1251 unsigned numDigits = (n + GIANT_BITS_PER_DIGIT - 1) >> 1252 GIANT_LOG2_BITS_PER_DIGIT; 1253 giantDigit lastWord = 0; 1254 giantDigit bits = 1; 1255 1256 if(g->sign >= 0) return; 1257 1258 /* 1259 * Cons up ((2**n)-1), add to g. 1260 */ 1261 scratch1 = borrowGiant(numDigits + 1); 1262 scratch1->sign = numDigits; 1263 for(j=0; j<(int)(numDigits-1); j++) { 1264 scratch1->n[j] = GIANT_DIGIT_MASK; 1265 } 1266 1267 /* 1268 * Last word has lower (n & (GIANT_BITS_PER_DIGIT-1)) bits set. 1269 */ 1270 for(j=0; j < (int)(n & (GIANT_BITS_PER_DIGIT-1)); j++) { 1271 lastWord |= bits; 1272 bits <<= 1; 1273 } 1274 scratch1->n[numDigits-1] = lastWord; 1275 addg(g, scratch1); /* One version. */ 1276 gtog(scratch1, g); 1277 returnGiant(scratch1); 1278 return; 1279 } 1280 if(b == n) { 1281 for(foundzero=0, j=0; j<b; j++) { 1282 if(bitval(g, j)==0) { 1283 foundzero = 1; 1284 break; 1285 } 1286 } 1287 if (!foundzero) { 1288 int_to_giant(0, g); 1289 return; 1290 } 1291 } 1292 1293 absg(g); 1294 scratch1 = borrowGiant(g->capacity); 1295 while ( ((unsigned)(g->sign) > digits) || 1296 ( ((unsigned)(g->sign)==digits) && (g->n[digits-1] > mask))) { 1297 extractbits(n, g, scratch1); 1298 gshiftright(n, g); 1299 addg(scratch1, g); 1300 } 1301 size = g->sign; 1302 1303/* Commence new negation routine - REC 17 Jan 1997. */ 1304 if (!isPositive) { /* Mersenne negation is just bitwise complement. */ 1305 for(j = digits-1; j >= size; j--) { 1306 g->n[j] = GIANT_DIGIT_MASK; 1307 } 1308 for(j = size-1; j >= 0; j--) { 1309 g->n[j] = ~g->n[j]; 1310 } 1311 g->n[digits-1] &= mask; 1312 j = digits-1; 1313 while((g->n[j] == 0) && (j > 0)) { 1314 --j; 1315 } 1316 size = j+1; 1317 } 1318/* End new negation routine. */ 1319 1320 g->sign = size; 1321 if (abs(g->sign) > g->capacity) { 1322 CKRaise("gmersennemod overflow"); 1323 } 1324 if (size < (int)digits) { 1325 goto bye; 1326 } 1327 if (g->n[size-1] != mask) { 1328 goto bye; 1329 } 1330 mask = GIANT_DIGIT_MASK; 1331 for(j=0; j<(size-1); j++) { 1332 if (g->n[j] != mask) { 1333 goto bye; 1334 } 1335 } 1336 g->sign = 0; 1337 bye: 1338 returnGiant(scratch1); 1339} 1340 1341#endif /* NEW_MERSENNE */ 1342 1343void mulg(giant a, giant b) { /* b becomes a*b. */ 1344 1345 int i; 1346 int asize, bsize; 1347 giantDigit *bptr = b->n; 1348 giantDigit mult; 1349 giant scratch1; 1350 giantDigit carry; 1351 giantDigit *scrPtr; 1352 1353 1354 if (isZero(b)) { 1355 return; 1356 } 1357 if (isZero(a)) { 1358 gtog(a, b); 1359 return; 1360 } 1361 if(a == b) { 1362 grammarSquare(b); 1363 return; 1364 } 1365 1366 bsize = abs(b->sign); 1367 asize = abs(a->sign); 1368 scratch1 = borrowGiant((asize+bsize)); 1369 scrPtr = scratch1->n; 1370 1371 for(i=0; i<asize+bsize; ++i) { 1372 scrPtr[i]=0; 1373 } 1374 1375 for(i=0; i<bsize; ++i, scrPtr++) { 1376 mult = bptr[i]; 1377 if (mult != 0) { 1378 carry = VectorMultiply(mult, 1379 a->n, 1380 asize, 1381 scrPtr); 1382 /* handle MSD carry */ 1383 scrPtr[asize] += carry; 1384 } 1385 } 1386 bsize+=asize; 1387 if(scratch1->n[bsize - 1] == 0) { 1388 --bsize; 1389 } 1390 scratch1->sign = gsign(a) * gsign(b) * bsize; 1391 if (abs(scratch1->sign) > scratch1->capacity) { 1392 CKRaise("GiantGrammarMul overflow"); 1393 } 1394 gtog(scratch1,b); 1395 returnGiant(scratch1); 1396 1397 #if FEE_DEBUG 1398 (void)bitlen(b); // Assertion.... 1399 #endif /* FEE_DEBUG */ 1400 PROF_INCR(numMulg); // for normal profiling 1401 INCR_MULGS; // for ellipticMeasure 1402} 1403 1404void grammarSquare(giant a) { 1405 /* 1406 * For now, we're going to match the old implementation line for 1407 * line by maintaining prod, carry, and temp as double precision 1408 * giantDigits. There is probably a much better implementation.... 1409 */ 1410 giantDigit prodLo; 1411 giantDigit prodHi; 1412 giantDigit carryLo = 0; 1413 giantDigit carryHi = 0; 1414 giantDigit tempLo; 1415 giantDigit tempHi; 1416 unsigned int cur_term; 1417 unsigned asize; 1418 unsigned max; 1419 giantDigit *ptr = a->n; 1420 giantDigit *ptr1; 1421 giantDigit *ptr2; 1422 giant scratch; 1423 1424 /* dmitch 11 Jan 1998 - special case for a == 0 */ 1425 if(a->sign == 0) { 1426 goto end; 1427 } 1428 /* end a == 0 case */ 1429 asize = abs(a->sign); 1430 max = asize * 2 - 1; 1431 scratch = borrowGiant(2 * asize); 1432 asize--; 1433 1434 /* 1435 * temp = *ptr; 1436 * temp *= temp; 1437 * scratch->n[0] = temp; 1438 * carry = temp >> 16; 1439 */ 1440 giantMulDigits(*ptr, *ptr, &tempLo, &tempHi); 1441 scratch->n[0] = tempLo; 1442 carryLo = tempHi; 1443 carryHi = 0; 1444 1445 for (cur_term = 1; cur_term < max; cur_term++) { 1446 ptr1 = ptr2 = ptr; 1447 if (cur_term <= asize) { 1448 ptr2 += cur_term; 1449 } else { 1450 ptr1 += cur_term - asize; 1451 ptr2 += asize; 1452 } 1453 1454 /* 1455 * prod = carry & 0xFFFF; 1456 * carry >>= 16; 1457 */ 1458 prodLo = carryLo; 1459 prodHi = 0; 1460 carryLo = carryHi; 1461 carryHi = 0; 1462 while(ptr1 < ptr2) { 1463 /* 1464 * temp = *ptr1++ * *ptr2--; 1465 */ 1466 giantMulDigits(*ptr1++, *ptr2--, &tempLo, &tempHi); 1467 1468 /* 1469 * prod += (temp << 1) & 0xFFFF; 1470 */ 1471 giantAddDouble(&prodLo, &prodHi, (tempLo << 1)); 1472 1473 /* 1474 * carry += (temp >> 15); 1475 * use bits from both product digits.. 1476 */ 1477 giantAddDouble(&carryLo, &carryHi, 1478 (tempLo >> (GIANT_BITS_PER_DIGIT - 1))); 1479 giantAddDouble(&carryLo, &carryHi, (tempHi << 1)); 1480 1481 /* snag the msb from that last shift */ 1482 carryHi += (tempHi >> (GIANT_BITS_PER_DIGIT - 1)); 1483 } 1484 if (ptr1 == ptr2) { 1485 /* 1486 * temp = *ptr1; 1487 * temp *= temp; 1488 */ 1489 giantMulDigits(*ptr1, *ptr1, &tempLo, &tempHi); 1490 1491 /* 1492 * prod += temp & 0xFFFF; 1493 */ 1494 giantAddDouble(&prodLo, &prodHi, tempLo); 1495 1496 /* 1497 * carry += (temp >> 16); 1498 */ 1499 giantAddDouble(&carryLo, &carryHi, tempHi); 1500 } 1501 1502 /* 1503 * carry += prod >> 16; 1504 */ 1505 giantAddDouble(&carryLo, &carryHi, prodHi); 1506 1507 scratch->n[cur_term] = prodLo; 1508 } 1509 if (carryLo) { 1510 scratch->n[cur_term] = carryLo; 1511 scratch->sign = cur_term+1; 1512 } else scratch->sign = cur_term; 1513 1514 gtog(scratch,a); 1515 returnGiant(scratch); 1516end: 1517 PROF_INCR(numGsquare); 1518} 1519 1520/* 1521 * Clear all of a giant's data fields, for secure erasure of sensitive data., 1522 */ 1523void clearGiant(giant g) 1524{ 1525 unsigned i; 1526 1527 for(i=0; i<g->capacity; i++) { 1528 g->n[i] = 0; 1529 } 1530 g->sign = 0; 1531} 1532 1533#if ENGINE_127_BITS 1534/* 1535 * only used by engineNSA127.c, which is obsolete as of 16 Jan 1997 1536 */ 1537int 1538scompg(int n, giant g) { 1539 if((g->sign == 1) && (g->n[0] == n)) return(1); 1540 return(0); 1541} 1542 1543#endif // ENGINE_127_BITS 1544 1545/* 1546 */ 1547 1548/* 1549 * Calculate the reciprocal of a demonimator used in divg_via_recip() and 1550 * modg_via_recip(). 1551 */ 1552void 1553make_recip(giant d, giant r) 1554/* r becomes the steady-state reciprocal 1555 2^(2b)/d, where b = bit-length of d-1. */ 1556{ 1557 int b; 1558 int giantSize = 4 * abs(d->sign); 1559 giant tmp = borrowGiant(giantSize); 1560 giant tmp2 = borrowGiant(giantSize); 1561 1562 if (isZero(d) || (d->sign < 0)) 1563 { 1564 CKRaise("illegal argument to make_recip"); 1565 } 1566 int_to_giant(1, r); subg(r, d); b = bitlen(d); addg(r, d); 1567 gshiftleft(b, r); gtog(r, tmp2); 1568 while(1) { 1569 gtog(r, tmp); 1570 gsquare(tmp); 1571 gshiftright(b, tmp); 1572 mulg(d, tmp); 1573 gshiftright(b, tmp); 1574 addg(r, r); subg(tmp, r); 1575 if(gcompg(r, tmp2) <= 0) break; 1576 gtog(r, tmp2); 1577 } 1578 int_to_giant(1, tmp); 1579 gshiftleft(2*b, tmp); 1580 gtog(r, tmp2); mulg(d, tmp2); 1581 subg(tmp2, tmp); 1582 int_to_giant(1, tmp2); 1583 while(tmp->sign < 0) { 1584 subg(tmp2, r); 1585 addg(d, tmp); 1586 } 1587 1588 returnGiant(tmp); 1589 returnGiant(tmp2); 1590 return; 1591} 1592 1593/* 1594 * Optimized divg, when reciprocal of denominator is known. 1595 */ 1596void 1597divg_via_recip(giant d, giant r, giant n) 1598/* n := n/d, where r is the precalculated 1599 steady-state reciprocal of d. */ 1600{ 1601 int s = 2*(bitlen(r)-1), sign = gsign(n); 1602 int giantSize = (4 * abs(d->sign)) + abs(n->sign); 1603 giant tmp = borrowGiant(giantSize); 1604 giant tmp2 = borrowGiant(giantSize); 1605 1606 if (isZero(d) || (d->sign < 0)) 1607 { 1608 CKRaise("illegal argument to divg_via_recip"); 1609 } 1610 n->sign = abs(n->sign); 1611 int_to_giant(0, tmp2); 1612 while(1) { 1613 gtog(n, tmp); 1614 mulg(r, tmp); 1615 gshiftright(s, tmp); 1616 addg(tmp, tmp2); 1617 mulg(d, tmp); 1618 subg(tmp, n); 1619 if(gcompg(n,d) >= 0) { 1620 subg(d,n); 1621 iaddg(1, tmp2); 1622 } 1623 if(gcompg(n,d) < 0) break; 1624 } 1625 gtog(tmp2, n); 1626 n->sign *= sign; 1627 returnGiant(tmp); 1628 returnGiant(tmp2); 1629 return; 1630} 1631 1632/* 1633 * Optimized modg, when reciprocal of denominator is known. 1634 */ 1635 1636/* New version, 24 Dec 1997. */ 1637 1638void 1639modg_via_recip( 1640 giant d, 1641 giant r, 1642 giant n 1643) 1644/* This is the fastest mod of the present collection. 1645 n := n % d, where r is the precalculated 1646 steady-state reciprocal of d. */ 1647 1648{ 1649 int s = (bitlen(r)-1), sign = n->sign; 1650 int giantSize = (4 * abs(d->sign)) + abs(n->sign); 1651 giant tmp, tmp2; 1652 1653 tmp = borrowGiant(giantSize); 1654 tmp2 = borrowGiant(giantSize); 1655 if (isZero(d) || (d->sign < 0)) 1656 { 1657 CKRaise("illegal argument to modg_via_recip"); 1658 } 1659 n->sign = abs(n->sign); 1660 while (1) 1661 { 1662 gtog(n, tmp); 1663 /* bug fix 13 Apr 1998 */ 1664 if(s == 0) { 1665 gshiftleft(1, tmp); 1666 } 1667 else { 1668 gshiftright(s-1, tmp); 1669 } 1670 /* end fix */ 1671 mulg(r, tmp); 1672 gshiftright(s+1, tmp); 1673 mulg(d, tmp); 1674 subg(tmp, n); 1675 if (gcompg(n,d) >= 0) 1676 subg(d,n); 1677 if (gcompg(n,d) < 0) 1678 break; 1679 } 1680 if (sign >= 0) 1681 goto done; 1682 if (isZero(n)) 1683 goto done; 1684 negg(n); 1685 addg(d,n); 1686done: 1687 returnGiant(tmp); 1688 returnGiant(tmp2); 1689 return; 1690} 1691 1692/* 1693 * Unoptimized, inefficient general modg, when reciprocal of denominator 1694 * is not known. 1695 */ 1696void 1697modg( 1698 giant d, 1699 giant n 1700) 1701{ 1702 /* n becomes n%d. n is arbitrary, but the denominator d must be 1703 * positive! */ 1704 1705 /* 1706 * 4/9/2001: seeing overflow on this recip. Alloc per 1707 * d->capacity, not d->sign. 1708 */ 1709 //giant recip = borrowGiant(2 * abs(d->sign)); 1710 giant recip = borrowGiant(2 * d->capacity); 1711 1712 #if WARN_UNOPTIMIZE 1713 dbgLog(("Warning: unoptimized modg!\n")); 1714 #endif // WARN_UNOPTIMIZE 1715 1716 make_recip(d, recip); 1717 modg_via_recip(d, recip, n); 1718 returnGiant(recip); 1719} 1720 1721/* 1722 * Unoptimized, inefficient general divg, when reciprocal of denominator 1723 * is not known. 1724 */ 1725void 1726divg( 1727 giant d, 1728 giant n 1729) 1730{ 1731 /* n becomes n/d. n is arbitrary, but the denominator d must be 1732 * positive! 1733 */ 1734 1735 giant recip = borrowGiant(2 * abs(d->sign)); 1736 1737 #if WARN_UNOPTIMIZE 1738 dbgLog(("Warning: unoptimized divg!\n")); 1739 #endif // WARN_UNOPTIMIZE 1740 1741 make_recip(d, recip); 1742 divg_via_recip(d, recip, n); 1743 returnGiant(recip); 1744} 1745