FORM v5.0.1-33-gdf7fc94
reken.c
Go to the documentation of this file.
1
11/* #[ License : */
12/*
13 * Copyright (C) 1984-2026 J.A.M. Vermaseren
14 * When using this file you are requested to refer to the publication
15 * J.A.M.Vermaseren "New features of FORM" math-ph/0010025
16 * This is considered a matter of courtesy as the development was paid
17 * for by FOM the Dutch physics granting agency and we would like to
18 * be able to track its scientific use to convince FOM of its value
19 * for the community.
20 *
21 * This file is part of FORM.
22 *
23 * FORM is free software: you can redistribute it and/or modify it under the
24 * terms of the GNU General Public License as published by the Free Software
25 * Foundation, either version 3 of the License, or (at your option) any later
26 * version.
27 *
28 * FORM is distributed in the hope that it will be useful, but WITHOUT ANY
29 * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
30 * FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
31 * details.
32 *
33 * You should have received a copy of the GNU General Public License along
34 * with FORM. If not, see <http://www.gnu.org/licenses/>.
35 */
36/* #] License : */
37/*
38 #[ Includes : reken.c
39*/
40
41#include "form3.h"
42#include <math.h>
43
44#ifdef WITHGMP
45#include <gmp.h>
46#define GMPSPREAD (GMP_LIMB_BITS/BITSINWORD)
47#endif
48
49#define GCDMAX 3
50
51#define NEWTRICK 1
52/*
53 #] Includes :
54 #[ RekenRational :
55 #[ Pack : void Pack(a,na,b,nb)
56
57 Packs the contents of the numerator a and the denominator b into
58 one normalized fraction a.
59
60*/
61
62void Pack(UWORD *a, WORD *na, UWORD *b, WORD nb)
63{
64 WORD c, sgn = 1, i;
65 UWORD *to,*from;
66 if ( (c = *na) == 0 ) {
67/* INTERNAL_ERROR_EXCL_START */
68 MLOCK(ErrorMessageLock);
69 MesPrint("!>Caught a zero in Pack");
70 MUNLOCK(ErrorMessageLock);
71 return;
72/* INTERNAL_ERROR_EXCL_STOP */
73 }
74 if ( nb == 0 ) {
75/* INTERNAL_ERROR_EXCL_START */
76 MLOCK(ErrorMessageLock);
77 MesPrint("!>Division by zero in Pack");
78 MUNLOCK(ErrorMessageLock);
79 return;
80/* INTERNAL_ERROR_EXCL_STOP */
81 }
82 if ( *na < 0 ) { sgn = -sgn; c = -c; }
83 if ( nb < 0 ) { sgn = -sgn; nb = -nb; }
84 *na = MaX(c,nb);
85 to = a + c;
86 i = *na - c;
87 while ( --i >= 0 ) *to++ = 0;
88 i = *na - nb;
89 from = b;
90 NCOPY(to,from,nb);
91 while ( --i >= 0 ) *to++ = 0;
92 if ( sgn < 0 ) *na = -*na;
93}
94
95/*
96 #] Pack :
97 #[ UnPack : void UnPack(a,na,denom,numer)
98
99 Determines the sizes of the numerator and the denominator in the
100 normalized fraction a with length na.
101
102*/
103
104void UnPack(UWORD *a, WORD na, WORD *denom, WORD *numer)
105{
106 UWORD *pos;
107 WORD i, sgn = na;
108 if ( na < 0 ) { na = -na; }
109 i = na;
110 if ( i > 1 ) { /* Find the respective leading words */
111 a += i;
112 a--;
113 pos = a + i;
114 while ( !(*a) ) { i--; a--; }
115 while ( !(*pos) ) { na--; pos--; }
116 }
117 *denom = na;
118 if ( sgn < 0 ) i = -i;
119 *numer = i;
120}
121
122/*
123 #] UnPack :
124 #[ Mully : WORD Mully(a,na,b,nb)
125
126 Multiplies the rational a by the Long b.
127
128*/
129
130int Mully(PHEAD UWORD *a, WORD *na, UWORD *b, WORD nb)
131{
132 GETBIDENTITY
133 UWORD *d, *e;
134 WORD i, sgn = 1;
135 WORD nd, ne, adenom, anumer;
136 if ( !nb ) { *na = 0; return(0); }
137 else if ( *b == 1 ) {
138 if ( nb == 1 ) return(0);
139 else if ( nb == -1 ) { *na = -*na; return(0); }
140 }
141 if ( *na < 0 ) { sgn = -sgn; *na = -*na; }
142 if ( nb < 0 ) { sgn = -sgn; nb = -nb; }
143 UnPack(a,*na,&adenom,&anumer);
144 d = NumberMalloc("Mully"); e = NumberMalloc("Mully");
145 for ( i = 0; i < nb; i++ ) { e[i] = *b++; }
146 ne = nb;
147 if ( Simplify(BHEAD a+*na,&adenom,e,&ne) ) goto MullyEr;
148 if ( MulLong(a,anumer,e,ne,d,&nd) ) goto MullyEr;
149 b = a+*na;
150 for ( i = 0; i < *na; i++ ) { e[i] = *b++; }
151 ne = adenom;
152 *na = nd;
153 b = d;
154 *na = nd;
155 for ( i = 0; i < *na; i++ ) { a[i] = *b++; }
156 Pack(a,na,e,ne);
157 if ( sgn < 0 ) *na = -*na;
158 NumberFree(d,"Mully"); NumberFree(e,"Mully");
159 return(0);
160MullyEr:
161 MLOCK(ErrorMessageLock);
162 MesCall("Mully");
163 MUNLOCK(ErrorMessageLock);
164 NumberFree(d,"Mully"); NumberFree(e,"Mully");
165 SETERROR(-1)
166}
167
168/*
169 #] Mully :
170 #[ Divvy : WORD Divvy(a,na,b,nb)
171
172 Divides the rational a by the Long b.
173
174*/
175
176int Divvy(PHEAD UWORD *a, WORD *na, UWORD *b, WORD nb)
177{
178 GETBIDENTITY
179 UWORD *d,*e;
180 WORD i, sgn = 1;
181 WORD nd, ne, adenom, anumer;
182 if ( !nb ) {
183/* INTERNAL_ERROR_EXCL_START */
184 MLOCK(ErrorMessageLock);
185 MesPrint("!>Division by zero in Divvy");
186 MUNLOCK(ErrorMessageLock);
187 return(-1);
188/* INTERNAL_ERROR_EXCL_STOP */
189 }
190 d = NumberMalloc("Divvy"); e = NumberMalloc("Divvy");
191 if ( nb < 0 ) { sgn = -sgn; nb = -nb; }
192 if ( *na < 0 ) { sgn = -sgn; *na = -*na; }
193 UnPack(a,*na,&adenom,&anumer);
194 for ( i = 0; i < nb; i++ ) { e[i] = *b++; }
195 ne = nb;
196 if ( Simplify(BHEAD a,&anumer,e,&ne) ) goto DivvyEr;
197 if ( MulLong(a+*na,adenom,e,ne,d,&nd) ) goto DivvyEr;
198 *na = anumer;
199 Pack(a,na,d,nd);
200 if ( sgn < 0 ) *na = -*na;
201 NumberFree(d,"Divvy"); NumberFree(e,"Divvy");
202 return(0);
203DivvyEr:
204 MLOCK(ErrorMessageLock);
205 MesCall("Divvy");
206 MUNLOCK(ErrorMessageLock);
207 NumberFree(d,"Divvy"); NumberFree(e,"Divvy");
208 SETERROR(-1)
209}
210
211/*
212 #] Divvy :
213 #[ AddRat : WORD AddRat(a,na,b,nb,c,nc)
214*/
215
216int AddRat(PHEAD UWORD *a, WORD na, UWORD *b, WORD nb, UWORD *c, WORD *nc)
217{
218 GETBIDENTITY
219 UWORD *d, *e, *f, *g;
220 WORD nd, ne, nf, ng, adenom, anumer, bdenom, bnumer;
221 if ( !na ) {
222 WORD i;
223 *nc = nb;
224 if ( nb < 0 ) nb = -nb;
225 nb *= 2;
226 for ( i = 0; i < nb; i++ ) *c++ = *b++;
227 return(0);
228 }
229 else if ( !nb ) {
230 WORD i;
231 *nc = na;
232 if ( na < 0 ) na = -na;
233 na *= 2;
234 for ( i = 0; i < na; i++ ) *c++ = *a++;
235 return(0);
236 }
237 else if ( b[1] == 1 && a[1] == 1 ) {
238 if ( na == 1 ) {
239 if ( nb == 1 ) {
240 *c = *a + *b;
241 c[1] = 1;
242 if ( *c < *a ) { c[2] = 1; c[3] = 0; *nc = 2; }
243 else { *nc = 1; }
244 return(0);
245 }
246 else if ( nb == -1 ) {
247 if ( *b > *a ) {
248 *c = *b - *a; *nc = -1;
249 }
250 else if ( *b < *a ) {
251 *c = *a - *b; *nc = 1;
252 }
253 else *nc = 0;
254 c[1] = 1;
255 return(0);
256 }
257 }
258 else if ( na == -1 ){
259 if ( nb == -1 ) {
260 c[1] = 1;
261 *c = *a + *b;
262 if ( *c < *a ) { c[2] = 1; c[3] = 0; *nc = -2; }
263 else { *nc = -1; }
264 return(0);
265 }
266 else if ( nb == 1 ) {
267 if ( *b > *a ) {
268 *c = *b - *a; *nc = 1;
269 }
270 else if ( *b < *a ) {
271 *c = *a - *b; *nc = -1;
272 }
273 else *nc = 0;
274 c[1] = 1;
275 return(0);
276 }
277 }
278 }
279 UnPack(a,na,&adenom,&anumer);
280 UnPack(b,nb,&bdenom,&bnumer);
281 if ( na < 0 ) na = -na;
282 if ( nb < 0 ) nb = -nb;
283 if ( na == 1 && nb == 1 ) {
284 RLONG t1, t2, t3;
285 t3 = ((RLONG)a[1])*((RLONG)b[1]);
286 t1 = ((RLONG)a[0])*((RLONG)b[1]);
287 t2 = ((RLONG)a[1])*((RLONG)b[0]);
288 if ( ( anumer > 0 && bnumer > 0 ) || ( anumer < 0 && bnumer < 0 ) ) {
289 if ( ( t1 = t1 + t2 ) < t2 ) {
290 c[2] = 1;
291 c[0] = (UWORD)t1;
292 c[1] = (UWORD)(t1 >> BITSINWORD);
293 *nc = 3;
294 }
295 else {
296 c[0] = (UWORD)t1;
297 if ( ( c[1] = (UWORD)(t1 >> BITSINWORD) ) != 0 ) *nc = 2;
298 else *nc = 1;
299 }
300 }
301 else {
302 if ( t1 == t2 ) { *nc = 0; return(0); }
303 if ( t1 > t2 ) {
304 t1 -= t2;
305 }
306 else {
307 t1 = t2 - t1;
308 anumer = -anumer;
309 }
310 c[0] = (UWORD)t1;
311 if ( ( c[1] = (UWORD)(t1 >> BITSINWORD) ) != 0 ) *nc = 2;
312 else *nc = 1;
313 }
314 if ( anumer < 0 ) *nc = -*nc;
315 d = NumberMalloc("AddRat");
316 d[0] = (UWORD)t3;
317 if ( ( d[1] = (UWORD)(t3 >> BITSINWORD) ) != 0 ) nd = 2;
318 else nd = 1;
319 if ( Simplify(BHEAD c,nc,d,&nd) ) goto AddRer1;
320 }
321/*
322 else if ( a[na] == 1 && b[nb] == 1 && adenom == 1 && bdenom == 1 ) {
323 if ( AddLong(a,na,b,nb,c,&nc) ) goto AddRer2;
324 i = ABS(nc); d = c + i; *d++ = 1;
325 while ( --i > 0 ) *d++ = 0 ;
326 return(0);
327 }
328*/
329 else {
330 d = NumberMalloc("AddRat"); e = NumberMalloc("AddRat");
331 f = NumberMalloc("AddRat"); g = NumberMalloc("AddRat");
332 if ( GcdLong(BHEAD a+na,adenom,b+nb,bdenom,d,&nd) ) goto AddRer;
333 if ( *d == 1 && nd == 1 ) nd = 0;
334 if ( nd ) {
335 if ( DivLong(a+na,adenom,d,nd,e,&ne,c,nc) ) goto AddRer;
336 if ( DivLong(b+nb,bdenom,d,nd,f,&nf,c,nc) ) goto AddRer;
337 if ( MulLong(a,anumer,f,nf,c,nc) ) goto AddRer;
338 if ( MulLong(b,bnumer,e,ne,g,&ng) ) goto AddRer;
339 }
340 else {
341 if ( MulLong(a+na,adenom,b,bnumer,c,nc) ) goto AddRer;
342 if ( MulLong(b+nb,bdenom,a,anumer,g,&ng) ) goto AddRer;
343 }
344 if ( AddLong(c,*nc,g,ng,c,nc) ) goto AddRer;
345 if ( !*nc ) {
346 NumberFree(g,"AddRat"); NumberFree(f,"AddRat");
347 NumberFree(e,"AddRat"); NumberFree(d,"AddRat");
348 return(0);
349 }
350 if ( nd ) {
351 if ( Simplify(BHEAD c,nc,d,&nd) ) goto AddRer;
352 if ( MulLong(e,ne,d,nd,g,&ng) ) goto AddRer;
353 if ( MulLong(g,ng,f,nf,d,&nd) ) goto AddRer;
354 }
355 else {
356 if ( MulLong(a+na,adenom,b+nb,bdenom,d,&nd) ) goto AddRer;
357 }
358 NumberFree(g,"AddRat"); NumberFree(f,"AddRat"); NumberFree(e,"AddRat");
359 }
360 Pack(c,nc,d,nd);
361 NumberFree(d,"AddRat");
362 return(0);
363AddRer:
364 NumberFree(g,"AddRat"); NumberFree(f,"AddRat"); NumberFree(e,"AddRat");
365AddRer1:
366 NumberFree(d,"AddRat");
367/* AddRer2: */
368 MLOCK(ErrorMessageLock);
369 MesCall("AddRat");
370 MUNLOCK(ErrorMessageLock);
371 SETERROR(-1)
372}
373
374/*
375 #] AddRat :
376 #[ MulRat : WORD MulRat(a,na,b,nb,c,nc)
377
378 Multiplies the rationals a and b. The Gcd of the individual
379 pieces is divided out first to minimize the chances of spurious
380 overflows.
381
382*/
383
384int MulRat(PHEAD UWORD *a, WORD na, UWORD *b, WORD nb, UWORD *c, WORD *nc)
385{
386 WORD i;
387 WORD sgn = 1;
388 if ( *b == 1 && b[1] == 1 ) {
389 if ( nb == 1 ) {
390 *nc = na;
391 i = ABS(na); i *= 2;
392 while ( --i >= 0 ) *c++ = *a++;
393 return(0);
394 }
395 else if ( nb == -1 ) {
396 *nc = - na;
397 i = ABS(na); i *= 2;
398 while ( --i >= 0 ) *c++ = *a++;
399 return(0);
400 }
401 }
402 if ( *a == 1 && a[1] == 1 ) {
403 if ( na == 1 ) {
404 *nc = nb;
405 i = ABS(nb); i *= 2;
406 while ( --i >= 0 ) *c++ = *b++;
407 return(0);
408 }
409 else if ( na == -1 ) {
410 *nc = - nb;
411 i = ABS(nb); i *= 2;
412 while ( --i >= 0 ) *c++ = *b++;
413 return(0);
414 }
415 }
416 if ( na < 0 ) { na = -na; sgn = -sgn; }
417 if ( nb < 0 ) { nb = -nb; sgn = -sgn; }
418 if ( !na || !nb ) { *nc = 0; return(0); }
419 if ( na != 1 || nb != 1 ) {
420 GETBIDENTITY
421 UWORD *xd,*xe, *xf,*xg;
422 WORD dden, dnumr, eden, enumr;
423 UnPack(a,na,&dden,&dnumr);
424 UnPack(b,nb,&eden,&enumr);
425 xd = NumberMalloc("MulRat"); xf = NumberMalloc("MulRat");
426 for ( i = 0; i < dnumr; i++ ) xd[i] = a[i];
427 a += na;
428 for ( i = 0; i < dden; i++ ) xf[i] = a[i];
429 xe = NumberMalloc("MulRat"); xg = NumberMalloc("MulRat");
430 for ( i = 0; i < enumr; i++ ) xe[i] = b[i];
431 b += nb;
432 for ( i = 0; i < eden; i++ ) xg[i] = b[i];
433 if ( Simplify(BHEAD xd,&dnumr,xg,&eden) ||
434 Simplify(BHEAD xe,&enumr,xf,&dden) ||
435 MulLong(xd,dnumr,xe,enumr,c,nc) ||
436 MulLong(xf,dden,xg,eden,xd,&dnumr) ) {
437 MLOCK(ErrorMessageLock);
438 MesCall("MulRat");
439 MUNLOCK(ErrorMessageLock);
440 NumberFree(xd,"MulRat"); NumberFree(xe,"MulRat"); NumberFree(xf,"MulRat"); NumberFree(xg,"MulRat");
441 SETERROR(-1)
442 }
443 Pack(c,nc,xd,dnumr);
444 NumberFree(xd,"MulRat"); NumberFree(xe,"MulRat"); NumberFree(xf,"MulRat"); NumberFree(xg,"MulRat");
445 }
446 else {
447 UWORD y;
448 UWORD a0,a1,b0,b1;
449 RLONG xx;
450 y = a[0]; b1=b[1];
451 do { a0 = y % b1; y = b1; } while ( ( b1 = a0 ) != 0 );
452 if ( y != 1 ) {
453 a0 = a[0] / y;
454 b1 = b[1] / y;
455 }
456 else {
457 a0 = a[0];
458 b1 = b[1];
459 }
460 y=b[0]; a1=a[1];
461 do { b0 = y % a1; y = a1; } while ( ( a1 = b0 ) != 0 );
462 if ( y != 1 ) {
463 a1 = a[1] / y;
464 b0 = b[0] / y;
465 }
466 else {
467 a1 = a[1];
468 b0 = b[0];
469 }
470 xx = ((RLONG)a0)*b0;
471 if ( xx & AWORDMASK ) {
472 *nc = 2;
473 c[0] = (UWORD)xx;
474 c[1] = (UWORD)(xx >> BITSINWORD);
475 xx = ((RLONG)a1)*b1;
476 c[2] = (UWORD)xx;
477 c[3] = (UWORD)(xx >> BITSINWORD);
478 }
479 else {
480 c[0] = (UWORD)xx;
481 xx = ((RLONG)a1)*b1;
482 if ( xx & AWORDMASK ) {
483 c[1] = 0;
484 c[2] = (UWORD)xx;
485 c[3] = (UWORD)(xx >> BITSINWORD);
486 *nc = 2;
487 }
488 else {
489 c[1] = (UWORD)xx;
490 *nc = 1;
491 }
492 }
493 }
494 if ( sgn < 0 ) *nc = -*nc;
495 return(0);
496}
497
498/*
499 #] MulRat :
500 #[ DivRat : WORD DivRat(a,na,b,nb,c,nc)
501
502 Divides the rational a by the rational b.
503
504*/
505
506int DivRat(PHEAD UWORD *a, WORD na, UWORD *b, WORD nb, UWORD *c, WORD *nc)
507{
508 GETBIDENTITY
509 WORD i, j;
510 int ret;
511 UWORD *xd,*xe,xx;
512 if ( !nb ) {
513/* INTERNAL_ERROR_EXCL_START */
514 MLOCK(ErrorMessageLock);
515 MesPrint("!>Rational division by zero");
516 MUNLOCK(ErrorMessageLock);
517 return(-1);
518/* INTERNAL_ERROR_EXCL_STOP */
519 }
520 j = i = (nb >= 0)? nb: -nb;
521 xd = b; xe = b + i;
522 do { xx = *xd; *xd++ = *xe; *xe++ = xx; } while ( --j > 0 );
523 ret = MulRat(BHEAD a,na,b,nb,c,nc);
524 xd = b; xe = b + i;
525 do { xx = *xd; *xd++ = *xe; *xe++ = xx; } while ( --i > 0 );
526 return(ret);
527}
528
529/*
530 #] DivRat :
531 #[ Simplify : WORD Simplify(a,na,b,nb)
532
533 Determines the greatest common denominator of a and b and
534 divides both by it. A possible sign is put in a. This is
535 the simplification of the fraction a/b.
536
537*/
538
539int Simplify(PHEAD UWORD *a, WORD *na, UWORD *b, WORD *nb)
540{
541 GETBIDENTITY
542 UWORD *x1,*x2,*x3;
543 UWORD *x4;
544 WORD n1,n2,n3,n4,sgn = 1;
545 WORD i;
546 UWORD *Siscrat5, *Siscrat6, *Siscrat7, *Siscrat8;
547 if ( *na < 0 ) { *na = -*na; sgn = -sgn; }
548 if ( *nb < 0 ) { *nb = -*nb; sgn = -sgn; }
549 Siscrat5 = NumberMalloc("Simplify"); Siscrat6 = NumberMalloc("Simplify");
550 Siscrat7 = NumberMalloc("Simplify"); Siscrat8 = NumberMalloc("Simplify");
551 x1 = Siscrat8; x2 = Siscrat7;
552 if ( *nb == 1 ) {
553 x3 = Siscrat6;
554 if ( DivLong(a,*na,b,*nb,x1,&n1,x2,&n2) ) goto SimpErr;
555 if ( !n2 ) {
556 for ( i = 0; i < n1; i++ ) *a++ = *x1++;
557 *na = n1;
558 *b = 1;
559 }
560 else {
561 UWORD y1, y2, y3;
562 y2 = *b;
563 y3 = *x2;
564 do { y1 = y2 % y3; y2 = y3; } while ( ( y3 = y1 ) != 0 );
565 if ( ( *x2 = y2 ) != 1 ) {
566 *b /= y2;
567 if ( DivLong(a,*na,x2,(WORD)1,x1,&n1,x3,&n3) ) goto SimpErr;
568 for ( i = 0; i < n1; i++ ) *a++ = *x1++;
569 *na = n1;
570 }
571 }
572 }
573#ifdef NEWTRICK
574 else if ( *na >= GCDMAX && *nb >= GCDMAX ) {
575 n1 = i = *na; x3 = a;
576 NCOPY(x1,x3,i);
577 x3 = b; n2 = i = *nb;
578 NCOPY(x2,x3,i);
579 x4 = Siscrat5;
580 x2 = Siscrat6;
581 x3 = Siscrat7;
582 if ( GcdLong(BHEAD Siscrat8,n1,Siscrat7,n2,x2,&n3) ) goto SimpErr;
583 n2 = n3;
584 if ( *x2 != 1 || n2 != 1 ) {
585 DivLong(a,*na,x2,n2,x1,&n1,x4,&n4);
586 *na = i = n1;
587 NCOPY(a,x1,i);
588 DivLong(b,*nb,x2,n2,x3,&n3,x4,&n4);
589 *nb = i = n3;
590 NCOPY(b,x3,i);
591 }
592 }
593#endif
594 else {
595 x4 = Siscrat5;
596 n1 = i = *na; x3 = a;
597 NCOPY(x1,x3,i);
598 x3 = b; n2 = i = *nb;
599 NCOPY(x2,x3,i);
600 x1 = Siscrat8; x2 = Siscrat7; x3 = Siscrat6;
601 for(;;){
602 if ( DivLong(x1,n1,x2,n2,x4,&n4,x3,&n3) ) goto SimpErr;
603 if ( !n3 ) break;
604 if ( n2 == 1 ) {
605 while ( ( *x1 = (*x2) % (*x3) ) != 0 ) { *x2 = *x3; *x3 = *x1; }
606 *x2 = *x3;
607 break;
608 }
609 if ( DivLong(x2,n2,x3,n3,x4,&n4,x1,&n1) ) goto SimpErr;
610 if ( !n1 ) { x2 = x3; n2 = n3; x3 = Siscrat7; break; }
611 if ( n3 == 1 ) {
612 while ( ( *x2 = (*x3) % (*x1) ) != 0 ) { *x3 = *x1; *x1 = *x2; }
613 *x2 = *x1;
614 n2 = 1;
615 break;
616 }
617 if ( DivLong(x3,n3,x1,n1,x4,&n4,x2,&n2) ) goto SimpErr;
618 if ( !n2 ) { x2 = x1; n2 = n1; x1 = Siscrat7; break; }
619 if ( n1 == 1 ) {
620 while ( ( *x3 = (*x1) % (*x2) ) != 0 ) { *x1 = *x2; *x2 = *x3; }
621 break;
622 }
623 }
624 if ( *x2 != 1 || n2 != 1 ) {
625 DivLong(a,*na,x2,n2,x1,&n1,x4,&n4);
626 *na = i = n1;
627 NCOPY(a,x1,i);
628 DivLong(b,*nb,x2,n2,x3,&n3,x4,&n4);
629 *nb = i = n3;
630 NCOPY(b,x3,i);
631 }
632 }
633 if ( sgn < 0 ) *na = -*na;
634 NumberFree(Siscrat5,"Simplify"); NumberFree(Siscrat6,"Simplify");
635 NumberFree(Siscrat7,"Simplify"); NumberFree(Siscrat8,"Simplify");
636 return(0);
637SimpErr:
638 MLOCK(ErrorMessageLock);
639 MesCall("Simplify");
640 MUNLOCK(ErrorMessageLock);
641 NumberFree(Siscrat5,"Simplify"); NumberFree(Siscrat6,"Simplify");
642 NumberFree(Siscrat7,"Simplify"); NumberFree(Siscrat8,"Simplify");
643 SETERROR(-1)
644}
645
646/*
647 #] Simplify :
648 #[ AccumGCD : WORD AccumGCD(PHEAD a,na,b,nb)
649
650 Routine takes the rational GCD of the fractions in a and b and
651 replaces a by the GCD of the two.
652 The rational GCD is defined as the rational that consists of
653 the GCD of the numerators divided by the GCD of the denominators
654*/
655
656int AccumGCD(PHEAD UWORD *a, WORD *na, UWORD *b, WORD nb)
657{
658 GETBIDENTITY
659 WORD nna,nnb,numa,numb,dena,denb,numc,denc;
660 UWORD *GCDbuffer = NumberMalloc("AccumGCD");
661 int i;
662 nna = *na; if ( nna < 0 ) nna = -nna; nna = (nna-1)/2;
663 nnb = nb; if ( nnb < 0 ) nnb = -nnb; nnb = (nnb-1)/2;
664 UnPack(a,nna,&dena,&numa);
665 UnPack(b,nnb,&denb,&numb);
666 if ( GcdLong(BHEAD a,numa,b,numb,GCDbuffer,&numc) ) goto AccErr;
667 numa = numc;
668 for ( i = 0; i < numa; i++ ) a[i] = GCDbuffer[i];
669 if ( GcdLong(BHEAD a+nna,dena,b+nnb,denb,GCDbuffer,&denc) ) goto AccErr;
670 dena = denc;
671 for ( i = 0; i < dena; i++ ) a[i+nna] = GCDbuffer[i];
672 Pack(a,&numa,a+nna,dena);
673 *na = INCLENG(numa);
674 NumberFree(GCDbuffer,"AccumGCD");
675 return(0);
676AccErr:
677 MLOCK(ErrorMessageLock);
678 MesCall("AccumGCD");
679 MUNLOCK(ErrorMessageLock);
680 NumberFree(GCDbuffer,"AccumGCD");
681 SETERROR(-1)
682}
683
684/*
685 #] AccumGCD :
686 #[ TakeRatRoot:
687*/
688
689int TakeRatRoot(UWORD *a, WORD *n, WORD power)
690{
691 WORD numer,denom, nn;
692 if ( ( power & 1 ) == 0 && *n < 0 ) return(1);
693 if ( ABS(*n) == 1 && a[0] == 1 && a[1] == 1 ) return(0);
694 nn = ABS(*n);
695 UnPack(a,nn,&denom,&numer);
696 if ( TakeLongRoot(a+nn,&denom,power) ) return(1);
697 if ( TakeLongRoot(a,&numer,power) ) return(1);
698 Pack(a,&numer,a+nn,denom);
699 if ( *n < 0 ) *n = -numer;
700 else *n = numer;
701 return(0);
702}
703
704/*
705 #] TakeRatRoot:
706 #] RekenRational :
707 #[ RekenLong :
708 #[ AddLong : WORD AddLong(a,na,b,nb,c,nc)
709
710 Long addition. Uses addition and subtraction of positive numbers.
711
712*/
713
714int AddLong(UWORD *a, WORD na, UWORD *b, WORD nb, UWORD *c, WORD *nc)
715{
716 WORD sgn, res;
717 if ( na < 0 ) {
718 if ( nb < 0 ) {
719 if ( AddPLon(a,-na,b,-nb,c,nc) ) return(-1);
720 *nc = -*nc;
721 return(0);
722 }
723 else {
724 na = -na;
725 sgn = -1;
726 }
727 }
728 else {
729 if ( nb < 0 ) {
730 nb = -nb;
731 sgn = 1;
732 }
733 else { return( AddPLon(a,na,b,nb,c,nc) ); }
734 }
735 if ( ( res = BigLong(a,na,b,nb) ) > 0 ) {
736 SubPLon(a,na,b,nb,c,nc);
737 if ( sgn < 0 ) *nc = -*nc;
738 }
739 else if ( res < 0 ) {
740 SubPLon(b,nb,a,na,c,nc);
741 if ( sgn > 0 ) *nc = -*nc;
742 }
743 else {
744 *nc = 0;
745 *c = 0;
746 }
747 return(0);
748}
749
750/*
751 #] AddLong :
752 #[ AddPLon : WORD AddPLon(a,na,b,nb,c,nc)
753
754 Adds two long integers a and b and puts the result in c.
755 The length of a and b are na and nb. The length of c is returned in nc.
756 c can be a or b.
757*/
758
759int AddPLon(UWORD *a, WORD na, UWORD *b, WORD nb, UWORD *c, WORD *nc)
760{
761 UWORD carry = 0, e, nd = 0;
762 while ( na && nb ) {
763 e = *a;
764 *c = e + *b + carry;
765 if ( carry ) {
766 if ( e < *c ) carry = 0;
767 }
768 else {
769 if ( e > *c ) carry = 1;
770 }
771 a++; b++; c++; nd++; na--; nb--;
772 }
773 while ( na ) {
774 if ( carry ) {
775 *c = *a++ + carry;
776 if ( *c++ ) carry = 0;
777 }
778 else *c++ = *a++;
779 nd++; na--;
780 }
781 while ( nb ) {
782 if ( carry ) {
783 *c = *b++ + carry;
784 if ( *c++ ) carry = 0;
785 }
786 else *c++ = *b++;
787 nd++; nb--;
788 }
789 if ( carry ) {
790 nd++;
791 if ( nd > (UWORD)AM.MaxTal ) {
792 MLOCK(ErrorMessageLock);
793 MesPrint("Overflow in addition");
794 MUNLOCK(ErrorMessageLock);
795 return(-1);
796 }
797 *c++ = carry;
798 }
799 *nc = nd;
800 return(0);
801}
802
803/*
804 #] AddPLon :
805 #[ SubPLon : void SubPLon(a,na,b,nb,c,nc)
806
807 Subtracts b from a. Assumes that a > b. Result in c.
808 c can be a or b.
809
810*/
811
812void SubPLon(UWORD *a, WORD na, UWORD *b, WORD nb, UWORD *c, WORD *nc)
813{
814 UWORD borrow = 0, e, nd = 0;
815 while ( nb ) {
816 e = *a;
817 if ( borrow ) {
818 *c = e - *b - borrow;
819 if ( *c < e ) borrow = 0;
820 }
821 else {
822 *c = e - *b;
823 if ( *c > e ) borrow = 1;
824 }
825 a++; b++; c++; na--; nb--; nd++;
826 }
827 while ( na ) {
828 if ( borrow ) {
829 if ( *a ) { *c++ = *a++ - 1; borrow = 0; }
830 else { *c++ = (UWORD)(-1); a++; }
831 }
832 else *c++ = *a++;
833 na--; nd++;
834 }
835 while ( nd && !*--c ) { nd--; }
836 *nc = (WORD)nd;
837}
838
839/*
840 #] SubPLon :
841 #[ MulLong : WORD MulLong(a,na,b,nb,c,nc)
842
843 Does a Long multiplication. Assumes that WORD is half the size
844 of a LONG to work out the scheme! The number of operations is
845 the canonical na*nm multiplications.
846 c should not overlap with a or b.
847
848*/
849
850int MulLong(UWORD *a, WORD na, UWORD *b, WORD nb, UWORD *c, WORD *nc)
851{
852 WORD sgn = 1;
853 UWORD i, *ic, *ia;
854 RLONG t, bb;
855 if ( !na || !nb ) { *nc = 0; return(0); }
856 if ( na < 0 ) { na = -na; sgn = -sgn; }
857 if ( nb < 0 ) { nb = -nb; sgn = -sgn; }
858 *nc = i = na + nb;
859 if ( i > (UWORD)(AM.MaxTal+1) ) goto MulLov;
860 ic = c;
861/*
862 #[ GMP stuff :
863*/
864#ifdef WITHGMP
865 if (na > 3 && nb > 3) {
866/* mp_limb_t res; */
867 UWORD *to, *from;
868 int j;
869 GETIDENTITY
870 UWORD *DLscrat9 = NumberMalloc("MulLong"), *DLscratA = NumberMalloc("MulLong"), *DLscratB = NumberMalloc("MulLong");
871#if ( GMPSPREAD != 1 )
872 if ( na & 1 ) {
873 from = a; a = to = DLscrat9; j = na; NCOPY(to, from, j);
874 a[na++] = 0;
875 ++*nc;
876 } else
877#endif
878 if ( (LONG)a & (sizeof(mp_limb_t)-1) ) {
879 from = a; a = to = DLscrat9; j = na; NCOPY(to, from, j);
880 }
881
882#if ( GMPSPREAD != 1 )
883 if ( nb & 1 ) {
884 from = b; b = to = DLscratA; j = nb; NCOPY(to, from, j);
885 b[nb++] = 0;
886 ++*nc;
887 } else
888#endif
889 if ( (LONG)b & (sizeof(mp_limb_t)-1) ) {
890 from = b; b = to = DLscratA; j = nb; NCOPY(to, from, j);
891 }
892
893 if ( ( *nc > (WORD)i ) || ( (LONG)c & (LONG)(sizeof(mp_limb_t)-1) ) ) {
894 ic = DLscratB;
895 }
896 if ( na < nb ) {
897 /* res = */
898 mpn_mul((mp_ptr)ic, (mp_srcptr)b, nb/GMPSPREAD, (mp_srcptr)a, na/GMPSPREAD);
899 } else {
900 /* res = */
901 mpn_mul((mp_ptr)ic, (mp_srcptr)a, na/GMPSPREAD, (mp_srcptr)b, nb/GMPSPREAD);
902 }
903 while ( ic[i-1] == 0 ) i--;
904 *nc = i;
905/*
906 if ( res == 0 ) *nc -= GMPSPREAD;
907 else if ( res <= WORDMASK ) --*nc;
908*/
909 if ( ic != c ) {
910 j = *nc; NCOPY(c, ic, j);
911 }
912 if ( sgn < 0 ) *nc = -(*nc);
913 NumberFree(DLscrat9,"MulLong"); NumberFree(DLscratA,"MulLong"); NumberFree(DLscratB,"MulLong");
914 return(0);
915 }
916#endif
917/*
918 #] GMP stuff :
919*/
920 do { *ic++ = 0; } while ( --i > 0 );
921 do {
922 ia = a;
923 ic = c++;
924 t = 0;
925 i = na;
926 bb = (RLONG)(*b++);
927 do {
928 t = (*ia++) * bb + t + *ic;
929 *ic++ = (WORD)t;
930 t >>= BITSINWORD; /* should actually be a swap */
931 } while ( --i > 0 );
932 if ( t ) *ic = (UWORD)t;
933 } while ( --nb > 0 );
934 if ( !*ic ) (*nc)--;
935 if ( *nc > AM.MaxTal ) goto MulLov;
936 if ( sgn < 0 ) *nc = -(*nc);
937 return(0);
938MulLov:
939 MLOCK(ErrorMessageLock);
940 MesPrint("Overflow in Multiplication");
941 MUNLOCK(ErrorMessageLock);
942 return(-1);
943}
944
945/*
946 #] MulLong :
947 #[ BigLong : WORD BigLong(a,na,b,nb)
948
949 Returns > 0 if a > b, < 0 if b > a and 0 if a == b
950
951*/
952
953int BigLong(UWORD *a, WORD na, UWORD *b, WORD nb)
954{
955 a += na;
956 b += nb;
957 while ( na && !*--a ) na--;
958 while ( nb && !*--b ) nb--;
959 if ( nb < na ) return(1);
960 if ( nb > na ) return(-1);
961 while ( --na >= 0 ) {
962 if ( *a > *b ) return(1);
963 else if ( *b > *a ) return(-1);
964 a--; b--;
965 }
966 return(0);
967}
968
969/*
970 #] BigLong :
971 #[ DivLong : WORD DivLong(a,na,b,nb,c,nc,d,nd)
972
973 This is the long division which knows a couple of exceptions.
974 It uses therefore a recursive call for the renormalization.
975 The quotient comes in c and the remainder in d.
976 d may be overlapping with b. It may also be identical to a.
977 c should not overlap with a, but it can overlap with b.
978
979*/
980
981int DivLong(UWORD *a, WORD na, UWORD *b, WORD nb, UWORD *c,
982 WORD *nc, UWORD *d, WORD *nd)
983{
984 WORD sgna = 1, sgnb = 1, ne, nf, ng, nh;
985 WORD i, ni;
986 UWORD *w1, *w2;
987 RLONG t, v;
988 UWORD *e, *f, *ff, *g, norm, estim;
989#ifdef WITHGMP
990 UWORD *DLscrat9, *DLscratA, *DLscratB, *DLscratC;
991#endif
992 RLONG esthelp;
993 if ( !nb ) {
994 MLOCK(ErrorMessageLock);
995 MesPrint("Division by zero");
996 MUNLOCK(ErrorMessageLock);
997 return(-1);
998 }
999 if ( !na ) { *nc = *nd = 0; return(0); }
1000 if ( na < 0 ) { sgna = -sgna; na = -na; }
1001 if ( nb < 0 ) { sgnb = -sgnb; nb = -nb; }
1002 if ( na < nb ) {
1003 for ( i = 0; i < na; i++ ) *d++ = *a++;
1004 *nd = na;
1005 *nc = 0;
1006 }
1007 else if ( nb == na && ( i = BigLong(b,nb,a,na) ) >= 0 ) {
1008 if ( i > 0 ) {
1009 for ( i = 0; i < na; i++ ) *d++ = *a++;
1010 *nd = na;
1011 *nc = 0;
1012 }
1013 else {
1014 *c = 1;
1015 *nc = 1;
1016 *nd = 0;
1017 }
1018 }
1019 else if ( nb == 1 ) {
1020 if ( *b == 1 ) {
1021 for ( i = 0; i < na; i++ ) *c++ = *a++;
1022 *nc = na;
1023 *nd = 0;
1024 }
1025 else {
1026 w1 = a+na;
1027 *nc = ni = na;
1028 *nd = 1;
1029 w2 = c+ni;
1030 v = (RLONG)(*b);
1031 t = (RLONG)(*--w1);
1032 while ( --ni >= 0 ) {
1033 *--w2 = t / v;
1034 t -= v * (*w2);
1035 if ( ni ) {
1036 t <<= BITSINWORD;
1037 t += *--w1;
1038 }
1039 }
1040 if ( ( *d = (UWORD)t ) == 0 ) *nd = 0;
1041 if ( !*(c+na-1) ) (*nc)--;
1042 }
1043 }
1044 else {
1045 GETIDENTITY
1046/*
1047 #[ GMP stuff :
1048
1049 We start with copying a and b.
1050 Then we make space for c and d.
1051 Next we call mpn_tdiv_qr
1052 We adjust sizes and copy to c and d if needed.
1053 Finally the signs are settled.
1054*/
1055#ifdef WITHGMP
1056 if ( na > 4 && nb > 3 ) {
1057 UWORD *ic, *id, *to, *from;
1058 int j = na - nb;
1059 DLscrat9 = NumberMalloc("DivLong"); DLscratA = NumberMalloc("DivLong");
1060 DLscratB = NumberMalloc("DivLong"); DLscratC = NumberMalloc("DivLong");
1061
1062#if ( GMPSPREAD != 1 )
1063 if ( na & 1 ) {
1064 from = a; a = to = DLscrat9; i = na; NCOPY(to, from, i);
1065 a[na++] = 0;
1066 } else
1067#endif
1068 if ( (LONG)a & (sizeof(mp_limb_t)-1) ) {
1069 from = a; a = to = DLscrat9; i = na; NCOPY(to, from, i);
1070 }
1071
1072#if ( GMPSPREAD != 1 )
1073 if ( nb & 1 ) {
1074 from = b; b = to = DLscratA; i = nb; NCOPY(to, from, i);
1075 b[nb++] = 0;
1076 } else
1077#endif
1078 if ( ( (LONG)b & (sizeof(mp_limb_t)-1) ) != 0 ) {
1079 from = b; b = to = DLscratA; i = nb; NCOPY(to, from, i);
1080 }
1081#if ( GMPSPREAD != 1 )
1082/*
1083 Not recognizing this case was a nasty and extremely rare bug.
1084 In the case that c == a and na is odd and nc == na and b
1085 still needs to be used, the least significant UWORD of b got
1086 overwritten by zero. (22-mar-2023)
1087*/
1088 ic = DLscratB; id = DLscratC;
1089#else
1090 if ( ( (LONG)c & (sizeof(mp_limb_t)-1) ) != 0 ) ic = DLscratB;
1091 else ic = c;
1092
1093 if ( ( (LONG)d & (sizeof(mp_limb_t)-1) ) != 0 ) id = DLscratC;
1094 else id = d;
1095#endif
1096 mpn_tdiv_qr((mp_limb_t *)ic,(mp_limb_t *)id,(mp_size_t)0,
1097 (const mp_limb_t *)a,(mp_size_t)(na/GMPSPREAD),
1098 (const mp_limb_t *)b,(mp_size_t)(nb/GMPSPREAD));
1099 while ( j >= 0 && ic[j] == 0 ) j--;
1100 j++; *nc = j;
1101 if ( c != ic ) { NCOPY(c,ic,j); }
1102 j = nb-1;
1103 while ( j >= 0 && id[j] == 0 ) j--;
1104 j++; *nd = j;
1105 if ( d != id ) { NCOPY(d,id,j); }
1106 if ( sgna < 0 ) { *nc = -(*nc); *nd = -(*nd); }
1107 if ( sgnb < 0 ) { *nc = -(*nc); }
1108 NumberFree(DLscrat9,"DivLong"); NumberFree(DLscratA,"DivLong");
1109 NumberFree(DLscratB,"DivLong"); NumberFree(DLscratC,"DivLong");
1110 return(0);
1111 }
1112#endif
1113/*
1114 #] GMP stuff :
1115*/
1116 /* Start with normalization operation */
1117
1118 e = NumberMalloc("DivLong"); f = NumberMalloc("DivLong"); g = NumberMalloc("DivLong");
1119 if ( b[nb-1] == (FULLMAX-1) ) norm = 1;
1120 else {
1121 norm = (UWORD)(((ULONG)FULLMAX) / (ULONG)((b[nb-1]+1L)));
1122 }
1123 f[na] = 0;
1124 if ( MulLong(b,nb,&norm,1,e,&ne) ||
1125 MulLong(a,na,&norm,1,f,&nf) ) {
1126 NumberFree(e,"DivLong"); NumberFree(f,"DivLong"); NumberFree(g,"DivLong");
1127 return(-1);
1128 }
1129 if ( BigLong(f+nf-ne,ne,e,ne) >= 0 ) {
1130 SubPLon(f+nf-ne,ne,e,ne,f+nf-ne,&nh);
1131 w1 = c + (nf-ne);
1132 *nc = nf-ne+1;
1133 }
1134 else {
1135 nh = ne;
1136 *nc = nf-ne;
1137 w1 = 0;
1138 }
1139 w2 = c; i = *nc; do { *w2++ = 0; } while ( --i > 0 );
1140 nf = na;
1141 ni = nf-ne;
1142 esthelp = (RLONG)(e[ne-1]) + 1L;
1143 while ( nf >= ne ) {
1144 if ( (WORD)esthelp == 0 ) {
1145 estim = (WORD)(((((RLONG)(f[nf]))<<BITSINWORD)+f[nf-1])>>BITSINWORD);
1146 }
1147 else {
1148 estim = (WORD)(((((RLONG)(f[nf]))<<BITSINWORD)+f[nf-1])/esthelp);
1149 }
1150 /* This estimate may be up to two too small */
1151 if ( estim ) {
1152 MulLong(e,ne,&estim,1,g,&ng);
1153 nh = ne + 1; if ( !f[ni+ne] ) nh--;
1154 SubPLon(f+ni,nh,g,ng,f+ni,&nh);
1155 }
1156 else {
1157 w2 = f+ni+ne; nh = ne+1;
1158 while ( ( nh > 0 ) && !*w2 ) { nh--; w2--; }
1159 }
1160 if ( BigLong(f+ni,nh,e,ne) >= 0 ) {
1161 estim++;
1162 SubPLon(f+ni,nh,e,ne,f+ni,&nh);
1163 if ( BigLong(f+ni,nh,e,ne) >= 0 ) {
1164 estim++;
1165 SubPLon(f+ni,nh,e,ne,f+ni,&nh);
1166 if ( BigLong(f+ni,nh,e,ne) >= 0 ) {
1167/* INTERNAL_ERROR_EXCL_START */
1168 MLOCK(ErrorMessageLock);
1169 MesPrint("!>Problems in DivLong");
1170 AO.OutSkip = 3;
1171 FiniLine();
1172 i = na;
1173 while ( --i >= 0 ) { TalToLine((UWORD)(*a++)); TokenToLine((UBYTE *)" "); }
1174 FiniLine();
1175 i = nb;
1176 while ( --i >= 0 ) { TalToLine((UWORD)(*b++)); TokenToLine((UBYTE *)" "); }
1177 AO.OutSkip = 0;
1178 FiniLine();
1179 MUNLOCK(ErrorMessageLock);
1180 NumberFree(e,"DivLong"); NumberFree(f,"DivLong"); NumberFree(g,"DivLong");
1181 return(-1);
1182/* INTERNAL_ERROR_EXCL_STOP */
1183 }
1184 }
1185 }
1186 c[ni] = estim;
1187 nf--;
1188 ni--;
1189 }
1190 if ( w1 ) *w1 = 1;
1191
1192 /* Finish with the renormalization operation */
1193
1194 if ( nh > 0 ) {
1195 if ( norm == 1 ) {
1196 *nd = i = nh; ff = f;
1197 NCOPY(d,ff,i);
1198 }
1199 else {
1200 w1 = f+nh;
1201 *nd = ni = nh;
1202 w2 = d+ni;
1203 v = norm;
1204 t = (RLONG)(*--w1);
1205 while ( --ni >= 0 ) {
1206 *--w2 = t / v;
1207 t -= v * (*w2);
1208 if ( ni ) {
1209 t <<= BITSINWORD;
1210 t += *--w1;
1211 }
1212 }
1213 if ( t ) {
1214/* INTERNAL_ERROR_EXCL_START */
1215 MLOCK(ErrorMessageLock);
1216 MesPrint("!>Error in DivLong");
1217 MUNLOCK(ErrorMessageLock);
1218 NumberFree(e,"DivLong"); NumberFree(f,"DivLong"); NumberFree(g,"DivLong");
1219 return(-1);
1220/* INTERNAL_ERROR_EXCL_STOP */
1221 }
1222 if ( !*(d+nh-1) ) (*nd)--;
1223 }
1224 }
1225 else { *nd = 0; }
1226 NumberFree(e,"DivLong"); NumberFree(f,"DivLong"); NumberFree(g,"DivLong");
1227 }
1228 if ( sgna < 0 ) { *nc = -(*nc); *nd = -(*nd); }
1229 if ( sgnb < 0 ) { *nc = -(*nc); }
1230 return(0);
1231}
1232
1233/*
1234 #] DivLong :
1235 #[ RaisPow : WORD RaisPow(a,na,b)
1236
1237 Raises a to the power b. a is a Long integer and b >= 0.
1238 The method that is used works with a bitdecomposition of b.
1239*/
1240
1241int RaisPow(PHEAD UWORD *a, WORD *na, UWORD b)
1242{
1243 GETBIDENTITY
1244 WORD i, nu;
1245 UWORD *it, *iu, c;
1246 UWORD *is, *iss;
1247 WORD ns, nt, nmod;
1248 nmod = ABS(AN.ncmod);
1249 if ( !*na || ( ( *na == 1 ) && ( *a == 1 ) ) ) return(0);
1250 if ( !b ) { *na=1; *a=1; return(0); }
1251 is = NumberMalloc("RaisPow");
1252 it = NumberMalloc("RaisPow");
1253 for ( i = 0; i < ABS(*na); i++ ) is[i] = a[i];
1254 ns = *na;
1255 c = b;
1256 for ( i = 0; i < BITSINWORD; i++ ) {
1257 if ( !c ) break;
1258 c /= 2;
1259 }
1260 i--;
1261 c = 1u << i;
1262 while ( --i >= 0 ) {
1263 c /= 2;
1264 if(MulLong(is,ns,is,ns,it,&nt)) goto RaisOvl;
1265 if ( b & c ) {
1266 if ( MulLong(it,nt,a,*na,is,&ns) ) goto RaisOvl;
1267 }
1268 else {
1269 iu = is; is = it; it = iu;
1270 nu = ns; ns = nt; nt = nu;
1271 }
1272 if ( nmod != 0 ) {
1273 if ( DivLong(is,ns,(UWORD *)AC.cmod,nmod,it,&nt,is,&ns) ) goto RaisOvl;
1274 }
1275 }
1276 if ( ( nmod != 0 ) && ( ( AC.modmode & POSNEG ) != 0 ) ) {
1277 NormalModulus(is,&ns);
1278 }
1279 if ( ( *na = i = ns ) != 0 ) { iss = is; i=ABS(i); NCOPY(a,iss,i); }
1280 NumberFree(is,"RaisPow"); NumberFree(it,"RaisPow");
1281 return(0);
1282RaisOvl:
1283 MLOCK(ErrorMessageLock);
1284 MesCall("RaisPow");
1285 MUNLOCK(ErrorMessageLock);
1286 NumberFree(is,"RaisPow"); NumberFree(it,"RaisPow");
1287 SETERROR(-1)
1288}
1289
1290/*
1291 #] RaisPow :
1292 #[ RaisPowCached :
1293*/
1294
1309void RaisPowCached (PHEAD WORD x, WORD n, UWORD **c, WORD *nc) {
1310
1311 int i,j;
1312 WORD new_small_power_maxx, new_small_power_maxn, ID;
1313 WORD *new_small_power_n;
1314 UWORD **new_small_power;
1315
1316 /* check whether to extend the array */
1317 if (x>=AT.small_power_maxx || n>=AT.small_power_maxn) {
1318
1319 new_small_power_maxx = AT.small_power_maxx;
1320 if (x>=AT.small_power_maxx)
1321 new_small_power_maxx = MaX(2*AT.small_power_maxx, x+1);
1322
1323 new_small_power_maxn = AT.small_power_maxn;
1324 if (n>=AT.small_power_maxn)
1325 new_small_power_maxn = MaX(2*AT.small_power_maxn, n+1);
1326
1327 new_small_power_n = (WORD*) Malloc1(new_small_power_maxx*new_small_power_maxn*sizeof(WORD),"RaisPowCached");
1328 new_small_power = (UWORD **) Malloc1(new_small_power_maxx*new_small_power_maxn*sizeof(UWORD *),"RaisPowCached");
1329
1330 for (i=0; i<new_small_power_maxx * new_small_power_maxn; i++) {
1331 new_small_power_n[i] = 0;
1332 new_small_power [i] = NULL;
1333 }
1334
1335 for (i=0; i<AT.small_power_maxx; i++)
1336 for (j=0; j<AT.small_power_maxn; j++) {
1337 new_small_power_n[i*new_small_power_maxn+j] = AT.small_power_n[i*AT.small_power_maxn+j];
1338 new_small_power [i*new_small_power_maxn+j] = AT.small_power [i*AT.small_power_maxn+j];
1339 }
1340
1341 if (AT.small_power_n != NULL) {
1342 M_free(AT.small_power_n,"RaisPowCached");
1343 M_free(AT.small_power,"RaisPowCached");
1344 }
1345
1346 AT.small_power_maxx = new_small_power_maxx;
1347 AT.small_power_maxn = new_small_power_maxn;
1348 AT.small_power_n = new_small_power_n;
1349 AT.small_power = new_small_power;
1350 }
1351
1352 /* check whether the results is already calculated */
1353 ID = x * AT.small_power_maxn + n;
1354
1355 if (AT.small_power[ID] == NULL) {
1356#ifdef OLDRAISPOWCACHED
1357 AT.small_power[ID] = NumberMalloc("RaisPowCached");
1358 AT.small_power_n[ID] = 1;
1359 AT.small_power[ID][0] = x;
1360 RaisPow(BHEAD AT.small_power[ID],&AT.small_power_n[ID],n);
1361#else
1362 UWORD *c = NumberMalloc("RaisPowCached");
1363 WORD i, k = 1;
1364 c[0] = x;
1365 RaisPow(BHEAD c,&k,n);
1366/*
1367 And now get the proper amount.
1368*/
1369 if ( AT.InNumMem < k ) { /* We should start a new buffer */
1370 AT.InNumMem = 5*AM.MaxTal;
1371 AT.NumMem = (UWORD *)Malloc1(AT.InNumMem*sizeof(UWORD),"RaisPowCached");
1372/*
1373 MesPrint(" Got an extra %l UWORDS in RaisPowCached",AT.InNumMem);
1374*/
1375 }
1376 for ( i = 0; i < k; i++ ) AT.NumMem[i] = c[i];
1377 AT.small_power[ID] = AT.NumMem;
1378 AT.small_power_n[ID] = k;
1379 AT.NumMem += k;
1380 AT.InNumMem -= k;
1381 NumberFree(c,"RaisPowCached");
1382#endif
1383 }
1384
1385 /* return the result */
1386 *c = AT.small_power[ID];
1387 *nc = AT.small_power_n[ID];
1388}
1389
1390/*
1391 #] RaisPowCached :
1392 #[ RaisPowMod :
1393
1394 Computes the power x^n mod m
1395 */
1396WORD RaisPowMod (WORD x, WORD n, WORD m) {
1397 LONG y=1, z=x;
1398 while (n) {
1399 if (n&1) { y*=z; y%=m; }
1400 z*=z; z%=m;
1401 n /= 2;
1402 }
1403 return (WORD)y;
1404}
1405
1406/*
1407 #] RaisPowMod :
1408 #[ NormalModulus : int NormalModulus(UWORD *a,WORD *na)
1409*/
1416int NormalModulus(UWORD *a,WORD *na)
1417{
1418 WORD n;
1419 if ( AC.halfmod == 0 ) {
1420 LOCK(AC.halfmodlock);
1421 if ( AC.halfmod == 0 ) {
1422 UWORD two[1],remain[1];
1423 WORD dummy;
1424 two[0] = 2;
1425 AC.halfmod = (UWORD *)Malloc1((ABS(AC.ncmod))*sizeof(UWORD),"halfmod");
1426 DivLong((UWORD *)AC.cmod,(ABS(AC.ncmod)),two,1
1427 ,(UWORD *)AC.halfmod,&(AC.nhalfmod),remain,&dummy);
1428 }
1429 UNLOCK(AC.halfmodlock);
1430 }
1431 n = ABS(*na);
1432 if ( BigLong(a,n,AC.halfmod,AC.nhalfmod) > 0 ) {
1433 SubPLon((UWORD *)AC.cmod,(ABS(AC.ncmod)),a,n,a,&n);
1434 if ( *na > 0 ) { *na = -n; }
1435 else { *na = n; }
1436 return(1);
1437 }
1438 return(0);
1439}
1440
1441/*
1442 #] NormalModulus :
1443 #[ MakeInverses :
1444*/
1453int MakeInverses(void)
1454{
1455 WORD n = AC.cmod[0], i, inv2;
1456 if ( AC.ncmod != 1 ) return(1);
1457 if ( AC.modinverses == 0 ) {
1458 LOCK(AC.halfmodlock);
1459 if ( AC.modinverses == 0 ) {
1460 AC.modinverses = (UWORD *)Malloc1(n*sizeof(UWORD),"modinverses");
1461 AC.modinverses[0] = 0;
1462 AC.modinverses[1] = 1;
1463 for ( i = 2; i < n; i++ ) {
1464 if ( GetModInverses(i,n,
1465 (WORD *)(&(AC.modinverses[i])),&inv2) ) {
1466 SETERROR(-1)
1467 }
1468 }
1469 }
1470 UNLOCK(AC.halfmodlock);
1471 }
1472 return(0);
1473}
1474
1475/*
1476 #] MakeInverses :
1477 #[ GetModInverses :
1478*/
1489int GetModInverses(WORD m1, WORD m2, WORD *im1, WORD *im2)
1490{
1491 WORD a1, a2, a3;
1492 WORD b1, b2, b3;
1493 WORD x = m1, y, c, d = m2;
1494 if ( x < 1 || d <= 1 ) goto somethingwrong;
1495 a1 = 0; a2 = 1;
1496 b1 = 1; b2 = 0;
1497 for(;;) {
1498 c = d/x; y = d%x; /* a good compiler makes this faster than y=d-c*x */
1499 if ( y == 0 ) break;
1500 a3 = a1-c*a2; a1 = a2; a2 = a3;
1501 b3 = b1-c*b2; b1 = b2; b2 = b3;
1502 d = x; x = y;
1503 }
1504 if ( x != 1 ) goto somethingwrong;
1505 if ( a2 < 0 ) a2 += m2;
1506 if ( b2 < 0 ) b2 += m1;
1507 if (im1!=NULL) *im1 = a2;
1508 if (im2!=NULL) *im2 = b2;
1509 return(0);
1510somethingwrong:
1511/* INTERNAL_ERROR_EXCL_START */
1512 MLOCK(ErrorMessageLock);
1513 MesPrint("!>Error trying to determine inverses in GetModInverses");
1514 MUNLOCK(ErrorMessageLock);
1515 return(-1);
1516/* INTERNAL_ERROR_EXCL_STOP */
1517}
1518/*
1519 #] GetModInverses :
1520 #[ GetLongModInverses :
1521*/
1522
1523int GetLongModInverses(PHEAD UWORD *a, WORD na, UWORD *b, WORD nb, UWORD *ia, WORD *nia, UWORD *ib, WORD *nib) {
1524
1525 UWORD *s, *t, *sa, *sb, *ta, *tb, *x, *y, *swap1;
1526 WORD ns, nt, nsa, nsb, nta, ntb, nx, ny, swap2;
1527
1528 s = NumberMalloc("GetLongModInverses");
1529 ns = na;
1530 WCOPY(s, a, ABS(ns));
1531
1532 t = NumberMalloc("GetLongModInverses");
1533 nt = nb;
1534 WCOPY(t, b, ABS(nt));
1535
1536 sa = NumberMalloc("GetLongModInverses");
1537 nsa = 1;
1538 sa[0] = 1;
1539
1540 sb = NumberMalloc("GetLongModInverses");
1541 nsb = 0;
1542
1543 ta = NumberMalloc("GetLongModInverses");
1544 nta = 0;
1545
1546 tb = NumberMalloc("GetLongModInverses");
1547 ntb = 1;
1548 tb[0] = 1;
1549
1550 x = NumberMalloc("GetLongModInverses");
1551 y = NumberMalloc("GetLongModInverses");
1552
1553 while (nt != 0) {
1554 DivLong(s,ns,t,nt,x,&nx,y,&ny);
1555 swap1=s; s=y; y=swap1;
1556 ns=ny;
1557 MulLong(x,nx,ta,nta,y,&ny);
1558 AddLong(sa,nsa,y,-ny,sa,&nsa);
1559 MulLong(x,nx,tb,ntb,y,&ny);
1560 AddLong(sb,nsb,y,-ny,sb,&nsb);
1561
1562 swap1=s; s=t; t=swap1;
1563 swap2=ns; ns=nt; nt=swap2;
1564 swap1=sa; sa=ta; ta=swap1;
1565 swap2=nsa; nsa=nta; nta=swap2;
1566 swap1=sb; sb=tb; tb=swap1;
1567 swap2=nsb; nsb=ntb; ntb=swap2;
1568 }
1569
1570 if (ia!=NULL) {
1571 *nia = nsa*ns;
1572 WCOPY(ia,sa,ABS(*nia));
1573 }
1574
1575 if (ib!=NULL) {
1576 *nib = nsb*ns;
1577 WCOPY(ib,sb,ABS(*nib));
1578 }
1579
1580 NumberFree(s,"GetLongModInverses");
1581 NumberFree(t,"GetLongModInverses");
1582 NumberFree(sa,"GetLongModInverses");
1583 NumberFree(sb,"GetLongModInverses");
1584 NumberFree(ta,"GetLongModInverses");
1585 NumberFree(tb,"GetLongModInverses");
1586 NumberFree(x,"GetLongModInverses");
1587 NumberFree(y,"GetLongModInverses");
1588
1589 return 0;
1590}
1591
1592/*
1593 #] GetLongModInverses :
1594 #[ Product : WORD Product(a,na,b)
1595
1596 Multiplies the Long number in a with the WORD b.
1597
1598*/
1599
1600int Product(UWORD *a, WORD *na, WORD b)
1601{
1602 WORD i, sgn = 1;
1603 RLONG t, u;
1604 if ( *na < 0 ) { *na = -(*na); sgn = -sgn; }
1605 if ( b < 0 ) { b = -b; sgn = -sgn; }
1606 t = 0;
1607 u = (RLONG)b;
1608 for ( i = 0; i < *na; i++ ) {
1609 t += *a * u;
1610 *a++ = (UWORD)t;
1611 t >>= BITSINWORD;
1612 }
1613 if ( t > 0 ) {
1614 if ( ++(*na) > AM.MaxTal ) {
1615 MLOCK(ErrorMessageLock);
1616 MesPrint("Overflow in Product");
1617 MUNLOCK(ErrorMessageLock);
1618 return(-1);
1619 }
1620 *a = (UWORD)t;
1621 }
1622 if ( sgn < 0 ) *na = -(*na);
1623 return(0);
1624}
1625
1626/*
1627 #] Product :
1628 #[ Quotient : UWORD Quotient(a,na,b)
1629
1630 Routine divides the long number a by b with the assumption that
1631 there is no remainder (like while computing binomials).
1632
1633*/
1634
1635UWORD Quotient(UWORD *a, WORD *na, WORD b)
1636{
1637 RLONG v, t;
1638 WORD i, j, sgn = 1;
1639 if ( ( i = *na ) < 0 ) { sgn = -1; i = -i; }
1640 if ( b < 0 ) { b = -b; sgn = -sgn; }
1641 if ( i == 1 ) {
1642 if ( ( *a /= (UWORD)b ) == 0 ) *na = 0;
1643 if ( sgn < 0 ) *na = -*na;
1644 return(0);
1645 }
1646 a += i;
1647 j = i;
1648 v = (RLONG)b;
1649 t = (RLONG)(*--a);
1650 while ( --i >= 0 ) {
1651 *a = t / v;
1652 t -= v * (*a);
1653 if ( i ) {
1654 t <<= BITSINWORD;
1655 t += *--a;
1656 }
1657 }
1658 a += j - 1;
1659 if ( !*a ) j--;
1660 if ( sgn < 0 ) j = -j;
1661 *na = j;
1662 return(0);
1663}
1664
1665/*
1666 #] Quotient :
1667 #[ Remain10 : WORD Remain10(a,na)
1668
1669 Routine divides a by 10 and gives the remainder as return value.
1670 The value of a will be the quotient! a must be positive.
1671
1672*/
1673
1674WORD Remain10(UWORD *a, WORD *na)
1675{
1676 WORD i;
1677 RLONG t, u;
1678 UWORD *b;
1679 i = *na;
1680 t = 0;
1681 b = a + i - 1;
1682 while ( --i >= 0 ) {
1683 t += *b;
1684 *b-- = u = t / 10;
1685 t -= u * 10;
1686 if ( i > 0 ) t <<= BITSINWORD;
1687 }
1688 if ( ( *na > 0 ) && !a[*na-1] ) (*na)--;
1689 return((WORD)t);
1690}
1691
1692/*
1693 #] Remain10 :
1694 #[ Remain4 : WORD Remain4(a,na)
1695
1696 Routine divides a by 10000 and gives the remainder as return value.
1697 The value of a will be the quotient! a must be positive.
1698
1699*/
1700
1701WORD Remain4(UWORD *a, WORD *na)
1702{
1703 WORD i;
1704 RLONG t, u;
1705 UWORD *b;
1706 i = *na;
1707 t = 0;
1708 b = a + i - 1;
1709 while ( --i >= 0 ) {
1710 t += *b;
1711 *b-- = u = t / 10000;
1712 t -= u * 10000;
1713 if ( i > 0 ) t <<= BITSINWORD;
1714 }
1715 if ( ( *na > 0 ) && !a[*na-1] ) (*na)--;
1716 return((WORD)t);
1717}
1718
1719/*
1720 #] Remain4 :
1721 #[ PrtLong : void PrtLong(a,na,s)
1722
1723 Puts the long number a in string s.
1724
1725*/
1726
1727void PrtLong(UWORD *a, WORD na, UBYTE *s)
1728{
1729 GETIDENTITY
1730 WORD q, i;
1731 UBYTE *sa, *sb;
1732 UBYTE c;
1733 UWORD *bb, *b;
1734
1735 if ( na < 0 ) {
1736 *s++ = '-';
1737 na = -na;
1738 }
1739
1740 b = NumberMalloc("PrtLong");
1741 bb = b;
1742 i = na; while ( --i >= 0 ) *bb++ = *a++;
1743 a = b;
1744 if ( na > 2 ) {
1745 sa = s;
1746 do {
1747 q = Remain4(a,&na);
1748 *sa++ = (UBYTE)('0' + (q%10));
1749 q /= 10;
1750 *sa++ = (UBYTE)('0' + (q%10));
1751 q /= 10;
1752 *sa++ = (UBYTE)('0' + (q%10));
1753 q /= 10;
1754 *sa++ = (UBYTE)('0' + (q%10));
1755 } while ( na );
1756 while ( sa[-1] == '0' ) sa--;
1757 sb = s;
1758 s = sa;
1759 sa--;
1760 while ( sa > sb ) { c = *sa; *sa = *sb; *sb = c; sa--; sb++; }
1761 }
1762 else if ( na ) {
1763 sa = s;
1764 do {
1765 q = Remain10(a,&na);
1766 *sa++ = (UBYTE)('0' + q);
1767 } while ( na );
1768 sb = s;
1769 s = sa;
1770 sa--;
1771 while ( sa > sb ) { c = *sa; *sa = *sb; *sb = c; sa--; sb++; }
1772 }
1773 else *s++ = '0';
1774 *s = '\0';
1775 NumberFree(b,"PrtLong");
1776}
1777
1778/*
1779 #] PrtLong :
1780 #[ GetLong : WORD GetLong(s,a,na)
1781
1782 Reads a long number from a string.
1783 The string is zero terminated and contains only digits!
1784
1785 New algorithm: try to read 4 digits together before the result
1786 is accumulated.
1787*/
1788
1789int GetLong(UBYTE *s, UWORD *a, WORD *na)
1790{
1791/*
1792 UWORD digit;
1793 *a = 0;
1794 *na = 0;
1795 while ( FG.cTable[*s] == 1 ) {
1796 digit = *s++ - '0';
1797 if ( *na && Product(a,na,(WORD)10) ) return(-1);
1798 if ( digit && AddLong(a,*na,&digit,(WORD)1,a,na) ) return(-1);
1799 }
1800 return(0);
1801*/
1802 UWORD digit, x = 0, y = 0;
1803 *a = 0;
1804 *na = 0;
1805 while ( FG.cTable[*s] == 1 ) {
1806 x = *s++ - '0';
1807 if ( FG.cTable[*s] != 1 ) { y = 10; break; }
1808 x = 10*x + *s++ - '0';
1809 if ( FG.cTable[*s] != 1 ) { y = 100; break; }
1810 x = 10*x + *s++ - '0';
1811 if ( FG.cTable[*s] != 1 ) { y = 1000; break; }
1812 x = 10*x + *s++ - '0';
1813 if ( *na && Product(a,na,(WORD)10000) ) return(-1);
1814 if ( ( digit = x ) != 0 && AddLong(a,*na,&digit,(WORD)1,a,na) )
1815 return(-1);
1816 y = 0;
1817 }
1818 if ( y ) {
1819 if ( *na && Product(a,na,(WORD)y) ) return(-1);
1820 if ( ( digit = x ) != 0 && AddLong(a,*na,&digit,(WORD)1,a,na) )
1821 return(-1);
1822 }
1823 return(0);
1824}
1825
1826/*
1827 #] GetLong :
1828 #[ GCD : WORD GCD(a,na,b,nb,c,nc)
1829
1830 Algorithm to compute the GCD of two long numbers.
1831 See Knuth, sec 4.5.2 algorithm L.
1832
1833 We assume that both numbers are positive
1834
1835 NOTE!!!!!. NumberMalloc gets called and it may not be freed
1836*/
1837
1838#ifdef EXTRAGCD
1839
1840#define Convert(ia,aa,naa) \
1841 if ( (LONG)ia < 0 ) { \
1842 ia = (ULONG)(-(LONG)ia); \
1843 aa[0] = ia; \
1844 if ( ( aa[1] = ia >> BITSINWORD ) != 0 ) naa = -2; \
1845 else naa = -1; \
1846 } \
1847 else if ( ia == 0 ) { aa[0] = 0; naa = 0; } \
1848 else { \
1849 aa[0] = ia; \
1850 if ( ( aa[1] = ia >> BITSINWORD ) != 0 ) naa = 2; \
1851 else naa = 1; \
1852 }
1853
1854void GCD(UWORD *a, WORD na, UWORD *b, WORD nb, UWORD *c, WORD *nc)
1855{
1856 int ja = 0, jb = 0, j;
1857 UWORD *r,*t;
1858 UWORD *x1, *x2, *x3;
1859 WORD nd,naa,nbb;
1860 ULONG ia,ib,ic,id,u,v,w,q,T;
1861 UWORD aa[2], bb[2];
1862/*
1863 First eliminate easy powers of 2^...
1864*/
1865 while ( a[0] == 0 ) { na--; ja++; a++; }
1866 while ( b[0] == 0 ) { nb--; jb++; b++; }
1867 if ( ja > jb ) ja = jb;
1868 if ( ja > 0 ) {
1869 j = ja;
1870 do { *c++ = 0; } while ( --j > 0 );
1871 }
1872/*
1873 Now arrange things such that a >= b
1874*/
1875 if ( na < nb ) {
1876 jb = na; na = nb; nb = jb;
1877exch:
1878 r = a; a = b; b = r;
1879 }
1880 else if ( na == nb ) {
1881 r = a+na;
1882 t = b+nb;
1883 j = na;
1884 while ( --j >= 0 ) {
1885 if ( *--r > *--t ) break;
1886 if ( *r < *t ) goto exch;
1887 }
1888 if ( j < 0 ) {
1889out:
1890 j = nb;
1891 NCOPY(c,b,j);
1892 *nc = nb+ja;
1893 return;
1894 }
1895 }
1896/*
1897 {
1898 MLOCK(ErrorMessageLock);
1899 MesPrint("Ordered input, ja = %d",(WORD)ja);
1900 AO.OutSkip = 3;
1901 FiniLine();
1902 j = na; r = a;
1903 while ( --j >= 0 ) { TalToLine((UWORD)(*r++)); TokenToLine((UBYTE *)" "); }
1904 FiniLine();
1905 j = nb; r = b;
1906 while ( --j >= 0 ) { TalToLine((UWORD)(*r++)); TokenToLine((UBYTE *)" "); }
1907 AO.OutSkip = 0;
1908 FiniLine();
1909 MUNLOCK(ErrorMessageLock);
1910 }
1911*/
1912/*
1913 We have now that A > B
1914 The loop recognizes the case that na-nb >= 1
1915 In that case we just have to divide!
1916*/
1917 r = x1 = NumberMalloc("GCD"); t = x2 = NumberMalloc("GCD"); x3 = NumberMalloc("GCD");
1918 j = na;
1919 NCOPY(r,a,j);
1920 j = nb;
1921 NCOPY(t,b,j);
1922
1923 for(;;) {
1924
1925 while ( na > nb ) {
1926toobad:
1927 DivLong(x1,na,x2,nb,c,nc,x3,&nd);
1928 if ( nd == 0 ) { b = x2; goto out; }
1929 t = x1; x1 = x2; x2 = x3; x3 = t; na = nb; nb = nd;
1930 if ( na == 2 ) break;
1931 }
1932/*
1933 Here we can use the shortcut.
1934*/
1935 if ( na == 2 ) {
1936 v = x1[0] + ( ((ULONG)x1[1]) << BITSINWORD );
1937 w = x2[0];
1938 if ( nb == 2 ) w += ((ULONG)x2[1]) << BITSINWORD;
1939#ifdef EXTRAGCD2
1940 v = GCD2(v,w);
1941#else
1942 do { u = v%w; v = w; w = u; } while ( w );
1943#endif
1944 c[0] = (UWORD)v;
1945 if ( ( c[1] = (UWORD)(v >> BITSINWORD) ) != 0 ) *nc = 2+ja;
1946 else *nc = 1+ja;
1947 NumberFree(x1,"GCD"); NumberFree(x2,"GCD"); NumberFree(x3,"GCD");
1948 return;
1949 }
1950 if ( na == 1 ) {
1951 UWORD ui, uj;
1952 ui = x1[0]; uj = x2[0];
1953#ifdef EXTRAGCD2
1954 ui = (UWORD)GCD2((ULONG)ui,(ULONG)uj);
1955#else
1956 do { nd = ui%uj; ui = uj; uj = nd; } while ( nd );
1957#endif
1958 c[0] = ui;
1959 *nc = 1 + ja;
1960 NumberFree(x1,"GCD"); NumberFree(x2,"GCD"); NumberFree(x3,"GCD");
1961 return;
1962 }
1963 ia = 1; ib = 0; ic = 0; id = 1;
1964 u = ( ((ULONG)x1[na-1]) << BITSINWORD ) + x1[na-2];
1965 v = ( ((ULONG)x2[nb-1]) << BITSINWORD ) + x2[nb-2];
1966
1967 while ( v+ic != 0 && v+id != 0 &&
1968 ( q = (u+ia)/(v+ic) ) == (u+ib)/(v+id) ) {
1969 T = ia-q*ic; ia = ic; ic = T;
1970 T = ib-q*id; ib = id; id = T;
1971 T = u - q*v; u = v; v = T;
1972 }
1973 if ( ib == 0 ) goto toobad;
1974 Convert(ia,aa,naa);
1975 Convert(ib,bb,nbb);
1976 MulLong(x1,na,aa,naa,x3,&nd);
1977 MulLong(x2,nb,bb,nbb,c,nc);
1978 AddLong(x3,nd,c,*nc,c,nc);
1979 Convert(ic,aa,naa);
1980 Convert(id,bb,nbb);
1981 MulLong(x1,na,aa,naa,x3,&nd);
1982 t = c; na = j = *nc; r = x1;
1983 NCOPY(r,t,j);
1984 MulLong(x2,nb,bb,nbb,c,nc);
1985 AddLong(x3,nd,c,*nc,x2,&nb);
1986 }
1987}
1988
1989#endif
1990
1991/*
1992 #] GCD :
1993 #[ GcdLong : WORD GcdLong(a,na,b,nb,c,nc)
1994
1995 Returns the Greatest Common Divider of a and b in c.
1996 If a and or b are zero an error message will be returned.
1997 The answer is always positive.
1998 In principle a and c can be the same.
1999*/
2000
2001#ifndef NEWTRICK
2002/*
2003 #[ Old Routine :
2004*/
2005
2006int GcdLong(PHEAD UWORD *a, WORD na, UWORD *b, WORD nb, UWORD *c, WORD *nc)
2007{
2008 GETBIDENTITY
2009 if ( !na || !nb ) {
2010 if ( !na && !nb ) {
2011/* INTERNAL_ERROR_EXCL_START */
2012 MLOCK(ErrorMessageLock);
2013 MesPrint("!>Cannot take gcd");
2014 MUNLOCK(ErrorMessageLock);
2015 return(-1);
2016/* INTERNAL_ERROR_EXCL_STOP */
2017 }
2018
2019 if ( !na ) {
2020 *nc = abs(nb);
2021 NCOPY(c,b,*nc);
2022 *nc = abs(nb);
2023 return(0);
2024 }
2025
2026 *nc = abs(na);
2027 NCOPY(c,a,*nc);
2028 *nc = abs(na);
2029 return(0);
2030 }
2031 if ( na < 0 ) na = -na;
2032 if ( nb < 0 ) nb = -nb;
2033 if ( na == 1 && nb == 1 ) {
2034#ifdef EXTRAGCD2
2035 *c = (UWORD)GCD2((ULONG)*a,(ULONG)*b);
2036#else
2037 UWORD x,y,z;
2038 x = *a;
2039 y = *b;
2040 do { z = x % y; x = y; } while ( ( y = z ) != 0 );
2041 *c = x;
2042#endif
2043 *nc = 1;
2044 }
2045 else if ( na <= 2 && nb <= 2 ) {
2046 RLONG lx,ly,lz;
2047 if ( na == 2 ) { lx = (((RLONG)(a[1]))<<BITSINWORD) + *a; }
2048 else { lx = *a; }
2049 if ( nb == 2 ) { ly = (((RLONG)(b[1]))<<BITSINWORD) + *b; }
2050 else { ly = *b; }
2051#ifdef LONGWORD
2052#ifdef EXTRAGCD2
2053 lx = GCD2(lx,ly);
2054#else
2055 do { lz = lx % ly; lx = ly; } while ( ( ly = lz ) != 0 );
2056#endif
2057#else
2058 if ( lx < ly ) { lz = lx; lx = ly; ly = lz; }
2059 do {
2060 lz = lx % ly; lx = ly;
2061 } while ( ( ly = lz ) != 0 && ( lx & AWORDMASK ) != 0 );
2062 if ( ly ) {
2063 do { *c = ((UWORD)lx)%((UWORD)ly); lx = ly; } while ( ( ly = *c ) != 0 );
2064 *c = (UWORD)lx;
2065 *nc = 1;
2066 }
2067 else
2068#endif
2069 {
2070 *c++ = (UWORD)lx;
2071 if ( ( *c = (UWORD)(lx >> BITSINWORD) ) != 0 ) *nc = 2;
2072 else *nc = 1;
2073 }
2074 }
2075 else {
2076#ifdef EXTRAGCD
2077 GCD(a,na,b,nb,c,nc);
2078#else
2079#ifdef NEWGCD
2080 UWORD *x3,*x1,*x2, *GLscrat7, *GLscrat8;
2081 WORD n1,n2,n3,n4;
2082 WORD i, j;
2083 x1 = c; x3 = a; n1 = i = na;
2084 NCOPY(x1,x3,i);
2085 GLscrat7 = NumberMalloc("GcdLong"); GLscrat8 = NumberMalloc("GcdLong");
2086 x2 = GLscrat8; x3 = b; n2 = i = nb;
2087 NCOPY(x2,x3,i);
2088 x1 = c; i = 0;
2089 while ( x1[0] == 0 ) { i += BITSINWORD; x1++; n1--; }
2090 while ( ( x1[0] & 1 ) == 0 ) { i++; SCHUIF(x1,n1) }
2091 x2 = GLscrat8; j = 0;
2092 while ( x2[0] == 0 ) { j += BITSINWORD; x2++; n2--; }
2093 while ( ( x2[0] & 1 ) == 0 ) { j++; SCHUIF(x2,n2) }
2094 if ( j > i ) j = i; /* powers of two in GCD */
2095 for(;;){
2096 if ( n1 > n2 ) {
2097firstbig:
2098 SubPLon(x1,n1,x2,n2,x1,&n3);
2099 n1 = n3;
2100 if ( n1 == 0 ) {
2101 x1 = c;
2102 n1 = i = n2; NCOPY(x1,x2,i);
2103 break;
2104 }
2105 while ( ( x1[0] & 1 ) == 0 ) SCHUIF(x1,n1)
2106 if ( n1 == 1 ) {
2107 if ( DivLong(x2,n2,x1,n1,GLscrat7,&n3,x2,&n4) ) goto GcdErr;
2108 n2 = n4;
2109 if ( n2 == 0 ) {
2110 i = n1; x2 = c; NCOPY(x2,x1,i);
2111 break;
2112 }
2113#ifdef EXTRAGCD2
2114 *c = (UWORD)GCD2((ULONG)x1[0],(ULONG)x2[0]);
2115#else
2116 {
2117 UWORD x,y,z;
2118 x = x1[0];
2119 y = x2[0];
2120 do { z = x % y; x = y; } while ( ( y = z ) != 0 );
2121 *c = x;
2122 }
2123#endif
2124 n1 = 1;
2125 break;
2126 }
2127 }
2128 else if ( n1 < n2 ) {
2129lastbig:
2130 SubPLon(x2,n2,x1,n1,x2,&n3);
2131 n2 = n3;
2132 if ( n2 == 0 ) {
2133 i = n1; x2 = c; NCOPY(x2,x1,i);
2134 break;
2135 }
2136 while ( ( x2[0] & 1 ) == 0 ) SCHUIF(x2,n2)
2137 if ( n2 == 1 ) {
2138 if ( DivLong(x1,n1,x2,n2,GLscrat7,&n3,x1,&n4) ) goto GcdErr;
2139 n1 = n4;
2140 if ( n1 == 0 ) {
2141 x1 = c;
2142 n1 = i = n2; NCOPY(x1,x2,i);
2143 break;
2144 }
2145#ifdef EXTRAGCD2
2146 *c = (UWORD)GCD2((ULONG)x2[0],(ULONG)x1[0]);
2147#else
2148 {
2149 UWORD x,y,z;
2150 x = x2[0];
2151 y = x1[0];
2152 do { z = x % y; x = y; } while ( ( y = z ) != 0 );
2153 *c = x;
2154 }
2155#endif
2156 n1 = 1;
2157 break;
2158 }
2159 }
2160 else {
2161 for ( i = n1-1; i >= 0; i-- ) {
2162 if ( x1[i] > x2[i] ) goto firstbig;
2163 else if ( x1[i] < x2[i] ) goto lastbig;
2164 }
2165 i = n1; x2 = c; NCOPY(x2,x1,i);
2166 break;
2167 }
2168 }
2169/*
2170 Now the GCD is in c but still needs j powers of 2.
2171*/
2172 x1 = c;
2173 while ( j >= BITSINWORD ) {
2174 for ( i = n1; i > 0; i-- ) x1[i] = x1[i-1];
2175 x1[0] = 0; n1++;
2176 j -= BITSINWORD;
2177 }
2178 if ( j > 0 ) {
2179 ULONG a1,a2 = 0;
2180 for ( i = 0; i < n1; i++ ) {
2181 a1 = x1[i]; a1 <<= j;
2182 a2 += a1;
2183 x1[i] = a2;
2184 a2 >>= BITSINWORD;
2185 }
2186 if ( a2 != 0 ) {
2187 x1[n1++] = a2;
2188 }
2189 }
2190 *nc = n1;
2191 NumberFree(GLscrat7,"GcdLong"); NumberFree(GLscrat8,"GcdLong");
2192#else
2193 UWORD *x1,*x2,*x3,*x4,*c1,*c2;
2194 WORD n1,n2,n3,n4,i;
2195 x1 = c; x3 = a; n1 = i = na;
2196 NCOPY(x1,x3,i);
2197 x1 = c; c1 = x2 = NumberMalloc("GcdLong"); x3 = NumberMalloc("GcdLong"); x4 = NumberMalloc("GcdLong");
2198 c2 = b; n2 = i = nb;
2199 NCOPY(c1,c2,i);
2200 for(;;){
2201 if ( DivLong(x1,n1,x2,n2,x4,&n4,x3,&n3) ) goto GcdErr;
2202 if ( !n3 ) { x1 = x2; n1 = n2; break; }
2203 if ( DivLong(x2,n2,x3,n3,x4,&n4,x1,&n1) ) goto GcdErr;
2204 if ( !n1 ) { x1 = x3; n1 = n3; break; }
2205 if ( DivLong(x3,n3,x1,n1,x4,&n4,x2,&n2) ) goto GcdErr;
2206 if ( !n2 ) {
2207 *nc = n1;
2208 NumberFree(x2,"GcdLong"); NumberFree(x3,"GcdLong"); NumberFree(x4,"GcdLong");
2209 return(0);
2210 }
2211 }
2212 *nc = i = n1;
2213 NCOPY(c,x1,i);
2214 NumberFree(x2,"GcdLong"); NumberFree(x3,"GcdLong"); NumberFree(x4,"GcdLong");
2215#endif
2216#endif
2217 }
2218 return(0);
2219GcdErr:
2220 MLOCK(ErrorMessageLock);
2221 MesCall("GcdLong");
2222 MUNLOCK(ErrorMessageLock);
2223 SETERROR(-1)
2224}
2225/*
2226 #] Old Routine :
2227*/
2228#else
2229
2230/*
2231 New routine for GcdLong that uses smart shortcuts.
2232 Algorithm by J. Vermaseren 15-nov-2006.
2233 It runs faster for very big numbers but only by a fixed factor.
2234 There is no improvement in the power behaviour.
2235 Improvement on the whole of hf9 (multiple zeta values at weight 9):
2236 Better than a factor 2 on a 32 bits architecture and 2.76 on a
2237 64 bits architecture.
2238 On hf10 (MZV's at weight 10), 64 bits architecture: factor 7.
2239
2240 If we have two long numbers (na,nb > GCDMAX) we will work in a
2241 truncated way. At the moment of writing (15-nov-2006) it isn't
2242 clear whether this algorithm is an invention or a reinvention.
2243 A short search on the web didn't show anything.
2244
2245 31-jul-2007:
2246 A better search shows that this is an adaptation of the Lehmer-Euclid
2247 algorithm, already described in Knuth. Here we can work without upper
2248 and lower limit because we are only interested in the GCD, not the
2249 extra numbers. Also it takes already some features of the double
2250 digit Lehmer-Euclid algorithm of Jebelean it seems.
2251
2252 Maybe this can be programmed slightly better and we can get another
2253 few percent speed increase. Further improvements for the asymptotic
2254 case come from splitting the calculation as in Karatsuba and working
2255 with FFT divisions and multiplications etc. But this is when hundreds
2256 of words are involved at the least.
2257
2258 Algorithm
2259
2260 1: while ( na > nb || nb < GCDMAX ) {
2261 if ( nb == 0 ) { result in a }
2262 c = a % b;
2263 a = b;
2264 b = c;
2265 }
2266 2: Make the truncated values in which a and b are the combinations
2267 of the top two words of a and b. The whole numbers are aa and bb now.
2268 3: ma1 = 1; ma2 = 0; mb1 = 0; mb2 = 1;
2269 4: A = a; B = b; m = a/b; c = a - m*b;
2270 c = ma1*a+ma2*b-m*(mb1*a+mb2*b) = (ma1-m*mb1)*a+(ma2-m*mb2)*b
2271 mc1 = ma1-m*mb1; mc2 = ma2-m*mb2;
2272 5: a = b; ma1 = mb1; ma2 = mb2;
2273 b = c; mb1 = mc1; mb2 = mc2;
2274 6: if ( b != 0 && nb >= FULLMAX ) goto 4;
2275 7: Now construct the new quantities
2276 ma1*aa+ma2*bb and mb1*aa+mb2*bb
2277 8: goto 1;
2278
2279 The essence of the above algorithm is that we do the divisions only
2280 on relatively short numbers. Also usually there are many steps 4&5
2281 for each step 7. This eliminates many operations.
2282 The termination at FULLMAX is that we make errors by not considering
2283 the tail of the number. If we run b down all the way, the errors combine
2284 in such a way that the new numbers may be of the same order as the old
2285 numbers. By stopping halfway we don't get the error beyond halfway
2286 either. Unfortunately this means that a >= FULLMAX and hence na > nb
2287 which means that next we will have a complete division. But just once.
2288 Running the steps 4-6 till a < FULLMAX runs already into problems.
2289 It may be necessary to experiment a bit to obtain the optimum value
2290 of GCDMAX.
2291*/
2292
2293int GcdLong(PHEAD UWORD *a, WORD na, UWORD *b, WORD nb, UWORD *c, WORD *nc)
2294{
2295 GETBIDENTITY
2296 UWORD x,y,z;
2297 UWORD *x1,*x2,*x3,*x4,*x5,*d;
2298 UWORD *GLscrat6, *GLscrat7, *GLscrat8, *GLscrat9, *GLscrat10;
2299 WORD n1,n2,n3,n4,n5,i;
2300 RLONG lx,ly,lz;
2301 LONG ma1, ma2, mb1, mb2, mc1, mc2, m;
2302 if ( !na || !nb ) {
2303 if ( !na && !nb ) {
2304/* INTERNAL_ERROR_EXCL_START */
2305 MLOCK(ErrorMessageLock);
2306 MesPrint("!>Cannot take gcd");
2307 MUNLOCK(ErrorMessageLock);
2308 return(-1);
2309/* INTERNAL_ERROR_EXCL_STOP */
2310 }
2311
2312 if ( !na ) {
2313 *nc = abs(nb);
2314 NCOPY(c,b,*nc);
2315 *nc = abs(nb);
2316 return(0);
2317 }
2318
2319 *nc = abs(na);
2320 NCOPY(c,a,*nc);
2321 *nc = abs(na);
2322 return(0);
2323 }
2324 if ( na < 0 ) na = -na;
2325 if ( nb < 0 ) nb = -nb;
2326/*
2327 #[ GMP stuff :
2328*/
2329#ifdef WITHGMP
2330 if ( na > 3 && nb > 3 ) {
2331 int ii;
2332 mp_limb_t *upa, *upb, *upc, xx;
2333 UWORD *uw, *u1, *u2;
2334 unsigned int tcounta, tcountb, tcounta1, tcountb1;
2335 mp_size_t ana, anb, anc;
2336
2337 u1 = uw = NumberMalloc("GcdLong");
2338 upa = (mp_limb_t *)u1;
2339 ana = na; tcounta1 = 0;
2340 while ( a[0] == 0 ) { a++; ana--; tcounta1++; }
2341 for ( ii = 0; ii < ana; ii++ ) { *uw++ = *a++; }
2342 if ( ( ana & 1 ) != 0 ) { *uw = 0; ana++; }
2343 ana /= 2;
2344
2345 u2 = uw = NumberMalloc("GcdLong");
2346 upb = (mp_limb_t *)u2;
2347 anb = nb; tcountb1 = 0;
2348 while ( b[0] == 0 ) { b++; anb--; tcountb1++; }
2349 for ( ii = 0; ii < anb; ii++ ) { *uw++ = *b++; }
2350 if ( ( anb & 1 ) != 0 ) { *uw = 0; anb++; }
2351 anb /= 2;
2352
2353 xx = upa[0]; tcounta = 0;
2354 while ( ( xx & 15 ) == 0 ) { tcounta += 4; xx >>= 4; }
2355 while ( ( xx & 1 ) == 0 ) { tcounta += 1; xx >>= 1; }
2356 xx = upb[0]; tcountb = 0;
2357 while ( ( xx & 15 ) == 0 ) { tcountb += 4; xx >>= 4; }
2358 while ( ( xx & 1 ) == 0 ) { tcountb += 1; xx >>= 1; }
2359
2360 if ( tcounta ) {
2361 mpn_rshift(upa,upa,ana,tcounta);
2362 if ( upa[ana-1] == 0 ) ana--;
2363 }
2364 if ( tcountb ) {
2365 mpn_rshift(upb,upb,anb,tcountb);
2366 if ( upb[anb-1] == 0 ) anb--;
2367 }
2368
2369 upc = (mp_limb_t *)(NumberMalloc("GcdLong"));
2370 if ( ( ana > anb ) || ( ( ana == anb ) && ( upa[ana-1] >= upb[ana-1] ) ) ) {
2371 anc = mpn_gcd(upc,upa,ana,upb,anb);
2372 }
2373 else {
2374 anc = mpn_gcd(upc,upb,anb,upa,ana);
2375 }
2376
2377 tcounta = tcounta1*BITSINWORD + tcounta;
2378 tcountb = tcountb1*BITSINWORD + tcountb;
2379 if ( tcountb > tcounta ) tcountb = tcounta;
2380 tcounta = tcountb/BITSINWORD;
2381 tcountb = tcountb%BITSINWORD;
2382
2383 if ( tcountb ) {
2384 xx = mpn_lshift(upc,upc,anc,tcountb);
2385 if ( xx ) { upc[anc] = xx; anc++; }
2386 }
2387
2388 uw = (UWORD *)upc; anc *= 2;
2389 while ( uw[anc-1] == 0 ) anc--;
2390 for ( ii = 0; ii < (int)tcounta; ii++ ) *c++ = 0;
2391 for ( ii = 0; ii < anc; ii++ ) *c++ = *uw++;
2392 *nc = anc + tcounta;
2393 NumberFree(u1,"GcdLong"); NumberFree(u2,"GcdLong"); NumberFree((UWORD *)(upc),"GcdLong");
2394 return(0);
2395 }
2396#endif
2397/*
2398 #] GMP stuff :
2399*/
2400/*
2401 #[ Easy cases :
2402*/
2403 if ( na == 1 && nb == 1 ) {
2404 x = *a;
2405 y = *b;
2406 do { z = x % y; x = y; } while ( ( y = z ) != 0 );
2407 *c = x;
2408 *nc = 1;
2409 return(0);
2410 }
2411 else if ( na <= 2 && nb <= 2 ) {
2412 if ( na == 2 ) { lx = (((RLONG)(a[1]))<<BITSINWORD) + *a; }
2413 else { lx = *a; }
2414 if ( nb == 2 ) { ly = (((RLONG)(b[1]))<<BITSINWORD) + *b; }
2415 else { ly = *b; }
2416 if ( lx < ly ) { lz = lx; lx = ly; ly = lz; }
2417#if ( BITSINWORD == 16 )
2418 do {
2419 lz = lx % ly; lx = ly;
2420 } while ( ( ly = lz ) != 0 );
2421#else
2422 do {
2423 lz = lx % ly; lx = ly;
2424 } while ( ( ly = lz ) != 0 && ( lx & AWORDMASK ) != 0 );
2425 if ( ly ) {
2426 x = (UWORD)lx; y = (UWORD)ly;
2427 do { *c = x % y; x = y; } while ( ( y = *c ) != 0 );
2428 *c = x;
2429 *nc = 1;
2430 }
2431 else
2432#endif
2433 {
2434 *c++ = (UWORD)lx;
2435 if ( ( *c = (UWORD)(lx >> BITSINWORD) ) != 0 ) *nc = 2;
2436 else *nc = 1;
2437 }
2438 return(0);
2439 }
2440/*
2441 #] Easy cases :
2442*/
2443 GLscrat6 = NumberMalloc("GcdLong"); GLscrat7 = NumberMalloc("GcdLong");
2444 GLscrat8 = NumberMalloc("GcdLong");
2445 GLscrat9 = NumberMalloc("GcdLong"); GLscrat10 = NumberMalloc("GcdLong");
2446restart:;
2447/*
2448 #[ Easy cases :
2449*/
2450 if ( na == 1 && nb == 1 ) {
2451 x = *a;
2452 y = *b;
2453 do { z = x % y; x = y; } while ( ( y = z ) != 0 );
2454 *c = x;
2455 *nc = 1;
2456 }
2457 else if ( na <= 2 && nb <= 2 ) {
2458 if ( na == 2 ) { lx = (((RLONG)(a[1]))<<BITSINWORD) + *a; }
2459 else { lx = *a; }
2460 if ( nb == 2 ) { ly = (((RLONG)(b[1]))<<BITSINWORD) + *b; }
2461 else { ly = *b; }
2462 if ( lx < ly ) { lz = lx; lx = ly; ly = lz; }
2463#if ( BITSINWORD == 16 )
2464 do {
2465 lz = lx % ly; lx = ly;
2466 } while ( ( ly = lz ) != 0 );
2467#else
2468 do {
2469 lz = lx % ly; lx = ly;
2470 } while ( ( ly = lz ) != 0 && ( lx & AWORDMASK ) != 0 );
2471 if ( ly ) {
2472 x = (UWORD)lx; y = (UWORD)ly;
2473 do { *c = x % y; x = y; } while ( ( y = *c ) != 0 );
2474 *c = x;
2475 *nc = 1;
2476 }
2477 else
2478#endif
2479 {
2480 *c++ = (UWORD)lx;
2481 if ( ( *c = (UWORD)(lx >> BITSINWORD) ) != 0 ) *nc = 2;
2482 else *nc = 1;
2483 }
2484 }
2485/*
2486 #] Easy cases :
2487 #[ Original code :
2488*/
2489 else if ( na < GCDMAX || nb < GCDMAX || na != nb ) {
2490 if ( na < nb ) {
2491 x2 = GLscrat8; x3 = a; n2 = i = na;
2492 NCOPY(x2,x3,i);
2493 x1 = c; x3 = b; n1 = i = nb;
2494 NCOPY(x1,x3,i);
2495 }
2496 else {
2497 x1 = c; x3 = a; n1 = i = na;
2498 NCOPY(x1,x3,i);
2499 x2 = GLscrat8; x3 = b; n2 = i = nb;
2500 NCOPY(x2,x3,i);
2501 }
2502 x1 = c; x2 = GLscrat8; x3 = GLscrat7; x4 = GLscrat6;
2503 for(;;){
2504 if ( DivLong(x1,n1,x2,n2,x4,&n4,x3,&n3) ) goto GcdErr;
2505 if ( !n3 ) { x1 = x2; n1 = n2; break; }
2506 if ( n2 <= 2 ) { a = x2; b = x3; na = n2; nb = n3; goto restart; }
2507 if ( n3 >= GCDMAX && n2 == n3 ) {
2508 a = GLscrat9; b = GLscrat10; na = n2; nb = n3;
2509 for ( i = 0; i < na; i++ ) a[i] = x2[i];
2510 for ( i = 0; i < nb; i++ ) b[i] = x3[i];
2511 goto newtrick;
2512 }
2513 if ( DivLong(x2,n2,x3,n3,x4,&n4,x1,&n1) ) goto GcdErr;
2514 if ( !n1 ) { x1 = x3; n1 = n3; break; }
2515 if ( n3 <= 2 ) { a = x3; b = x1; na = n3; nb = n1; goto restart; }
2516 if ( n1 >= GCDMAX && n1 == n3 ) {
2517 a = GLscrat9; b = GLscrat10; na = n3; nb = n1;
2518 for ( i = 0; i < na; i++ ) a[i] = x3[i];
2519 for ( i = 0; i < nb; i++ ) b[i] = x1[i];
2520 goto newtrick;
2521 }
2522 if ( DivLong(x3,n3,x1,n1,x4,&n4,x2,&n2) ) goto GcdErr;
2523 if ( !n2 ) { *nc = n1; goto normalend; }
2524 if ( n1 <= 2 ) { a = x1; b = x2; na = n1; nb = n2; goto restart; }
2525 if ( n2 >= GCDMAX && n2 == n1 ) {
2526 a = GLscrat9; b = GLscrat10; na = n1; nb = n2;
2527 for ( i = 0; i < na; i++ ) a[i] = x1[i];
2528 for ( i = 0; i < nb; i++ ) b[i] = x2[i];
2529 goto newtrick;
2530 }
2531 }
2532 *nc = i = n1;
2533 NCOPY(c,x1,i);
2534 }
2535/*
2536 #] Original code :
2537 #[ New code :
2538*/
2539 else {
2540/*
2541 This is the new algorithm starting at step 3.
2542
2543 3: ma1 = 1; ma2 = 0; mb1 = 0; mb2 = 1;
2544 4: A = a; B = b; m = a/b; c = a - m*b;
2545 c = ma1*a+ma2*b-m*(mb1*a+mb2*b) = (ma1-m*mb1)*a+(ma2-m*mb2)*b
2546 mc1 = ma1-m*mb1; mc2 = ma2-m*mb2;
2547 5: a = b; ma1 = mb1; ma2 = mb2;
2548 b = c; mb1 = mc1; mb2 = mc2;
2549 6: if ( b != 0 ) goto 4;
2550*/
2551newtrick:;
2552 ma1 = 1; ma2 = 0; mb1 = 0; mb2 = 1;
2553 lx = (((RLONG)(a[na-1]))<<BITSINWORD) + a[na-2];
2554 ly = (((RLONG)(b[nb-1]))<<BITSINWORD) + b[nb-2];
2555 if ( ly > lx ) { lz = lx; lx = ly; ly = lz; d = a; a = b; b = d; }
2556 do {
2557 m = lx/ly;
2558 mc1 = ma1-m*mb1; mc2 = ma2-m*mb2;
2559 ma1 = mb1; ma2 = mb2; mb1 = mc1; mb2 = mc2;
2560 lz = lx - m*ly; lx = ly; ly = lz;
2561 } while ( ly >= FULLMAX );
2562/*
2563 Next the construction of the two new numbers
2564
2565 7: Now construct the new quantities
2566 a = ma1*aa+ma2*bb and b = mb1*aa+mb2*bb
2567*/
2568 x1 = GLscrat6;
2569 x2 = GLscrat7;
2570 x3 = GLscrat8;
2571 x5 = GLscrat10;
2572 if ( ma1 < 0 ) {
2573 ma1 = -ma1;
2574 x1[0] = (UWORD)ma1;
2575 x1[1] = (UWORD)(ma1 >> BITSINWORD);
2576 if ( x1[1] ) n1 = -2;
2577 else n1 = -1;
2578 }
2579 else {
2580 x1[0] = (UWORD)ma1;
2581 x1[1] = (UWORD)(ma1 >> BITSINWORD);
2582 if ( x1[1] ) n1 = 2;
2583 else n1 = 1;
2584 }
2585 if ( MulLong(a,na,x1,n1,x2,&n2) ) goto GcdErr;
2586 if ( ma2 < 0 ) {
2587 ma2 = -ma2;
2588 x1[0] = (UWORD)ma2;
2589 x1[1] = (UWORD)(ma2 >> BITSINWORD);
2590 if ( x1[1] ) n1 = -2;
2591 else n1 = -1;
2592 }
2593 else {
2594 x1[0] = (UWORD)ma2;
2595 x1[1] = (UWORD)(ma2 >> BITSINWORD);
2596 if ( x1[1] ) n1 = 2;
2597 else n1 = 1;
2598 }
2599 if ( MulLong(b,nb,x1,n1,x3,&n3) ) goto GcdErr;
2600 if ( AddLong(x2,n2,x3,n3,c,&n4) ) goto GcdErr;
2601 if ( mb1 < 0 ) {
2602 mb1 = -mb1;
2603 x1[0] = (UWORD)mb1;
2604 x1[1] = (UWORD)(mb1 >> BITSINWORD);
2605 if ( x1[1] ) n1 = -2;
2606 else n1 = -1;
2607 }
2608 else {
2609 x1[0] = (UWORD)mb1;
2610 x1[1] = (UWORD)(mb1 >> BITSINWORD);
2611 if ( x1[1] ) n1 = 2;
2612 else n1 = 1;
2613 }
2614 if ( MulLong(a,na,x1,n1,x2,&n2) ) goto GcdErr;
2615 if ( mb2 < 0 ) {
2616 mb2 = -mb2;
2617 x1[0] = (UWORD)mb2;
2618 x1[1] = (UWORD)(mb2 >> BITSINWORD);
2619 if ( x1[1] ) n1 = -2;
2620 else n1 = -1;
2621 }
2622 else {
2623 x1[0] = (UWORD)mb2;
2624 x1[1] = (UWORD)(mb2 >> BITSINWORD);
2625 if ( x1[1] ) n1 = 2;
2626 else n1 = 1;
2627 }
2628 if ( MulLong(b,nb,x1,n1,x3,&n3) ) goto GcdErr;
2629 if ( AddLong(x2,n2,x3,n3,x5,&n5) ) goto GcdErr;
2630 a = c; na = n4; b = x5; nb = n5;
2631 if ( nb == 0 ) { *nc = n4; goto normalend; }
2632 x4 = GLscrat9;
2633 for ( i = 0; i < na; i++ ) x4[i] = a[i];
2634 a = x4;
2635 if ( na < 0 ) na = -na;
2636 if ( nb < 0 ) nb = -nb;
2637/*
2638 The typical case now is that in a we have the last step to go
2639 to loose the leading word, while in b we have lost the leading word.
2640 We could go to DivLong now but we can also add an extra step that
2641 is less wasteful.
2642 In the case that the new leading word of b is extrememly short (like 1)
2643 we make a rather large error of course. In the worst case the whole
2644 will be intercepted by DivLong after all, but that is so rare that
2645 it shouldn't influence any timing in a measurable way.
2646*/
2647 if ( nb >= GCDMAX && na == nb+1 && b[nb-1] >= HALFMAX && b[nb-1] > a[na-1] ) {
2648 lx = (((RLONG)(a[na-1]))<<BITSINWORD) + a[na-2];
2649 x1[0] = lx/b[nb-1]; n1 = 1;
2650 MulLong(b,nb,x1,n1,x2,&n2);
2651 n2 = -n2;
2652 AddLong(a,na,x2,n2,x4,&n4);
2653 if ( n4 == 0 ) {
2654 *nc = nb;
2655 for ( i = 0; i < nb; i++ ) c[i] = b[i];
2656 goto normalend;
2657 }
2658 if ( n4 < 0 ) n4 = -n4;
2659 a = b; na = nb; b = x4; nb = n4;
2660 }
2661 goto restart;
2662/*
2663 #] New code :
2664*/
2665 }
2666normalend:
2667 NumberFree(GLscrat6,"GcdLong"); NumberFree(GLscrat7,"GcdLong"); NumberFree(GLscrat8,"GcdLong");
2668 NumberFree(GLscrat9,"GcdLong"); NumberFree(GLscrat10,"GcdLong");
2669 return(0);
2670GcdErr:
2671 MLOCK(ErrorMessageLock);
2672 MesCall("GcdLong");
2673 MUNLOCK(ErrorMessageLock);
2674 NumberFree(GLscrat6,"GcdLong"); NumberFree(GLscrat7,"GcdLong"); NumberFree(GLscrat8,"GcdLong");
2675 NumberFree(GLscrat9,"GcdLong"); NumberFree(GLscrat10,"GcdLong");
2676 SETERROR(-1)
2677}
2678
2679#endif
2680
2681/*
2682 #] GcdLong :
2683 #[ GetBinom : WORD GetBinom(a,na,i1,i2)
2684*/
2685
2686int GetBinom(UWORD *a, WORD *na, WORD i1, WORD i2)
2687{
2688 GETIDENTITY
2689 WORD j, k, l;
2690 UWORD *GBscrat3, *GBscrat4;
2691 if ( i1-i2 < i2 ) i2 = i1-i2;
2692 if ( i2 == 0 ) { *a = 1; *na = 1; return(0); }
2693 if ( i2 > i1 ) { *a = 0; *na = 0; return(0); }
2694 *a = i1; *na = 1;
2695 GBscrat3 = NumberMalloc("GetBinom"); GBscrat4 = NumberMalloc("GetBinom");
2696 for ( j = 2; j <= i2; j++ ) {
2697 GBscrat3[0] = i1+1-j;
2698 if ( MulLong(a,*na,GBscrat3,(WORD)1,GBscrat4,&k) ) goto CalledFrom;
2699 GBscrat3[0] = j;
2700 if ( DivLong(GBscrat4,k,GBscrat3,(WORD)1,a,na,GBscrat3,&l) ) goto CalledFrom;
2701 }
2702 NumberFree(GBscrat3,"GetBinom"); NumberFree(GBscrat4,"GetBinom");
2703 return(0);
2704CalledFrom:
2705 MLOCK(ErrorMessageLock);
2706 MesCall("GetBinom");
2707 MUNLOCK(ErrorMessageLock);
2708 NumberFree(GBscrat3,"GetBinom"); NumberFree(GBscrat4,"GetBinom");
2709 SETERROR(-1)
2710}
2711
2712/*
2713 #] GetBinom :
2714 #[ LcmLong : WORD LcmLong(a,na,b,nb)
2715
2716 Computes the LCM of the long numbers a and b and puts the result
2717 in c. c is allowed to be equal to a.
2718*/
2719
2720int LcmLong(PHEAD UWORD *a, WORD na, UWORD *b, WORD nb, UWORD *c, WORD *nc)
2721{
2722 int error = 0;
2723 UWORD *d = NumberMalloc("LcmLong");
2724 UWORD *e = NumberMalloc("LcmLong");
2725 UWORD *f = NumberMalloc("LcmLong");
2726 WORD nd, ne, nf;
2727 GcdLong(BHEAD a, na, b, nb, d, &nd);
2728 DivLong(a,na,d,nd,e,&ne,f,&nf);
2729 if ( MulLong(b,nb,e,ne,c,nc) ) {
2730 MLOCK(ErrorMessageLock);
2731 MesCall("LcmLong");
2732 MUNLOCK(ErrorMessageLock);
2733 error = -1;
2734 }
2735 NumberFree(f,"LcmLong");
2736 NumberFree(e,"LcmLong");
2737 NumberFree(d,"LcmLong");
2738 return(error);
2739}
2740
2741/*
2742 #] LcmLong :
2743 #[ TakeLongRoot: int TakeLongRoot(a,n,power)
2744
2745 Takes the 'power'-root of the long number in a.
2746 If the root could be taken the return value is zero.
2747 If the root could not be taken, the return value is 1.
2748 The root will be in a if it could be taken, otherwise there will be garbage
2749 Algorithm: (assume b is guess of root, b' better guess)
2750 b' = (a-(power-1)*b^power)/(n*b^(power-1))
2751 Note: power should be positive!
2752*/
2753
2754int TakeLongRoot(UWORD *a, WORD *n, WORD power)
2755{
2756 GETIDENTITY
2757 int numbits, guessbits, i, retval = 0;
2758 UWORD x, *b, *c, *d, *e;
2759 WORD na, nb, nc, nd, ne;
2760 if ( *n < 0 && ( power & 1 ) == 0 ) return(1);
2761 if ( power == 1 ) return(0);
2762 if ( *n < 0 ) { na = -*n; }
2763 else { na = *n; }
2764 if ( na == 1 ) {
2765/* Special cases that are the most frequent */
2766 if ( a[0] == 1 ) return(0);
2767 if ( power < BITSINWORD && na == 1 && a[0] == (UWORD)(1<<power) ) {
2768 a[0] = 2; return(0);
2769 }
2770 if ( 2*power < BITSINWORD && na == 1 && a[0] == (UWORD)(1<<(2*power)) ) {
2771 a[0] = 4; return(0);
2772 }
2773 }
2774/*
2775 Step 1: make a guess. We count the number of bits.
2776 numbits will be the 1+2log(a)
2777*/
2778 numbits = BITSINWORD*(na-1);
2779 x = a[na-1];
2780 while ( ( x >> 8 ) != 0 ) { numbits += 8; x >>= 8; }
2781 if ( ( x >> 4 ) != 0 ) { numbits += 4; x >>= 4; }
2782 if ( ( x >> 2 ) != 0 ) { numbits += 2; x >>= 2; }
2783 if ( ( x >> 1 ) != 0 ) numbits++;
2784 guessbits = numbits / power;
2785 if ( guessbits <= 0 ) return(1); /* root < 2 and 1 we did already */
2786 nb = guessbits/BITSINWORD;
2787/*
2788 The recursion is:
2789 (b'-b) = (a/b^(power-1)-b)/n
2790 = (a/c-b)/n
2791 = (d-b)/n (remainder of a/c is e)
2792 = c/n (we reuse the scratch array c)
2793 Termination can be tricky. When a/c has no remainder and = b we have a root.
2794 When d = b but the remainder of a/c != 0, there is definitely no root.
2795*/
2796 b = NumberMalloc("TakeLongRoot"); c = NumberMalloc("TakeLongRoot");
2797 d = NumberMalloc("TakeLongRoot"); e = NumberMalloc("TakeLongRoot");
2798 for ( i = 0; i < nb; i++ ) { b[i] = 0; }
2799 b[nb] = 1 << (guessbits%BITSINWORD);
2800 nb++;
2801 for(;;) {
2802 nc = nb;
2803 for ( i = 0; i < nb; i++ ) c[i] = b[i];
2804 if ( RaisPow(BHEAD c,&nc,power-1) ) goto TLcall;
2805 if ( DivLong(a,na,c,nc,d,&nd,e,&ne) ) goto TLcall;
2806 nb = -nb;
2807 if ( AddLong(d,nd,b,nb,c,&nc) ) goto TLcall;
2808 nb = -nb;
2809 if ( nc == 0 ) {
2810 if ( ne == 0 ) break;
2811 retval = 1; break;
2812/*
2813 else {
2814 NumberFree(b,"TakeLongRoot"); NumberFree(c,"TakeLongRoot");
2815 NumberFree(d,"TakeLongRoot"); NumberFree(e,"TakeLongRoot");
2816 return(1);
2817 }
2818*/
2819 }
2820 DivLong(c,nc,(UWORD *)(&power),1,d,&nd,e,&ne);
2821 if ( nd == 0 ) {
2822 retval = 1;
2823 break;
2824/*
2825 NumberFree(b,"TakeLongRoot"); NumberFree(c,"TakeLongRoot");
2826 NumberFree(d,"TakeLongRoot"); NumberFree(e,"TakeLongRoot");
2827 return(1);
2828*/
2829/*
2830 This code tries b+1 as a final possibility.
2831 We believe this is not needed
2832 UWORD one = 1;
2833 if ( AddLong(b,nb,&one,1,c,&nc) ) goto TLcall;
2834 if ( RaisPow(BHEAD c,&nc,power-1) ) goto TLcall;
2835 if ( DivLong(a,na,c,nc,d,&nd,e,&ne) ) goto TLcall;
2836 if ( ne != 0 ) return(1);
2837 nb = -nb;
2838 if ( SubLong(d,nd,b,nb,c,&nc) ) goto TLcall;
2839 nb = -nb;
2840 if ( nc != 0 ) {
2841 NumberFree(b,"TakeLongRoot"); NumberFree(c,"TakeLongRoot");
2842 NumberFree(d,"TakeLongRoot"); NumberFree(e,"TakeLongRoot");
2843 return(1);
2844 }
2845 break;
2846*/
2847 }
2848 if ( AddLong(b,nb,d,nd,b,&nb) ) goto TLcall;
2849 }
2850 for ( i = 0; i < nb; i++ ) a[i] = b[i];
2851 if ( *n < 0 ) *n = -nb;
2852 else *n = nb;
2853 NumberFree(b,"TakeLongRoot"); NumberFree(c,"TakeLongRoot");
2854 NumberFree(d,"TakeLongRoot"); NumberFree(e,"TakeLongRoot");
2855 return(retval);
2856TLcall:
2857 MLOCK(ErrorMessageLock);
2858 MesCall("TakeLongRoot");
2859 MUNLOCK(ErrorMessageLock);
2860 NumberFree(b,"TakeLongRoot"); NumberFree(c,"TakeLongRoot");
2861 NumberFree(d,"TakeLongRoot"); NumberFree(e,"TakeLongRoot");
2862 Terminate(-1);
2863 return(-1);
2864}
2865
2866/*
2867 #] TakeLongRoot:
2868 #[ MakeRational:
2869
2870 Makes the integer a mod m into a traction b/c with |b|,|c| < sqrt(m)
2871 For the algorithm, see MakeLongRational.
2872*/
2873
2874int MakeRational(WORD a,WORD m, WORD *b, WORD *c)
2875{
2876 LONG x1,x2,x3,x4,y1,y2;
2877 if ( a < 0 ) { a = a+m; }
2878 if ( a <= 1 ) {
2879 if ( a > m/2 ) a = a-m;
2880 *b = a; *c = 1; return(0);
2881 }
2882 x1 = m; x2 = a;
2883 if ( x2*x2 >= m ) {
2884 y1 = x1/x2; y2 = x1%x2; x3 = 1; x4 = -y1; x1 = x2; x2 = y2;
2885 while ( x2*x2 >= m ) {
2886 y1 = x1/x2; y2 = x1%x2; x1 = x2; x2 = y2; y2 = x3-y1*x4; x3 = x4; x4 = y2;
2887 }
2888 }
2889 else x4 = 1;
2890 if ( x2 == 0 ) { return(1); }
2891 if ( x2 > m/2 ) *b = x2-m;
2892 else *b = x2;
2893 if ( x4 > m/2 ) { *c = x4-m; *c = -*c; *b = -*b; }
2894 else if ( x4 <= -m/2 ) { x4 += m; *c = x4; }
2895 else if ( x4 < 0 ) { x4 = -x4; *c = x4; *b = -*b; }
2896 else *c = x4;
2897 return(0);
2898}
2899
2900/*
2901 #] MakeRational:
2902 #[ MakeLongRational:
2903
2904 Converts the long number a mod m into the fraction b
2905 One of the properties of b is that num,den < sqrt(m)
2906 The algorithm: Start with: m 0
2907 a 1
2908 Make now c=m%a, c1=m/a c c2=0-c1*1
2909 Make now d=a%c d1=a/c d d2=1-d1*c2
2910 Make now e=c%d e1=c/d e e2=1-e1*d2
2911 etc till in the first column we get a number < sqrt(m)
2912 We have then f,f2 and the fraction is f/f2.
2913 If at any moment we get a zero, m contained an unlucky prime.
2914
2915 Note that this can be made a lot faster when we make the same
2916 improvements as in the GCD routine. That is something for later.
2917#ifdef WITHMAKERATIONAL
2918*/
2919
2920#define COPYLONG(x1,nx1,x2,nx2) { int i; for(i=0;i<ABS(nx2);i++)x1[i]=x2[i];nx1=nx2; }
2921
2922int MakeLongRational(PHEAD UWORD *a, WORD na, UWORD *m, WORD nm, UWORD *b, WORD *nb)
2923{
2924 UWORD *root = NumberMalloc("MakeRational");
2925 UWORD *x1 = NumberMalloc("MakeRational");
2926 UWORD *x2 = NumberMalloc("MakeRational");
2927 UWORD *x3 = NumberMalloc("MakeRational");
2928 UWORD *x4 = NumberMalloc("MakeRational");
2929 UWORD *y1 = NumberMalloc("MakeRational");
2930 UWORD *y2 = NumberMalloc("MakeRational");
2931 WORD nroot,nx1,nx2,nx3,nx4,ny1,ny2 = 0,retval = 0;
2932 WORD sign = 1;
2933/*
2934 Step 1: Take the square root of m
2935*/
2936 COPYLONG(root,nroot,m,nm)
2937 TakeLongRoot(root,&nroot,2);
2938/*
2939 Step 2: Set the start values
2940*/
2941 if ( na < 0 ) { na = -na; sign = -sign; }
2942 COPYLONG(x1,nx1,m,nm)
2943 COPYLONG(x2,nx2,a,na)
2944/*
2945 x3[0] = 0, nx3 = 0;
2946 x4[0] = 1, nx4 = 1;
2947*/
2948/*
2949 The start operation needs some special attention because of the zero.
2950*/
2951 if ( BigLong(x2,nx2,root,nroot) <= 0 ) {
2952 x4[0] = 1, nx4 = 1;
2953 goto gottheanswer;
2954 }
2955 DivLong(x1,nx1,x2,nx2,y1,&ny1,y2,&ny2);
2956 if ( ny2 == 0 ) { retval = 1; goto cleanup; }
2957 COPYLONG(x1,nx1,x2,nx2)
2958 COPYLONG(x2,nx2,y2,ny2)
2959 x3[0] = 1; nx3 = 1;
2960 COPYLONG(x4,nx4,y1,ny1)
2961 nx4 = -nx4;
2962/*
2963 Now the loop.
2964*/
2965 while ( BigLong(x2,nx2,root,nroot) > 0 ) {
2966 DivLong(x1,nx1,x2,nx2,y1,&ny1,y2,&ny2);
2967 if ( ny2 == 0 ) { retval = 1; goto cleanup; }
2968 COPYLONG(x1,nx1,x2,nx2)
2969 COPYLONG(x2,nx2,y2,ny2)
2970 MulLong(y1,ny1,x4,nx4,y2,&ny2);
2971 ny2 = -ny2;
2972 AddLong(x3,nx3,y2,ny2,y1,&ny1);
2973 COPYLONG(x3,nx3,x4,nx4)
2974 COPYLONG(x4,nx4,y1,ny1)
2975 }
2976/*
2977 Now we have the answer. It is x2/x4. It has to be packed into b.
2978*/
2979gottheanswer:
2980 if ( nx4 < 0 ) { sign = -sign; nx4 = -nx4; }
2981 COPYLONG(b,*nb,x2,nx2)
2982 Pack(b,nb,x4,nx4);
2983 if ( sign < 0 ) *nb = -*nb;
2984cleanup:
2985 NumberFree(y2,"MakeRational");
2986 NumberFree(y1,"MakeRational");
2987 NumberFree(x4,"MakeRational");
2988 NumberFree(x3,"MakeRational");
2989 NumberFree(x2,"MakeRational");
2990 NumberFree(x1,"MakeRational");
2991 NumberFree(root,"MakeRational");
2992 return(retval);
2993}
2994
2995/*
2996#endif
2997 #] MakeLongRational:
2998 #[ ChineseRemainder:
2999*/
3011#ifdef WITHCHINESEREMAINDER
3012
3013int ChineseRemainder(PHEAD MODNUM *a1, MODNUM *a2, MODNUM *a)
3014{
3015 UWORD *inv1 = NumberMalloc("ChineseRemainder");
3016 UWORD *inv2 = NumberMalloc("ChineseRemainder");
3017 UWORD *fac1 = NumberMalloc("ChineseRemainder");
3018 UWORD *fac2 = NumberMalloc("ChineseRemainder");
3019 UWORD two[1];
3020 WORD ninv1, ninv2, nfac1, nfac2;
3021 if ( a1->na < 0 ) {
3022 AddLong(a1->a,a1->na,a1->m,a1->nm,a1->a,&(a1->na));
3023 }
3024 if ( a2->na < 0 ) {
3025 AddLong(a2->a,a2->na,a2->m,a2->nm,a2->a,&(a2->na));
3026 }
3027 MulLong(a1->m,a1->nm,a2->m,a2->nm,a->m,&(a->nm));
3028
3029 GetLongModInverses(BHEAD a1->m,a1->nm,a2->m,a2->nm,inv1,&ninv1,inv2,&ninv2);
3030 MulLong(inv1,ninv1,a1->m,a1->nm,fac1,&nfac1);
3031 MulLong(inv2,ninv2,a2->m,a2->nm,fac2,&nfac2);
3032
3033 MulLong(fac1,nfac1,a2->a,a2->na,inv1,&ninv1);
3034 MulLong(fac2,nfac2,a1->a,a1->na,inv2,&ninv2);
3035 AddLong(inv1,ninv1,inv2,ninv2,a->a,&(a->na));
3036
3037 two[0] = 2;
3038 MulLong(a->a,a->na,two,1,fac1,&nfac1);
3039 if ( BigLong(fac1,nfac1,a->m,a->nm) > 0 ) {
3040 a->nm = -a->nm;
3041 AddLong(a->a,a->na,a->m,a->nm,a->a,&(a->na));
3042 a->nm = -a->nm;
3043 }
3044 NumberFree(fac2,"ChineseRemainder");
3045 NumberFree(fac1,"ChineseRemainder");
3046 NumberFree(inv2,"ChineseRemainder");
3047 NumberFree(inv1,"ChineseRemainder");
3048 return(0);
3049}
3050
3051#endif
3052
3053/*
3054 #] ChineseRemainder:
3055 #] RekenLong :
3056 #[ RekenTerms :
3057 #[ CompCoef : WORD CompCoef(term1,term2)
3058
3059 Compares the coefficients of term1 and term2 by subtracting them.
3060 This does more work than needed but this routine is only called
3061 when sorting functions and function arguments.
3062 (and comparing values
3063*/
3064/* #define 64SAVE */
3065
3066WORD CompCoef(WORD *term1, WORD *term2)
3067{
3068 GETIDENTITY
3069 UWORD *c;
3070 WORD n1,n2,n3,*a;
3071 GETCOEF(term1,n1);
3072 GETCOEF(term2,n2);
3073 if ( term1[1] == 0 && n1 == 1 ) {
3074 if ( term2[1] == 0 && n2 == 1 ) return(0);
3075 if ( n2 < 0 ) return(1);
3076 return(-1);
3077 }
3078 else if ( term2[1] == 0 && n2 == 1 ) {
3079 if ( n1 < 0 ) return(-1);
3080 return(1);
3081 }
3082 if ( n1 > 0 ) {
3083 if ( n2 < 0 ) return(1);
3084 }
3085 else {
3086 if ( n2 > 0 ) return(-1);
3087 a = term1; term1 = term2; term2 = a;
3088 n3 = -n1; n1 = -n2; n2 = n3;
3089 }
3090 if ( term1[1] == 1 && term2[1] == 1 && n1 == 1 && n2 == 1 ) {
3091 if ( (UWORD)*term1 > (UWORD)*term2 ) return(1);
3092 else if ( (UWORD)*term1 < (UWORD)*term2 ) return(-1);
3093 else return(0);
3094 }
3095
3096/*
3097 The next call should get dedicated code, as AddRat does more than
3098 strictly needed. Also more attention should be given to overflow.
3099*/
3100 c = NumberMalloc("CompCoef");
3101 if ( AddRat(BHEAD (UWORD *)term1,n1,(UWORD *)term2,-n2,c,&n3) ) {
3102 MLOCK(ErrorMessageLock);
3103 MesCall("CompCoef");
3104 MUNLOCK(ErrorMessageLock);
3105 NumberFree(c,"CompCoef");
3106 SETERROR(-1)
3107 }
3108 NumberFree(c,"CompCoef");
3109 return(n3);
3110}
3111
3112/*
3113 #] CompCoef :
3114 #[ Modulus : WORD Modulus(term)
3115
3116 Routine takes the coefficient of term modulus b. The answer
3117 is in term again and the length of term is adjusted.
3118
3119*/
3120
3121int Modulus(WORD *term)
3122{
3123 WORD *t;
3124 WORD n1;
3125 t = term;
3126 GETCOEF(t,n1);
3127 if ( TakeModulus((UWORD *)t,&n1,AC.cmod,AC.ncmod,UNPACK) ) {
3128 MLOCK(ErrorMessageLock);
3129 MesCall("Modulus");
3130 MUNLOCK(ErrorMessageLock);
3131 SETERROR(-1)
3132 }
3133 if ( !n1 ) {
3134 *term = 0;
3135 return(0);
3136 }
3137 else if ( n1 > 0 ) {
3138 n1 *= 2;
3139 t += n1; /* Note that n1 >= 0 */
3140 n1++;
3141 }
3142 else if ( n1 < 0 ) {
3143 n1 *= 2;
3144 t += -n1;
3145 n1--;
3146 }
3147 *t++ = n1;
3148 *term = WORDDIF(t,term);
3149 return(0);
3150}
3151
3152/*
3153 #] Modulus :
3154 #[ TakeModulus : WORD TakeModulus(a,na,cmodvec,ncmod,par)
3155
3156 Routine gets the rational number in a with reduced length na.
3157 It is called when AC.ncmod != 0 and the number in AC.cmod is the
3158 number wrt which we need the modulus.
3159 The result is returned in a and na again.
3160
3161 If par == NOUNPACK we only do a single number, not a fraction.
3162 In addition we don't do fancy. We want a positive number and
3163 the input was supposed to be positive.
3164 We don't pack the result. The calling routine is responsible for that.
3165 This may not be a good idea. To be checked.
3166*/
3167
3168int TakeModulus(UWORD *a, WORD *na, UWORD *cmodvec, WORD ncmod, WORD par)
3169{
3170 GETIDENTITY
3171 UWORD *c, *d, *e, *f, *g, *h;
3172 UWORD *x4,*x2;
3173 UWORD *x3,*x1,*x5,*x6,*x7,*x8;
3174 WORD y3,y1,y5,y6;
3175 WORD n1, i, y2, y4;
3176 WORD nh, tdenom, tnumer, nmod;
3177 LONG x;
3178 if ( ncmod == 0 ) return(0); /* No modulus operation */
3179 nmod = ABS(ncmod);
3180 n1 = *na;
3181 if ( ( par & UNPACK ) != 0 ) UnPack(a,n1,&tdenom,&tnumer);
3182 else { tnumer = n1; }
3183/*
3184 We fish out the special case that the coefficient is short as well.
3185 There is no need to make lots of calls etc
3186*/
3187 if ( ( ( par & UNPACK ) == 0 ) && nmod == 1 && ( n1 == 1 || n1 == -1 ) ) {
3188 goto simplecase;
3189 }
3190 else if ( nmod == 1 && ( n1 == 1 || n1 == -1 ) ) {
3191 if ( a[1] != 1 ) {
3192 a[1] = a[1] % cmodvec[0];
3193 if ( a[1] == 0 ) {
3194 MesPrint("Division by zero in short modulus arithmetic");
3195 return(-1);
3196 }
3197 y1 = 0;
3198 if ( ( AC.modinverses != 0 ) && ( ( par & NOINVERSES ) == 0 ) ) {
3199 y1 = AC.modinverses[a[1]];
3200 }
3201 else {
3202 GetModInverses(a[1],cmodvec[0],&y1,&y2);
3203 }
3204 x = a[0];
3205 a[0] = (x*y1) % cmodvec[0];
3206 a[1] = 1;
3207 }
3208 else {
3209simplecase:
3210 a[0] = a[0] % cmodvec[0];
3211 }
3212 if ( a[0] == 0 ) { *na = 0; return(0); }
3213 if ( ( AC.modmode & POSNEG ) != 0 ) {
3214 if ( a[0] > (UWORD)(cmodvec[0]/2) ) {
3215 a[0] = cmodvec[0] - a[0];
3216 *na = -*na;
3217 }
3218 }
3219 else if ( *na < 0 ) {
3220 *na = 1; a[0] = cmodvec[0] - a[0];
3221 }
3222 return(0);
3223 }
3224 c = NumberMalloc("TakeModulus"); d = NumberMalloc("TakeModulus"); e = NumberMalloc("TakeModulus");
3225 f = NumberMalloc("TakeModulus"); g = NumberMalloc("TakeModulus"); h = NumberMalloc("TakeModulus");
3226 n1 = ABS(n1);
3227 if ( DivLong(a,tnumer,(UWORD *)cmodvec,nmod,
3228 c,&nh,a,&tnumer) ) goto ModErr;
3229 if ( tnumer == 0 ) { *na = 0; goto normalreturn; }
3230 if ( ( par & UNPACK ) == 0 ) {
3231 if ( ( AC.modmode & POSNEG ) != 0 ) {
3232 NormalModulus(a,&tnumer);
3233 }
3234 else if ( tnumer < 0 ) {
3235 SubPLon((UWORD *)cmodvec,nmod,a,-tnumer,a,&tnumer);
3236 }
3237 *na = tnumer;
3238 goto normalreturn;
3239 }
3240 if ( tdenom == 1 && a[n1] == 1 ) {
3241 if ( ( AC.modmode & POSNEG ) != 0 ) {
3242 NormalModulus(a,&tnumer);
3243 }
3244 else if ( tnumer < 0 ) {
3245 SubPLon((UWORD *)cmodvec,nmod,a,-tnumer,a,&tnumer);
3246 }
3247 *na = tnumer;
3248 i = ABS(tnumer);
3249 a += i;
3250 *a++ = 1;
3251 while ( --i > 0 ) *a++ = 0;
3252 goto normalreturn;
3253 }
3254 if ( DivLong(a+n1,tdenom,(UWORD *)cmodvec,nmod,c,&nh,a+n1,&tdenom) ) goto ModErr;
3255 if ( !tdenom ) {
3256 MLOCK(ErrorMessageLock);
3257 MesPrint("Division by zero in modulus arithmetic");
3258 if ( AP.DebugFlag ) {
3259 AO.OutSkip = 3;
3260 FiniLine();
3261 i = *na;
3262 if ( i < 0 ) i = -i;
3263 while ( --i >= 0 ) { TalToLine((UWORD)(*a++)); TokenToLine((UBYTE *)" "); }
3264 i = *na;
3265 if ( i < 0 ) i = -i;
3266 while ( --i >= 0 ) { TalToLine((UWORD)(*a++)); TokenToLine((UBYTE *)" "); }
3267 TalToLine((UWORD)(*na));
3268 AO.OutSkip = 0;
3269 FiniLine();
3270 }
3271 MUNLOCK(ErrorMessageLock);
3272 NumberFree(c,"TakeModulus"); NumberFree(d,"TakeModulus"); NumberFree(e,"TakeModulus");
3273 NumberFree(f,"TakeModulus"); NumberFree(g,"TakeModulus"); NumberFree(h,"TakeModulus");
3274 return(-1);
3275 }
3276 if ( ( AC.modinverses != 0 ) && ( ( par & NOINVERSES ) == 0 )
3277 && ( tdenom == 1 || tdenom == -1 ) ) {
3278 *d = AC.modinverses[a[n1]]; y1 = 1; y2 = tdenom;
3279 if ( MulLong(a,tnumer,d,y1,c,&y3) ) goto ModErr;
3280 if ( DivLong(c,y3,(UWORD *)cmodvec,nmod,d,&y5,a,&tdenom) ) goto ModErr;
3281 if ( y2 < 0 ) tdenom = -tdenom;
3282 }
3283 else {
3284 x2 = (UWORD *)cmodvec; x1 = c; i = nmod; while ( --i >= 0 ) *x1++ = *x2++;
3285 x1 = c; x2 = a+n1; x3 = d; x4 = e; x5 = f; x6 = g;
3286 y1 = nmod; y2 = tdenom; y4 = 0; y5 = 1; *x5 = 1;
3287 for(;;) {
3288 if ( DivLong(x1,y1,x2,y2,h,&nh,x3,&y3) ) goto ModErr;
3289 if ( MulLong(x5,y5,h,nh,x6,&y6) ) goto ModErr;
3290 if ( AddLong(x4,y4,x6,-y6,x6,&y6) ) goto ModErr;
3291 if ( !y3 ) {
3292 if ( y2 != 1 || *x2 != 1 ) {
3293 MLOCK(ErrorMessageLock);
3294 MesPrint("Inverse in modulus arithmetic doesn't exist");
3295 MesPrint("Denominator and modulus are not relative prime");
3296 MUNLOCK(ErrorMessageLock);
3297 goto ModErr;
3298 }
3299 break;
3300 }
3301 x7 = x1; x1 = x2; y1 = y2; x2 = x3; y2 = y3; x3 = x7;
3302 x8 = x4; x4 = x5; y4 = y5; x5 = x6; y5 = y6; x6 = x8;
3303 }
3304 if ( y5 < 0 && AddLong((UWORD *)cmodvec,nmod,x5,y5,x5,&y5) ) goto ModErr;
3305 if ( MulLong(a,tnumer,x5,y5,c,&y3) ) goto ModErr;
3306 if ( DivLong(c,y3,(UWORD *)cmodvec,nmod,d,&y5,a,&tdenom) ) goto ModErr;
3307 }
3308 if ( !tdenom ) { *na = 0; goto normalreturn; }
3309 if ( ( ( AC.modmode & POSNEG ) != 0 ) && ( ( par & FROMFUNCTION ) == 0 ) ) {
3310 NormalModulus(a,&tdenom);
3311 }
3312 else if ( tdenom < 0 ) {
3313 SubPLon((UWORD *)cmodvec,nmod,a,-tdenom,a,&tdenom);
3314 }
3315 *na = tdenom;
3316 i = ABS(tdenom);
3317 a += i;
3318 *a++ = 1;
3319 while ( --i > 0 ) *a++ = 0;
3320normalreturn:
3321 NumberFree(c,"TakeModulus"); NumberFree(d,"TakeModulus"); NumberFree(e,"TakeModulus");
3322 NumberFree(f,"TakeModulus"); NumberFree(g,"TakeModulus"); NumberFree(h,"TakeModulus");
3323 return(0);
3324ModErr:
3325 MLOCK(ErrorMessageLock);
3326 MesCall("TakeModulus");
3327 MUNLOCK(ErrorMessageLock);
3328 NumberFree(c,"TakeModulus"); NumberFree(d,"TakeModulus"); NumberFree(e,"TakeModulus");
3329 NumberFree(f,"TakeModulus"); NumberFree(g,"TakeModulus"); NumberFree(h,"TakeModulus");
3330 SETERROR(-1)
3331}
3332
3333/*
3334 #] TakeModulus :
3335 #[ TakeNormalModulus : WORD TakeNormalModulus(a,na,par)
3336
3337 added by Jan [01-09-2010]
3338*/
3339
3340int TakeNormalModulus (UWORD *a, WORD *na, UWORD *c, WORD nc, WORD par)
3341{
3342 WORD n;
3343 WORD nhalfc;
3344 UWORD *halfc;
3345
3346 GETIDENTITY;
3347
3348 /* determine c/2 by right shifting */
3349 halfc = NumberMalloc("TakeNormalModulus");
3350 nhalfc=nc;
3351 WCOPY(halfc,c,nc);
3352
3353 for (n=0; n<nhalfc; n++) {
3354 halfc[n] /= 2;
3355 if (n+1<nc) halfc[n] |= c[n+1] << (BITSINWORD-1);
3356 }
3357
3358 if (halfc[nhalfc-1]==0)
3359 nhalfc--;
3360
3361 /* takes care of the number never expanding, e.g., -1(mod 100) -> 99 -> -1 */
3362 if (BigLong(a,ABS(*na),halfc,nhalfc) > 0) {
3363
3364 TakeModulus(a,na,c,nc,par);
3365
3366 n = ABS(*na);
3367 if (BigLong(a,n,halfc,nhalfc) > 0) {
3368 SubPLon(c,nc,a,n,a,&n);
3369 *na = (*na > 0 ? -n : n);
3370 }
3371 }
3372
3373 NumberFree(halfc,"TakeNormalModulus");
3374 return(0);
3375}
3376
3377/*
3378 #] TakeNormalModulus :
3379 #[ MakeModTable : WORD MakeModTable()
3380*/
3381
3382int MakeModTable(void)
3383{
3384 LONG size, i, j, n;
3385 n = ABS(AC.ncmod);
3386 if ( AC.modpowers ) {
3387 M_free(AC.modpowers,"AC.modpowers");
3388 AC.modpowers = NULL;
3389 }
3390 if ( n > 2 ) {
3391/* INTERNAL_ERROR_EXCL_START */
3392 MLOCK(ErrorMessageLock);
3393 MesPrint("!>No memory for modulus generator power table");
3394 MUNLOCK(ErrorMessageLock);
3395 Terminate(-1);
3396/* INTERNAL_ERROR_EXCL_STOP */
3397 }
3398 if ( n == 0 ) return(0);
3399 size = (LONG)(*AC.cmod);
3400 if ( n == 2 ) size += (((LONG)AC.cmod[1])<<BITSINWORD);
3401 AC.modpowers = (UWORD *)Malloc1(size*n*sizeof(UWORD),"table for powers of modulus");
3402 if ( n == 1 ) {
3403 j = 1;
3404 for ( i = 0; i < size; i++ ) AC.modpowers[i] = 0;
3405 for ( i = 0; i < size; i++ ) {
3406 AC.modpowers[j] = (WORD)i;
3407 j *= *AC.powmod;
3408 j %= *AC.cmod;
3409 }
3410 for ( i = 2; i < size; i++ ) {
3411 if ( AC.modpowers[i] == 0 ) {
3412 MLOCK(ErrorMessageLock);
3413 MesPrint("&improper generator for this modulus");
3414 MUNLOCK(ErrorMessageLock);
3415 M_free(AC.modpowers,"AC.modpowers");
3416 return(-1);
3417 }
3418 }
3419 AC.modpowers[1] = 0;
3420 }
3421 else {
3422 GETIDENTITY
3423 WORD nScrat, n2;
3424 UWORD *MMscrat7 = NumberMalloc("MakeModTable"), *MMscratC = NumberMalloc("MakeModTable");
3425 *MMscratC = 1;
3426 nScrat = 1;
3427 j = size * 2;
3428 for ( i = 0; i < j; i+=2 ) { AC.modpowers[i] = 0; AC.modpowers[i+1] = 0; }
3429 for ( i = 0; i < size; i++ ) {
3430 j = *MMscratC + (((LONG)MMscratC[1])<<BITSINWORD);
3431 j *= 2;
3432 AC.modpowers[j] = (WORD)(i & WORDMASK);
3433 AC.modpowers[j+1] = (WORD)(i >> BITSINWORD);
3434 MulLong((UWORD *)MMscratC,nScrat,(UWORD *)AC.powmod,
3435 AC.npowmod,(UWORD *)MMscrat7,&n2);
3436 TakeModulus(MMscrat7,&n2,AC.cmod,AC.ncmod,NOUNPACK);
3437 *MMscratC = *MMscrat7; MMscratC[1] = MMscrat7[1]; nScrat = n2;
3438 }
3439 NumberFree(MMscrat7,"MakeModTable"); NumberFree(MMscratC,"MakeModTable");
3440 j = size * 2;
3441 for ( i = 4; i < j; i+=2 ) {
3442 if ( AC.modpowers[i] == 0 && AC.modpowers[i+1] == 0 ) {
3443 MLOCK(ErrorMessageLock);
3444 MesPrint("&improper generator for this modulus");
3445 MUNLOCK(ErrorMessageLock);
3446 M_free(AC.modpowers,"AC.modpowers");
3447 return(-1);
3448 }
3449 }
3450 AC.modpowers[2] = AC.modpowers[3] = 0;
3451 }
3452 return(0);
3453}
3454
3455/*
3456 #] MakeModTable :
3457 #] RekenTerms :
3458 #[ Functions :
3459 #[ Factorial : WORD Factorial(n,a,na)
3460
3461 Starts with only the value of fac_(0).
3462 Builds up what is needed and remembers it for the next time.
3463
3464 We have:
3465 AT.nfac: the number of the highest stored factorial
3466 AT.pfac: the array of locations in the array of stored factorials
3467 AT.factorials: the array with stored factorials
3468*/
3469
3470int Factorial(PHEAD WORD n, UWORD *a, WORD *na)
3471{
3472 GETBIDENTITY
3473 UWORD *b, *c;
3474 WORD nc;
3475 int i, j;
3476 LONG ii;
3477 if ( n > AT.nfac ) {
3478 if ( AT.factorials == 0 ) {
3479 AT.nfac = 0; AT.mfac = 50; AT.sfact = 400;
3480 AT.pfac = (LONG *)Malloc1((AT.mfac+2)*sizeof(LONG),"factorials");
3481 AT.factorials = (UWORD *)Malloc1(AT.sfact*sizeof(UWORD),"factorials");
3482 AT.factorials[0] = 1; AT.pfac[0] = 0; AT.pfac[1] = 1;
3483 }
3484 b = a;
3485 c = AT.factorials+AT.pfac[AT.nfac];
3486 nc = i = AT.pfac[AT.nfac+1] - AT.pfac[AT.nfac];
3487 while ( --i >= 0 ) *b++ = *c++;
3488 for ( j = AT.nfac+1; j <= n; j++ ) {
3489 Product(a,&nc,j);
3490 if ( nc > AM.MaxTal ) {
3491 MLOCK(ErrorMessageLock);
3492 MesPrint("Overflow in factorial. MaxTal = %d",AM.MaxTal);
3493 MesPrint("Increase MaxTerm in %s",setupfilename);
3494 MUNLOCK(ErrorMessageLock);
3495 return(-1);
3496 }
3497 if ( j > AT.mfac ) { /* double the pfac buffer */
3498 LONG *p;
3499 p = (LONG *)Malloc1((AT.mfac*2+2)*sizeof(LONG),"factorials");
3500 i = AT.mfac;
3501 for ( i = AT.mfac+1; i >= 0; i-- ) p[i] = AT.pfac[i];
3502 M_free(AT.pfac,"factorial offsets"); AT.pfac = p; AT.mfac *= 2;
3503 }
3504 if ( AT.pfac[j] + nc >= AT.sfact ) { /* double the factorials buffer */
3505 UWORD *f;
3506 f = (UWORD *)Malloc1(AT.sfact*2*sizeof(UWORD),"factorials");
3507 ii = AT.sfact;
3508 c = AT.factorials; b = f;
3509 while ( --ii >= 0 ) *b++ = *c++;
3510 M_free(AT.factorials,"factorials");
3511 AT.factorials = f;
3512 AT.sfact *= 2;
3513 }
3514 b = a; c = AT.factorials + AT.pfac[j]; i = nc;
3515 while ( --i >= 0 ) *c++ = *b++;
3516 AT.pfac[j+1] = AT.pfac[j] + nc;
3517 }
3518 *na = nc;
3519 AT.nfac = n;
3520 }
3521 else if ( n == 0 ) {
3522 *a = 1; *na = 1;
3523 }
3524 else {
3525 *na = i = AT.pfac[n+1] - AT.pfac[n];
3526 b = AT.factorials + AT.pfac[n];
3527 while ( --i >= 0 ) *a++ = *b++;
3528 }
3529 return(0);
3530}
3531
3532/*
3533 #] Factorial :
3534 #[ Bernoulli : WORD Bernoulli(n,a,na)
3535
3536 Starts with only the value of bernoulli_(0).
3537 Builds up what is needed and remembers it for the next time.
3538 b_0 = 1
3539 (n+1)*b_n = -b_{n-1}-sum_(i,1,n-1,b_i*b_{n-i})
3540 The n-1 plays only a role for b_2.
3541 We have hard coded b_0,b_1,b_2 and b_odd. After that:
3542 (2n+1)*b_2n = -sum_(i,1,n-1,b_2i*b_{2n-2i})
3543
3544 We have:
3545 AT.nBer: the number of the highest stored Bernoulli number
3546 AT.pBer: the array of locations in the array of stored Bernoulli numbers
3547 AT.bernoullis: the array with stored Bernoulli numbers
3548*/
3549
3550int Bernoulli(WORD n, UWORD *a, WORD *na)
3551{
3552 GETIDENTITY
3553 UWORD *b, *c, *scrib, *ntop, *ntop1;
3554 WORD i, i1, i2, nhalf, nqua, nscrib, nntop, nntop1, *oldworkpointer;
3555 UWORD twee = 2, twonplus1;
3556 int j;
3557 LONG ii;
3558 if ( n <= 1 ) {
3559 if ( n == 0 ) { a[0] = a[1] = 1; *na = 3; }
3560 else if ( n == 1 ) { a[0] = 1; a[1] = 2; *na = 3; }
3561 return(0);
3562 }
3563 if ( ( n & 1 ) != 0 ) { a[0] = a[1] = 0; *na = 0; return(0); }
3564 nhalf = n/2;
3565 if ( nhalf > AT.nBer ) {
3566 oldworkpointer = AT.WorkPointer;
3567 if ( AT.bernoullis == 0 ) {
3568 AT.nBer = 1; AT.mBer = 50; AT.sBer = 400;
3569 AT.pBer = (LONG *)Malloc1((AT.mBer+2)*sizeof(LONG),"bernoullis");
3570 AT.bernoullis = (UWORD *)Malloc1(AT.sBer*sizeof(UWORD),"bernoullis");
3571 AT.pBer[1] = 0; AT.pBer[2] = 3;
3572 AT.bernoullis[0] = 3; AT.bernoullis[1] = 1; AT.bernoullis[2] = 12;
3573 if ( nhalf == 1 ) {
3574 a[0] = 1; a[1] = 12; *na = 3; return(0);
3575 }
3576 }
3577 while ( nhalf > AT.mBer ) {
3578 LONG *p;
3579 p = (LONG *)Malloc1((AT.mBer*2+1)*sizeof(LONG),"bernoullis");
3580 i = AT.mBer;
3581 for ( i = AT.mBer; i >= 0; i-- ) p[i] = AT.pBer[i];
3582 M_free(AT.pBer,"factorial pointers"); AT.pBer = p; AT.mBer *= 2;
3583 }
3584 for ( n = AT.nBer+1; n <= nhalf; n++ ) {
3585 scrib = (UWORD *)(AT.WorkPointer);
3586 nqua = n/2;
3587 if ( ( n & 1 ) == 1 ) {
3588 nscrib = 0; ntop = scrib;
3589 }
3590 else {
3591 b = AT.bernoullis + AT.pBer[nqua];
3592 nscrib = *b++;
3593 i = (WORD)(REDLENG(nscrib));
3594 MulRat(BHEAD b,i,b,i,scrib,&nscrib);
3595 ntop = scrib + 2*nscrib;
3596 nqua--;
3597 }
3598 for ( j = 1; j <= nqua; j++ ) {
3599 b = AT.bernoullis + AT.pBer[j];
3600 c = AT.bernoullis + AT.pBer[n-j];
3601 i1 = (WORD)(*b); i2 = (WORD)(*c);
3602 i1 = REDLENG(i1);
3603 i2 = REDLENG(i2);
3604 MulRat(BHEAD b+1,i1,c+1,i2,ntop,&nntop);
3605 Mully(BHEAD ntop,&nntop,&twee,1);
3606 if ( nscrib ) {
3607 i = (WORD)nntop; if ( i < 0 ) i = -i;
3608 ntop1 = ntop + 2*i;
3609 AddRat(BHEAD ntop,nntop,scrib,nscrib,ntop1,&nntop1);
3610 }
3611 else {
3612 ntop1 = ntop; nntop1 = nntop;
3613 }
3614 nscrib = i1 = (WORD)nntop1;
3615 if ( i1 < 0 ) i1 = - i1;
3616 i1 = 2*i1;
3617 for ( i = 0; i < i1; i++ ) scrib[i] = ntop1[i];
3618 ntop = scrib + i1;
3619 }
3620 twonplus1 = 2*n+1;
3621 Divvy(BHEAD scrib,&nscrib,&twonplus1,-1);
3622 i1 = INCLENG(nscrib);
3623 i2 = i1; if ( i2 < 0 ) i2 = -i2;
3624 i = (WORD)(AT.bernoullis[AT.pBer[n-1]]);
3625 if ( i < 0 ) i = -i;
3626 AT.pBer[n] = AT.pBer[n-1]+i;
3627 if ( AT.pBer[n] + i2 >= AT.sBer ) {
3628 UWORD *f;
3629 f = (UWORD *)Malloc1(AT.sBer*2*sizeof(UWORD),"bernoullis");
3630 ii = AT.sBer;
3631 c = AT.bernoullis; b = f;
3632 while ( --ii >= 0 ) *b++ = *c++;
3633 M_free(AT.bernoullis,"bernoullis");
3634 AT.bernoullis = f;
3635 AT.sBer *= 2;
3636 }
3637 c = AT.bernoullis + AT.pBer[n]; b = scrib;
3638 *c++ = i1;
3639 for ( i = 1; i < i2; i++ ) *c++ = *b++;
3640 }
3641 AT.nBer = nhalf;
3642 AT.WorkPointer = oldworkpointer;
3643 }
3644 b = AT.bernoullis + AT.pBer[nhalf];
3645 *na = i = (WORD)(*b++);
3646 if ( i < 0 ) i = -i;
3647 i--;
3648 while ( --i >= 0 ) *a++ = *b++;
3649 return(0);
3650}
3651
3652/*
3653 #] Bernoulli :
3654 #[ NextPrime :
3655*/
3667#if ( BITSINWORD == 32 )
3668
3669void StartPrimeList(PHEAD0)
3670{
3671 int i, j;
3672 AR.PrimeList[AR.numinprimelist++] = 3;
3673 for ( i = 5; i < 46340; i += 2 ) {
3674 for ( j = 0; j < AR.numinprimelist && AR.PrimeList[j]*AR.PrimeList[j] <= i; j++ ) {
3675 if ( i % AR.PrimeList[j] == 0 ) goto nexti;
3676 }
3677 AR.PrimeList[AR.numinprimelist++] = i;
3678nexti:;
3679 }
3680 AR.notfirstprime = 1;
3681}
3682
3683#endif
3684
3685WORD NextPrime(PHEAD WORD num)
3686{
3687 int i, j;
3688 WORD *newpl;
3689 LONG newsize, x;
3690#if ( BITSINWORD == 32 )
3691 if ( AR.notfirstprime == 0 ) StartPrimeList(BHEAD0);
3692#endif
3693 if ( num > AT.inprimelist ) {
3694 while ( AT.inprimelist < num ) {
3695 if ( num >= AT.sizeprimelist ) {
3696 if ( AT.sizeprimelist == 0 ) newsize = 32;
3697 else newsize = 2*AT.sizeprimelist;
3698 while ( num >= newsize ) newsize = newsize*2;
3699 newpl = (WORD *)Malloc1(newsize*sizeof(WORD),"NextPrime");
3700 for ( i = 0; i < AT.sizeprimelist; i++ ) {
3701 newpl[i] = AT.primelist[i];
3702 }
3703 if ( AT.sizeprimelist > 0 ) {
3704 M_free(AT.primelist,"NextPrime");
3705 }
3706 AT.sizeprimelist = newsize;
3707 AT.primelist = newpl;
3708 }
3709 if ( AT.inprimelist < 0 ) { i = MAXPOSITIVE; }
3710 else { i = AT.primelist[AT.inprimelist]; }
3711 while ( i > MAXPOWER ) {
3712 i -= 2; x = i;
3713#if ( BITSINWORD == 32 )
3714 for ( j = 0; j < AR.numinprimelist && AR.PrimeList[j]*(LONG)(AR.PrimeList[j]) <= x; j++ ) {
3715 if ( x % AR.PrimeList[j] == 0 ) goto nexti;
3716 }
3717#else
3718 for ( j = 3; j*((LONG)j) <= x; j += 2 ) {
3719 if ( x % j == 0 ) goto nexti;
3720 }
3721#endif
3722 AT.inprimelist++;
3723 AT.primelist[AT.inprimelist] = i;
3724 break;
3725nexti:;
3726 }
3727 if ( i < MAXPOWER ) {
3728/* INTERNAL_ERROR_EXCL_START */
3729 MLOCK(ErrorMessageLock);
3730 MesPrint("!>There are not enough short prime numbers for this calculation");
3731 MesPrint("Try to use a computer with a %d-bits architecture",
3732 (int)(BITSINWORD*4));
3733 MUNLOCK(ErrorMessageLock);
3734 Terminate(-1);
3735/* INTERNAL_ERROR_EXCL_STOP */
3736 }
3737 }
3738 }
3739 return(AT.primelist[num]);
3740}
3741
3742/*
3743 #] NextPrime :
3744 #[ Moebius :
3745
3746 Returns the value of the Moebius function if n fits inside
3747 a WORD.
3748 The method we use is a bit like the sieve of Erathostenes.
3749*/
3750
3751WORD Moebius(PHEAD WORD nn)
3752{
3753 WORD i,n = nn, x;
3754 LONG newsize;
3755 SBYTE *newtable, mu;
3756#if ( BITSINWORD == 32 )
3757 if ( AR.notfirstprime == 0 ) StartPrimeList(BHEAD0);
3758#endif
3759/*
3760 First we make sure that:
3761 a: the table is big enough.
3762 b: the number is not already in the table.
3763*/
3764 if ( nn >= AR.moebiustablesize ) {
3765 if ( AR.moebiustablesize <= 0 ) { newsize = (LONG)nn + 20; }
3766 else { newsize = (LONG)nn*2; }
3767 if ( newsize > MAXPOSITIVE ) newsize = MAXPOSITIVE;
3768 newtable = (SBYTE *)Malloc1(newsize*sizeof(SBYTE),"Moebius");
3769 for ( i = 0; i < AR.moebiustablesize; i++ ) newtable[i] = AR.moebiustable[i];
3770 for ( ; i < newsize; i++ ) newtable[i] = 2;
3771 if ( AR.moebiustablesize > 0 ) M_free(AR.moebiustable,"Moebius");
3772 AR.moebiustable = newtable;
3773 AR.moebiustablesize = newsize;
3774 }
3775 /* NOTE: nn == MAXPOSITIVE never fits in moebiustable. */
3776 if ( nn != MAXPOSITIVE && AR.moebiustable[nn] != 2 ) return((WORD)AR.moebiustable[nn]);
3777 mu = 1;
3778 if ( n == 1 ) goto putvalue;
3779 if ( n % 2 == 0 ) {
3780 n /= 2;
3781 if ( n % 2 == 0 ) { mu = 0; goto putvalue; }
3782 if ( AR.moebiustable[n] != 2 ) { mu = -AR.moebiustable[n]; goto putvalue; }
3783 mu = -mu;
3784 if ( n == 1 ) goto putvalue;
3785 }
3786#if ( BITSINWORD == 32 )
3787 for ( i = 0; i < AR.numinprimelist; i++ ) {
3788 x = AR.PrimeList[i];
3789#else
3790 for ( x = 3; x < MAXPOSITIVE; x += 2 ) {
3791#endif
3792 if ( n % x == 0 ) {
3793 n /= x;
3794 if ( n % x == 0 ) { mu = 0; goto putvalue; }
3795 if ( AR.moebiustable[n] != 2 ) { mu = -AR.moebiustable[n]; goto putvalue; }
3796 mu = -mu;
3797 if ( n == 1 ) goto putvalue;
3798 }
3799 if ( n < x*x ) break; /* notice that x*x always fits inside a WORD */
3800 }
3801 mu = -mu;
3802putvalue:
3803 if ( nn != MAXPOSITIVE ) AR.moebiustable[nn] = mu;
3804 return((WORD)mu);
3805}
3806
3807/*
3808 #] Moebius :
3809 #[ wranf :
3810
3811 A random number generator that generates random WORDs with a very
3812 long sequence. It is based on the Knuth generator.
3813
3814 We take some care that each thread can run its own, but each
3815 uses its own startup. Hence the seed includes the identity of
3816 the thread.
3817
3818 For NPAIR1, NPAIR2 we can use any pair from the table on page 28.
3819 Candidates are 24,55 (the example on the pages 171,172)
3820 or (33,97) or (38,89)
3821 These values are defined in fsizes.h and used in startup.c and threads.c
3822*/
3823
3824#define WARMUP 6
3825
3826static void wranfnew(PHEAD0)
3827{
3828 int i;
3829 LONG j;
3830 for ( i = 0; i < AR.wranfnpair1; i++ ) {
3831 j = AR.wranfia[i] - AR.wranfia[i+(AR.wranfnpair2-AR.wranfnpair1)];
3832 if ( j < 0 ) j += (LONG)1 << (2*BITSINWORD-2);
3833 AR.wranfia[i] = j;
3834 }
3835 for ( i = AR.wranfnpair1; i < AR.wranfnpair2; i++ ) {
3836 j = AR.wranfia[i] - AR.wranfia[i-AR.wranfnpair1];
3837 if ( j < 0 ) j += (LONG)1 << (2*BITSINWORD-2);
3838 AR.wranfia[i] = j;
3839 }
3840}
3841
3842void iniwranf(PHEAD0)
3843{
3844 int imax = AR.wranfnpair2-1;
3845 ULONG i, ii, seed = AR.wranfseed;
3846 LONG j, k;
3847 ULONG offset = 12345;
3848#ifdef PARALLELCODE
3849 int id;
3850#if defined(WITHPTHREADS)
3851 id = AT.identity;
3852#elif defined(WITHMPI)
3853 id = PF.me;
3854#endif
3855 seed += id;
3856 i = id + 1;
3857 if ( i > 1 ) {
3858 ULONG pow, accu;
3859 pow = offset; accu = 1;
3860 while ( i ) {
3861 if ( ( i & 1 ) != 0 ) accu *= pow;
3862 i /= 2; pow = pow*pow;
3863 }
3864 offset = accu;
3865 }
3866#endif
3867 if ( seed < ((LONG)1<<(BITSINWORD-1)) ) {
3868 j = ( (seed+31459L) << (BITSINWORD-2))+offset;
3869 }
3870 else if ( seed < ((LONG)1<<(BITSINWORD+10-1)) ) {
3871 j = ( (seed+31459L) << (BITSINWORD-10-2))+offset;
3872 }
3873 else {
3874 j = ( (seed+31459L) << 1)+offset;
3875 }
3876 if ( ( seed & 1 ) == 1 ) seed++;
3877 j += seed;
3878 AR.wranfia[imax] = j;
3879 k = 1;
3880 for ( i = 0; i <= (ULONG)(imax); i++ ) {
3881 ii = (AR.wranfnpair1*i)%AR.wranfnpair2;
3882 AR.wranfia[ii] = k;
3883 k = ULongToLong((ULONG)j - (ULONG)k);
3884 if ( k < 0 ) k += (LONG)1 << (2*BITSINWORD-2);
3885 j = AR.wranfia[ii];
3886 }
3887 for ( i = 0; i < WARMUP; i++ ) wranfnew(BHEAD0);
3888 AR.wranfcall = 0;
3889}
3890
3891UWORD wranf(PHEAD0)
3892{
3893 UWORD wval;
3894 if ( AR.wranfia == 0 ) {
3895 AR.wranfia = (ULONG *)Malloc1(AR.wranfnpair2*sizeof(ULONG),"wranf");
3896 iniwranf(BHEAD0);
3897 }
3898 if ( AR.wranfcall >= AR.wranfnpair2) {
3899 wranfnew(BHEAD0);
3900 AR.wranfcall = 0;
3901 }
3902 wval = (UWORD)(AR.wranfia[AR.wranfcall++]>>(BITSINWORD-1));
3903 return(wval);
3904}
3905
3906/*
3907 Returns a random UWORD in the range (0,...,imax-1)
3908*/
3909
3910UWORD iranf(PHEAD UWORD imax)
3911{
3912 UWORD i;
3913 ULONG x, xmax;
3914 if (imax < 2) return 0;
3915 x = (LONG)1 << BITSINWORD;
3916 xmax = x - x%imax;
3917 while ( ( i = wranf(BHEAD0) ) >= xmax ) {}
3918 return(i%imax);
3919}
3920
3921/*
3922 #] wranf :
3923 #[ PreRandom :
3924
3925 The random number generator of the preprocessor.
3926 This one is completely different from the execution time generator
3927 random_(number). In the preprocessor we generate a floating point
3928 number in a string according to a distribution.
3929 Currently allowed are:
3930 RANDOM_(log,min,max)
3931 RANDOM_(lin,min,max)
3932 The return value is a string with the floating point number.
3933*/
3934
3935UBYTE *PreRandom(UBYTE *s)
3936{
3937 GETIDENTITY
3938 UBYTE *mode,*mins = 0,*maxs = 0, *outval;
3939 float num;
3940 double minval, maxval, value = 0;
3941 int linlog = -1;
3942 mode = s;
3943 while ( FG.cTable[*s] <= 1 ) s++;
3944 if ( *s == ',' ) { *s = 0; s++; }
3945 mins = s;
3946 while ( *s && *s != ',' ) s++;
3947 if ( *s == ',' ) { *s = 0; s++; }
3948 maxs = s;
3949 while ( *s && *s != ',' ) s++;
3950 if ( *s || *maxs == 0 || *mins == 0 ) {
3951 MesPrint("@Illegal arguments in macro RANDOM_");
3952 Terminate(-1);
3953 }
3954 if ( StrICmp(mode,(UBYTE *)"lin") == 0 ) {
3955 linlog = 0;
3956 }
3957 else if ( StrICmp(mode,(UBYTE *)"log") == 0 ) {
3958 linlog = 1;
3959 }
3960 else {
3961 MesPrint("@Illegal mode argument in macro RANDOM_");
3962 Terminate(-1);
3963 }
3964
3965 sscanf((char *)mins,"%f",&num); minval = num;
3966 sscanf((char *)maxs,"%f",&num); maxval = num;
3967
3968 /*
3969 * Note on ParFORM: we should use the same random number on all the
3970 * processes in the complication phase. The random number is generated
3971 * on the master and broadcast to the other processes.
3972 */
3973 {
3974 UWORD x;
3975 double xx;
3976#ifdef WITHMPI
3977 x = 0;
3978 if ( PF.me == MASTER ) {
3979 x = wranf(BHEAD0);
3980 }
3981 x = (UWORD)PF_BroadcastNumber((LONG)x);
3982#else
3983 x = wranf(BHEAD0);
3984#endif
3985 xx = x/pow(2.0,(double)(BITSINWORD-1));
3986 if ( linlog == 0 ) {
3987 value = minval + (maxval-minval)*xx;
3988 }
3989 else if ( linlog == 1 ) {
3990 value = minval * pow(maxval/minval,xx);
3991 }
3992 }
3993
3994 outval = (UBYTE *)Malloc1(64,"PreRandom");
3995 if ( ABS(value) < 0.00001 || ABS(value) > 1000000. ) {
3996 snprintf((char *)outval,64,"%e",value);
3997 }
3998 else if ( ABS(value) < 0.0001 ) { snprintf((char *)outval,64,"%10f",value); }
3999 else if ( ABS(value) < 0.001 ) { snprintf((char *)outval,64,"%9f",value); }
4000 else if ( ABS(value) < 0.01 ) { snprintf((char *)outval,64,"%8f",value); }
4001 else if ( ABS(value) < 0.1 ) { snprintf((char *)outval,64,"%7f",value); }
4002 else if ( ABS(value) < 1. ) { snprintf((char *)outval,64,"%6f",value); }
4003 else if ( ABS(value) < 10. ) { snprintf((char *)outval,64,"%5f",value); }
4004 else if ( ABS(value) < 100. ) { snprintf((char *)outval,64,"%4f",value); }
4005 else if ( ABS(value) < 1000. ) { snprintf((char *)outval,64,"%3f",value); }
4006 else if ( ABS(value) < 10000. ) { snprintf((char *)outval,64,"%2f",value); }
4007 else { snprintf((char *)outval,64,"%1f",value); }
4008 return(outval);
4009}
4010
4011/*
4012 #] PreRandom :
4013 #] Functions :
4014*/
PF_BroadcastNumber
LONG PF_BroadcastNumber(LONG x)
Definition parallel.c:2098
GetModInverses
int GetModInverses(WORD m1, WORD m2, WORD *im1, WORD *im2)
Definition reken.c:1489
RaisPowCached
void RaisPowCached(PHEAD WORD x, WORD n, UWORD **c, WORD *nc)
Definition reken.c:1309
NormalModulus
int NormalModulus(UWORD *a, WORD *na)
Definition reken.c:1416
MakeInverses
int MakeInverses(void)
Definition reken.c:1453
NextPrime
WORD NextPrime(PHEAD WORD num)
Definition reken.c:3685
CompCoef
WORD CompCoef(WORD *term1, WORD *term2)
Definition reken.c:3066
Generated by doxygen 1.9.8