46 #define GMPSPREAD (GMP_LIMB_BITS/BITSINWORD) 62 VOID Pack(UWORD *a, WORD *na, UWORD *b, WORD nb)
66 if ( (c = *na) == 0 ) {
67 MLOCK(ErrorMessageLock);
68 MesPrint(
"Caught a zero in Pack");
69 MUNLOCK(ErrorMessageLock);
73 MLOCK(ErrorMessageLock);
74 MesPrint(
"Division by zero in Pack");
75 MUNLOCK(ErrorMessageLock);
78 if ( *na < 0 ) { sgn = -sgn; c = -c; }
79 if ( nb < 0 ) { sgn = -sgn; nb = -nb; }
83 while ( --i >= 0 ) *to++ = 0;
87 while ( --i >= 0 ) *to++ = 0;
88 if ( sgn < 0 ) *na = -*na;
100 VOID UnPack(UWORD *a, WORD na, WORD *denom, WORD *numer)
104 if ( na < 0 ) { na = -na; }
110 while ( !(*a) ) { i--; a--; }
111 while ( !(*pos) ) { na--; pos--; }
114 if ( sgn < 0 ) i = -i;
126 WORD Mully(
PHEAD UWORD *a, WORD *na, UWORD *b, WORD nb)
131 WORD nd, ne, adenom, anumer;
132 if ( !nb ) { *na = 0;
return(0); }
133 else if ( *b == 1 ) {
134 if ( nb == 1 )
return(0);
135 else if ( nb == -1 ) { *na = -*na;
return(0); }
137 if ( *na < 0 ) { sgn = -sgn; *na = -*na; }
138 if ( nb < 0 ) { sgn = -sgn; nb = -nb; }
139 UnPack(a,*na,&adenom,&anumer);
140 d = NumberMalloc(
"Mully"); e = NumberMalloc(
"Mully");
141 for ( i = 0; i < nb; i++ ) { e[i] = *b++; }
143 if ( Simplify(BHEAD a+*na,&adenom,e,&ne) )
goto MullyEr;
144 if ( MulLong(a,anumer,e,ne,d,&nd) )
goto MullyEr;
146 for ( i = 0; i < *na; i++ ) { e[i] = *b++; }
151 for ( i = 0; i < *na; i++ ) { a[i] = *b++; }
153 if ( sgn < 0 ) *na = -*na;
154 NumberFree(d,
"Mully"); NumberFree(e,
"Mully");
157 MLOCK(ErrorMessageLock);
159 MUNLOCK(ErrorMessageLock);
160 NumberFree(d,
"Mully"); NumberFree(e,
"Mully");
172 WORD Divvy(
PHEAD UWORD *a, WORD *na, UWORD *b, WORD nb)
177 WORD nd, ne, adenom, anumer;
179 MLOCK(ErrorMessageLock);
180 MesPrint(
"Division by zero in Divvy");
181 MUNLOCK(ErrorMessageLock);
184 d = NumberMalloc(
"Divvy"); e = NumberMalloc(
"Divvy");
185 if ( nb < 0 ) { sgn = -sgn; nb = -nb; }
186 if ( *na < 0 ) { sgn = -sgn; *na = -*na; }
187 UnPack(a,*na,&adenom,&anumer);
188 for ( i = 0; i < nb; i++ ) { e[i] = *b++; }
190 if ( Simplify(BHEAD a,&anumer,e,&ne) )
goto DivvyEr;
191 if ( MulLong(a+*na,adenom,e,ne,d,&nd) )
goto DivvyEr;
194 if ( sgn < 0 ) *na = -*na;
195 NumberFree(d,
"Divvy"); NumberFree(e,
"Divvy");
198 MLOCK(ErrorMessageLock);
200 MUNLOCK(ErrorMessageLock);
201 NumberFree(d,
"Divvy"); NumberFree(e,
"Divvy");
210 WORD AddRat(
PHEAD UWORD *a, WORD na, UWORD *b, WORD nb, UWORD *c, WORD *nc)
213 UWORD *d, *e, *f, *g;
214 WORD nd, ne, nf, ng, adenom, anumer, bdenom, bnumer;
218 if ( nb < 0 ) nb = -nb;
220 for ( i = 0; i < nb; i++ ) *c++ = *b++;
226 if ( na < 0 ) na = -na;
228 for ( i = 0; i < na; i++ ) *c++ = *a++;
231 else if ( b[1] == 1 && a[1] == 1 ) {
236 if ( *c < *a ) { c[2] = 1; c[3] = 0; *nc = 2; }
240 else if ( nb == -1 ) {
242 *c = *b - *a; *nc = -1;
244 else if ( *b < *a ) {
245 *c = *a - *b; *nc = 1;
252 else if ( na == -1 ){
256 if ( *c < *a ) { c[2] = 1; c[3] = 0; *nc = -2; }
260 else if ( nb == 1 ) {
262 *c = *b - *a; *nc = 1;
264 else if ( *b < *a ) {
265 *c = *a - *b; *nc = -1;
273 UnPack(a,na,&adenom,&anumer);
274 UnPack(b,nb,&bdenom,&bnumer);
275 if ( na < 0 ) na = -na;
276 if ( nb < 0 ) nb = -nb;
277 if ( na == 1 && nb == 1 ) {
279 t3 = ((RLONG)a[1])*((RLONG)b[1]);
280 t1 = ((RLONG)a[0])*((RLONG)b[1]);
281 t2 = ((RLONG)a[1])*((RLONG)b[0]);
282 if ( ( anumer > 0 && bnumer > 0 ) || ( anumer < 0 && bnumer < 0 ) ) {
283 if ( ( t1 = t1 + t2 ) < t2 ) {
286 c[1] = (UWORD)(t1 >> BITSINWORD);
291 if ( ( c[1] = (UWORD)(t1 >> BITSINWORD) ) != 0 ) *nc = 2;
296 if ( t1 == t2 ) { *nc = 0;
return(0); }
305 if ( ( c[1] = (UWORD)(t1 >> BITSINWORD) ) != 0 ) *nc = 2;
308 if ( anumer < 0 ) *nc = -*nc;
309 d = NumberMalloc(
"AddRat");
311 if ( ( d[1] = (UWORD)(t3 >> BITSINWORD) ) != 0 ) nd = 2;
313 if ( Simplify(BHEAD c,nc,d,&nd) )
goto AddRer1;
324 d = NumberMalloc(
"AddRat"); e = NumberMalloc(
"AddRat");
325 f = NumberMalloc(
"AddRat"); g = NumberMalloc(
"AddRat");
326 if ( GcdLong(BHEAD a+na,adenom,b+nb,bdenom,d,&nd) )
goto AddRer;
327 if ( *d == 1 && nd == 1 ) nd = 0;
329 if ( DivLong(a+na,adenom,d,nd,e,&ne,c,nc) )
goto AddRer;
330 if ( DivLong(b+nb,bdenom,d,nd,f,&nf,c,nc) )
goto AddRer;
331 if ( MulLong(a,anumer,f,nf,c,nc) )
goto AddRer;
332 if ( MulLong(b,bnumer,e,ne,g,&ng) )
goto AddRer;
335 if ( MulLong(a+na,adenom,b,bnumer,c,nc) )
goto AddRer;
336 if ( MulLong(b+nb,bdenom,a,anumer,g,&ng) )
goto AddRer;
338 if ( AddLong(c,*nc,g,ng,c,nc) )
goto AddRer;
340 NumberFree(g,
"AddRat"); NumberFree(f,
"AddRat");
341 NumberFree(e,
"AddRat"); NumberFree(d,
"AddRat");
345 if ( Simplify(BHEAD c,nc,d,&nd) )
goto AddRer;
346 if ( MulLong(e,ne,d,nd,g,&ng) )
goto AddRer;
347 if ( MulLong(g,ng,f,nf,d,&nd) )
goto AddRer;
350 if ( MulLong(a+na,adenom,b+nb,bdenom,d,&nd) )
goto AddRer;
352 NumberFree(g,
"AddRat"); NumberFree(f,
"AddRat"); NumberFree(e,
"AddRat");
355 NumberFree(d,
"AddRat");
358 NumberFree(g,
"AddRat"); NumberFree(f,
"AddRat"); NumberFree(e,
"AddRat");
360 NumberFree(d,
"AddRat");
362 MLOCK(ErrorMessageLock);
364 MUNLOCK(ErrorMessageLock);
378 WORD MulRat(
PHEAD UWORD *a, WORD na, UWORD *b, WORD nb, UWORD *c, WORD *nc)
382 if ( *b == 1 && b[1] == 1 ) {
385 i = ABS(na); i <<= 1;
386 while ( --i >= 0 ) *c++ = *a++;
389 else if ( nb == -1 ) {
391 i = ABS(na); i <<= 1;
392 while ( --i >= 0 ) *c++ = *a++;
396 if ( *a == 1 && a[1] == 1 ) {
399 i = ABS(nb); i <<= 1;
400 while ( --i >= 0 ) *c++ = *b++;
403 else if ( na == -1 ) {
405 i = ABS(nb); i <<= 1;
406 while ( --i >= 0 ) *c++ = *b++;
410 if ( na < 0 ) { na = -na; sgn = -sgn; }
411 if ( nb < 0 ) { nb = -nb; sgn = -sgn; }
412 if ( !na || !nb ) { *nc = 0;
return(0); }
413 if ( na != 1 || nb != 1 ) {
415 UWORD *xd,*xe, *xf,*xg;
416 WORD dden, dnumr, eden, enumr;
417 UnPack(a,na,&dden,&dnumr);
418 UnPack(b,nb,&eden,&enumr);
419 xd = NumberMalloc(
"MulRat"); xf = NumberMalloc(
"MulRat");
420 for ( i = 0; i < dnumr; i++ ) xd[i] = a[i];
422 for ( i = 0; i < dden; i++ ) xf[i] = a[i];
423 xe = NumberMalloc(
"MulRat"); xg = NumberMalloc(
"MulRat");
424 for ( i = 0; i < enumr; i++ ) xe[i] = b[i];
426 for ( i = 0; i < eden; i++ ) xg[i] = b[i];
427 if ( Simplify(BHEAD xd,&dnumr,xg,&eden) ||
428 Simplify(BHEAD xe,&enumr,xf,&dden) ||
429 MulLong(xd,dnumr,xe,enumr,c,nc) ||
430 MulLong(xf,dden,xg,eden,xd,&dnumr) ) {
431 MLOCK(ErrorMessageLock);
433 MUNLOCK(ErrorMessageLock);
434 NumberFree(xd,
"MulRat"); NumberFree(xe,
"MulRat"); NumberFree(xf,
"MulRat"); NumberFree(xg,
"MulRat");
438 NumberFree(xd,
"MulRat"); NumberFree(xe,
"MulRat"); NumberFree(xf,
"MulRat"); NumberFree(xg,
"MulRat");
445 do { a0 = y % b1; y = b1; }
while ( ( b1 = a0 ) != 0 );
455 do { b0 = y % a1; y = a1; }
while ( ( a1 = b0 ) != 0 );
465 if ( xx & AWORDMASK ) {
468 c[1] = (UWORD)(xx >> BITSINWORD);
471 c[3] = (UWORD)(xx >> BITSINWORD);
476 if ( xx & AWORDMASK ) {
479 c[3] = (UWORD)(xx >> BITSINWORD);
488 if ( sgn < 0 ) *nc = -*nc;
500 WORD DivRat(
PHEAD UWORD *a, WORD na, UWORD *b, WORD nb, UWORD *c, WORD *nc)
506 MLOCK(ErrorMessageLock);
507 MesPrint(
"Rational division by zero");
508 MUNLOCK(ErrorMessageLock);
511 j = i = (nb >= 0)? nb: -nb;
513 do { xx = *xd; *xd++ = *xe; *xe++ = xx; }
while ( --j > 0 );
514 j = MulRat(BHEAD a,na,b,nb,c,nc);
516 do { xx = *xd; *xd++ = *xe; *xe++ = xx; }
while ( --i > 0 );
530 WORD Simplify(
PHEAD UWORD *a, WORD *na, UWORD *b, WORD *nb)
535 WORD n1,n2,n3,n4,sgn = 1;
537 UWORD *Siscrat5, *Siscrat6, *Siscrat7, *Siscrat8;
538 if ( *na < 0 ) { *na = -*na; sgn = -sgn; }
539 if ( *nb < 0 ) { *nb = -*nb; sgn = -sgn; }
540 Siscrat5 = NumberMalloc(
"Simplify"); Siscrat6 = NumberMalloc(
"Simplify");
541 Siscrat7 = NumberMalloc(
"Simplify"); Siscrat8 = NumberMalloc(
"Simplify");
542 x1 = Siscrat8; x2 = Siscrat7;
545 if ( DivLong(a,*na,b,*nb,x1,&n1,x2,&n2) )
goto SimpErr;
547 for ( i = 0; i < n1; i++ ) *a++ = *x1++;
555 do { y1 = y2 % y3; y2 = y3; }
while ( ( y3 = y1 ) != 0 );
556 if ( ( *x2 = y2 ) != 1 ) {
558 if ( DivLong(a,*na,x2,(WORD)1,x1,&n1,x3,&n3) )
goto SimpErr;
559 for ( i = 0; i < n1; i++ ) *a++ = *x1++;
565 else if ( *na >= GCDMAX && *nb >= GCDMAX ) {
566 n1 = i = *na; x3 = a;
568 x3 = b; n2 = i = *nb;
573 if ( GcdLong(BHEAD Siscrat8,n1,Siscrat7,n2,x2,&n3) )
goto SimpErr;
575 if ( *x2 != 1 || n2 != 1 ) {
576 DivLong(a,*na,x2,n2,x1,&n1,x4,&n4);
579 DivLong(b,*nb,x2,n2,x3,&n3,x4,&n4);
587 n1 = i = *na; x3 = a;
589 x3 = b; n2 = i = *nb;
591 x1 = Siscrat8; x2 = Siscrat7; x3 = Siscrat6;
593 if ( DivLong(x1,n1,x2,n2,x4,&n4,x3,&n3) )
goto SimpErr;
596 while ( ( *x1 = (*x2) % (*x3) ) != 0 ) { *x2 = *x3; *x3 = *x1; }
600 if ( DivLong(x2,n2,x3,n3,x4,&n4,x1,&n1) )
goto SimpErr;
601 if ( !n1 ) { x2 = x3; n2 = n3; x3 = Siscrat7;
break; }
603 while ( ( *x2 = (*x3) % (*x1) ) != 0 ) { *x3 = *x1; *x1 = *x2; }
608 if ( DivLong(x3,n3,x1,n1,x4,&n4,x2,&n2) )
goto SimpErr;
609 if ( !n2 ) { x2 = x1; n2 = n1; x1 = Siscrat7;
break; }
611 while ( ( *x3 = (*x1) % (*x2) ) != 0 ) { *x1 = *x2; *x2 = *x3; }
615 if ( *x2 != 1 || n2 != 1 ) {
616 DivLong(a,*na,x2,n2,x1,&n1,x4,&n4);
619 DivLong(b,*nb,x2,n2,x3,&n3,x4,&n4);
624 if ( sgn < 0 ) *na = -*na;
625 NumberFree(Siscrat5,
"Simplify"); NumberFree(Siscrat6,
"Simplify");
626 NumberFree(Siscrat7,
"Simplify"); NumberFree(Siscrat8,
"Simplify");
629 MLOCK(ErrorMessageLock);
631 MUNLOCK(ErrorMessageLock);
632 NumberFree(Siscrat5,
"Simplify"); NumberFree(Siscrat6,
"Simplify");
633 NumberFree(Siscrat7,
"Simplify"); NumberFree(Siscrat8,
"Simplify");
647 WORD AccumGCD(
PHEAD UWORD *a, WORD *na, UWORD *b, WORD nb)
650 WORD nna,nnb,numa,numb,dena,denb,numc,denc;
651 UWORD *GCDbuffer = NumberMalloc(
"AccumGCD");
653 nna = *na;
if ( nna < 0 ) nna = -nna; nna = (nna-1)/2;
654 nnb = nb;
if ( nnb < 0 ) nnb = -nnb; nnb = (nnb-1)/2;
655 UnPack(a,nna,&dena,&numa);
656 UnPack(b,nnb,&denb,&numb);
657 if ( GcdLong(BHEAD a,numa,b,numb,GCDbuffer,&numc) )
goto AccErr;
659 for ( i = 0; i < numa; i++ ) a[i] = GCDbuffer[i];
660 if ( GcdLong(BHEAD a+nna,dena,b+nnb,denb,GCDbuffer,&denc) )
goto AccErr;
662 for ( i = 0; i < dena; i++ ) a[i+nna] = GCDbuffer[i];
663 Pack(a,&numa,a+nna,dena);
665 NumberFree(GCDbuffer,
"AccumGCD");
668 MLOCK(ErrorMessageLock);
670 MUNLOCK(ErrorMessageLock);
671 NumberFree(GCDbuffer,
"AccumGCD");
680 int TakeRatRoot(UWORD *a, WORD *n, WORD power)
682 WORD numer,denom, nn;
683 if ( ( power & 1 ) == 0 && *n < 0 )
return(1);
684 if ( ABS(*n) == 1 && a[0] == 1 && a[1] == 1 )
return(0);
686 UnPack(a,nn,&denom,&numer);
687 if ( TakeLongRoot(a+nn,&denom,power) )
return(1);
688 if ( TakeLongRoot(a,&numer,power) )
return(1);
689 Pack(a,&numer,a+nn,denom);
690 if ( *n < 0 ) *n = -numer;
705 WORD AddLong(UWORD *a, WORD na, UWORD *b, WORD nb, UWORD *c, WORD *nc)
710 if ( AddPLon(a,-na,b,-nb,c,nc) )
return(-1);
724 else {
return( AddPLon(a,na,b,nb,c,nc) ); }
726 if ( ( res = BigLong(a,na,b,nb) ) > 0 ) {
727 SubPLon(a,na,b,nb,c,nc);
728 if ( sgn < 0 ) *nc = -*nc;
730 else if ( res < 0 ) {
731 SubPLon(b,nb,a,na,c,nc);
732 if ( sgn > 0 ) *nc = -*nc;
750 WORD AddPLon(UWORD *a, WORD na, UWORD *b, WORD nb, UWORD *c, WORD *nc)
752 UWORD carry = 0, e, nd = 0;
757 if ( e < *c ) carry = 0;
760 if ( e > *c ) carry = 1;
762 a++; b++; c++; nd++; na--; nb--;
767 if ( *c++ ) carry = 0;
775 if ( *c++ ) carry = 0;
782 if ( nd > (UWORD)AM.MaxTal ) {
783 MLOCK(ErrorMessageLock);
784 MesPrint(
"Overflow in addition");
785 MUNLOCK(ErrorMessageLock);
803 VOID SubPLon(UWORD *a, WORD na, UWORD *b, WORD nb, UWORD *c, WORD *nc)
805 UWORD borrow = 0, e, nd = 0;
809 *c = e - *b - borrow;
810 if ( *c < e ) borrow = 0;
814 if ( *c > e ) borrow = 1;
816 a++; b++; c++; na--; nb--; nd++;
820 if ( *a ) { *c++ = *a++ - 1; borrow = 0; }
821 else { *c++ = (UWORD)(-1); a++; }
826 while ( nd && !*--c ) { nd--; }
841 WORD MulLong(UWORD *a, WORD na, UWORD *b, WORD nb, UWORD *c, WORD *nc)
846 if ( !na || !nb ) { *nc = 0;
return(0); }
847 if ( na < 0 ) { na = -na; sgn = -sgn; }
848 if ( nb < 0 ) { nb = -nb; sgn = -sgn; }
850 if ( i > (UWORD)(AM.MaxTal+1) )
goto MulLov;
856 if (na > 3 && nb > 3) {
861 UWORD *DLscrat9 = NumberMalloc(
"MulLong"), *DLscratA = NumberMalloc(
"MulLong"), *DLscratB = NumberMalloc(
"MulLong");
862 #if ( GMPSPREAD != 1 ) 864 from = a; a = to = DLscrat9; j = na; NCOPY(to, from, j);
869 if ( (LONG)a & (
sizeof(mp_limb_t)-1) ) {
870 from = a; a = to = DLscrat9; j = na; NCOPY(to, from, j);
873 #if ( GMPSPREAD != 1 ) 875 from = b; b = to = DLscratA; j = nb; NCOPY(to, from, j);
880 if ( (LONG)b & (
sizeof(mp_limb_t)-1) ) {
881 from = b; b = to = DLscratA; j = nb; NCOPY(to, from, j);
884 if ( ( *nc > (WORD)i ) || ( (LONG)c & (LONG)(
sizeof(mp_limb_t)-1) ) ) {
889 mpn_mul((mp_ptr)ic, (mp_srcptr)b, nb/GMPSPREAD, (mp_srcptr)a, na/GMPSPREAD);
892 mpn_mul((mp_ptr)ic, (mp_srcptr)a, na/GMPSPREAD, (mp_srcptr)b, nb/GMPSPREAD);
894 while ( ic[i-1] == 0 ) i--;
901 j = *nc; NCOPY(c, ic, j);
903 if ( sgn < 0 ) *nc = -(*nc);
904 NumberFree(DLscrat9,
"MulLong"); NumberFree(DLscratA,
"MulLong"); NumberFree(DLscratB,
"MulLong");
911 do { *ic++ = 0; }
while ( --i > 0 );
919 t = (*ia++) * bb + t + *ic;
923 if ( t ) *ic = (UWORD)t;
924 }
while ( --nb > 0 );
926 if ( *nc > AM.MaxTal )
goto MulLov;
927 if ( sgn < 0 ) *nc = -(*nc);
930 MLOCK(ErrorMessageLock);
931 MesPrint(
"Overflow in Multiplication");
932 MUNLOCK(ErrorMessageLock);
944 WORD BigLong(UWORD *a, WORD na, UWORD *b, WORD nb)
948 while ( na && !*--a ) na--;
949 while ( nb && !*--b ) nb--;
950 if ( nb < na )
return(1);
951 if ( nb > na )
return(-1);
952 while ( --na >= 0 ) {
953 if ( *a > *b )
return(1);
954 else if ( *b > *a )
return(-1);
972 WORD DivLong(UWORD *a, WORD na, UWORD *b, WORD nb, UWORD *c,
973 WORD *nc, UWORD *d, WORD *nd)
975 WORD sgna = 1, sgnb = 1, ne, nf, ng, nh;
979 UWORD *e, *f, *ff, *g, norm, estim;
981 UWORD *DLscrat9, *DLscratA, *DLscratB, *DLscratC;
985 MLOCK(ErrorMessageLock);
986 MesPrint(
"Division by zero");
987 MUNLOCK(ErrorMessageLock);
990 if ( !na ) { *nc = *nd = 0;
return(0); }
991 if ( na < 0 ) { sgna = -sgna; na = -na; }
992 if ( nb < 0 ) { sgnb = -sgnb; nb = -nb; }
994 for ( i = 0; i < na; i++ ) *d++ = *a++;
998 else if ( nb == na && ( i = BigLong(b,nb,a,na) ) >= 0 ) {
1000 for ( i = 0; i < na; i++ ) *d++ = *a++;
1010 else if ( nb == 1 ) {
1012 for ( i = 0; i < na; i++ ) *c++ = *a++;
1023 while ( --ni >= 0 ) {
1031 if ( ( *d = (UWORD)t ) == 0 ) *nd = 0;
1032 if ( !*(c+na-1) ) (*nc)--;
1047 if ( na > 4 && nb > 3 ) {
1048 UWORD *ic, *id, *to, *from;
1050 DLscrat9 = NumberMalloc(
"DivLong"); DLscratA = NumberMalloc(
"DivLong");
1051 DLscratB = NumberMalloc(
"DivLong"); DLscratC = NumberMalloc(
"DivLong");
1053 #if ( GMPSPREAD != 1 ) 1055 from = a; a = to = DLscrat9; i = na; NCOPY(to, from, i);
1059 if ( (LONG)a & (
sizeof(mp_limb_t)-1) ) {
1060 from = a; a = to = DLscrat9; i = na; NCOPY(to, from, i);
1063 #if ( GMPSPREAD != 1 ) 1065 from = b; b = to = DLscratA; i = nb; NCOPY(to, from, i);
1069 if ( ( (LONG)b & (
sizeof(mp_limb_t)-1) ) != 0 ) {
1070 from = b; b = to = DLscratA; i = nb; NCOPY(to, from, i);
1072 if ( ( (LONG)c & (
sizeof(mp_limb_t)-1) ) != 0 ) ic = DLscratB;
1075 if ( ( (LONG)d & (
sizeof(mp_limb_t)-1) ) != 0 )
id = DLscratC;
1077 mpn_tdiv_qr((mp_limb_t *)ic,(mp_limb_t *)
id,(mp_size_t)0,
1078 (
const mp_limb_t *)a,(mp_size_t)(na/GMPSPREAD),
1079 (
const mp_limb_t *)b,(mp_size_t)(nb/GMPSPREAD));
1080 while ( j >= 0 && ic[j] == 0 ) j--;
1082 if ( c != ic ) { NCOPY(c,ic,j); }
1084 while ( j >= 0 &&
id[j] == 0 ) j--;
1086 if ( d !=
id ) { NCOPY(d,
id,j); }
1087 if ( sgna < 0 ) { *nc = -(*nc); *nd = -(*nd); }
1088 if ( sgnb < 0 ) { *nc = -(*nc); }
1089 NumberFree(DLscrat9,
"DivLong"); NumberFree(DLscratA,
"DivLong");
1090 NumberFree(DLscratB,
"DivLong"); NumberFree(DLscratC,
"DivLong");
1099 e = NumberMalloc(
"DivLong"); f = NumberMalloc(
"DivLong"); g = NumberMalloc(
"DivLong");
1100 if ( b[nb-1] == (FULLMAX-1) ) norm = 1;
1102 norm = (UWORD)(((ULONG)FULLMAX) / (ULONG)((b[nb-1]+1L)));
1105 if ( MulLong(b,nb,&norm,1,e,&ne) ||
1106 MulLong(a,na,&norm,1,f,&nf) ) {
1107 NumberFree(e,
"DivLong"); NumberFree(f,
"DivLong"); NumberFree(g,
"DivLong");
1110 if ( BigLong(f+nf-ne,ne,e,ne) >= 0 ) {
1111 SubPLon(f+nf-ne,ne,e,ne,f+nf-ne,&nh);
1120 w2 = c; i = *nc;
do { *w2++ = 0; }
while ( --i > 0 );
1123 esthelp = (RLONG)(e[ne-1]) + 1L;
1124 while ( nf >= ne ) {
1125 if ( (WORD)esthelp == 0 ) {
1126 estim = (WORD)(((((RLONG)(f[nf]))<<BITSINWORD)+f[nf-1])>>BITSINWORD);
1129 estim = (WORD)(((((RLONG)(f[nf]))<<BITSINWORD)+f[nf-1])/esthelp);
1133 MulLong(e,ne,&estim,1,g,&ng);
1134 nh = ne + 1;
if ( !f[ni+ne] ) nh--;
1135 SubPLon(f+ni,nh,g,ng,f+ni,&nh);
1138 w2 = f+ni+ne; nh = ne+1;
1139 while ( ( nh > 0 ) && !*w2 ) { nh--; w2--; }
1141 if ( BigLong(f+ni,nh,e,ne) >= 0 ) {
1143 SubPLon(f+ni,nh,e,ne,f+ni,&nh);
1144 if ( BigLong(f+ni,nh,e,ne) >= 0 ) {
1146 SubPLon(f+ni,nh,e,ne,f+ni,&nh);
1147 if ( BigLong(f+ni,nh,e,ne) >= 0 ) {
1148 MLOCK(ErrorMessageLock);
1149 MesPrint(
"Problems in DivLong");
1153 while ( --i >= 0 ) { TalToLine((UWORD)(*a++)); TokenToLine((UBYTE *)
" "); }
1156 while ( --i >= 0 ) { TalToLine((UWORD)(*b++)); TokenToLine((UBYTE *)
" "); }
1159 MUNLOCK(ErrorMessageLock);
1160 NumberFree(e,
"DivLong"); NumberFree(f,
"DivLong"); NumberFree(g,
"DivLong");
1175 *nd = i = nh; ff = f;
1184 while ( --ni >= 0 ) {
1193 MLOCK(ErrorMessageLock);
1194 MesPrint(
"Error in DivLong");
1195 MUNLOCK(ErrorMessageLock);
1196 NumberFree(e,
"DivLong"); NumberFree(f,
"DivLong"); NumberFree(g,
"DivLong");
1199 if ( !*(d+nh-1) ) (*nd)--;
1203 NumberFree(e,
"DivLong"); NumberFree(f,
"DivLong"); NumberFree(g,
"DivLong");
1205 if ( sgna < 0 ) { *nc = -(*nc); *nd = -(*nd); }
1206 if ( sgnb < 0 ) { *nc = -(*nc); }
1218 WORD RaisPow(
PHEAD UWORD *a, WORD *na, UWORD b)
1225 nmod = ABS(AN.ncmod);
1226 if ( !*na || ( ( *na == 1 ) && ( *a == 1 ) ) )
return(0);
1227 if ( !b ) { *na=1; *a=1;
return(0); }
1228 is = NumberMalloc(
"RaisPow");
1229 it = NumberMalloc(
"RaisPow");
1230 for ( i = 0; i < ABS(*na); i++ ) is[i] = a[i];
1233 for ( i = 0; i < BITSINWORD; i++ ) {
1239 while ( --i >= 0 ) {
1241 if(MulLong(is,ns,is,ns,it,&nt))
goto RaisOvl;
1243 if ( MulLong(it,nt,a,*na,is,&ns) )
goto RaisOvl;
1246 iu = is; is = it; it = iu;
1247 nu = ns; ns = nt; nt = nu;
1250 if ( DivLong(is,ns,(UWORD *)AC.cmod,nmod,it,&nt,is,&ns) )
goto RaisOvl;
1253 if ( ( nmod != 0 ) && ( ( AC.modmode & POSNEG ) != 0 ) ) {
1256 if ( ( *na = i = ns ) != 0 ) { iss = is; i=ABS(i); NCOPY(a,iss,i); }
1257 NumberFree(is,
"RaisPow"); NumberFree(it,
"RaisPow");
1260 MLOCK(ErrorMessageLock);
1262 MUNLOCK(ErrorMessageLock);
1263 NumberFree(is,
"RaisPow"); NumberFree(it,
"RaisPow");
1289 WORD new_small_power_maxx, new_small_power_maxn, ID;
1290 WORD *new_small_power_n;
1291 UWORD **new_small_power;
1294 if (x>=AT.small_power_maxx || n>=AT.small_power_maxn) {
1296 new_small_power_maxx = AT.small_power_maxx;
1297 if (x>=AT.small_power_maxx)
1298 new_small_power_maxx = MaX(2*AT.small_power_maxx, x+1);
1300 new_small_power_maxn = AT.small_power_maxn;
1301 if (n>=AT.small_power_maxn)
1302 new_small_power_maxn = MaX(2*AT.small_power_maxn, n+1);
1304 new_small_power_n = (WORD*) Malloc1(new_small_power_maxx*new_small_power_maxn*
sizeof(WORD),
"RaisPowCached");
1305 new_small_power = (UWORD **) Malloc1(new_small_power_maxx*new_small_power_maxn*
sizeof(UWORD *),
"RaisPowCached");
1307 for (i=0; i<new_small_power_maxx * new_small_power_maxn; i++) {
1308 new_small_power_n[i] = 0;
1309 new_small_power [i] = NULL;
1312 for (i=0; i<AT.small_power_maxx; i++)
1313 for (j=0; j<AT.small_power_maxn; j++) {
1314 new_small_power_n[i*new_small_power_maxn+j] = AT.small_power_n[i*AT.small_power_maxn+j];
1315 new_small_power [i*new_small_power_maxn+j] = AT.small_power [i*AT.small_power_maxn+j];
1318 if (AT.small_power_n != NULL) {
1319 M_free(AT.small_power_n,
"RaisPowCached");
1320 M_free(AT.small_power,
"RaisPowCached");
1323 AT.small_power_maxx = new_small_power_maxx;
1324 AT.small_power_maxn = new_small_power_maxn;
1325 AT.small_power_n = new_small_power_n;
1326 AT.small_power = new_small_power;
1330 ID = x * AT.small_power_maxn + n;
1332 if (AT.small_power[ID] == NULL) {
1333 AT.small_power[ID] = NumberMalloc(
"RaisPowCached");
1334 AT.small_power_n[ID] = 1;
1335 AT.small_power[ID][0] = x;
1336 RaisPow(BHEAD AT.small_power[ID],&AT.small_power_n[ID],n);
1340 *c = AT.small_power[ID];
1341 *nc = AT.small_power_n[ID];
1350 WORD RaisPowMod (WORD x, WORD n, WORD m) {
1353 if (n&1) { y*=z; y%=m; }
1373 if ( AC.halfmod == 0 ) {
1374 LOCK(AC.halfmodlock);
1375 if ( AC.halfmod == 0 ) {
1376 UWORD two[1],remain[1];
1379 AC.halfmod = (UWORD *)Malloc1((ABS(AC.ncmod))*
sizeof(UWORD),
"halfmod");
1380 DivLong((UWORD *)AC.cmod,(ABS(AC.ncmod)),two,1
1381 ,(UWORD *)AC.halfmod,&(AC.nhalfmod),remain,&dummy);
1383 UNLOCK(AC.halfmodlock);
1386 if ( BigLong(a,n,AC.halfmod,AC.nhalfmod) > 0 ) {
1387 SubPLon((UWORD *)AC.cmod,(ABS(AC.ncmod)),a,n,a,&n);
1388 if ( *na > 0 ) { *na = -n; }
1409 WORD n = AC.cmod[0], i, inv2;
1410 if ( AC.ncmod != 1 )
return(1);
1411 if ( AC.modinverses == 0 ) {
1412 LOCK(AC.halfmodlock);
1413 if ( AC.modinverses == 0 ) {
1414 AC.modinverses = (UWORD *)Malloc1(n*
sizeof(UWORD),
"modinverses");
1415 AC.modinverses[0] = 0;
1416 AC.modinverses[1] = 1;
1417 for ( i = 2; i < n; i++ ) {
1419 (WORD *)(&(AC.modinverses[i])),&inv2) ) {
1424 UNLOCK(AC.halfmodlock);
1447 WORD x = m1, y, c, d = m2;
1448 if ( x < 1 || d <= 1 )
goto somethingwrong;
1453 if ( y == 0 )
break;
1454 a3 = a1-c*a2; a1 = a2; a2 = a3;
1455 b3 = b1-c*b2; b1 = b2; b2 = b3;
1458 if ( x != 1 )
goto somethingwrong;
1459 if ( a2 < 0 ) a2 += m2;
1460 if ( b2 < 0 ) b2 += m1;
1461 if (im1!=NULL) *im1 = a2;
1462 if (im2!=NULL) *im2 = b2;
1465 MLOCK(ErrorMessageLock);
1466 MesPrint(
"Error trying to determine inverses in GetModInverses");
1467 MUNLOCK(ErrorMessageLock);
1475 int GetLongModInverses(
PHEAD UWORD *a, WORD na, UWORD *b, WORD nb, UWORD *ia, WORD *nia, UWORD *ib, WORD *nib) {
1477 UWORD *s, *t, *sa, *sb, *ta, *tb, *x, *y, *swap1;
1478 WORD ns, nt, nsa, nsb, nta, ntb, nx, ny, swap2;
1480 s = NumberMalloc(
"GetLongModInverses");
1482 WCOPY(s, a, ABS(ns));
1484 t = NumberMalloc(
"GetLongModInverses");
1486 WCOPY(t, b, ABS(nt));
1488 sa = NumberMalloc(
"GetLongModInverses");
1492 sb = NumberMalloc(
"GetLongModInverses");
1495 ta = NumberMalloc(
"GetLongModInverses");
1498 tb = NumberMalloc(
"GetLongModInverses");
1502 x = NumberMalloc(
"GetLongModInverses");
1503 y = NumberMalloc(
"GetLongModInverses");
1506 DivLong(s,ns,t,nt,x,&nx,y,&ny);
1507 swap1=s; s=y; y=swap1;
1509 MulLong(x,nx,ta,nta,y,&ny);
1510 AddLong(sa,nsa,y,-ny,sa,&nsa);
1511 MulLong(x,nx,tb,ntb,y,&ny);
1512 AddLong(sb,nsb,y,-ny,sb,&nsb);
1514 swap1=s; s=t; t=swap1;
1515 swap2=ns; ns=nt; nt=swap2;
1516 swap1=sa; sa=ta; ta=swap1;
1517 swap2=nsa; nsa=nta; nta=swap2;
1518 swap1=sb; sb=tb; tb=swap1;
1519 swap2=nsb; nsb=ntb; ntb=swap2;
1524 WCOPY(ia,sa,ABS(*nia));
1529 WCOPY(ib,sb,ABS(*nib));
1532 NumberFree(s,
"GetLongModInverses");
1533 NumberFree(t,
"GetLongModInverses");
1534 NumberFree(sa,
"GetLongModInverses");
1535 NumberFree(sb,
"GetLongModInverses");
1536 NumberFree(ta,
"GetLongModInverses");
1537 NumberFree(tb,
"GetLongModInverses");
1538 NumberFree(x,
"GetLongModInverses");
1539 NumberFree(y,
"GetLongModInverses");
1552 WORD Product(UWORD *a, WORD *na, WORD b)
1556 if ( *na < 0 ) { *na = -(*na); sgn = -sgn; }
1557 if ( b < 0 ) { b = -b; sgn = -sgn; }
1560 for ( i = 0; i < *na; i++ ) {
1566 if ( ++(*na) > AM.MaxTal ) {
1567 MLOCK(ErrorMessageLock);
1568 MesPrint(
"Overflow in Product");
1569 MUNLOCK(ErrorMessageLock);
1574 if ( sgn < 0 ) *na = -(*na);
1587 UWORD Quotient(UWORD *a, WORD *na, WORD b)
1591 if ( ( i = *na ) < 0 ) { sgn = -1; i = -i; }
1592 if ( b < 0 ) { b = -b; sgn = -sgn; }
1594 if ( ( *a /= (UWORD)b ) == 0 ) *na = 0;
1595 if ( sgn < 0 ) *na = -*na;
1602 while ( --i >= 0 ) {
1612 if ( sgn < 0 ) j = -j;
1626 WORD Remain10(UWORD *a, WORD *na)
1634 while ( --i >= 0 ) {
1638 if ( i > 0 ) t <<= BITSINWORD;
1640 if ( ( *na > 0 ) && !a[*na-1] ) (*na)--;
1653 WORD Remain4(UWORD *a, WORD *na)
1661 while ( --i >= 0 ) {
1663 *b-- = u = t / 10000;
1665 if ( i > 0 ) t <<= BITSINWORD;
1667 if ( ( *na > 0 ) && !a[*na-1] ) (*na)--;
1679 VOID PrtLong(UWORD *a, WORD na, UBYTE *s)
1692 b = NumberMalloc(
"PrtLong");
1694 i = na;
while ( --i >= 0 ) *bb++ = *a++;
1700 *sa++ = (UBYTE)(
'0' + (q%10));
1702 *sa++ = (UBYTE)(
'0' + (q%10));
1704 *sa++ = (UBYTE)(
'0' + (q%10));
1706 *sa++ = (UBYTE)(
'0' + (q%10));
1708 while ( sa[-1] ==
'0' ) sa--;
1712 while ( sa > sb ) { c = *sa; *sa = *sb; *sb = c; sa--; sb++; }
1717 q = Remain10(a,&na);
1718 *sa++ = (UBYTE)(
'0' + q);
1723 while ( sa > sb ) { c = *sa; *sa = *sb; *sb = c; sa--; sb++; }
1727 NumberFree(b,
"PrtLong");
1741 WORD GetLong(UBYTE *s, UWORD *a, WORD *na)
1754 UWORD digit, x = 0, y = 0;
1757 while ( FG.cTable[*s] == 1 ) {
1759 if ( FG.cTable[*s] != 1 ) { y = 10;
break; }
1760 x = 10*x + *s++ -
'0';
1761 if ( FG.cTable[*s] != 1 ) { y = 100;
break; }
1762 x = 10*x + *s++ -
'0';
1763 if ( FG.cTable[*s] != 1 ) { y = 1000;
break; }
1764 x = 10*x + *s++ -
'0';
1765 if ( *na && Product(a,na,(WORD)10000) )
return(-1);
1766 if ( ( digit = x ) != 0 && AddLong(a,*na,&digit,(WORD)1,a,na) )
1771 if ( *na && Product(a,na,(WORD)y) )
return(-1);
1772 if ( ( digit = x ) != 0 && AddLong(a,*na,&digit,(WORD)1,a,na) )
1790 #define Convert(ia,aa,naa) \ 1791 if ( (LONG)ia < 0 ) { \ 1792 ia = (ULONG)(-(LONG)ia); \ 1794 if ( ( aa[1] = ia >> BITSINWORD ) != 0 ) naa = -2; \ 1797 else if ( ia == 0 ) { aa[0] = 0; naa = 0; } \ 1800 if ( ( aa[1] = ia >> BITSINWORD ) != 0 ) naa = 2; \ 1804 VOID GCD(UWORD *a, WORD na, UWORD *b, WORD nb, UWORD *c, WORD *nc)
1806 int ja = 0, jb = 0, j;
1808 UWORD *x1, *x2, *x3;
1810 ULONG ia,ib,ic,id,u,v,w,q,T;
1815 while ( a[0] == 0 ) { na--; ja++; a++; }
1816 while ( b[0] == 0 ) { nb--; jb++; b++; }
1817 if ( ja > jb ) ja = jb;
1820 do { *c++ = 0; }
while ( --j > 0 );
1826 jb = na; na = nb; nb = jb;
1828 r = a; a = b; b = r;
1830 else if ( na == nb ) {
1834 while ( --j >= 0 ) {
1835 if ( *--r > *--t )
break;
1836 if ( *r < *t )
goto exch;
1867 r = x1 = NumberMalloc(
"GCD"); t = x2 = NumberMalloc(
"GCD"); x3 = NumberMalloc(
"GCD");
1877 DivLong(x1,na,x2,nb,c,nc,x3,&nd);
1878 if ( nd == 0 ) { b = x2;
goto out; }
1879 t = x1; x1 = x2; x2 = x3; x3 = t; na = nb; nb = nd;
1880 if ( na == 2 )
break;
1886 v = x1[0] + ( ((ULONG)x1[1]) << BITSINWORD );
1888 if ( nb == 2 ) w += ((ULONG)x2[1]) << BITSINWORD;
1892 do { u = v%w; v = w; w = u; }
while ( w );
1895 if ( ( c[1] = (UWORD)(v >> BITSINWORD) ) != 0 ) *nc = 2+ja;
1897 NumberFree(x1,
"GCD"); NumberFree(x2,
"GCD"); NumberFree(x3,
"GCD");
1902 ui = x1[0]; uj = x2[0];
1904 ui = (UWORD)GCD2((ULONG)ui,(ULONG)uj);
1906 do { nd = ui%uj; ui = uj; uj = nd; }
while ( nd );
1910 NumberFree(x1,
"GCD"); NumberFree(x2,
"GCD"); NumberFree(x3,
"GCD");
1913 ia = 1; ib = 0; ic = 0;
id = 1;
1914 u = ( ((ULONG)x1[na-1]) << BITSINWORD ) + x1[na-2];
1915 v = ( ((ULONG)x2[nb-1]) << BITSINWORD ) + x2[nb-2];
1917 while ( v+ic != 0 && v+
id != 0 &&
1918 ( q = (u+ia)/(v+ic) ) == (u+ib)/(v+
id) ) {
1919 T = ia-q*ic; ia = ic; ic = T;
1920 T = ib-q*id; ib = id;
id = T;
1921 T = u - q*v; u = v; v = T;
1923 if ( ib == 0 )
goto toobad;
1926 MulLong(x1,na,aa,naa,x3,&nd);
1927 MulLong(x2,nb,bb,nbb,c,nc);
1928 AddLong(x3,nd,c,*nc,c,nc);
1931 MulLong(x1,na,aa,naa,x3,&nd);
1932 t = c; na = j = *nc; r = x1;
1934 MulLong(x2,nb,bb,nbb,c,nc);
1935 AddLong(x3,nd,c,*nc,x2,&nb);
1956 WORD GcdLong(
PHEAD UWORD *a, WORD na, UWORD *b, WORD nb, UWORD *c, WORD *nc)
1961 MLOCK(ErrorMessageLock);
1962 MesPrint(
"Cannot take gcd");
1963 MUNLOCK(ErrorMessageLock);
1979 if ( na < 0 ) na = -na;
1980 if ( nb < 0 ) nb = -nb;
1981 if ( na == 1 && nb == 1 ) {
1983 *c = (UWORD)GCD2((ULONG)*a,(ULONG)*b);
1988 do { z = x % y; x = y; }
while ( ( y = z ) != 0 );
1993 else if ( na <= 2 && nb <= 2 ) {
1995 if ( na == 2 ) { lx = (((RLONG)(a[1]))<<BITSINWORD) + *a; }
1997 if ( nb == 2 ) { ly = (((RLONG)(b[1]))<<BITSINWORD) + *b; }
2003 do { lz = lx % ly; lx = ly; }
while ( ( ly = lz ) != 0 );
2006 if ( lx < ly ) { lz = lx; lx = ly; ly = lz; }
2008 lz = lx % ly; lx = ly;
2009 }
while ( ( ly = lz ) != 0 && ( lx & AWORDMASK ) != 0 );
2011 do { *c = ((UWORD)lx)%((UWORD)ly); lx = ly; }
while ( ( ly = *c ) != 0 );
2019 if ( ( *c = (UWORD)(lx >> BITSINWORD) ) != 0 ) *nc = 2;
2025 GCD(a,na,b,nb,c,nc);
2028 UWORD *x3,*x1,*x2, *GLscrat7, *GLscrat8;
2031 x1 = c; x3 = a; n1 = i = na;
2033 GLscrat7 = NumberMalloc(
"GcdLong"); GLscrat8 = NumberMalloc(
"GcdLong");
2034 x2 = GLscrat8; x3 = b; n2 = i = nb;
2037 while ( x1[0] == 0 ) { i += BITSINWORD; x1++; n1--; }
2038 while ( ( x1[0] & 1 ) == 0 ) { i++; SCHUIF(x1,n1) }
2039 x2 = GLscrat8; j = 0;
2040 while ( x2[0] == 0 ) { j += BITSINWORD; x2++; n2--; }
2041 while ( ( x2[0] & 1 ) == 0 ) { j++; SCHUIF(x2,n2) }
2046 SubPLon(x1,n1,x2,n2,x1,&n3);
2050 n1 = i = n2; NCOPY(x1,x2,i);
2053 while ( ( x1[0] & 1 ) == 0 ) SCHUIF(x1,n1)
2055 if ( DivLong(x2,n2,x1,n1,GLscrat7,&n3,x2,&n4) )
goto GcdErr;
2058 i = n1; x2 = c; NCOPY(x2,x1,i);
2062 *c = (UWORD)GCD2((ULONG)x1[0],(ULONG)x2[0]);
2068 do { z = x % y; x = y; }
while ( ( y = z ) != 0 );
2076 else if ( n1 < n2 ) {
2078 SubPLon(x2,n2,x1,n1,x2,&n3);
2081 i = n1; x2 = c; NCOPY(x2,x1,i);
2084 while ( ( x2[0] & 1 ) == 0 ) SCHUIF(x2,n2)
2086 if ( DivLong(x1,n1,x2,n2,GLscrat7,&n3,x1,&n4) )
goto GcdErr;
2090 n1 = i = n2; NCOPY(x1,x2,i);
2094 *c = (UWORD)GCD2((ULONG)x2[0],(ULONG)x1[0]);
2100 do { z = x % y; x = y; }
while ( ( y = z ) != 0 );
2109 for ( i = n1-1; i >= 0; i-- ) {
2110 if ( x1[i] > x2[i] )
goto firstbig;
2111 else if ( x1[i] < x2[i] )
goto lastbig;
2113 i = n1; x2 = c; NCOPY(x2,x1,i);
2121 while ( j >= BITSINWORD ) {
2122 for ( i = n1; i > 0; i-- ) x1[i] = x1[i-1];
2128 for ( i = 0; i < n1; i++ ) {
2129 a1 = x1[i]; a1 <<= j;
2139 NumberFree(GLscrat7,
"GcdLong"); NumberFree(GLscrat8,
"GcdLong");
2141 UWORD *x1,*x2,*x3,*x4,*c1,*c2;
2143 x1 = c; x3 = a; n1 = i = na;
2145 x1 = c; c1 = x2 = NumberMalloc(
"GcdLong"); x3 = NumberMalloc(
"GcdLong"); x4 = NumberMalloc(
"GcdLong");
2146 c2 = b; n2 = i = nb;
2149 if ( DivLong(x1,n1,x2,n2,x4,&n4,x3,&n3) )
goto GcdErr;
2150 if ( !n3 ) { x1 = x2; n1 = n2;
break; }
2151 if ( DivLong(x2,n2,x3,n3,x4,&n4,x1,&n1) )
goto GcdErr;
2152 if ( !n1 ) { x1 = x3; n1 = n3;
break; }
2153 if ( DivLong(x3,n3,x1,n1,x4,&n4,x2,&n2) )
goto GcdErr;
2156 NumberFree(x2,
"GcdLong"); NumberFree(x3,
"GcdLong"); NumberFree(x4,
"GcdLong");
2162 NumberFree(x2,
"GcdLong"); NumberFree(x3,
"GcdLong"); NumberFree(x4,
"GcdLong");
2168 MLOCK(ErrorMessageLock);
2170 MUNLOCK(ErrorMessageLock);
2241 WORD GcdLong(
PHEAD UWORD *a, WORD na, UWORD *b, WORD nb, UWORD *c, WORD *nc)
2245 UWORD *x1,*x2,*x3,*x4,*x5,*d;
2246 UWORD *GLscrat6, *GLscrat7, *GLscrat8, *GLscrat9, *GLscrat10;
2247 WORD n1,n2,n3,n4,n5,i;
2249 LONG ma1, ma2, mb1, mb2, mc1, mc2, m;
2252 MLOCK(ErrorMessageLock);
2253 MesPrint(
"Cannot take gcd");
2254 MUNLOCK(ErrorMessageLock);
2270 if ( na < 0 ) na = -na;
2271 if ( nb < 0 ) nb = -nb;
2276 if ( na > 3 && nb > 3 ) {
2278 mp_limb_t *upa, *upb, *upc, xx;
2279 UWORD *uw, *u1, *u2;
2280 unsigned int tcounta, tcountb, tcounta1, tcountb1;
2281 mp_size_t ana, anb, anc;
2283 u1 = uw = NumberMalloc(
"GcdLong");
2284 upa = (mp_limb_t *)u1;
2285 ana = na; tcounta1 = 0;
2286 while ( a[0] == 0 ) { a++; ana--; tcounta1++; }
2287 for ( ii = 0; ii < ana; ii++ ) { *uw++ = *a++; }
2288 if ( ( ana & 1 ) != 0 ) { *uw = 0; ana++; }
2291 u2 = uw = NumberMalloc(
"GcdLong");
2292 upb = (mp_limb_t *)u2;
2293 anb = nb; tcountb1 = 0;
2294 while ( b[0] == 0 ) { b++; anb--; tcountb1++; }
2295 for ( ii = 0; ii < anb; ii++ ) { *uw++ = *b++; }
2296 if ( ( anb & 1 ) != 0 ) { *uw = 0; anb++; }
2299 xx = upa[0]; tcounta = 0;
2300 while ( ( xx & 15 ) == 0 ) { tcounta += 4; xx >>= 4; }
2301 while ( ( xx & 1 ) == 0 ) { tcounta += 1; xx >>= 1; }
2302 xx = upb[0]; tcountb = 0;
2303 while ( ( xx & 15 ) == 0 ) { tcountb += 4; xx >>= 4; }
2304 while ( ( xx & 1 ) == 0 ) { tcountb += 1; xx >>= 1; }
2307 mpn_rshift(upa,upa,ana,tcounta);
2308 if ( upa[ana-1] == 0 ) ana--;
2311 mpn_rshift(upb,upb,anb,tcountb);
2312 if ( upb[anb-1] == 0 ) anb--;
2315 upc = (mp_limb_t *)(NumberMalloc(
"GcdLong"));
2316 if ( ( ana > anb ) || ( ( ana == anb ) && ( upa[ana-1] >= upb[ana-1] ) ) ) {
2317 anc = mpn_gcd(upc,upa,ana,upb,anb);
2320 anc = mpn_gcd(upc,upb,anb,upa,ana);
2323 tcounta = tcounta1*BITSINWORD + tcounta;
2324 tcountb = tcountb1*BITSINWORD + tcountb;
2325 if ( tcountb > tcounta ) tcountb = tcounta;
2326 tcounta = tcountb/BITSINWORD;
2327 tcountb = tcountb%BITSINWORD;
2330 xx = mpn_lshift(upc,upc,anc,tcountb);
2331 if ( xx ) { upc[anc] = xx; anc++; }
2334 uw = (UWORD *)upc; anc *= 2;
2335 while ( uw[anc-1] == 0 ) anc--;
2336 for ( ii = 0; ii < (int)tcounta; ii++ ) *c++ = 0;
2337 for ( ii = 0; ii < anc; ii++ ) *c++ = *uw++;
2338 *nc = anc + tcounta;
2339 NumberFree(u1,
"GcdLong"); NumberFree(u2,
"GcdLong"); NumberFree((UWORD *)(upc),
"GcdLong");
2349 if ( na == 1 && nb == 1 ) {
2352 do { z = x % y; x = y; }
while ( ( y = z ) != 0 );
2357 else if ( na <= 2 && nb <= 2 ) {
2358 if ( na == 2 ) { lx = (((RLONG)(a[1]))<<BITSINWORD) + *a; }
2360 if ( nb == 2 ) { ly = (((RLONG)(b[1]))<<BITSINWORD) + *b; }
2362 if ( lx < ly ) { lz = lx; lx = ly; ly = lz; }
2363 #if ( BITSINWORD == 16 ) 2365 lz = lx % ly; lx = ly;
2366 }
while ( ( ly = lz ) != 0 );
2369 lz = lx % ly; lx = ly;
2370 }
while ( ( ly = lz ) != 0 && ( lx & AWORDMASK ) != 0 );
2372 x = (UWORD)lx; y = (UWORD)ly;
2373 do { *c = x % y; x = y; }
while ( ( y = *c ) != 0 );
2381 if ( ( *c = (UWORD)(lx >> BITSINWORD) ) != 0 ) *nc = 2;
2389 GLscrat6 = NumberMalloc(
"GcdLong"); GLscrat7 = NumberMalloc(
"GcdLong");
2390 GLscrat8 = NumberMalloc(
"GcdLong");
2391 GLscrat9 = NumberMalloc(
"GcdLong"); GLscrat10 = NumberMalloc(
"GcdLong");
2396 if ( na == 1 && nb == 1 ) {
2399 do { z = x % y; x = y; }
while ( ( y = z ) != 0 );
2403 else if ( na <= 2 && nb <= 2 ) {
2404 if ( na == 2 ) { lx = (((RLONG)(a[1]))<<BITSINWORD) + *a; }
2406 if ( nb == 2 ) { ly = (((RLONG)(b[1]))<<BITSINWORD) + *b; }
2408 if ( lx < ly ) { lz = lx; lx = ly; ly = lz; }
2409 #if ( BITSINWORD == 16 ) 2411 lz = lx % ly; lx = ly;
2412 }
while ( ( ly = lz ) != 0 );
2415 lz = lx % ly; lx = ly;
2416 }
while ( ( ly = lz ) != 0 && ( lx & AWORDMASK ) != 0 );
2418 x = (UWORD)lx; y = (UWORD)ly;
2419 do { *c = x % y; x = y; }
while ( ( y = *c ) != 0 );
2427 if ( ( *c = (UWORD)(lx >> BITSINWORD) ) != 0 ) *nc = 2;
2435 else if ( na < GCDMAX || nb < GCDMAX || na != nb ) {
2437 x2 = GLscrat8; x3 = a; n2 = i = na;
2439 x1 = c; x3 = b; n1 = i = nb;
2443 x1 = c; x3 = a; n1 = i = na;
2445 x2 = GLscrat8; x3 = b; n2 = i = nb;
2448 x1 = c; x2 = GLscrat8; x3 = GLscrat7; x4 = GLscrat6;
2450 if ( DivLong(x1,n1,x2,n2,x4,&n4,x3,&n3) )
goto GcdErr;
2451 if ( !n3 ) { x1 = x2; n1 = n2;
break; }
2452 if ( n2 <= 2 ) { a = x2; b = x3; na = n2; nb = n3;
goto restart; }
2453 if ( n3 >= GCDMAX && n2 == n3 ) {
2454 a = GLscrat9; b = GLscrat10; na = n2; nb = n3;
2455 for ( i = 0; i < na; i++ ) a[i] = x2[i];
2456 for ( i = 0; i < nb; i++ ) b[i] = x3[i];
2459 if ( DivLong(x2,n2,x3,n3,x4,&n4,x1,&n1) )
goto GcdErr;
2460 if ( !n1 ) { x1 = x3; n1 = n3;
break; }
2461 if ( n3 <= 2 ) { a = x3; b = x1; na = n3; nb = n1;
goto restart; }
2462 if ( n1 >= GCDMAX && n1 == n3 ) {
2463 a = GLscrat9; b = GLscrat10; na = n3; nb = n1;
2464 for ( i = 0; i < na; i++ ) a[i] = x3[i];
2465 for ( i = 0; i < nb; i++ ) b[i] = x1[i];
2468 if ( DivLong(x3,n3,x1,n1,x4,&n4,x2,&n2) )
goto GcdErr;
2469 if ( !n2 ) { *nc = n1;
goto normalend; }
2470 if ( n1 <= 2 ) { a = x1; b = x2; na = n1; nb = n2;
goto restart; }
2471 if ( n2 >= GCDMAX && n2 == n1 ) {
2472 a = GLscrat9; b = GLscrat10; na = n1; nb = n2;
2473 for ( i = 0; i < na; i++ ) a[i] = x1[i];
2474 for ( i = 0; i < nb; i++ ) b[i] = x2[i];
2498 ma1 = 1; ma2 = 0; mb1 = 0; mb2 = 1;
2499 lx = (((RLONG)(a[na-1]))<<BITSINWORD) + a[na-2];
2500 ly = (((RLONG)(b[nb-1]))<<BITSINWORD) + b[nb-2];
2501 if ( ly > lx ) { lz = lx; lx = ly; ly = lz; d = a; a = b; b = d; }
2504 mc1 = ma1-m*mb1; mc2 = ma2-m*mb2;
2505 ma1 = mb1; ma2 = mb2; mb1 = mc1; mb2 = mc2;
2506 lz = lx - m*ly; lx = ly; ly = lz;
2507 }
while ( ly >= FULLMAX );
2521 x1[1] = (UWORD)(ma1 >> BITSINWORD);
2522 if ( x1[1] ) n1 = -2;
2527 x1[1] = (UWORD)(ma1 >> BITSINWORD);
2528 if ( x1[1] ) n1 = 2;
2531 if ( MulLong(a,na,x1,n1,x2,&n2) )
goto GcdErr;
2535 x1[1] = (UWORD)(ma2 >> BITSINWORD);
2536 if ( x1[1] ) n1 = -2;
2541 x1[1] = (UWORD)(ma2 >> BITSINWORD);
2542 if ( x1[1] ) n1 = 2;
2545 if ( MulLong(b,nb,x1,n1,x3,&n3) )
goto GcdErr;
2546 if ( AddLong(x2,n2,x3,n3,c,&n4) )
goto GcdErr;
2550 x1[1] = (UWORD)(mb1 >> BITSINWORD);
2551 if ( x1[1] ) n1 = -2;
2556 x1[1] = (UWORD)(mb1 >> BITSINWORD);
2557 if ( x1[1] ) n1 = 2;
2560 if ( MulLong(a,na,x1,n1,x2,&n2) )
goto GcdErr;
2564 x1[1] = (UWORD)(mb2 >> BITSINWORD);
2565 if ( x1[1] ) n1 = -2;
2570 x1[1] = (UWORD)(mb2 >> BITSINWORD);
2571 if ( x1[1] ) n1 = 2;
2574 if ( MulLong(b,nb,x1,n1,x3,&n3) )
goto GcdErr;
2575 if ( AddLong(x2,n2,x3,n3,x5,&n5) )
goto GcdErr;
2576 a = c; na = n4; b = x5; nb = n5;
2577 if ( nb == 0 ) { *nc = n4;
goto normalend; }
2579 for ( i = 0; i < na; i++ ) x4[i] = a[i];
2581 if ( na < 0 ) na = -na;
2582 if ( nb < 0 ) nb = -nb;
2593 if ( nb >= GCDMAX && na == nb+1 && b[nb-1] >= HALFMAX && b[nb-1] > a[na-1] ) {
2594 lx = (((RLONG)(a[na-1]))<<BITSINWORD) + a[na-2];
2595 x1[0] = lx/b[nb-1]; n1 = 1;
2596 MulLong(b,nb,x1,n1,x2,&n2);
2598 AddLong(a,na,x2,n2,x4,&n4);
2601 for ( i = 0; i < nb; i++ ) c[i] = b[i];
2604 if ( n4 < 0 ) n4 = -n4;
2605 a = b; na = nb; b = x4; nb = n4;
2613 NumberFree(GLscrat6,
"GcdLong"); NumberFree(GLscrat7,
"GcdLong"); NumberFree(GLscrat8,
"GcdLong");
2614 NumberFree(GLscrat9,
"GcdLong"); NumberFree(GLscrat10,
"GcdLong");
2617 MLOCK(ErrorMessageLock);
2619 MUNLOCK(ErrorMessageLock);
2620 NumberFree(GLscrat6,
"GcdLong"); NumberFree(GLscrat7,
"GcdLong"); NumberFree(GLscrat8,
"GcdLong");
2621 NumberFree(GLscrat9,
"GcdLong"); NumberFree(GLscrat10,
"GcdLong");
2632 WORD GetBinom(UWORD *a, WORD *na, WORD i1, WORD i2)
2636 UWORD *GBscrat3, *GBscrat4;
2637 if ( i1-i2 < i2 ) i2 = i1-i2;
2638 if ( i2 == 0 ) { *a = 1; *na = 1;
return(0); }
2639 if ( i2 > i1 ) { *a = 0; *na = 0;
return(0); }
2641 GBscrat3 = NumberMalloc(
"GetBinom"); GBscrat4 = NumberMalloc(
"GetBinom");
2642 for ( j = 2; j <= i2; j++ ) {
2643 GBscrat3[0] = i1+1-j;
2644 if ( MulLong(a,*na,GBscrat3,(WORD)1,GBscrat4,&k) )
goto CalledFrom;
2646 if ( DivLong(GBscrat4,k,GBscrat3,(WORD)1,a,na,GBscrat3,&l) )
goto CalledFrom;
2648 NumberFree(GBscrat3,
"GetBinom"); NumberFree(GBscrat4,
"GetBinom");
2651 MLOCK(ErrorMessageLock);
2652 MesCall(
"GetBinom");
2653 MUNLOCK(ErrorMessageLock);
2654 NumberFree(GBscrat3,
"GetBinom"); NumberFree(GBscrat4,
"GetBinom");
2666 WORD LcmLong(
PHEAD UWORD *a, WORD na, UWORD *b, WORD nb, UWORD *c, WORD *nc)
2669 UWORD *d = NumberMalloc(
"LcmLong");
2670 UWORD *e = NumberMalloc(
"LcmLong");
2671 UWORD *f = NumberMalloc(
"LcmLong");
2673 GcdLong(BHEAD a, na, b, nb, d, &nd);
2674 DivLong(a,na,d,nd,e,&ne,f,&nf);
2675 if ( MulLong(b,nb,e,ne,c,nc) ) {
2676 MLOCK(ErrorMessageLock);
2678 MUNLOCK(ErrorMessageLock);
2681 NumberFree(f,
"LcmLong");
2682 NumberFree(e,
"LcmLong");
2683 NumberFree(d,
"LcmLong");
2700 int TakeLongRoot(UWORD *a, WORD *n, WORD power)
2703 int numbits, guessbits, i, retval = 0;
2704 UWORD x, *b, *c, *d, *e;
2705 WORD na, nb, nc, nd, ne;
2706 if ( *n < 0 && ( power & 1 ) == 0 )
return(1);
2707 if ( power == 1 )
return(0);
2708 if ( *n < 0 ) { na = -*n; }
2712 if ( a[0] == 1 )
return(0);
2713 if ( power < BITSINWORD && na == 1 && a[0] == (UWORD)(1<<power) ) {
2714 a[0] = 2;
return(0);
2716 if ( 2*power < BITSINWORD && na == 1 && a[0] == (UWORD)(1<<(2*power)) ) {
2717 a[0] = 4;
return(0);
2724 numbits = BITSINWORD*(na-1);
2726 while ( ( x >> 8 ) != 0 ) { numbits += 8; x >>= 8; }
2727 if ( ( x >> 4 ) != 0 ) { numbits += 4; x >>= 4; }
2728 if ( ( x >> 2 ) != 0 ) { numbits += 2; x >>= 2; }
2729 if ( ( x >> 1 ) != 0 ) numbits++;
2730 guessbits = numbits / power;
2731 if ( guessbits <= 0 )
return(1);
2732 nb = guessbits/BITSINWORD;
2742 b = NumberMalloc(
"TakeLongRoot"); c = NumberMalloc(
"TakeLongRoot");
2743 d = NumberMalloc(
"TakeLongRoot"); e = NumberMalloc(
"TakeLongRoot");
2744 for ( i = 0; i < nb; i++ ) { b[i] = 0; }
2745 b[nb] = 1 << (guessbits%BITSINWORD);
2749 for ( i = 0; i < nb; i++ ) c[i] = b[i];
2750 if ( RaisPow(BHEAD c,&nc,power-1) )
goto TLcall;
2751 if ( DivLong(a,na,c,nc,d,&nd,e,&ne) )
goto TLcall;
2753 if ( AddLong(d,nd,b,nb,c,&nc) )
goto TLcall;
2756 if ( ne == 0 )
break;
2766 DivLong(c,nc,(UWORD *)(&power),1,d,&nd,e,&ne);
2794 if ( AddLong(b,nb,d,nd,b,&nb) )
goto TLcall;
2796 for ( i = 0; i < nb; i++ ) a[i] = b[i];
2797 if ( *n < 0 ) *n = -nb;
2799 NumberFree(b,
"TakeLongRoot"); NumberFree(c,
"TakeLongRoot");
2800 NumberFree(d,
"TakeLongRoot"); NumberFree(e,
"TakeLongRoot");
2803 MLOCK(ErrorMessageLock);
2804 MesCall(
"TakeLongRoot");
2805 MUNLOCK(ErrorMessageLock);
2806 NumberFree(b,
"TakeLongRoot"); NumberFree(c,
"TakeLongRoot");
2807 NumberFree(d,
"TakeLongRoot"); NumberFree(e,
"TakeLongRoot");
2820 int MakeRational(WORD a,WORD m, WORD *b, WORD *c)
2822 LONG x1,x2,x3,x4,y1,y2;
2823 if ( a < 0 ) { a = a+m; }
2825 if ( a > m/2 ) a = a-m;
2826 *b = a; *c = 1;
return(0);
2830 y1 = x1/x2; y2 = x1%x2; x3 = 1; x4 = -y1; x1 = x2; x2 = y2;
2831 while ( x2*x2 >= m ) {
2832 y1 = x1/x2; y2 = x1%x2; x1 = x2; x2 = y2; y2 = x3-y1*x4; x3 = x4; x4 = y2;
2836 if ( x2 == 0 ) {
return(1); }
2837 if ( x2 > m/2 ) *b = x2-m;
2839 if ( x4 > m/2 ) { *c = x4-m; *c = -*c; *b = -*b; }
2840 else if ( x4 <= -m/2 ) { x4 += m; *c = x4; }
2841 else if ( x4 < 0 ) { x4 = -x4; *c = x4; *b = -*b; }
2866 #define COPYLONG(x1,nx1,x2,nx2) { int i; for(i=0;i<ABS(nx2);i++)x1[i]=x2[i];nx1=nx2; } 2868 int MakeLongRational(
PHEAD UWORD *a, WORD na, UWORD *m, WORD nm, UWORD *b, WORD *nb)
2870 UWORD *root = NumberMalloc(
"MakeRational");
2871 UWORD *x1 = NumberMalloc(
"MakeRational");
2872 UWORD *x2 = NumberMalloc(
"MakeRational");
2873 UWORD *x3 = NumberMalloc(
"MakeRational");
2874 UWORD *x4 = NumberMalloc(
"MakeRational");
2875 UWORD *y1 = NumberMalloc(
"MakeRational");
2876 UWORD *y2 = NumberMalloc(
"MakeRational");
2877 WORD nroot,nx1,nx2,nx3,nx4,ny1,ny2,retval = 0;
2882 COPYLONG(root,nroot,m,nm)
2883 TakeLongRoot(root,&nroot,2);
2887 if ( na < 0 ) { na = -na; sign = -sign; }
2888 COPYLONG(x1,nx1,m,nm)
2889 COPYLONG(x2,nx2,a,na)
2897 if ( BigLong(x2,nx2,root,nroot) <= 0 ) {
2901 DivLong(x1,nx1,x2,nx2,y1,&ny1,y2,&ny2);
2902 if ( ny2 == 0 ) { retval = 1;
goto cleanup; }
2903 COPYLONG(x1,nx1,x2,nx2)
2904 COPYLONG(x2,nx2,y2,ny2)
2906 COPYLONG(x4,nx4,y1,ny1)
2911 while ( BigLong(x2,nx2,root,nroot) > 0 ) {
2912 DivLong(x1,nx1,x2,nx2,y1,&ny1,y2,&ny2);
2913 if ( ny2 == 0 ) { retval = 1;
goto cleanup; }
2914 COPYLONG(x1,nx1,x2,nx2)
2915 COPYLONG(x2,nx2,y2,ny2)
2916 MulLong(y1,ny1,x4,nx4,y2,&ny2);
2918 AddLong(x3,nx3,y2,ny2,y1,&ny1);
2919 COPYLONG(x3,nx3,x4,nx4)
2920 COPYLONG(x4,nx4,y1,ny1)
2926 if ( nx4 < 0 ) { sign = -sign; nx4 = -nx4; }
2927 COPYLONG(b,*nb,x2,nx2)
2929 if ( sign < 0 ) *nb = -*nb;
2931 NumberFree(y2,
"MakeRational");
2932 NumberFree(y1,
"MakeRational");
2933 NumberFree(x4,
"MakeRational");
2934 NumberFree(x3,
"MakeRational");
2935 NumberFree(x2,
"MakeRational");
2936 NumberFree(x1,
"MakeRational");
2937 NumberFree(root,
"MakeRational");
2957 #ifdef WITHCHINESEREMAINDER 2961 UWORD *inv1 = NumberMalloc(
"ChineseRemainder");
2962 UWORD *inv2 = NumberMalloc(
"ChineseRemainder");
2963 UWORD *fac1 = NumberMalloc(
"ChineseRemainder");
2964 UWORD *fac2 = NumberMalloc(
"ChineseRemainder");
2966 WORD ninv1, ninv2, nfac1, nfac2;
2968 AddLong(a1->a,a1->na,a1->m,a1->nm,a1->a,&(a1->na));
2971 AddLong(a2->a,a2->na,a2->m,a2->nm,a2->a,&(a2->na));
2973 MulLong(a1->m,a1->nm,a2->m,a2->nm,a->m,&(a->nm));
2975 GetLongModInverses(BHEAD a1->m,a1->nm,a2->m,a2->nm,inv1,&ninv1,inv2,&ninv2);
2976 MulLong(inv1,ninv1,a1->m,a1->nm,fac1,&nfac1);
2977 MulLong(inv2,ninv2,a2->m,a2->nm,fac2,&nfac2);
2979 MulLong(fac1,nfac1,a2->a,a2->na,inv1,&ninv1);
2980 MulLong(fac2,nfac2,a1->a,a1->na,inv2,&ninv2);
2981 AddLong(inv1,ninv1,inv2,ninv2,a->a,&(a->na));
2984 MulLong(a->a,a->na,two,1,fac1,&nfac1);
2985 if ( BigLong(fac1,nfac1,a->m,a->nm) > 0 ) {
2987 AddLong(a->a,a->na,a->m,a->nm,a->a,&(a->na));
2990 NumberFree(fac2,
"ChineseRemainder");
2991 NumberFree(fac1,
"ChineseRemainder");
2992 NumberFree(inv2,
"ChineseRemainder");
2993 NumberFree(inv1,
"ChineseRemainder");
3019 if ( term1[1] == 0 && n1 == 1 ) {
3020 if ( term2[1] == 0 && n2 == 1 )
return(0);
3021 if ( n2 < 0 )
return(1);
3024 else if ( term2[1] == 0 && n2 == 1 ) {
3025 if ( n1 < 0 )
return(-1);
3029 if ( n2 < 0 )
return(1);
3032 if ( n2 > 0 )
return(-1);
3033 a = term1; term1 = term2; term2 = a;
3034 n3 = -n1; n1 = -n2; n2 = n3;
3036 if ( term1[1] == 1 && term2[1] == 1 && n1 == 1 && n2 == 1 ) {
3037 if ( (UWORD)*term1 > (UWORD)*term2 )
return(1);
3038 else if ( (UWORD)*term1 < (UWORD)*term2 )
return(-1);
3046 c = NumberMalloc(
"CompCoef");
3047 if ( AddRat(BHEAD (UWORD *)term1,n1,(UWORD *)term2,-n2,c,&n3) ) {
3048 MLOCK(ErrorMessageLock);
3049 MesCall(
"CompCoef");
3050 MUNLOCK(ErrorMessageLock);
3051 NumberFree(c,
"CompCoef");
3054 NumberFree(c,
"CompCoef");
3067 WORD Modulus(WORD *term)
3073 if ( TakeModulus((UWORD *)t,&n1,AC.cmod,AC.ncmod,UNPACK) ) {
3074 MLOCK(ErrorMessageLock);
3076 MUNLOCK(ErrorMessageLock);
3083 else if ( n1 > 0 ) {
3088 else if ( n1 < 0 ) {
3094 *term = WORDDIF(t,term);
3114 WORD TakeModulus(UWORD *a, WORD *na, UWORD *cmodvec, WORD ncmod, WORD par)
3117 UWORD *c, *d, *e, *f, *g, *h;
3119 UWORD *x3,*x1,*x5,*x6,*x7,*x8;
3122 WORD nh, tdenom, tnumer, nmod;
3124 if ( ncmod == 0 )
return(0);
3127 if ( ( par & UNPACK ) != 0 ) UnPack(a,n1,&tdenom,&tnumer);
3128 else { tnumer = n1; }
3133 if ( ( ( par & UNPACK ) == 0 ) && nmod == 1 && ( n1 == 1 || n1 == -1 ) ) {
3136 else if ( nmod == 1 && ( n1 == 1 || n1 == -1 ) ) {
3138 a[1] = a[1] % cmodvec[0];
3140 MesPrint(
"Division by zero in short modulus arithmetic");
3144 if ( ( AC.modinverses != 0 ) && ( ( par & NOINVERSES ) == 0 ) ) {
3145 y1 = AC.modinverses[a[1]];
3151 a[0] = (x*y1) % cmodvec[0];
3156 a[0] = a[0] % cmodvec[0];
3158 if ( a[0] == 0 ) { *na = 0;
return(0); }
3159 if ( ( AC.modmode & POSNEG ) != 0 ) {
3160 if ( a[0] > (UWORD)(cmodvec[0]/2) ) {
3161 a[0] = cmodvec[0] - a[0];
3165 else if ( *na < 0 ) {
3166 *na = 1; a[0] = cmodvec[0] - a[0];
3170 c = NumberMalloc(
"TakeModulus"); d = NumberMalloc(
"TakeModulus"); e = NumberMalloc(
"TakeModulus");
3171 f = NumberMalloc(
"TakeModulus"); g = NumberMalloc(
"TakeModulus"); h = NumberMalloc(
"TakeModulus");
3173 if ( DivLong(a,tnumer,(UWORD *)cmodvec,nmod,
3174 c,&nh,a,&tnumer) )
goto ModErr;
3175 if ( tnumer == 0 ) { *na = 0;
goto normalreturn; }
3176 if ( ( par & UNPACK ) == 0 ) {
3177 if ( ( AC.modmode & POSNEG ) != 0 ) {
3180 else if ( tnumer < 0 ) {
3181 SubPLon((UWORD *)cmodvec,nmod,a,-tnumer,a,&tnumer);
3186 if ( tdenom == 1 && a[n1] == 1 ) {
3187 if ( ( AC.modmode & POSNEG ) != 0 ) {
3190 else if ( tnumer < 0 ) {
3191 SubPLon((UWORD *)cmodvec,nmod,a,-tnumer,a,&tnumer);
3197 while ( --i > 0 ) *a++ = 0;
3200 if ( DivLong(a+n1,tdenom,(UWORD *)cmodvec,nmod,c,&nh,a+n1,&tdenom) )
goto ModErr;
3202 MLOCK(ErrorMessageLock);
3203 MesPrint(
"Division by zero in modulus arithmetic");
3204 if ( AP.DebugFlag ) {
3208 if ( i < 0 ) i = -i;
3209 while ( --i >= 0 ) { TalToLine((UWORD)(*a++)); TokenToLine((UBYTE *)
" "); }
3211 if ( i < 0 ) i = -i;
3212 while ( --i >= 0 ) { TalToLine((UWORD)(*a++)); TokenToLine((UBYTE *)
" "); }
3213 TalToLine((UWORD)(*na));
3217 MUNLOCK(ErrorMessageLock);
3218 NumberFree(c,
"TakeModulus"); NumberFree(d,
"TakeModulus"); NumberFree(e,
"TakeModulus");
3219 NumberFree(f,
"TakeModulus"); NumberFree(g,
"TakeModulus"); NumberFree(h,
"TakeModulus");
3222 if ( ( AC.modinverses != 0 ) && ( ( par & NOINVERSES ) == 0 )
3223 && ( tdenom == 1 || tdenom == -1 ) ) {
3224 *d = AC.modinverses[a[n1]]; y1 = 1; y2 = tdenom;
3225 if ( MulLong(a,tnumer,d,y1,c,&y3) )
goto ModErr;
3226 if ( DivLong(c,y3,(UWORD *)cmodvec,nmod,d,&y5,a,&tdenom) )
goto ModErr;
3227 if ( y2 < 0 ) tdenom = -tdenom;
3230 x2 = (UWORD *)cmodvec; x1 = c; i = nmod;
while ( --i >= 0 ) *x1++ = *x2++;
3231 x1 = c; x2 = a+n1; x3 = d; x4 = e; x5 = f; x6 = g;
3232 y1 = nmod; y2 = tdenom; y4 = 0; y5 = 1; *x5 = 1;
3234 if ( DivLong(x1,y1,x2,y2,h,&nh,x3,&y3) )
goto ModErr;
3235 if ( MulLong(x5,y5,h,nh,x6,&y6) )
goto ModErr;
3236 if ( AddLong(x4,y4,x6,-y6,x6,&y6) )
goto ModErr;
3238 if ( y2 != 1 || *x2 != 1 ) {
3239 MLOCK(ErrorMessageLock);
3240 MesPrint(
"Inverse in modulus arithmetic doesn't exist");
3241 MesPrint(
"Denominator and modulus are not relative prime");
3242 MUNLOCK(ErrorMessageLock);
3247 x7 = x1; x1 = x2; y1 = y2; x2 = x3; y2 = y3; x3 = x7;
3248 x8 = x4; x4 = x5; y4 = y5; x5 = x6; y5 = y6; x6 = x8;
3250 if ( y5 < 0 && AddLong((UWORD *)cmodvec,nmod,x5,y5,x5,&y5) )
goto ModErr;
3251 if ( MulLong(a,tnumer,x5,y5,c,&y3) )
goto ModErr;
3252 if ( DivLong(c,y3,(UWORD *)cmodvec,nmod,d,&y5,a,&tdenom) )
goto ModErr;
3254 if ( !tdenom ) { *na = 0;
goto normalreturn; }
3255 if ( ( ( AC.modmode & POSNEG ) != 0 ) && ( ( par & FROMFUNCTION ) == 0 ) ) {
3258 else if ( tdenom < 0 ) {
3259 SubPLon((UWORD *)cmodvec,nmod,a,-tdenom,a,&tdenom);
3265 while ( --i > 0 ) *a++ = 0;
3267 NumberFree(c,
"TakeModulus"); NumberFree(d,
"TakeModulus"); NumberFree(e,
"TakeModulus");
3268 NumberFree(f,
"TakeModulus"); NumberFree(g,
"TakeModulus"); NumberFree(h,
"TakeModulus");
3271 MLOCK(ErrorMessageLock);
3272 MesCall(
"TakeModulus");
3273 MUNLOCK(ErrorMessageLock);
3274 NumberFree(c,
"TakeModulus"); NumberFree(d,
"TakeModulus"); NumberFree(e,
"TakeModulus");
3275 NumberFree(f,
"TakeModulus"); NumberFree(g,
"TakeModulus"); NumberFree(h,
"TakeModulus");
3286 WORD TakeNormalModulus (UWORD *a, WORD *na, UWORD *c, WORD nc, WORD par)
3295 halfc = NumberMalloc(
"TakeNormalModulus");
3299 for (n=0; n<nhalfc; n++) {
3301 if (n+1<nc) halfc[n] |= c[n+1] << (BITSINWORD-1);
3304 if (halfc[nhalfc-1]==0)
3308 if (BigLong(a,ABS(*na),halfc,nhalfc) > 0) {
3310 TakeModulus(a,na,c,nc,par);
3313 if (BigLong(a,n,halfc,nhalfc) > 0) {
3314 SubPLon(c,nc,a,n,a,&n);
3315 *na = (*na > 0 ? -n : n);
3319 NumberFree(halfc,
"TakeNormalModulus");
3332 if ( AC.modpowers ) {
3333 M_free(AC.modpowers,
"AC.modpowers");
3334 AC.modpowers = NULL;
3337 MLOCK(ErrorMessageLock);
3338 MesPrint(
"&No memory for modulus generator power table");
3339 MUNLOCK(ErrorMessageLock);
3342 if ( n == 0 )
return(0);
3343 size = (LONG)(*AC.cmod);
3344 if ( n == 2 ) size += (((LONG)AC.cmod[1])<<BITSINWORD);
3345 AC.modpowers = (UWORD *)Malloc1(size*n*
sizeof(UWORD),
"table for powers of modulus");
3348 for ( i = 0; i < size; i++ ) AC.modpowers[i] = 0;
3349 for ( i = 0; i < size; i++ ) {
3350 AC.modpowers[j] = (WORD)i;
3354 for ( i = 2; i < size; i++ ) {
3355 if ( AC.modpowers[i] == 0 ) {
3356 MLOCK(ErrorMessageLock);
3357 MesPrint(
"&improper generator for this modulus");
3358 MUNLOCK(ErrorMessageLock);
3359 M_free(AC.modpowers,
"AC.modpowers");
3363 AC.modpowers[1] = 0;
3368 UWORD *MMscrat7 = NumberMalloc(
"MakeModTable"), *MMscratC = NumberMalloc(
"MakeModTable");
3372 for ( i = 0; i < j; i+=2 ) { AC.modpowers[i] = 0; AC.modpowers[i+1] = 0; }
3373 for ( i = 0; i < size; i++ ) {
3374 j = *MMscratC + (((LONG)MMscratC[1])<<BITSINWORD);
3376 AC.modpowers[j] = (WORD)(i & WORDMASK);
3377 AC.modpowers[j+1] = (WORD)(i >> BITSINWORD);
3378 MulLong((UWORD *)MMscratC,nScrat,(UWORD *)AC.powmod,
3379 AC.npowmod,(UWORD *)MMscrat7,&n2);
3380 TakeModulus(MMscrat7,&n2,AC.cmod,AC.ncmod,NOUNPACK);
3381 *MMscratC = *MMscrat7; MMscratC[1] = MMscrat7[1]; nScrat = n2;
3383 NumberFree(MMscrat7,
"MakeModTable"); NumberFree(MMscratC,
"MakeModTable");
3385 for ( i = 4; i < j; i+=2 ) {
3386 if ( AC.modpowers[i] == 0 && AC.modpowers[i+1] == 0 ) {
3387 MLOCK(ErrorMessageLock);
3388 MesPrint(
"&improper generator for this modulus");
3389 MUNLOCK(ErrorMessageLock);
3390 M_free(AC.modpowers,
"AC.modpowers");
3394 AC.modpowers[2] = AC.modpowers[3] = 0;
3414 int Factorial(
PHEAD WORD n, UWORD *a, WORD *na)
3421 if ( n > AT.nfac ) {
3422 if ( AT.factorials == 0 ) {
3423 AT.nfac = 0; AT.mfac = 50; AT.sfact = 400;
3424 AT.pfac = (LONG *)Malloc1((AT.mfac+2)*
sizeof(LONG),
"factorials");
3425 AT.factorials = (UWORD *)Malloc1(AT.sfact*
sizeof(UWORD),
"factorials");
3426 AT.factorials[0] = 1; AT.pfac[0] = 0; AT.pfac[1] = 1;
3429 c = AT.factorials+AT.pfac[AT.nfac];
3430 nc = i = AT.pfac[AT.nfac+1] - AT.pfac[AT.nfac];
3431 while ( --i >= 0 ) *b++ = *c++;
3432 for ( j = AT.nfac+1; j <= n; j++ ) {
3434 if ( nc > AM.MaxTal ) {
3435 MLOCK(ErrorMessageLock);
3436 MesPrint(
"Overflow in factorial. MaxTal = %d",AM.MaxTal);
3437 MesPrint(
"Increase MaxTerm in %s",setupfilename);
3438 MUNLOCK(ErrorMessageLock);
3441 if ( j > AT.mfac ) {
3443 p = (LONG *)Malloc1((AT.mfac*2+2)*
sizeof(LONG),
"factorials");
3445 for ( i = AT.mfac+1; i >= 0; i-- ) p[i] = AT.pfac[i];
3446 M_free(AT.pfac,
"factorial pointers"); AT.pfac = p; AT.mfac *= 2;
3448 if ( AT.pfac[j] + nc >= AT.sfact ) {
3450 f = (UWORD *)Malloc1(AT.sfact*2*
sizeof(UWORD),
"factorials");
3452 c = AT.factorials; b = f;
3453 while ( --ii >= 0 ) *b++ = *c++;
3454 M_free(AT.factorials,
"factorials");
3458 b = a; c = AT.factorials + AT.pfac[j]; i = nc;
3459 while ( --i >= 0 ) *c++ = *b++;
3460 AT.pfac[j+1] = AT.pfac[j] + nc;
3465 else if ( n == 0 ) {
3469 *na = i = AT.pfac[n+1] - AT.pfac[n];
3470 b = AT.factorials + AT.pfac[n];
3471 while ( --i >= 0 ) *a++ = *b++;
3494 int Bernoulli(WORD n, UWORD *a, WORD *na)
3497 UWORD *b, *c, *scrib, *ntop, *ntop1;
3498 WORD i, i1, i2, nhalf, nqua, nscrib, nntop, nntop1, *oldworkpointer;
3499 UWORD twee = 2, twonplus1;
3503 if ( n == 0 ) { a[0] = a[1] = 1; *na = 3; }
3504 else if ( n == 1 ) { a[0] = 1; a[1] = 2; *na = 3; }
3507 if ( ( n & 1 ) != 0 ) { a[0] = a[1] = 0; *na = 0;
return(0); }
3509 if ( nhalf > AT.nBer ) {
3510 oldworkpointer = AT.WorkPointer;
3511 if ( AT.bernoullis == 0 ) {
3512 AT.nBer = 1; AT.mBer = 50; AT.sBer = 400;
3513 AT.pBer = (LONG *)Malloc1((AT.mBer+2)*
sizeof(LONG),
"bernoullis");
3514 AT.bernoullis = (UWORD *)Malloc1(AT.sBer*
sizeof(UWORD),
"bernoullis");
3515 AT.pBer[1] = 0; AT.pBer[2] = 3;
3516 AT.bernoullis[0] = 3; AT.bernoullis[1] = 1; AT.bernoullis[2] = 12;
3518 a[0] = 1; a[1] = 12; *na = 3;
return(0);
3521 while ( nhalf > AT.mBer ) {
3523 p = (LONG *)Malloc1((AT.mBer*2+1)*
sizeof(LONG),
"bernoullis");
3525 for ( i = AT.mBer; i >= 0; i-- ) p[i] = AT.pBer[i];
3526 M_free(AT.pBer,
"factorial pointers"); AT.pBer = p; AT.mBer *= 2;
3528 for ( n = AT.nBer+1; n <= nhalf; n++ ) {
3529 scrib = (UWORD *)(AT.WorkPointer);
3531 if ( ( n & 1 ) == 1 ) {
3532 nscrib = 0; ntop = scrib;
3535 b = AT.bernoullis + AT.pBer[nqua];
3537 i = (WORD)(REDLENG(nscrib));
3538 MulRat(BHEAD b,i,b,i,scrib,&nscrib);
3539 ntop = scrib + 2*nscrib;
3542 for ( j = 1; j <= nqua; j++ ) {
3543 b = AT.bernoullis + AT.pBer[j];
3544 c = AT.bernoullis + AT.pBer[n-j];
3545 i1 = (WORD)(*b); i2 = (WORD)(*c);
3548 MulRat(BHEAD b+1,i1,c+1,i2,ntop,&nntop);
3549 Mully(BHEAD ntop,&nntop,&twee,1);
3551 i = (WORD)nntop;
if ( i < 0 ) i = -i;
3553 AddRat(BHEAD ntop,nntop,scrib,nscrib,ntop1,&nntop1);
3556 ntop1 = ntop; nntop1 = nntop;
3558 nscrib = i1 = (WORD)nntop1;
3559 if ( i1 < 0 ) i1 = - i1;
3561 for ( i = 0; i < i1; i++ ) scrib[i] = ntop1[i];
3565 Divvy(BHEAD scrib,&nscrib,&twonplus1,-1);
3566 i1 = INCLENG(nscrib);
3567 i2 = i1;
if ( i2 < 0 ) i2 = -i2;
3568 i = (WORD)(AT.bernoullis[AT.pBer[n-1]]);
3569 if ( i < 0 ) i = -i;
3570 AT.pBer[n] = AT.pBer[n-1]+i;
3571 if ( AT.pBer[n] + i2 >= AT.sBer ) {
3573 f = (UWORD *)Malloc1(AT.sBer*2*
sizeof(UWORD),
"bernoullis");
3575 c = AT.bernoullis; b = f;
3576 while ( --ii >= 0 ) *b++ = *c++;
3577 M_free(AT.bernoullis,
"bernoullis");
3581 c = AT.bernoullis + AT.pBer[n]; b = scrib;
3583 for ( i = 1; i < i2; i++ ) *c++ = *b++;
3586 AT.WorkPointer = oldworkpointer;
3588 b = AT.bernoullis + AT.pBer[nhalf];
3589 *na = i = (WORD)(*b++);
3590 if ( i < 0 ) i = -i;
3592 while ( --i >= 0 ) *a++ = *b++;
3611 #if ( BITSINWORD == 32 ) 3613 void StartPrimeList(PHEAD0)
3616 AR.PrimeList[AR.numinprimelist++] = 3;
3617 for ( i = 5; i < 46340; i += 2 ) {
3618 for ( j = 0; j < AR.numinprimelist && AR.PrimeList[j]*AR.PrimeList[j] <= i; j++ ) {
3619 if ( i % AR.PrimeList[j] == 0 )
goto nexti;
3621 AR.PrimeList[AR.numinprimelist++] = i;
3624 AR.notfirstprime = 1;
3634 #if ( BITSINWORD == 32 ) 3635 if ( AR.notfirstprime == 0 ) StartPrimeList(BHEAD0);
3637 if ( num > AT.inprimelist ) {
3638 while ( AT.inprimelist < num ) {
3639 if ( num >= AT.sizeprimelist ) {
3640 if ( AT.sizeprimelist == 0 ) newsize = 32;
3641 else newsize = 2*AT.sizeprimelist;
3642 while ( num >= newsize ) newsize = newsize*2;
3643 newpl = (WORD *)Malloc1(newsize*
sizeof(WORD),
"NextPrime");
3644 for ( i = 0; i < AT.sizeprimelist; i++ ) {
3645 newpl[i] = AT.primelist[i];
3647 if ( AT.sizeprimelist > 0 ) {
3648 M_free(AT.primelist,
"NextPrime");
3650 AT.sizeprimelist = newsize;
3651 AT.primelist = newpl;
3653 if ( AT.inprimelist < 0 ) { i = MAXPOSITIVE; }
3654 else { i = AT.primelist[AT.inprimelist]; }
3655 while ( i > MAXPOWER ) {
3657 #if ( BITSINWORD == 32 ) 3658 for ( j = 0; j < AR.numinprimelist && AR.PrimeList[j]*(LONG)(AR.PrimeList[j]) <= x; j++ ) {
3659 if ( x % AR.PrimeList[j] == 0 )
goto nexti;
3662 for ( j = 3; j*((LONG)j) <= x; j += 2 ) {
3663 if ( x % j == 0 )
goto nexti;
3667 AT.primelist[AT.inprimelist] = i;
3671 if ( i < MAXPOWER ) {
3672 MLOCK(ErrorMessageLock);
3673 MesPrint(
"There are not enough short prime numbers for this calculation");
3674 MesPrint(
"Try to use a computer with a %d-bits architecture",
3675 (
int)(BITSINWORD*4));
3676 MUNLOCK(ErrorMessageLock);
3681 return(AT.primelist[num]);
3703 static void wranfnew(PHEAD0)
3707 for ( i = 0; i < AR.wranfnpair1; i++ ) {
3708 j = AR.wranfia[i] - AR.wranfia[i+(AR.wranfnpair2-AR.wranfnpair1)];
3709 if ( j < 0 ) j += (LONG)1 << (2*BITSINWORD-2);
3712 for ( i = AR.wranfnpair1; i < AR.wranfnpair2; i++ ) {
3713 j = AR.wranfia[i] - AR.wranfia[i-AR.wranfnpair1];
3714 if ( j < 0 ) j += (LONG)1 << (2*BITSINWORD-2);
3719 void iniwranf(PHEAD0)
3721 int imax = AR.wranfnpair2-1;
3722 ULONG i, ii, seed = AR.wranfseed;
3724 ULONG offset = 12345;
3727 #if defined(WITHPTHREADS) 3729 #elif defined(WITHMPI) 3736 pow = offset; accu = 1;
3738 if ( ( i & 1 ) != 0 ) accu *= pow;
3739 i >>= 1; pow = pow*pow;
3744 if ( seed < ((LONG)1<<(BITSINWORD-1)) ) {
3745 j = ( (seed+31459L) << (BITSINWORD-2))+offset;
3747 else if ( seed < ((LONG)1<<(BITSINWORD+10-1)) ) {
3748 j = ( (seed+31459L) << (BITSINWORD-10-2))+offset;
3751 j = ( (seed+31459L) << 1)+offset;
3753 if ( ( seed & 1 ) == 1 ) seed++;
3755 AR.wranfia[imax] = j;
3757 for ( i = 0; i <= (ULONG)(imax); i++ ) {
3758 ii = (AR.wranfnpair1*i)%AR.wranfnpair2;
3761 if ( k < 0 ) k += (LONG)1 << (2*BITSINWORD-2);
3764 for ( i = 0; i < WARMUP; i++ ) wranfnew(BHEAD0);
3771 if ( AR.wranfia == 0 ) {
3772 AR.wranfia = (ULONG *)Malloc1(AR.wranfnpair2*
sizeof(ULONG),
"wranf");
3775 if ( AR.wranfcall >= AR.wranfnpair2) {
3779 wval = (UWORD)(AR.wranfia[AR.wranfcall++]>>(BITSINWORD-1));
3787 UWORD iranf(
PHEAD UWORD imax)
3790 ULONG x = (LONG)1 << BITSINWORD, xmax = x - x%imax;
3791 while ( ( i = wranf(BHEAD0) ) >= xmax ) {}
3809 UBYTE *PreRandom(UBYTE *s)
3812 UBYTE *mode,*mins = 0,*maxs = 0, *outval;
3814 double minval, maxval, value = 0;
3817 while ( FG.cTable[*s] <= 1 ) s++;
3818 if ( *s ==
',' ) { *s = 0; s++; }
3820 while ( *s && *s !=
',' ) s++;
3821 if ( *s ==
',' ) { *s = 0; s++; }
3823 while ( *s && *s !=
',' ) s++;
3824 if ( *s || *maxs == 0 || *mins == 0 ) {
3825 MesPrint(
"@Illegal arguments in macro RANDOM_");
3828 if ( StrICmp(mode,(UBYTE *)
"lin") == 0 ) {
3831 else if ( StrICmp(mode,(UBYTE *)
"log") == 0 ) {
3835 MesPrint(
"@Illegal mode argument in macro RANDOM_");
3839 sscanf((
char *)mins,
"%f",&num); minval = num;
3840 sscanf((
char *)maxs,
"%f",&num); maxval = num;
3851 if ( PF.me == MASTER ) {
3858 xx = x/pow(2.0,(
double)(BITSINWORD-1));
3859 if ( linlog == 0 ) {
3860 value = minval + (maxval-minval)*xx;
3862 else if ( linlog == 1 ) {
3863 value = minval * pow(maxval/minval,xx);
3867 outval = (UBYTE *)Malloc1(64,
"PreRandom");
3868 if ( ABS(value) < 0.00001 || ABS(value) > 1000000. ) {
3869 sprintf((
char *)outval,
"%e",value);
3871 else if ( ABS(value) < 0.0001 ) { sprintf((
char *)outval,
"%10f",value); }
3872 else if ( ABS(value) < 0.001 ) { sprintf((
char *)outval,
"%9f",value); }
3873 else if ( ABS(value) < 0.01 ) { sprintf((
char *)outval,
"%8f",value); }
3874 else if ( ABS(value) < 0.1 ) { sprintf((
char *)outval,
"%7f",value); }
3875 else if ( ABS(value) < 1. ) { sprintf((
char *)outval,
"%6f",value); }
3876 else if ( ABS(value) < 10. ) { sprintf((
char *)outval,
"%5f",value); }
3877 else if ( ABS(value) < 100. ) { sprintf((
char *)outval,
"%4f",value); }
3878 else if ( ABS(value) < 1000. ) { sprintf((
char *)outval,
"%3f",value); }
3879 else if ( ABS(value) < 10000. ) { sprintf((
char *)outval,
"%2f",value); }
3880 else { sprintf((
char *)outval,
"%1f",value); }
int GetModInverses(WORD m1, WORD m2, WORD *im1, WORD *im2)
WORD CompCoef(WORD *term1, WORD *term2)
int NormalModulus(UWORD *a, WORD *na)
LONG PF_BroadcastNumber(LONG x)
WORD NextPrime(PHEAD WORD num)
VOID RaisPowCached(PHEAD WORD x, WORD n, UWORD **c, WORD *nc)