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