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