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 ) {
386 while ( --i >= 0 ) *c++ = *a++;
389 else if ( nb == -1 ) {
392 while ( --i >= 0 ) *c++ = *a++;
396 if ( *a == 1 && a[1] == 1 ) {
400 while ( --i >= 0 ) *c++ = *b++;
403 else if ( na == -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 #ifdef OLDRAISPOWCACHED 1334 AT.small_power[ID] = NumberMalloc(
"RaisPowCached");
1335 AT.small_power_n[ID] = 1;
1336 AT.small_power[ID][0] = x;
1337 RaisPow(BHEAD AT.small_power[ID],&AT.small_power_n[ID],n);
1339 UWORD *c = NumberMalloc(
"RaisPowCached");
1342 RaisPow(BHEAD c,&k,n);
1346 if ( AT.InNumMem < k ) {
1347 AT.InNumMem = 5*AM.MaxTal;
1348 AT.NumMem = (UWORD *)Malloc1(AT.InNumMem*
sizeof(UWORD),
"RaisPowCached");
1353 for ( i = 0; i < k; i++ ) AT.NumMem[i] = c[i];
1354 AT.small_power[ID] = AT.NumMem;
1355 AT.small_power_n[ID] = k;
1358 NumberFree(c,
"RaisPowCached");
1363 *c = AT.small_power[ID];
1364 *nc = AT.small_power_n[ID];
1373 WORD RaisPowMod (WORD x, WORD n, WORD m) {
1376 if (n&1) { y*=z; y%=m; }
1396 if ( AC.halfmod == 0 ) {
1397 LOCK(AC.halfmodlock);
1398 if ( AC.halfmod == 0 ) {
1399 UWORD two[1],remain[1];
1402 AC.halfmod = (UWORD *)Malloc1((ABS(AC.ncmod))*
sizeof(UWORD),
"halfmod");
1403 DivLong((UWORD *)AC.cmod,(ABS(AC.ncmod)),two,1
1404 ,(UWORD *)AC.halfmod,&(AC.nhalfmod),remain,&dummy);
1406 UNLOCK(AC.halfmodlock);
1409 if ( BigLong(a,n,AC.halfmod,AC.nhalfmod) > 0 ) {
1410 SubPLon((UWORD *)AC.cmod,(ABS(AC.ncmod)),a,n,a,&n);
1411 if ( *na > 0 ) { *na = -n; }
1432 WORD n = AC.cmod[0], i, inv2;
1433 if ( AC.ncmod != 1 )
return(1);
1434 if ( AC.modinverses == 0 ) {
1435 LOCK(AC.halfmodlock);
1436 if ( AC.modinverses == 0 ) {
1437 AC.modinverses = (UWORD *)Malloc1(n*
sizeof(UWORD),
"modinverses");
1438 AC.modinverses[0] = 0;
1439 AC.modinverses[1] = 1;
1440 for ( i = 2; i < n; i++ ) {
1442 (WORD *)(&(AC.modinverses[i])),&inv2) ) {
1447 UNLOCK(AC.halfmodlock);
1470 WORD x = m1, y, c, d = m2;
1471 if ( x < 1 || d <= 1 )
goto somethingwrong;
1476 if ( y == 0 )
break;
1477 a3 = a1-c*a2; a1 = a2; a2 = a3;
1478 b3 = b1-c*b2; b1 = b2; b2 = b3;
1481 if ( x != 1 )
goto somethingwrong;
1482 if ( a2 < 0 ) a2 += m2;
1483 if ( b2 < 0 ) b2 += m1;
1484 if (im1!=NULL) *im1 = a2;
1485 if (im2!=NULL) *im2 = b2;
1488 MLOCK(ErrorMessageLock);
1489 MesPrint(
"Error trying to determine inverses in GetModInverses");
1490 MUNLOCK(ErrorMessageLock);
1498 int GetLongModInverses(PHEAD UWORD *a, WORD na, UWORD *b, WORD nb, UWORD *ia, WORD *nia, UWORD *ib, WORD *nib) {
1500 UWORD *s, *t, *sa, *sb, *ta, *tb, *x, *y, *swap1;
1501 WORD ns, nt, nsa, nsb, nta, ntb, nx, ny, swap2;
1503 s = NumberMalloc(
"GetLongModInverses");
1505 WCOPY(s, a, ABS(ns));
1507 t = NumberMalloc(
"GetLongModInverses");
1509 WCOPY(t, b, ABS(nt));
1511 sa = NumberMalloc(
"GetLongModInverses");
1515 sb = NumberMalloc(
"GetLongModInverses");
1518 ta = NumberMalloc(
"GetLongModInverses");
1521 tb = NumberMalloc(
"GetLongModInverses");
1525 x = NumberMalloc(
"GetLongModInverses");
1526 y = NumberMalloc(
"GetLongModInverses");
1529 DivLong(s,ns,t,nt,x,&nx,y,&ny);
1530 swap1=s; s=y; y=swap1;
1532 MulLong(x,nx,ta,nta,y,&ny);
1533 AddLong(sa,nsa,y,-ny,sa,&nsa);
1534 MulLong(x,nx,tb,ntb,y,&ny);
1535 AddLong(sb,nsb,y,-ny,sb,&nsb);
1537 swap1=s; s=t; t=swap1;
1538 swap2=ns; ns=nt; nt=swap2;
1539 swap1=sa; sa=ta; ta=swap1;
1540 swap2=nsa; nsa=nta; nta=swap2;
1541 swap1=sb; sb=tb; tb=swap1;
1542 swap2=nsb; nsb=ntb; ntb=swap2;
1547 WCOPY(ia,sa,ABS(*nia));
1552 WCOPY(ib,sb,ABS(*nib));
1555 NumberFree(s,
"GetLongModInverses");
1556 NumberFree(t,
"GetLongModInverses");
1557 NumberFree(sa,
"GetLongModInverses");
1558 NumberFree(sb,
"GetLongModInverses");
1559 NumberFree(ta,
"GetLongModInverses");
1560 NumberFree(tb,
"GetLongModInverses");
1561 NumberFree(x,
"GetLongModInverses");
1562 NumberFree(y,
"GetLongModInverses");
1575 WORD Product(UWORD *a, WORD *na, WORD b)
1579 if ( *na < 0 ) { *na = -(*na); sgn = -sgn; }
1580 if ( b < 0 ) { b = -b; sgn = -sgn; }
1583 for ( i = 0; i < *na; i++ ) {
1589 if ( ++(*na) > AM.MaxTal ) {
1590 MLOCK(ErrorMessageLock);
1591 MesPrint(
"Overflow in Product");
1592 MUNLOCK(ErrorMessageLock);
1597 if ( sgn < 0 ) *na = -(*na);
1610 UWORD Quotient(UWORD *a, WORD *na, WORD b)
1614 if ( ( i = *na ) < 0 ) { sgn = -1; i = -i; }
1615 if ( b < 0 ) { b = -b; sgn = -sgn; }
1617 if ( ( *a /= (UWORD)b ) == 0 ) *na = 0;
1618 if ( sgn < 0 ) *na = -*na;
1625 while ( --i >= 0 ) {
1635 if ( sgn < 0 ) j = -j;
1649 WORD Remain10(UWORD *a, WORD *na)
1657 while ( --i >= 0 ) {
1661 if ( i > 0 ) t <<= BITSINWORD;
1663 if ( ( *na > 0 ) && !a[*na-1] ) (*na)--;
1676 WORD Remain4(UWORD *a, WORD *na)
1684 while ( --i >= 0 ) {
1686 *b-- = u = t / 10000;
1688 if ( i > 0 ) t <<= BITSINWORD;
1690 if ( ( *na > 0 ) && !a[*na-1] ) (*na)--;
1702 VOID PrtLong(UWORD *a, WORD na, UBYTE *s)
1715 b = NumberMalloc(
"PrtLong");
1717 i = na;
while ( --i >= 0 ) *bb++ = *a++;
1723 *sa++ = (UBYTE)(
'0' + (q%10));
1725 *sa++ = (UBYTE)(
'0' + (q%10));
1727 *sa++ = (UBYTE)(
'0' + (q%10));
1729 *sa++ = (UBYTE)(
'0' + (q%10));
1731 while ( sa[-1] ==
'0' ) sa--;
1735 while ( sa > sb ) { c = *sa; *sa = *sb; *sb = c; sa--; sb++; }
1740 q = Remain10(a,&na);
1741 *sa++ = (UBYTE)(
'0' + q);
1746 while ( sa > sb ) { c = *sa; *sa = *sb; *sb = c; sa--; sb++; }
1750 NumberFree(b,
"PrtLong");
1764 WORD GetLong(UBYTE *s, UWORD *a, WORD *na)
1777 UWORD digit, x = 0, y = 0;
1780 while ( FG.cTable[*s] == 1 ) {
1782 if ( FG.cTable[*s] != 1 ) { y = 10;
break; }
1783 x = 10*x + *s++ -
'0';
1784 if ( FG.cTable[*s] != 1 ) { y = 100;
break; }
1785 x = 10*x + *s++ -
'0';
1786 if ( FG.cTable[*s] != 1 ) { y = 1000;
break; }
1787 x = 10*x + *s++ -
'0';
1788 if ( *na && Product(a,na,(WORD)10000) )
return(-1);
1789 if ( ( digit = x ) != 0 && AddLong(a,*na,&digit,(WORD)1,a,na) )
1794 if ( *na && Product(a,na,(WORD)y) )
return(-1);
1795 if ( ( digit = x ) != 0 && AddLong(a,*na,&digit,(WORD)1,a,na) )
1815 #define Convert(ia,aa,naa) \ 1816 if ( (LONG)ia < 0 ) { \ 1817 ia = (ULONG)(-(LONG)ia); \ 1819 if ( ( aa[1] = ia >> BITSINWORD ) != 0 ) naa = -2; \ 1822 else if ( ia == 0 ) { aa[0] = 0; naa = 0; } \ 1825 if ( ( aa[1] = ia >> BITSINWORD ) != 0 ) naa = 2; \ 1829 VOID GCD(UWORD *a, WORD na, UWORD *b, WORD nb, UWORD *c, WORD *nc)
1831 int ja = 0, jb = 0, j;
1833 UWORD *x1, *x2, *x3;
1835 ULONG ia,ib,ic,id,u,v,w,q,T;
1840 while ( a[0] == 0 ) { na--; ja++; a++; }
1841 while ( b[0] == 0 ) { nb--; jb++; b++; }
1842 if ( ja > jb ) ja = jb;
1845 do { *c++ = 0; }
while ( --j > 0 );
1851 jb = na; na = nb; nb = jb;
1853 r = a; a = b; b = r;
1855 else if ( na == nb ) {
1859 while ( --j >= 0 ) {
1860 if ( *--r > *--t )
break;
1861 if ( *r < *t )
goto exch;
1892 r = x1 = NumberMalloc(
"GCD"); t = x2 = NumberMalloc(
"GCD"); x3 = NumberMalloc(
"GCD");
1902 DivLong(x1,na,x2,nb,c,nc,x3,&nd);
1903 if ( nd == 0 ) { b = x2;
goto out; }
1904 t = x1; x1 = x2; x2 = x3; x3 = t; na = nb; nb = nd;
1905 if ( na == 2 )
break;
1911 v = x1[0] + ( ((ULONG)x1[1]) << BITSINWORD );
1913 if ( nb == 2 ) w += ((ULONG)x2[1]) << BITSINWORD;
1917 do { u = v%w; v = w; w = u; }
while ( w );
1920 if ( ( c[1] = (UWORD)(v >> BITSINWORD) ) != 0 ) *nc = 2+ja;
1922 NumberFree(x1,
"GCD"); NumberFree(x2,
"GCD"); NumberFree(x3,
"GCD");
1927 ui = x1[0]; uj = x2[0];
1929 ui = (UWORD)GCD2((ULONG)ui,(ULONG)uj);
1931 do { nd = ui%uj; ui = uj; uj = nd; }
while ( nd );
1935 NumberFree(x1,
"GCD"); NumberFree(x2,
"GCD"); NumberFree(x3,
"GCD");
1938 ia = 1; ib = 0; ic = 0;
id = 1;
1939 u = ( ((ULONG)x1[na-1]) << BITSINWORD ) + x1[na-2];
1940 v = ( ((ULONG)x2[nb-1]) << BITSINWORD ) + x2[nb-2];
1942 while ( v+ic != 0 && v+
id != 0 &&
1943 ( q = (u+ia)/(v+ic) ) == (u+ib)/(v+
id) ) {
1944 T = ia-q*ic; ia = ic; ic = T;
1945 T = ib-q*id; ib = id;
id = T;
1946 T = u - q*v; u = v; v = T;
1948 if ( ib == 0 )
goto toobad;
1951 MulLong(x1,na,aa,naa,x3,&nd);
1952 MulLong(x2,nb,bb,nbb,c,nc);
1953 AddLong(x3,nd,c,*nc,c,nc);
1956 MulLong(x1,na,aa,naa,x3,&nd);
1957 t = c; na = j = *nc; r = x1;
1959 MulLong(x2,nb,bb,nbb,c,nc);
1960 AddLong(x3,nd,c,*nc,x2,&nb);
1981 WORD GcdLong(PHEAD UWORD *a, WORD na, UWORD *b, WORD nb, UWORD *c, WORD *nc)
1986 MLOCK(ErrorMessageLock);
1987 MesPrint(
"Cannot take gcd");
1988 MUNLOCK(ErrorMessageLock);
2004 if ( na < 0 ) na = -na;
2005 if ( nb < 0 ) nb = -nb;
2006 if ( na == 1 && nb == 1 ) {
2008 *c = (UWORD)GCD2((ULONG)*a,(ULONG)*b);
2013 do { z = x % y; x = y; }
while ( ( y = z ) != 0 );
2018 else if ( na <= 2 && nb <= 2 ) {
2020 if ( na == 2 ) { lx = (((RLONG)(a[1]))<<BITSINWORD) + *a; }
2022 if ( nb == 2 ) { ly = (((RLONG)(b[1]))<<BITSINWORD) + *b; }
2028 do { lz = lx % ly; lx = ly; }
while ( ( ly = lz ) != 0 );
2031 if ( lx < ly ) { lz = lx; lx = ly; ly = lz; }
2033 lz = lx % ly; lx = ly;
2034 }
while ( ( ly = lz ) != 0 && ( lx & AWORDMASK ) != 0 );
2036 do { *c = ((UWORD)lx)%((UWORD)ly); lx = ly; }
while ( ( ly = *c ) != 0 );
2044 if ( ( *c = (UWORD)(lx >> BITSINWORD) ) != 0 ) *nc = 2;
2050 GCD(a,na,b,nb,c,nc);
2053 UWORD *x3,*x1,*x2, *GLscrat7, *GLscrat8;
2056 x1 = c; x3 = a; n1 = i = na;
2058 GLscrat7 = NumberMalloc(
"GcdLong"); GLscrat8 = NumberMalloc(
"GcdLong");
2059 x2 = GLscrat8; x3 = b; n2 = i = nb;
2062 while ( x1[0] == 0 ) { i += BITSINWORD; x1++; n1--; }
2063 while ( ( x1[0] & 1 ) == 0 ) { i++; SCHUIF(x1,n1) }
2064 x2 = GLscrat8; j = 0;
2065 while ( x2[0] == 0 ) { j += BITSINWORD; x2++; n2--; }
2066 while ( ( x2[0] & 1 ) == 0 ) { j++; SCHUIF(x2,n2) }
2071 SubPLon(x1,n1,x2,n2,x1,&n3);
2075 n1 = i = n2; NCOPY(x1,x2,i);
2078 while ( ( x1[0] & 1 ) == 0 ) SCHUIF(x1,n1)
2080 if ( DivLong(x2,n2,x1,n1,GLscrat7,&n3,x2,&n4) )
goto GcdErr;
2083 i = n1; x2 = c; NCOPY(x2,x1,i);
2087 *c = (UWORD)GCD2((ULONG)x1[0],(ULONG)x2[0]);
2093 do { z = x % y; x = y; }
while ( ( y = z ) != 0 );
2101 else if ( n1 < n2 ) {
2103 SubPLon(x2,n2,x1,n1,x2,&n3);
2106 i = n1; x2 = c; NCOPY(x2,x1,i);
2109 while ( ( x2[0] & 1 ) == 0 ) SCHUIF(x2,n2)
2111 if ( DivLong(x1,n1,x2,n2,GLscrat7,&n3,x1,&n4) )
goto GcdErr;
2115 n1 = i = n2; NCOPY(x1,x2,i);
2119 *c = (UWORD)GCD2((ULONG)x2[0],(ULONG)x1[0]);
2125 do { z = x % y; x = y; }
while ( ( y = z ) != 0 );
2134 for ( i = n1-1; i >= 0; i-- ) {
2135 if ( x1[i] > x2[i] )
goto firstbig;
2136 else if ( x1[i] < x2[i] )
goto lastbig;
2138 i = n1; x2 = c; NCOPY(x2,x1,i);
2146 while ( j >= BITSINWORD ) {
2147 for ( i = n1; i > 0; i-- ) x1[i] = x1[i-1];
2153 for ( i = 0; i < n1; i++ ) {
2154 a1 = x1[i]; a1 <<= j;
2164 NumberFree(GLscrat7,
"GcdLong"); NumberFree(GLscrat8,
"GcdLong");
2166 UWORD *x1,*x2,*x3,*x4,*c1,*c2;
2168 x1 = c; x3 = a; n1 = i = na;
2170 x1 = c; c1 = x2 = NumberMalloc(
"GcdLong"); x3 = NumberMalloc(
"GcdLong"); x4 = NumberMalloc(
"GcdLong");
2171 c2 = b; n2 = i = nb;
2174 if ( DivLong(x1,n1,x2,n2,x4,&n4,x3,&n3) )
goto GcdErr;
2175 if ( !n3 ) { x1 = x2; n1 = n2;
break; }
2176 if ( DivLong(x2,n2,x3,n3,x4,&n4,x1,&n1) )
goto GcdErr;
2177 if ( !n1 ) { x1 = x3; n1 = n3;
break; }
2178 if ( DivLong(x3,n3,x1,n1,x4,&n4,x2,&n2) )
goto GcdErr;
2181 NumberFree(x2,
"GcdLong"); NumberFree(x3,
"GcdLong"); NumberFree(x4,
"GcdLong");
2187 NumberFree(x2,
"GcdLong"); NumberFree(x3,
"GcdLong"); NumberFree(x4,
"GcdLong");
2193 MLOCK(ErrorMessageLock);
2195 MUNLOCK(ErrorMessageLock);
2266 WORD GcdLong(PHEAD UWORD *a, WORD na, UWORD *b, WORD nb, UWORD *c, WORD *nc)
2270 UWORD *x1,*x2,*x3,*x4,*x5,*d;
2271 UWORD *GLscrat6, *GLscrat7, *GLscrat8, *GLscrat9, *GLscrat10;
2272 WORD n1,n2,n3,n4,n5,i;
2274 LONG ma1, ma2, mb1, mb2, mc1, mc2, m;
2277 MLOCK(ErrorMessageLock);
2278 MesPrint(
"Cannot take gcd");
2279 MUNLOCK(ErrorMessageLock);
2295 if ( na < 0 ) na = -na;
2296 if ( nb < 0 ) nb = -nb;
2301 if ( na > 3 && nb > 3 ) {
2303 mp_limb_t *upa, *upb, *upc, xx;
2304 UWORD *uw, *u1, *u2;
2305 unsigned int tcounta, tcountb, tcounta1, tcountb1;
2306 mp_size_t ana, anb, anc;
2308 u1 = uw = NumberMalloc(
"GcdLong");
2309 upa = (mp_limb_t *)u1;
2310 ana = na; tcounta1 = 0;
2311 while ( a[0] == 0 ) { a++; ana--; tcounta1++; }
2312 for ( ii = 0; ii < ana; ii++ ) { *uw++ = *a++; }
2313 if ( ( ana & 1 ) != 0 ) { *uw = 0; ana++; }
2316 u2 = uw = NumberMalloc(
"GcdLong");
2317 upb = (mp_limb_t *)u2;
2318 anb = nb; tcountb1 = 0;
2319 while ( b[0] == 0 ) { b++; anb--; tcountb1++; }
2320 for ( ii = 0; ii < anb; ii++ ) { *uw++ = *b++; }
2321 if ( ( anb & 1 ) != 0 ) { *uw = 0; anb++; }
2324 xx = upa[0]; tcounta = 0;
2325 while ( ( xx & 15 ) == 0 ) { tcounta += 4; xx >>= 4; }
2326 while ( ( xx & 1 ) == 0 ) { tcounta += 1; xx >>= 1; }
2327 xx = upb[0]; tcountb = 0;
2328 while ( ( xx & 15 ) == 0 ) { tcountb += 4; xx >>= 4; }
2329 while ( ( xx & 1 ) == 0 ) { tcountb += 1; xx >>= 1; }
2332 mpn_rshift(upa,upa,ana,tcounta);
2333 if ( upa[ana-1] == 0 ) ana--;
2336 mpn_rshift(upb,upb,anb,tcountb);
2337 if ( upb[anb-1] == 0 ) anb--;
2340 upc = (mp_limb_t *)(NumberMalloc(
"GcdLong"));
2341 if ( ( ana > anb ) || ( ( ana == anb ) && ( upa[ana-1] >= upb[ana-1] ) ) ) {
2342 anc = mpn_gcd(upc,upa,ana,upb,anb);
2345 anc = mpn_gcd(upc,upb,anb,upa,ana);
2348 tcounta = tcounta1*BITSINWORD + tcounta;
2349 tcountb = tcountb1*BITSINWORD + tcountb;
2350 if ( tcountb > tcounta ) tcountb = tcounta;
2351 tcounta = tcountb/BITSINWORD;
2352 tcountb = tcountb%BITSINWORD;
2355 xx = mpn_lshift(upc,upc,anc,tcountb);
2356 if ( xx ) { upc[anc] = xx; anc++; }
2359 uw = (UWORD *)upc; anc *= 2;
2360 while ( uw[anc-1] == 0 ) anc--;
2361 for ( ii = 0; ii < (int)tcounta; ii++ ) *c++ = 0;
2362 for ( ii = 0; ii < anc; ii++ ) *c++ = *uw++;
2363 *nc = anc + tcounta;
2364 NumberFree(u1,
"GcdLong"); NumberFree(u2,
"GcdLong"); NumberFree((UWORD *)(upc),
"GcdLong");
2374 if ( na == 1 && nb == 1 ) {
2377 do { z = x % y; x = y; }
while ( ( y = z ) != 0 );
2382 else if ( na <= 2 && nb <= 2 ) {
2383 if ( na == 2 ) { lx = (((RLONG)(a[1]))<<BITSINWORD) + *a; }
2385 if ( nb == 2 ) { ly = (((RLONG)(b[1]))<<BITSINWORD) + *b; }
2387 if ( lx < ly ) { lz = lx; lx = ly; ly = lz; }
2388 #if ( BITSINWORD == 16 ) 2390 lz = lx % ly; lx = ly;
2391 }
while ( ( ly = lz ) != 0 );
2394 lz = lx % ly; lx = ly;
2395 }
while ( ( ly = lz ) != 0 && ( lx & AWORDMASK ) != 0 );
2397 x = (UWORD)lx; y = (UWORD)ly;
2398 do { *c = x % y; x = y; }
while ( ( y = *c ) != 0 );
2406 if ( ( *c = (UWORD)(lx >> BITSINWORD) ) != 0 ) *nc = 2;
2414 GLscrat6 = NumberMalloc(
"GcdLong"); GLscrat7 = NumberMalloc(
"GcdLong");
2415 GLscrat8 = NumberMalloc(
"GcdLong");
2416 GLscrat9 = NumberMalloc(
"GcdLong"); GLscrat10 = NumberMalloc(
"GcdLong");
2421 if ( na == 1 && nb == 1 ) {
2424 do { z = x % y; x = y; }
while ( ( y = z ) != 0 );
2428 else if ( na <= 2 && nb <= 2 ) {
2429 if ( na == 2 ) { lx = (((RLONG)(a[1]))<<BITSINWORD) + *a; }
2431 if ( nb == 2 ) { ly = (((RLONG)(b[1]))<<BITSINWORD) + *b; }
2433 if ( lx < ly ) { lz = lx; lx = ly; ly = lz; }
2434 #if ( BITSINWORD == 16 ) 2436 lz = lx % ly; lx = ly;
2437 }
while ( ( ly = lz ) != 0 );
2440 lz = lx % ly; lx = ly;
2441 }
while ( ( ly = lz ) != 0 && ( lx & AWORDMASK ) != 0 );
2443 x = (UWORD)lx; y = (UWORD)ly;
2444 do { *c = x % y; x = y; }
while ( ( y = *c ) != 0 );
2452 if ( ( *c = (UWORD)(lx >> BITSINWORD) ) != 0 ) *nc = 2;
2460 else if ( na < GCDMAX || nb < GCDMAX || na != nb ) {
2462 x2 = GLscrat8; x3 = a; n2 = i = na;
2464 x1 = c; x3 = b; n1 = i = nb;
2468 x1 = c; x3 = a; n1 = i = na;
2470 x2 = GLscrat8; x3 = b; n2 = i = nb;
2473 x1 = c; x2 = GLscrat8; x3 = GLscrat7; x4 = GLscrat6;
2475 if ( DivLong(x1,n1,x2,n2,x4,&n4,x3,&n3) )
goto GcdErr;
2476 if ( !n3 ) { x1 = x2; n1 = n2;
break; }
2477 if ( n2 <= 2 ) { a = x2; b = x3; na = n2; nb = n3;
goto restart; }
2478 if ( n3 >= GCDMAX && n2 == n3 ) {
2479 a = GLscrat9; b = GLscrat10; na = n2; nb = n3;
2480 for ( i = 0; i < na; i++ ) a[i] = x2[i];
2481 for ( i = 0; i < nb; i++ ) b[i] = x3[i];
2484 if ( DivLong(x2,n2,x3,n3,x4,&n4,x1,&n1) )
goto GcdErr;
2485 if ( !n1 ) { x1 = x3; n1 = n3;
break; }
2486 if ( n3 <= 2 ) { a = x3; b = x1; na = n3; nb = n1;
goto restart; }
2487 if ( n1 >= GCDMAX && n1 == n3 ) {
2488 a = GLscrat9; b = GLscrat10; na = n3; nb = n1;
2489 for ( i = 0; i < na; i++ ) a[i] = x3[i];
2490 for ( i = 0; i < nb; i++ ) b[i] = x1[i];
2493 if ( DivLong(x3,n3,x1,n1,x4,&n4,x2,&n2) )
goto GcdErr;
2494 if ( !n2 ) { *nc = n1;
goto normalend; }
2495 if ( n1 <= 2 ) { a = x1; b = x2; na = n1; nb = n2;
goto restart; }
2496 if ( n2 >= GCDMAX && n2 == n1 ) {
2497 a = GLscrat9; b = GLscrat10; na = n1; nb = n2;
2498 for ( i = 0; i < na; i++ ) a[i] = x1[i];
2499 for ( i = 0; i < nb; i++ ) b[i] = x2[i];
2523 ma1 = 1; ma2 = 0; mb1 = 0; mb2 = 1;
2524 lx = (((RLONG)(a[na-1]))<<BITSINWORD) + a[na-2];
2525 ly = (((RLONG)(b[nb-1]))<<BITSINWORD) + b[nb-2];
2526 if ( ly > lx ) { lz = lx; lx = ly; ly = lz; d = a; a = b; b = d; }
2529 mc1 = ma1-m*mb1; mc2 = ma2-m*mb2;
2530 ma1 = mb1; ma2 = mb2; mb1 = mc1; mb2 = mc2;
2531 lz = lx - m*ly; lx = ly; ly = lz;
2532 }
while ( ly >= FULLMAX );
2546 x1[1] = (UWORD)(ma1 >> BITSINWORD);
2547 if ( x1[1] ) n1 = -2;
2552 x1[1] = (UWORD)(ma1 >> BITSINWORD);
2553 if ( x1[1] ) n1 = 2;
2556 if ( MulLong(a,na,x1,n1,x2,&n2) )
goto GcdErr;
2560 x1[1] = (UWORD)(ma2 >> BITSINWORD);
2561 if ( x1[1] ) n1 = -2;
2566 x1[1] = (UWORD)(ma2 >> BITSINWORD);
2567 if ( x1[1] ) n1 = 2;
2570 if ( MulLong(b,nb,x1,n1,x3,&n3) )
goto GcdErr;
2571 if ( AddLong(x2,n2,x3,n3,c,&n4) )
goto GcdErr;
2575 x1[1] = (UWORD)(mb1 >> BITSINWORD);
2576 if ( x1[1] ) n1 = -2;
2581 x1[1] = (UWORD)(mb1 >> BITSINWORD);
2582 if ( x1[1] ) n1 = 2;
2585 if ( MulLong(a,na,x1,n1,x2,&n2) )
goto GcdErr;
2589 x1[1] = (UWORD)(mb2 >> BITSINWORD);
2590 if ( x1[1] ) n1 = -2;
2595 x1[1] = (UWORD)(mb2 >> BITSINWORD);
2596 if ( x1[1] ) n1 = 2;
2599 if ( MulLong(b,nb,x1,n1,x3,&n3) )
goto GcdErr;
2600 if ( AddLong(x2,n2,x3,n3,x5,&n5) )
goto GcdErr;
2601 a = c; na = n4; b = x5; nb = n5;
2602 if ( nb == 0 ) { *nc = n4;
goto normalend; }
2604 for ( i = 0; i < na; i++ ) x4[i] = a[i];
2606 if ( na < 0 ) na = -na;
2607 if ( nb < 0 ) nb = -nb;
2618 if ( nb >= GCDMAX && na == nb+1 && b[nb-1] >= HALFMAX && b[nb-1] > a[na-1] ) {
2619 lx = (((RLONG)(a[na-1]))<<BITSINWORD) + a[na-2];
2620 x1[0] = lx/b[nb-1]; n1 = 1;
2621 MulLong(b,nb,x1,n1,x2,&n2);
2623 AddLong(a,na,x2,n2,x4,&n4);
2626 for ( i = 0; i < nb; i++ ) c[i] = b[i];
2629 if ( n4 < 0 ) n4 = -n4;
2630 a = b; na = nb; b = x4; nb = n4;
2638 NumberFree(GLscrat6,
"GcdLong"); NumberFree(GLscrat7,
"GcdLong"); NumberFree(GLscrat8,
"GcdLong");
2639 NumberFree(GLscrat9,
"GcdLong"); NumberFree(GLscrat10,
"GcdLong");
2642 MLOCK(ErrorMessageLock);
2644 MUNLOCK(ErrorMessageLock);
2645 NumberFree(GLscrat6,
"GcdLong"); NumberFree(GLscrat7,
"GcdLong"); NumberFree(GLscrat8,
"GcdLong");
2646 NumberFree(GLscrat9,
"GcdLong"); NumberFree(GLscrat10,
"GcdLong");
2657 WORD GetBinom(UWORD *a, WORD *na, WORD i1, WORD i2)
2661 UWORD *GBscrat3, *GBscrat4;
2662 if ( i1-i2 < i2 ) i2 = i1-i2;
2663 if ( i2 == 0 ) { *a = 1; *na = 1;
return(0); }
2664 if ( i2 > i1 ) { *a = 0; *na = 0;
return(0); }
2666 GBscrat3 = NumberMalloc(
"GetBinom"); GBscrat4 = NumberMalloc(
"GetBinom");
2667 for ( j = 2; j <= i2; j++ ) {
2668 GBscrat3[0] = i1+1-j;
2669 if ( MulLong(a,*na,GBscrat3,(WORD)1,GBscrat4,&k) )
goto CalledFrom;
2671 if ( DivLong(GBscrat4,k,GBscrat3,(WORD)1,a,na,GBscrat3,&l) )
goto CalledFrom;
2673 NumberFree(GBscrat3,
"GetBinom"); NumberFree(GBscrat4,
"GetBinom");
2676 MLOCK(ErrorMessageLock);
2677 MesCall(
"GetBinom");
2678 MUNLOCK(ErrorMessageLock);
2679 NumberFree(GBscrat3,
"GetBinom"); NumberFree(GBscrat4,
"GetBinom");
2691 WORD LcmLong(PHEAD UWORD *a, WORD na, UWORD *b, WORD nb, UWORD *c, WORD *nc)
2694 UWORD *d = NumberMalloc(
"LcmLong");
2695 UWORD *e = NumberMalloc(
"LcmLong");
2696 UWORD *f = NumberMalloc(
"LcmLong");
2698 GcdLong(BHEAD a, na, b, nb, d, &nd);
2699 DivLong(a,na,d,nd,e,&ne,f,&nf);
2700 if ( MulLong(b,nb,e,ne,c,nc) ) {
2701 MLOCK(ErrorMessageLock);
2703 MUNLOCK(ErrorMessageLock);
2706 NumberFree(f,
"LcmLong");
2707 NumberFree(e,
"LcmLong");
2708 NumberFree(d,
"LcmLong");
2725 int TakeLongRoot(UWORD *a, WORD *n, WORD power)
2728 int numbits, guessbits, i, retval = 0;
2729 UWORD x, *b, *c, *d, *e;
2730 WORD na, nb, nc, nd, ne;
2731 if ( *n < 0 && ( power & 1 ) == 0 )
return(1);
2732 if ( power == 1 )
return(0);
2733 if ( *n < 0 ) { na = -*n; }
2737 if ( a[0] == 1 )
return(0);
2738 if ( power < BITSINWORD && na == 1 && a[0] == (UWORD)(1<<power) ) {
2739 a[0] = 2;
return(0);
2741 if ( 2*power < BITSINWORD && na == 1 && a[0] == (UWORD)(1<<(2*power)) ) {
2742 a[0] = 4;
return(0);
2749 numbits = BITSINWORD*(na-1);
2751 while ( ( x >> 8 ) != 0 ) { numbits += 8; x >>= 8; }
2752 if ( ( x >> 4 ) != 0 ) { numbits += 4; x >>= 4; }
2753 if ( ( x >> 2 ) != 0 ) { numbits += 2; x >>= 2; }
2754 if ( ( x >> 1 ) != 0 ) numbits++;
2755 guessbits = numbits / power;
2756 if ( guessbits <= 0 )
return(1);
2757 nb = guessbits/BITSINWORD;
2767 b = NumberMalloc(
"TakeLongRoot"); c = NumberMalloc(
"TakeLongRoot");
2768 d = NumberMalloc(
"TakeLongRoot"); e = NumberMalloc(
"TakeLongRoot");
2769 for ( i = 0; i < nb; i++ ) { b[i] = 0; }
2770 b[nb] = 1 << (guessbits%BITSINWORD);
2774 for ( i = 0; i < nb; i++ ) c[i] = b[i];
2775 if ( RaisPow(BHEAD c,&nc,power-1) )
goto TLcall;
2776 if ( DivLong(a,na,c,nc,d,&nd,e,&ne) )
goto TLcall;
2778 if ( AddLong(d,nd,b,nb,c,&nc) )
goto TLcall;
2781 if ( ne == 0 )
break;
2791 DivLong(c,nc,(UWORD *)(&power),1,d,&nd,e,&ne);
2819 if ( AddLong(b,nb,d,nd,b,&nb) )
goto TLcall;
2821 for ( i = 0; i < nb; i++ ) a[i] = b[i];
2822 if ( *n < 0 ) *n = -nb;
2824 NumberFree(b,
"TakeLongRoot"); NumberFree(c,
"TakeLongRoot");
2825 NumberFree(d,
"TakeLongRoot"); NumberFree(e,
"TakeLongRoot");
2828 MLOCK(ErrorMessageLock);
2829 MesCall(
"TakeLongRoot");
2830 MUNLOCK(ErrorMessageLock);
2831 NumberFree(b,
"TakeLongRoot"); NumberFree(c,
"TakeLongRoot");
2832 NumberFree(d,
"TakeLongRoot"); NumberFree(e,
"TakeLongRoot");
2845 int MakeRational(WORD a,WORD m, WORD *b, WORD *c)
2847 LONG x1,x2,x3,x4,y1,y2;
2848 if ( a < 0 ) { a = a+m; }
2850 if ( a > m/2 ) a = a-m;
2851 *b = a; *c = 1;
return(0);
2855 y1 = x1/x2; y2 = x1%x2; x3 = 1; x4 = -y1; x1 = x2; x2 = y2;
2856 while ( x2*x2 >= m ) {
2857 y1 = x1/x2; y2 = x1%x2; x1 = x2; x2 = y2; y2 = x3-y1*x4; x3 = x4; x4 = y2;
2861 if ( x2 == 0 ) {
return(1); }
2862 if ( x2 > m/2 ) *b = x2-m;
2864 if ( x4 > m/2 ) { *c = x4-m; *c = -*c; *b = -*b; }
2865 else if ( x4 <= -m/2 ) { x4 += m; *c = x4; }
2866 else if ( x4 < 0 ) { x4 = -x4; *c = x4; *b = -*b; }
2891 #define COPYLONG(x1,nx1,x2,nx2) { int i; for(i=0;i<ABS(nx2);i++)x1[i]=x2[i];nx1=nx2; } 2893 int MakeLongRational(PHEAD UWORD *a, WORD na, UWORD *m, WORD nm, UWORD *b, WORD *nb)
2895 UWORD *root = NumberMalloc(
"MakeRational");
2896 UWORD *x1 = NumberMalloc(
"MakeRational");
2897 UWORD *x2 = NumberMalloc(
"MakeRational");
2898 UWORD *x3 = NumberMalloc(
"MakeRational");
2899 UWORD *x4 = NumberMalloc(
"MakeRational");
2900 UWORD *y1 = NumberMalloc(
"MakeRational");
2901 UWORD *y2 = NumberMalloc(
"MakeRational");
2902 WORD nroot,nx1,nx2,nx3,nx4,ny1,ny2,retval = 0;
2907 COPYLONG(root,nroot,m,nm)
2908 TakeLongRoot(root,&nroot,2);
2912 if ( na < 0 ) { na = -na; sign = -sign; }
2913 COPYLONG(x1,nx1,m,nm)
2914 COPYLONG(x2,nx2,a,na)
2922 if ( BigLong(x2,nx2,root,nroot) <= 0 ) {
2926 DivLong(x1,nx1,x2,nx2,y1,&ny1,y2,&ny2);
2927 if ( ny2 == 0 ) { retval = 1;
goto cleanup; }
2928 COPYLONG(x1,nx1,x2,nx2)
2929 COPYLONG(x2,nx2,y2,ny2)
2931 COPYLONG(x4,nx4,y1,ny1)
2936 while ( BigLong(x2,nx2,root,nroot) > 0 ) {
2937 DivLong(x1,nx1,x2,nx2,y1,&ny1,y2,&ny2);
2938 if ( ny2 == 0 ) { retval = 1;
goto cleanup; }
2939 COPYLONG(x1,nx1,x2,nx2)
2940 COPYLONG(x2,nx2,y2,ny2)
2941 MulLong(y1,ny1,x4,nx4,y2,&ny2);
2943 AddLong(x3,nx3,y2,ny2,y1,&ny1);
2944 COPYLONG(x3,nx3,x4,nx4)
2945 COPYLONG(x4,nx4,y1,ny1)
2951 if ( nx4 < 0 ) { sign = -sign; nx4 = -nx4; }
2952 COPYLONG(b,*nb,x2,nx2)
2954 if ( sign < 0 ) *nb = -*nb;
2956 NumberFree(y2,
"MakeRational");
2957 NumberFree(y1,
"MakeRational");
2958 NumberFree(x4,
"MakeRational");
2959 NumberFree(x3,
"MakeRational");
2960 NumberFree(x2,
"MakeRational");
2961 NumberFree(x1,
"MakeRational");
2962 NumberFree(root,
"MakeRational");
2982 #ifdef WITHCHINESEREMAINDER 2986 UWORD *inv1 = NumberMalloc(
"ChineseRemainder");
2987 UWORD *inv2 = NumberMalloc(
"ChineseRemainder");
2988 UWORD *fac1 = NumberMalloc(
"ChineseRemainder");
2989 UWORD *fac2 = NumberMalloc(
"ChineseRemainder");
2991 WORD ninv1, ninv2, nfac1, nfac2;
2993 AddLong(a1->a,a1->na,a1->m,a1->nm,a1->a,&(a1->na));
2996 AddLong(a2->a,a2->na,a2->m,a2->nm,a2->a,&(a2->na));
2998 MulLong(a1->m,a1->nm,a2->m,a2->nm,a->m,&(a->nm));
3000 GetLongModInverses(BHEAD a1->m,a1->nm,a2->m,a2->nm,inv1,&ninv1,inv2,&ninv2);
3001 MulLong(inv1,ninv1,a1->m,a1->nm,fac1,&nfac1);
3002 MulLong(inv2,ninv2,a2->m,a2->nm,fac2,&nfac2);
3004 MulLong(fac1,nfac1,a2->a,a2->na,inv1,&ninv1);
3005 MulLong(fac2,nfac2,a1->a,a1->na,inv2,&ninv2);
3006 AddLong(inv1,ninv1,inv2,ninv2,a->a,&(a->na));
3009 MulLong(a->a,a->na,two,1,fac1,&nfac1);
3010 if ( BigLong(fac1,nfac1,a->m,a->nm) > 0 ) {
3012 AddLong(a->a,a->na,a->m,a->nm,a->a,&(a->na));
3015 NumberFree(fac2,
"ChineseRemainder");
3016 NumberFree(fac1,
"ChineseRemainder");
3017 NumberFree(inv2,
"ChineseRemainder");
3018 NumberFree(inv1,
"ChineseRemainder");
3044 if ( term1[1] == 0 && n1 == 1 ) {
3045 if ( term2[1] == 0 && n2 == 1 )
return(0);
3046 if ( n2 < 0 )
return(1);
3049 else if ( term2[1] == 0 && n2 == 1 ) {
3050 if ( n1 < 0 )
return(-1);
3054 if ( n2 < 0 )
return(1);
3057 if ( n2 > 0 )
return(-1);
3058 a = term1; term1 = term2; term2 = a;
3059 n3 = -n1; n1 = -n2; n2 = n3;
3061 if ( term1[1] == 1 && term2[1] == 1 && n1 == 1 && n2 == 1 ) {
3062 if ( (UWORD)*term1 > (UWORD)*term2 )
return(1);
3063 else if ( (UWORD)*term1 < (UWORD)*term2 )
return(-1);
3071 c = NumberMalloc(
"CompCoef");
3072 if ( AddRat(BHEAD (UWORD *)term1,n1,(UWORD *)term2,-n2,c,&n3) ) {
3073 MLOCK(ErrorMessageLock);
3074 MesCall(
"CompCoef");
3075 MUNLOCK(ErrorMessageLock);
3076 NumberFree(c,
"CompCoef");
3079 NumberFree(c,
"CompCoef");
3092 WORD Modulus(WORD *term)
3098 if ( TakeModulus((UWORD *)t,&n1,AC.cmod,AC.ncmod,UNPACK) ) {
3099 MLOCK(ErrorMessageLock);
3101 MUNLOCK(ErrorMessageLock);
3108 else if ( n1 > 0 ) {
3113 else if ( n1 < 0 ) {
3119 *term = WORDDIF(t,term);
3139 WORD TakeModulus(UWORD *a, WORD *na, UWORD *cmodvec, WORD ncmod, WORD par)
3142 UWORD *c, *d, *e, *f, *g, *h;
3144 UWORD *x3,*x1,*x5,*x6,*x7,*x8;
3147 WORD nh, tdenom, tnumer, nmod;
3149 if ( ncmod == 0 )
return(0);
3152 if ( ( par & UNPACK ) != 0 ) UnPack(a,n1,&tdenom,&tnumer);
3153 else { tnumer = n1; }
3158 if ( ( ( par & UNPACK ) == 0 ) && nmod == 1 && ( n1 == 1 || n1 == -1 ) ) {
3161 else if ( nmod == 1 && ( n1 == 1 || n1 == -1 ) ) {
3163 a[1] = a[1] % cmodvec[0];
3165 MesPrint(
"Division by zero in short modulus arithmetic");
3169 if ( ( AC.modinverses != 0 ) && ( ( par & NOINVERSES ) == 0 ) ) {
3170 y1 = AC.modinverses[a[1]];
3176 a[0] = (x*y1) % cmodvec[0];
3181 a[0] = a[0] % cmodvec[0];
3183 if ( a[0] == 0 ) { *na = 0;
return(0); }
3184 if ( ( AC.modmode & POSNEG ) != 0 ) {
3185 if ( a[0] > (UWORD)(cmodvec[0]/2) ) {
3186 a[0] = cmodvec[0] - a[0];
3190 else if ( *na < 0 ) {
3191 *na = 1; a[0] = cmodvec[0] - a[0];
3195 c = NumberMalloc(
"TakeModulus"); d = NumberMalloc(
"TakeModulus"); e = NumberMalloc(
"TakeModulus");
3196 f = NumberMalloc(
"TakeModulus"); g = NumberMalloc(
"TakeModulus"); h = NumberMalloc(
"TakeModulus");
3198 if ( DivLong(a,tnumer,(UWORD *)cmodvec,nmod,
3199 c,&nh,a,&tnumer) )
goto ModErr;
3200 if ( tnumer == 0 ) { *na = 0;
goto normalreturn; }
3201 if ( ( par & UNPACK ) == 0 ) {
3202 if ( ( AC.modmode & POSNEG ) != 0 ) {
3205 else if ( tnumer < 0 ) {
3206 SubPLon((UWORD *)cmodvec,nmod,a,-tnumer,a,&tnumer);
3211 if ( tdenom == 1 && a[n1] == 1 ) {
3212 if ( ( AC.modmode & POSNEG ) != 0 ) {
3215 else if ( tnumer < 0 ) {
3216 SubPLon((UWORD *)cmodvec,nmod,a,-tnumer,a,&tnumer);
3222 while ( --i > 0 ) *a++ = 0;
3225 if ( DivLong(a+n1,tdenom,(UWORD *)cmodvec,nmod,c,&nh,a+n1,&tdenom) )
goto ModErr;
3227 MLOCK(ErrorMessageLock);
3228 MesPrint(
"Division by zero in modulus arithmetic");
3229 if ( AP.DebugFlag ) {
3233 if ( i < 0 ) i = -i;
3234 while ( --i >= 0 ) { TalToLine((UWORD)(*a++)); TokenToLine((UBYTE *)
" "); }
3236 if ( i < 0 ) i = -i;
3237 while ( --i >= 0 ) { TalToLine((UWORD)(*a++)); TokenToLine((UBYTE *)
" "); }
3238 TalToLine((UWORD)(*na));
3242 MUNLOCK(ErrorMessageLock);
3243 NumberFree(c,
"TakeModulus"); NumberFree(d,
"TakeModulus"); NumberFree(e,
"TakeModulus");
3244 NumberFree(f,
"TakeModulus"); NumberFree(g,
"TakeModulus"); NumberFree(h,
"TakeModulus");
3247 if ( ( AC.modinverses != 0 ) && ( ( par & NOINVERSES ) == 0 )
3248 && ( tdenom == 1 || tdenom == -1 ) ) {
3249 *d = AC.modinverses[a[n1]]; y1 = 1; y2 = tdenom;
3250 if ( MulLong(a,tnumer,d,y1,c,&y3) )
goto ModErr;
3251 if ( DivLong(c,y3,(UWORD *)cmodvec,nmod,d,&y5,a,&tdenom) )
goto ModErr;
3252 if ( y2 < 0 ) tdenom = -tdenom;
3255 x2 = (UWORD *)cmodvec; x1 = c; i = nmod;
while ( --i >= 0 ) *x1++ = *x2++;
3256 x1 = c; x2 = a+n1; x3 = d; x4 = e; x5 = f; x6 = g;
3257 y1 = nmod; y2 = tdenom; y4 = 0; y5 = 1; *x5 = 1;
3259 if ( DivLong(x1,y1,x2,y2,h,&nh,x3,&y3) )
goto ModErr;
3260 if ( MulLong(x5,y5,h,nh,x6,&y6) )
goto ModErr;
3261 if ( AddLong(x4,y4,x6,-y6,x6,&y6) )
goto ModErr;
3263 if ( y2 != 1 || *x2 != 1 ) {
3264 MLOCK(ErrorMessageLock);
3265 MesPrint(
"Inverse in modulus arithmetic doesn't exist");
3266 MesPrint(
"Denominator and modulus are not relative prime");
3267 MUNLOCK(ErrorMessageLock);
3272 x7 = x1; x1 = x2; y1 = y2; x2 = x3; y2 = y3; x3 = x7;
3273 x8 = x4; x4 = x5; y4 = y5; x5 = x6; y5 = y6; x6 = x8;
3275 if ( y5 < 0 && AddLong((UWORD *)cmodvec,nmod,x5,y5,x5,&y5) )
goto ModErr;
3276 if ( MulLong(a,tnumer,x5,y5,c,&y3) )
goto ModErr;
3277 if ( DivLong(c,y3,(UWORD *)cmodvec,nmod,d,&y5,a,&tdenom) )
goto ModErr;
3279 if ( !tdenom ) { *na = 0;
goto normalreturn; }
3280 if ( ( ( AC.modmode & POSNEG ) != 0 ) && ( ( par & FROMFUNCTION ) == 0 ) ) {
3283 else if ( tdenom < 0 ) {
3284 SubPLon((UWORD *)cmodvec,nmod,a,-tdenom,a,&tdenom);
3290 while ( --i > 0 ) *a++ = 0;
3292 NumberFree(c,
"TakeModulus"); NumberFree(d,
"TakeModulus"); NumberFree(e,
"TakeModulus");
3293 NumberFree(f,
"TakeModulus"); NumberFree(g,
"TakeModulus"); NumberFree(h,
"TakeModulus");
3296 MLOCK(ErrorMessageLock);
3297 MesCall(
"TakeModulus");
3298 MUNLOCK(ErrorMessageLock);
3299 NumberFree(c,
"TakeModulus"); NumberFree(d,
"TakeModulus"); NumberFree(e,
"TakeModulus");
3300 NumberFree(f,
"TakeModulus"); NumberFree(g,
"TakeModulus"); NumberFree(h,
"TakeModulus");
3311 WORD TakeNormalModulus (UWORD *a, WORD *na, UWORD *c, WORD nc, WORD par)
3320 halfc = NumberMalloc(
"TakeNormalModulus");
3324 for (n=0; n<nhalfc; n++) {
3326 if (n+1<nc) halfc[n] |= c[n+1] << (BITSINWORD-1);
3329 if (halfc[nhalfc-1]==0)
3333 if (BigLong(a,ABS(*na),halfc,nhalfc) > 0) {
3335 TakeModulus(a,na,c,nc,par);
3338 if (BigLong(a,n,halfc,nhalfc) > 0) {
3339 SubPLon(c,nc,a,n,a,&n);
3340 *na = (*na > 0 ? -n : n);
3344 NumberFree(halfc,
"TakeNormalModulus");
3357 if ( AC.modpowers ) {
3358 M_free(AC.modpowers,
"AC.modpowers");
3359 AC.modpowers = NULL;
3362 MLOCK(ErrorMessageLock);
3363 MesPrint(
"&No memory for modulus generator power table");
3364 MUNLOCK(ErrorMessageLock);
3367 if ( n == 0 )
return(0);
3368 size = (LONG)(*AC.cmod);
3369 if ( n == 2 ) size += (((LONG)AC.cmod[1])<<BITSINWORD);
3370 AC.modpowers = (UWORD *)Malloc1(size*n*
sizeof(UWORD),
"table for powers of modulus");
3373 for ( i = 0; i < size; i++ ) AC.modpowers[i] = 0;
3374 for ( i = 0; i < size; i++ ) {
3375 AC.modpowers[j] = (WORD)i;
3379 for ( i = 2; i < size; i++ ) {
3380 if ( AC.modpowers[i] == 0 ) {
3381 MLOCK(ErrorMessageLock);
3382 MesPrint(
"&improper generator for this modulus");
3383 MUNLOCK(ErrorMessageLock);
3384 M_free(AC.modpowers,
"AC.modpowers");
3388 AC.modpowers[1] = 0;
3393 UWORD *MMscrat7 = NumberMalloc(
"MakeModTable"), *MMscratC = NumberMalloc(
"MakeModTable");
3397 for ( i = 0; i < j; i+=2 ) { AC.modpowers[i] = 0; AC.modpowers[i+1] = 0; }
3398 for ( i = 0; i < size; i++ ) {
3399 j = *MMscratC + (((LONG)MMscratC[1])<<BITSINWORD);
3401 AC.modpowers[j] = (WORD)(i & WORDMASK);
3402 AC.modpowers[j+1] = (WORD)(i >> BITSINWORD);
3403 MulLong((UWORD *)MMscratC,nScrat,(UWORD *)AC.powmod,
3404 AC.npowmod,(UWORD *)MMscrat7,&n2);
3405 TakeModulus(MMscrat7,&n2,AC.cmod,AC.ncmod,NOUNPACK);
3406 *MMscratC = *MMscrat7; MMscratC[1] = MMscrat7[1]; nScrat = n2;
3408 NumberFree(MMscrat7,
"MakeModTable"); NumberFree(MMscratC,
"MakeModTable");
3410 for ( i = 4; i < j; i+=2 ) {
3411 if ( AC.modpowers[i] == 0 && AC.modpowers[i+1] == 0 ) {
3412 MLOCK(ErrorMessageLock);
3413 MesPrint(
"&improper generator for this modulus");
3414 MUNLOCK(ErrorMessageLock);
3415 M_free(AC.modpowers,
"AC.modpowers");
3419 AC.modpowers[2] = AC.modpowers[3] = 0;
3439 int Factorial(PHEAD WORD n, UWORD *a, WORD *na)
3446 if ( n > AT.nfac ) {
3447 if ( AT.factorials == 0 ) {
3448 AT.nfac = 0; AT.mfac = 50; AT.sfact = 400;
3449 AT.pfac = (LONG *)Malloc1((AT.mfac+2)*
sizeof(LONG),
"factorials");
3450 AT.factorials = (UWORD *)Malloc1(AT.sfact*
sizeof(UWORD),
"factorials");
3451 AT.factorials[0] = 1; AT.pfac[0] = 0; AT.pfac[1] = 1;
3454 c = AT.factorials+AT.pfac[AT.nfac];
3455 nc = i = AT.pfac[AT.nfac+1] - AT.pfac[AT.nfac];
3456 while ( --i >= 0 ) *b++ = *c++;
3457 for ( j = AT.nfac+1; j <= n; j++ ) {
3459 if ( nc > AM.MaxTal ) {
3460 MLOCK(ErrorMessageLock);
3461 MesPrint(
"Overflow in factorial. MaxTal = %d",AM.MaxTal);
3462 MesPrint(
"Increase MaxTerm in %s",setupfilename);
3463 MUNLOCK(ErrorMessageLock);
3466 if ( j > AT.mfac ) {
3468 p = (LONG *)Malloc1((AT.mfac*2+2)*
sizeof(LONG),
"factorials");
3470 for ( i = AT.mfac+1; i >= 0; i-- ) p[i] = AT.pfac[i];
3471 M_free(AT.pfac,
"factorial pointers"); AT.pfac = p; AT.mfac *= 2;
3473 if ( AT.pfac[j] + nc >= AT.sfact ) {
3475 f = (UWORD *)Malloc1(AT.sfact*2*
sizeof(UWORD),
"factorials");
3477 c = AT.factorials; b = f;
3478 while ( --ii >= 0 ) *b++ = *c++;
3479 M_free(AT.factorials,
"factorials");
3483 b = a; c = AT.factorials + AT.pfac[j]; i = nc;
3484 while ( --i >= 0 ) *c++ = *b++;
3485 AT.pfac[j+1] = AT.pfac[j] + nc;
3490 else if ( n == 0 ) {
3494 *na = i = AT.pfac[n+1] - AT.pfac[n];
3495 b = AT.factorials + AT.pfac[n];
3496 while ( --i >= 0 ) *a++ = *b++;
3519 int Bernoulli(WORD n, UWORD *a, WORD *na)
3522 UWORD *b, *c, *scrib, *ntop, *ntop1;
3523 WORD i, i1, i2, nhalf, nqua, nscrib, nntop, nntop1, *oldworkpointer;
3524 UWORD twee = 2, twonplus1;
3528 if ( n == 0 ) { a[0] = a[1] = 1; *na = 3; }
3529 else if ( n == 1 ) { a[0] = 1; a[1] = 2; *na = 3; }
3532 if ( ( n & 1 ) != 0 ) { a[0] = a[1] = 0; *na = 0;
return(0); }
3534 if ( nhalf > AT.nBer ) {
3535 oldworkpointer = AT.WorkPointer;
3536 if ( AT.bernoullis == 0 ) {
3537 AT.nBer = 1; AT.mBer = 50; AT.sBer = 400;
3538 AT.pBer = (LONG *)Malloc1((AT.mBer+2)*
sizeof(LONG),
"bernoullis");
3539 AT.bernoullis = (UWORD *)Malloc1(AT.sBer*
sizeof(UWORD),
"bernoullis");
3540 AT.pBer[1] = 0; AT.pBer[2] = 3;
3541 AT.bernoullis[0] = 3; AT.bernoullis[1] = 1; AT.bernoullis[2] = 12;
3543 a[0] = 1; a[1] = 12; *na = 3;
return(0);
3546 while ( nhalf > AT.mBer ) {
3548 p = (LONG *)Malloc1((AT.mBer*2+1)*
sizeof(LONG),
"bernoullis");
3550 for ( i = AT.mBer; i >= 0; i-- ) p[i] = AT.pBer[i];
3551 M_free(AT.pBer,
"factorial pointers"); AT.pBer = p; AT.mBer *= 2;
3553 for ( n = AT.nBer+1; n <= nhalf; n++ ) {
3554 scrib = (UWORD *)(AT.WorkPointer);
3556 if ( ( n & 1 ) == 1 ) {
3557 nscrib = 0; ntop = scrib;
3560 b = AT.bernoullis + AT.pBer[nqua];
3562 i = (WORD)(REDLENG(nscrib));
3563 MulRat(BHEAD b,i,b,i,scrib,&nscrib);
3564 ntop = scrib + 2*nscrib;
3567 for ( j = 1; j <= nqua; j++ ) {
3568 b = AT.bernoullis + AT.pBer[j];
3569 c = AT.bernoullis + AT.pBer[n-j];
3570 i1 = (WORD)(*b); i2 = (WORD)(*c);
3573 MulRat(BHEAD b+1,i1,c+1,i2,ntop,&nntop);
3574 Mully(BHEAD ntop,&nntop,&twee,1);
3576 i = (WORD)nntop;
if ( i < 0 ) i = -i;
3578 AddRat(BHEAD ntop,nntop,scrib,nscrib,ntop1,&nntop1);
3581 ntop1 = ntop; nntop1 = nntop;
3583 nscrib = i1 = (WORD)nntop1;
3584 if ( i1 < 0 ) i1 = - i1;
3586 for ( i = 0; i < i1; i++ ) scrib[i] = ntop1[i];
3590 Divvy(BHEAD scrib,&nscrib,&twonplus1,-1);
3591 i1 = INCLENG(nscrib);
3592 i2 = i1;
if ( i2 < 0 ) i2 = -i2;
3593 i = (WORD)(AT.bernoullis[AT.pBer[n-1]]);
3594 if ( i < 0 ) i = -i;
3595 AT.pBer[n] = AT.pBer[n-1]+i;
3596 if ( AT.pBer[n] + i2 >= AT.sBer ) {
3598 f = (UWORD *)Malloc1(AT.sBer*2*
sizeof(UWORD),
"bernoullis");
3600 c = AT.bernoullis; b = f;
3601 while ( --ii >= 0 ) *b++ = *c++;
3602 M_free(AT.bernoullis,
"bernoullis");
3606 c = AT.bernoullis + AT.pBer[n]; b = scrib;
3608 for ( i = 1; i < i2; i++ ) *c++ = *b++;
3611 AT.WorkPointer = oldworkpointer;
3613 b = AT.bernoullis + AT.pBer[nhalf];
3614 *na = i = (WORD)(*b++);
3615 if ( i < 0 ) i = -i;
3617 while ( --i >= 0 ) *a++ = *b++;
3636 #if ( BITSINWORD == 32 ) 3638 void StartPrimeList(PHEAD0)
3641 AR.PrimeList[AR.numinprimelist++] = 3;
3642 for ( i = 5; i < 46340; i += 2 ) {
3643 for ( j = 0; j < AR.numinprimelist && AR.PrimeList[j]*AR.PrimeList[j] <= i; j++ ) {
3644 if ( i % AR.PrimeList[j] == 0 )
goto nexti;
3646 AR.PrimeList[AR.numinprimelist++] = i;
3649 AR.notfirstprime = 1;
3659 #if ( BITSINWORD == 32 ) 3660 if ( AR.notfirstprime == 0 ) StartPrimeList(BHEAD0);
3662 if ( num > AT.inprimelist ) {
3663 while ( AT.inprimelist < num ) {
3664 if ( num >= AT.sizeprimelist ) {
3665 if ( AT.sizeprimelist == 0 ) newsize = 32;
3666 else newsize = 2*AT.sizeprimelist;
3667 while ( num >= newsize ) newsize = newsize*2;
3668 newpl = (WORD *)Malloc1(newsize*
sizeof(WORD),
"NextPrime");
3669 for ( i = 0; i < AT.sizeprimelist; i++ ) {
3670 newpl[i] = AT.primelist[i];
3672 if ( AT.sizeprimelist > 0 ) {
3673 M_free(AT.primelist,
"NextPrime");
3675 AT.sizeprimelist = newsize;
3676 AT.primelist = newpl;
3678 if ( AT.inprimelist < 0 ) { i = MAXPOSITIVE; }
3679 else { i = AT.primelist[AT.inprimelist]; }
3680 while ( i > MAXPOWER ) {
3682 #if ( BITSINWORD == 32 ) 3683 for ( j = 0; j < AR.numinprimelist && AR.PrimeList[j]*(LONG)(AR.PrimeList[j]) <= x; j++ ) {
3684 if ( x % AR.PrimeList[j] == 0 )
goto nexti;
3687 for ( j = 3; j*((LONG)j) <= x; j += 2 ) {
3688 if ( x % j == 0 )
goto nexti;
3692 AT.primelist[AT.inprimelist] = i;
3696 if ( i < MAXPOWER ) {
3697 MLOCK(ErrorMessageLock);
3698 MesPrint(
"There are not enough short prime numbers for this calculation");
3699 MesPrint(
"Try to use a computer with a %d-bits architecture",
3700 (
int)(BITSINWORD*4));
3701 MUNLOCK(ErrorMessageLock);
3706 return(AT.primelist[num]);
3718 WORD Moebius(PHEAD WORD nn)
3723 #if ( BITSINWORD == 32 ) 3724 if ( AR.notfirstprime == 0 ) StartPrimeList(BHEAD0);
3731 if ( nn >= AR.moebiustablesize ) {
3732 if ( AR.moebiustablesize <= 0 ) { newsize = nn + 20; }
3733 else { newsize = nn*2; }
3734 if ( newsize > MAXPOSITIVE ) newsize = MAXPOSITIVE;
3735 newtable = (
char *)Malloc1(newsize*
sizeof(
char),
"Moebius");
3736 for ( i = 0; i < AR.moebiustablesize; i++ ) newtable[i] = AR.moebiustable[i];
3737 for ( ; i < newsize; i++ ) newtable[i] = 2;
3738 if ( AR.moebiustablesize > 0 ) M_free(AR.moebiustable,
"Moebius");
3739 AR.moebiustable = newtable;
3740 AR.moebiustablesize = newsize;
3742 if ( AR.moebiustable[nn] != 2 )
return((WORD)AR.moebiustable[nn]);
3744 if ( n == 1 )
goto putvalue;
3747 if ( n % 2 == 0 ) { mu = 0;
goto putvalue; }
3748 if ( AR.moebiustable[n] != 2 ) { mu = -AR.moebiustable[n];
goto putvalue; }
3750 if ( n == 1 )
goto putvalue;
3752 for ( i = 0; i < AR.numinprimelist; i++ ) {
3753 x = AR.PrimeList[i];
3756 if ( n % x == 0 ) { mu = 0;
goto putvalue; }
3757 if ( AR.moebiustable[n] != 2 ) { mu = -AR.moebiustable[n];
goto putvalue; }
3759 if ( n == 1 )
goto putvalue;
3761 if ( n < x*x )
break;
3765 AR.moebiustable[nn] = mu;
3788 static void wranfnew(PHEAD0)
3792 for ( i = 0; i < AR.wranfnpair1; i++ ) {
3793 j = AR.wranfia[i] - AR.wranfia[i+(AR.wranfnpair2-AR.wranfnpair1)];
3794 if ( j < 0 ) j += (LONG)1 << (2*BITSINWORD-2);
3797 for ( i = AR.wranfnpair1; i < AR.wranfnpair2; i++ ) {
3798 j = AR.wranfia[i] - AR.wranfia[i-AR.wranfnpair1];
3799 if ( j < 0 ) j += (LONG)1 << (2*BITSINWORD-2);
3804 void iniwranf(PHEAD0)
3806 int imax = AR.wranfnpair2-1;
3807 ULONG i, ii, seed = AR.wranfseed;
3809 ULONG offset = 12345;
3812 #if defined(WITHPTHREADS) 3814 #elif defined(WITHMPI) 3821 pow = offset; accu = 1;
3823 if ( ( i & 1 ) != 0 ) accu *= pow;
3824 i /= 2; pow = pow*pow;
3829 if ( seed < ((LONG)1<<(BITSINWORD-1)) ) {
3830 j = ( (seed+31459L) << (BITSINWORD-2))+offset;
3832 else if ( seed < ((LONG)1<<(BITSINWORD+10-1)) ) {
3833 j = ( (seed+31459L) << (BITSINWORD-10-2))+offset;
3836 j = ( (seed+31459L) << 1)+offset;
3838 if ( ( seed & 1 ) == 1 ) seed++;
3840 AR.wranfia[imax] = j;
3842 for ( i = 0; i <= (ULONG)(imax); i++ ) {
3843 ii = (AR.wranfnpair1*i)%AR.wranfnpair2;
3845 k = ULongToLong((ULONG)j - (ULONG)k);
3846 if ( k < 0 ) k += (LONG)1 << (2*BITSINWORD-2);
3849 for ( i = 0; i < WARMUP; i++ ) wranfnew(BHEAD0);
3856 if ( AR.wranfia == 0 ) {
3857 AR.wranfia = (ULONG *)Malloc1(AR.wranfnpair2*
sizeof(ULONG),
"wranf");
3860 if ( AR.wranfcall >= AR.wranfnpair2) {
3864 wval = (UWORD)(AR.wranfia[AR.wranfcall++]>>(BITSINWORD-1));
3872 UWORD iranf(PHEAD UWORD imax)
3876 if (imax < 2)
return 0;
3877 x = (LONG)1 << BITSINWORD;
3879 while ( ( i = wranf(BHEAD0) ) >= xmax ) {}
3897 UBYTE *PreRandom(UBYTE *s)
3900 UBYTE *mode,*mins = 0,*maxs = 0, *outval;
3902 double minval, maxval, value = 0;
3905 while ( FG.cTable[*s] <= 1 ) s++;
3906 if ( *s ==
',' ) { *s = 0; s++; }
3908 while ( *s && *s !=
',' ) s++;
3909 if ( *s ==
',' ) { *s = 0; s++; }
3911 while ( *s && *s !=
',' ) s++;
3912 if ( *s || *maxs == 0 || *mins == 0 ) {
3913 MesPrint(
"@Illegal arguments in macro RANDOM_");
3916 if ( StrICmp(mode,(UBYTE *)
"lin") == 0 ) {
3919 else if ( StrICmp(mode,(UBYTE *)
"log") == 0 ) {
3923 MesPrint(
"@Illegal mode argument in macro RANDOM_");
3927 sscanf((
char *)mins,
"%f",&num); minval = num;
3928 sscanf((
char *)maxs,
"%f",&num); maxval = num;
3940 if ( PF.me == MASTER ) {
3947 xx = x/pow(2.0,(
double)(BITSINWORD-1));
3948 if ( linlog == 0 ) {
3949 value = minval + (maxval-minval)*xx;
3951 else if ( linlog == 1 ) {
3952 value = minval * pow(maxval/minval,xx);
3956 outval = (UBYTE *)Malloc1(64,
"PreRandom");
3957 if ( ABS(value) < 0.00001 || ABS(value) > 1000000. ) {
3958 sprintf((
char *)outval,
"%e",value);
3960 else if ( ABS(value) < 0.0001 ) { sprintf((
char *)outval,
"%10f",value); }
3961 else if ( ABS(value) < 0.001 ) { sprintf((
char *)outval,
"%9f",value); }
3962 else if ( ABS(value) < 0.01 ) { sprintf((
char *)outval,
"%8f",value); }
3963 else if ( ABS(value) < 0.1 ) { sprintf((
char *)outval,
"%7f",value); }
3964 else if ( ABS(value) < 1. ) { sprintf((
char *)outval,
"%6f",value); }
3965 else if ( ABS(value) < 10. ) { sprintf((
char *)outval,
"%5f",value); }
3966 else if ( ABS(value) < 100. ) { sprintf((
char *)outval,
"%4f",value); }
3967 else if ( ABS(value) < 1000. ) { sprintf((
char *)outval,
"%3f",value); }
3968 else if ( ABS(value) < 10000. ) { sprintf((
char *)outval,
"%2f",value); }
3969 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)