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