FORM  4.1
reken.c
Go to the documentation of this file.
1 
11 /* #[ License : */
12 /*
13  * Copyright (C) 1984-2013 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 <<= 1;
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 <<= 1;
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 <<= 1;
386  while ( --i >= 0 ) *c++ = *a++;
387  return(0);
388  }
389  else if ( nb == -1 ) {
390  *nc = - na;
391  i = ABS(na); i <<= 1;
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 <<= 1;
400  while ( --i >= 0 ) *c++ = *b++;
401  return(0);
402  }
403  else if ( na == -1 ) {
404  *nc = - nb;
405  i = ABS(nb); i <<= 1;
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 >>= 1;
1236  }
1237  i--;
1238  c = 1 << i;
1239  while ( --i >= 0 ) {
1240  c >>= 1;
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  AT.small_power[ID] = NumberMalloc("RaisPowCached");
1334  AT.small_power_n[ID] = 1;
1335  AT.small_power[ID][0] = x;
1336  RaisPow(BHEAD AT.small_power[ID],&AT.small_power_n[ID],n);
1337  }
1338 
1339  /* return the result */
1340  *c = AT.small_power[ID];
1341  *nc = AT.small_power_n[ID];
1342 }
1343 
1344 /*
1345  #] RaisPowCached :
1346  #[ RaisPowMod :
1347 
1348  Computes the power x^n mod m
1349  */
1350 WORD RaisPowMod (WORD x, WORD n, WORD m) {
1351  LONG y=1, z=x;
1352  while (n) {
1353  if (n&1) { y*=z; y%=m; }
1354  z*=z; z%=m;
1355  n>>=1;
1356  }
1357  return (WORD)y;
1358 }
1359 
1360 /*
1361  #] RaisPowMod :
1362  #[ NormalModulus : int NormalModulus(UWORD *a,WORD *na)
1363 */
1370 int NormalModulus(UWORD *a,WORD *na)
1371 {
1372  WORD n;
1373  if ( AC.halfmod == 0 ) {
1374  LOCK(AC.halfmodlock);
1375  if ( AC.halfmod == 0 ) {
1376  UWORD two[1],remain[1];
1377  WORD dummy;
1378  two[0] = 2;
1379  AC.halfmod = (UWORD *)Malloc1((ABS(AC.ncmod))*sizeof(UWORD),"halfmod");
1380  DivLong((UWORD *)AC.cmod,(ABS(AC.ncmod)),two,1
1381  ,(UWORD *)AC.halfmod,&(AC.nhalfmod),remain,&dummy);
1382  }
1383  UNLOCK(AC.halfmodlock);
1384  }
1385  n = ABS(*na);
1386  if ( BigLong(a,n,AC.halfmod,AC.nhalfmod) > 0 ) {
1387  SubPLon((UWORD *)AC.cmod,(ABS(AC.ncmod)),a,n,a,&n);
1388  if ( *na > 0 ) { *na = -n; }
1389  else { *na = n; }
1390  return(1);
1391  }
1392  return(0);
1393 }
1394 
1395 /*
1396  #] NormalModulus :
1397  #[ MakeInverses :
1398 */
1408 {
1409  WORD n = AC.cmod[0], i, inv2;
1410  if ( AC.ncmod != 1 ) return(1);
1411  if ( AC.modinverses == 0 ) {
1412  LOCK(AC.halfmodlock);
1413  if ( AC.modinverses == 0 ) {
1414  AC.modinverses = (UWORD *)Malloc1(n*sizeof(UWORD),"modinverses");
1415  AC.modinverses[0] = 0;
1416  AC.modinverses[1] = 1;
1417  for ( i = 2; i < n; i++ ) {
1418  if ( GetModInverses(i,n,
1419  (WORD *)(&(AC.modinverses[i])),&inv2) ) {
1420  SETERROR(-1)
1421  }
1422  }
1423  }
1424  UNLOCK(AC.halfmodlock);
1425  }
1426  return(0);
1427 }
1428 
1429 /*
1430  #] MakeInverses :
1431  #[ GetModInverses :
1432 */
1443 int GetModInverses(WORD m1, WORD m2, WORD *im1, WORD *im2)
1444 {
1445  WORD a1, a2, a3;
1446  WORD b1, b2, b3;
1447  WORD x = m1, y, c, d = m2;
1448  if ( x < 1 || d <= 1 ) goto somethingwrong;
1449  a1 = 0; a2 = 1;
1450  b1 = 1; b2 = 0;
1451  for(;;) {
1452  c = d/x; y = d%x; /* a good compiler makes this faster than y=d-c*x */
1453  if ( y == 0 ) break;
1454  a3 = a1-c*a2; a1 = a2; a2 = a3;
1455  b3 = b1-c*b2; b1 = b2; b2 = b3;
1456  d = x; x = y;
1457  }
1458  if ( x != 1 ) goto somethingwrong;
1459  if ( a2 < 0 ) a2 += m2;
1460  if ( b2 < 0 ) b2 += m1;
1461  if (im1!=NULL) *im1 = a2;
1462  if (im2!=NULL) *im2 = b2;
1463  return(0);
1464 somethingwrong:
1465  MLOCK(ErrorMessageLock);
1466  MesPrint("Error trying to determine inverses in GetModInverses");
1467  MUNLOCK(ErrorMessageLock);
1468  return(-1);
1469 }
1470 /*
1471  #] GetModInverses :
1472  #[ GetLongModInverses :
1473 */
1474 
1475 int GetLongModInverses(PHEAD UWORD *a, WORD na, UWORD *b, WORD nb, UWORD *ia, WORD *nia, UWORD *ib, WORD *nib) {
1476 
1477  UWORD *s, *t, *sa, *sb, *ta, *tb, *x, *y, *swap1;
1478  WORD ns, nt, nsa, nsb, nta, ntb, nx, ny, swap2;
1479 
1480  s = NumberMalloc("GetLongModInverses");
1481  ns = na;
1482  WCOPY(s, a, ABS(ns));
1483 
1484  t = NumberMalloc("GetLongModInverses");
1485  nt = nb;
1486  WCOPY(t, b, ABS(nt));
1487 
1488  sa = NumberMalloc("GetLongModInverses");
1489  nsa = 1;
1490  sa[0] = 1;
1491 
1492  sb = NumberMalloc("GetLongModInverses");
1493  nsb = 0;
1494 
1495  ta = NumberMalloc("GetLongModInverses");
1496  nta = 0;
1497 
1498  tb = NumberMalloc("GetLongModInverses");
1499  ntb = 1;
1500  tb[0] = 1;
1501 
1502  x = NumberMalloc("GetLongModInverses");
1503  y = NumberMalloc("GetLongModInverses");
1504 
1505  while (nt != 0) {
1506  DivLong(s,ns,t,nt,x,&nx,y,&ny);
1507  swap1=s; s=y; y=swap1;
1508  ns=ny;
1509  MulLong(x,nx,ta,nta,y,&ny);
1510  AddLong(sa,nsa,y,-ny,sa,&nsa);
1511  MulLong(x,nx,tb,ntb,y,&ny);
1512  AddLong(sb,nsb,y,-ny,sb,&nsb);
1513 
1514  swap1=s; s=t; t=swap1;
1515  swap2=ns; ns=nt; nt=swap2;
1516  swap1=sa; sa=ta; ta=swap1;
1517  swap2=nsa; nsa=nta; nta=swap2;
1518  swap1=sb; sb=tb; tb=swap1;
1519  swap2=nsb; nsb=ntb; ntb=swap2;
1520  }
1521 
1522  if (ia!=NULL) {
1523  *nia = nsa*ns;
1524  WCOPY(ia,sa,ABS(*nia));
1525  }
1526 
1527  if (ib!=NULL) {
1528  *nib = nsb*ns;
1529  WCOPY(ib,sb,ABS(*nib));
1530  }
1531 
1532  NumberFree(s,"GetLongModInverses");
1533  NumberFree(t,"GetLongModInverses");
1534  NumberFree(sa,"GetLongModInverses");
1535  NumberFree(sb,"GetLongModInverses");
1536  NumberFree(ta,"GetLongModInverses");
1537  NumberFree(tb,"GetLongModInverses");
1538  NumberFree(x,"GetLongModInverses");
1539  NumberFree(y,"GetLongModInverses");
1540 
1541  return 0;
1542 }
1543 
1544 /*
1545  #] GetLongModInverses :
1546  #[ Product : WORD Product(a,na,b)
1547 
1548  Multiplies the Long number in a with the WORD b.
1549 
1550 */
1551 
1552 WORD Product(UWORD *a, WORD *na, WORD b)
1553 {
1554  WORD i, sgn = 1;
1555  RLONG t, u;
1556  if ( *na < 0 ) { *na = -(*na); sgn = -sgn; }
1557  if ( b < 0 ) { b = -b; sgn = -sgn; }
1558  t = 0;
1559  u = (RLONG)b;
1560  for ( i = 0; i < *na; i++ ) {
1561  t += *a * u;
1562  *a++ = (UWORD)t;
1563  t >>= BITSINWORD;
1564  }
1565  if ( t > 0 ) {
1566  if ( ++(*na) > AM.MaxTal ) {
1567  MLOCK(ErrorMessageLock);
1568  MesPrint("Overflow in Product");
1569  MUNLOCK(ErrorMessageLock);
1570  return(-1);
1571  }
1572  *a = (UWORD)t;
1573  }
1574  if ( sgn < 0 ) *na = -(*na);
1575  return(0);
1576 }
1577 
1578 /*
1579  #] Product :
1580  #[ Quotient : UWORD Quotient(a,na,b)
1581 
1582  Routine divides the long number a by b with the assumption that
1583  there is no remainder (like while computing binomials).
1584 
1585 */
1586 
1587 UWORD Quotient(UWORD *a, WORD *na, WORD b)
1588 {
1589  RLONG v, t;
1590  WORD i, j, sgn = 1;
1591  if ( ( i = *na ) < 0 ) { sgn = -1; i = -i; }
1592  if ( b < 0 ) { b = -b; sgn = -sgn; }
1593  if ( i == 1 ) {
1594  if ( ( *a /= (UWORD)b ) == 0 ) *na = 0;
1595  if ( sgn < 0 ) *na = -*na;
1596  return(0);
1597  }
1598  a += i;
1599  j = i;
1600  v = (RLONG)b;
1601  t = (RLONG)(*--a);
1602  while ( --i >= 0 ) {
1603  *a = t / v;
1604  t -= v * (*a);
1605  if ( i ) {
1606  t <<= BITSINWORD;
1607  t += *--a;
1608  }
1609  }
1610  a += j - 1;
1611  if ( !*a ) j--;
1612  if ( sgn < 0 ) j = -j;
1613  *na = j;
1614  return(0);
1615 }
1616 
1617 /*
1618  #] Quotient :
1619  #[ Remain10 : WORD Remain10(a,na)
1620 
1621  Routine devides a by 10 and gives the remainder as return value.
1622  The value of a will be the quotient! a must be positive.
1623 
1624 */
1625 
1626 WORD Remain10(UWORD *a, WORD *na)
1627 {
1628  WORD i;
1629  RLONG t, u;
1630  UWORD *b;
1631  i = *na;
1632  t = 0;
1633  b = a + i - 1;
1634  while ( --i >= 0 ) {
1635  t += *b;
1636  *b-- = u = t / 10;
1637  t -= u * 10;
1638  if ( i > 0 ) t <<= BITSINWORD;
1639  }
1640  if ( ( *na > 0 ) && !a[*na-1] ) (*na)--;
1641  return((WORD)t);
1642 }
1643 
1644 /*
1645  #] Remain10 :
1646  #[ Remain4 : WORD Remain4(a,na)
1647 
1648  Routine devides a by 10000 and gives the remainder as return value.
1649  The value of a will be the quotient! a must be positive.
1650 
1651 */
1652 
1653 WORD Remain4(UWORD *a, WORD *na)
1654 {
1655  WORD i;
1656  RLONG t, u;
1657  UWORD *b;
1658  i = *na;
1659  t = 0;
1660  b = a + i - 1;
1661  while ( --i >= 0 ) {
1662  t += *b;
1663  *b-- = u = t / 10000;
1664  t -= u * 10000;
1665  if ( i > 0 ) t <<= BITSINWORD;
1666  }
1667  if ( ( *na > 0 ) && !a[*na-1] ) (*na)--;
1668  return((WORD)t);
1669 }
1670 
1671 /*
1672  #] Remain4 :
1673  #[ PrtLong : VOID PrtLong(a,na,s)
1674 
1675  Puts the long number a in string s.
1676 
1677 */
1678 
1679 VOID PrtLong(UWORD *a, WORD na, UBYTE *s)
1680 {
1681  GETIDENTITY
1682  WORD q, i;
1683  UBYTE *sa, *sb;
1684  UBYTE c;
1685  UWORD *bb, *b;
1686 
1687  if ( na < 0 ) {
1688  *s++ = '-';
1689  na = -na;
1690  }
1691 
1692  b = NumberMalloc("PrtLong");
1693  bb = b;
1694  i = na; while ( --i >= 0 ) *bb++ = *a++;
1695  a = b;
1696  if ( na > 2 ) {
1697  sa = s;
1698  do {
1699  q = Remain4(a,&na);
1700  *sa++ = (UBYTE)('0' + (q%10));
1701  q /= 10;
1702  *sa++ = (UBYTE)('0' + (q%10));
1703  q /= 10;
1704  *sa++ = (UBYTE)('0' + (q%10));
1705  q /= 10;
1706  *sa++ = (UBYTE)('0' + (q%10));
1707  } while ( na );
1708  while ( sa[-1] == '0' ) sa--;
1709  sb = s;
1710  s = sa;
1711  sa--;
1712  while ( sa > sb ) { c = *sa; *sa = *sb; *sb = c; sa--; sb++; }
1713  }
1714  else if ( na ) {
1715  sa = s;
1716  do {
1717  q = Remain10(a,&na);
1718  *sa++ = (UBYTE)('0' + q);
1719  } while ( na );
1720  sb = s;
1721  s = sa;
1722  sa--;
1723  while ( sa > sb ) { c = *sa; *sa = *sb; *sb = c; sa--; sb++; }
1724  }
1725  else *s++ = '0';
1726  *s = '\0';
1727  NumberFree(b,"PrtLong");
1728 }
1729 
1730 /*
1731  #] PrtLong :
1732  #[ GetLong : WORD GetLong(s,a,na)
1733 
1734  Reads a long number from a string.
1735  The string is zero terminated and contains only digits!
1736 
1737  New algorithm: try to read 4 digits together before the result
1738  is accumulated.
1739 */
1740 
1741 WORD GetLong(UBYTE *s, UWORD *a, WORD *na)
1742 {
1743 /*
1744  UWORD digit;
1745  *a = 0;
1746  *na = 0;
1747  while ( FG.cTable[*s] == 1 ) {
1748  digit = *s++ - '0';
1749  if ( *na && Product(a,na,(WORD)10) ) return(-1);
1750  if ( digit && AddLong(a,*na,&digit,(WORD)1,a,na) ) return(-1);
1751  }
1752  return(0);
1753 */
1754  UWORD digit, x = 0, y = 0;
1755  *a = 0;
1756  *na = 0;
1757  while ( FG.cTable[*s] == 1 ) {
1758  x = *s++ - '0';
1759  if ( FG.cTable[*s] != 1 ) { y = 10; break; }
1760  x = 10*x + *s++ - '0';
1761  if ( FG.cTable[*s] != 1 ) { y = 100; break; }
1762  x = 10*x + *s++ - '0';
1763  if ( FG.cTable[*s] != 1 ) { y = 1000; break; }
1764  x = 10*x + *s++ - '0';
1765  if ( *na && Product(a,na,(WORD)10000) ) return(-1);
1766  if ( ( digit = x ) != 0 && AddLong(a,*na,&digit,(WORD)1,a,na) )
1767  return(-1);
1768  y = 0;
1769  }
1770  if ( y ) {
1771  if ( *na && Product(a,na,(WORD)y) ) return(-1);
1772  if ( ( digit = x ) != 0 && AddLong(a,*na,&digit,(WORD)1,a,na) )
1773  return(-1);
1774  }
1775  return(0);
1776 }
1777 
1778 /*
1779  #] GetLong :
1780  #[ GCD : WORD GCD(a,na,b,nb,c,nc)
1781 
1782  Algorithm to compute the GCD of two long numbers.
1783  See Knuth, sec 4.5.2 algorithm L.
1784 
1785  We assume that both numbers are positive
1786 */
1787 
1788 #ifdef EXTRAGCD
1789 
1790 #define Convert(ia,aa,naa) \
1791  if ( (LONG)ia < 0 ) { \
1792  ia = (ULONG)(-(LONG)ia); \
1793  aa[0] = ia; \
1794  if ( ( aa[1] = ia >> BITSINWORD ) != 0 ) naa = -2; \
1795  else naa = -1; \
1796  } \
1797  else if ( ia == 0 ) { aa[0] = 0; naa = 0; } \
1798  else { \
1799  aa[0] = ia; \
1800  if ( ( aa[1] = ia >> BITSINWORD ) != 0 ) naa = 2; \
1801  else naa = 1; \
1802  }
1803 
1804 VOID GCD(UWORD *a, WORD na, UWORD *b, WORD nb, UWORD *c, WORD *nc)
1805 {
1806  int ja = 0, jb = 0, j;
1807  UWORD *r,*t;
1808  UWORD *x1, *x2, *x3;
1809  WORD nd,naa,nbb;
1810  ULONG ia,ib,ic,id,u,v,w,q,T;
1811  UWORD aa[2], bb[2];
1812 /*
1813  First eliminate easy powers of 2^...
1814 */
1815  while ( a[0] == 0 ) { na--; ja++; a++; }
1816  while ( b[0] == 0 ) { nb--; jb++; b++; }
1817  if ( ja > jb ) ja = jb;
1818  if ( ja > 0 ) {
1819  j = ja;
1820  do { *c++ = 0; } while ( --j > 0 );
1821  }
1822 /*
1823  Now arrange things such that a >= b
1824 */
1825  if ( na < nb ) {
1826  jb = na; na = nb; nb = jb;
1827 exch:
1828  r = a; a = b; b = r;
1829  }
1830  else if ( na == nb ) {
1831  r = a+na;
1832  t = b+nb;
1833  j = na;
1834  while ( --j >= 0 ) {
1835  if ( *--r > *--t ) break;
1836  if ( *r < *t ) goto exch;
1837  }
1838  if ( j < 0 ) {
1839 out:
1840  j = nb;
1841  NCOPY(c,b,j);
1842  *nc = nb+ja;
1843  return;
1844  }
1845  }
1846 /*
1847  {
1848  MLOCK(ErrorMessageLock);
1849  MesPrint("Ordered input, ja = %d",(WORD)ja);
1850  AO.OutSkip = 3;
1851  FiniLine();
1852  j = na; r = a;
1853  while ( --j >= 0 ) { TalToLine((UWORD)(*r++)); TokenToLine((UBYTE *)" "); }
1854  FiniLine();
1855  j = nb; r = b;
1856  while ( --j >= 0 ) { TalToLine((UWORD)(*r++)); TokenToLine((UBYTE *)" "); }
1857  AO.OutSkip = 0;
1858  FiniLine();
1859  MUNLOCK(ErrorMessageLock);
1860  }
1861 */
1862 /*
1863  We have now that A > B
1864  The loop recognizes the case that na-nb >= 1
1865  In that case we just have to divide!
1866 */
1867  r = x1 = NumberMalloc("GCD"); t = x2 = NumberMalloc("GCD"); x3 = NumberMalloc("GCD");
1868  j = na;
1869  NCOPY(r,a,j);
1870  j = nb;
1871  NCOPY(t,b,j);
1872 
1873  for(;;) {
1874 
1875  while ( na > nb ) {
1876 toobad:
1877  DivLong(x1,na,x2,nb,c,nc,x3,&nd);
1878  if ( nd == 0 ) { b = x2; goto out; }
1879  t = x1; x1 = x2; x2 = x3; x3 = t; na = nb; nb = nd;
1880  if ( na == 2 ) break;
1881  }
1882 /*
1883  Here we can use the shortcut.
1884 */
1885  if ( na == 2 ) {
1886  v = x1[0] + ( ((ULONG)x1[1]) << BITSINWORD );
1887  w = x2[0];
1888  if ( nb == 2 ) w += ((ULONG)x2[1]) << BITSINWORD;
1889 #ifdef EXTRAGCD2
1890  v = GCD2(v,w);
1891 #else
1892  do { u = v%w; v = w; w = u; } while ( w );
1893 #endif
1894  c[0] = (UWORD)v;
1895  if ( ( c[1] = (UWORD)(v >> BITSINWORD) ) != 0 ) *nc = 2+ja;
1896  else *nc = 1+ja;
1897  NumberFree(x1,"GCD"); NumberFree(x2,"GCD"); NumberFree(x3,"GCD");
1898  return;
1899  }
1900  if ( na == 1 ) {
1901  UWORD ui, uj;
1902  ui = x1[0]; uj = x2[0];
1903 #ifdef EXTRAGCD2
1904  ui = (UWORD)GCD2((ULONG)ui,(ULONG)uj);
1905 #else
1906  do { nd = ui%uj; ui = uj; uj = nd; } while ( nd );
1907 #endif
1908  c[0] = ui;
1909  *nc = 1 + ja;
1910  NumberFree(x1,"GCD"); NumberFree(x2,"GCD"); NumberFree(x3,"GCD");
1911  return;
1912  }
1913  ia = 1; ib = 0; ic = 0; id = 1;
1914  u = ( ((ULONG)x1[na-1]) << BITSINWORD ) + x1[na-2];
1915  v = ( ((ULONG)x2[nb-1]) << BITSINWORD ) + x2[nb-2];
1916 
1917  while ( v+ic != 0 && v+id != 0 &&
1918  ( q = (u+ia)/(v+ic) ) == (u+ib)/(v+id) ) {
1919  T = ia-q*ic; ia = ic; ic = T;
1920  T = ib-q*id; ib = id; id = T;
1921  T = u - q*v; u = v; v = T;
1922  }
1923  if ( ib == 0 ) goto toobad;
1924  Convert(ia,aa,naa);
1925  Convert(ib,bb,nbb);
1926  MulLong(x1,na,aa,naa,x3,&nd);
1927  MulLong(x2,nb,bb,nbb,c,nc);
1928  AddLong(x3,nd,c,*nc,c,nc);
1929  Convert(ic,aa,naa);
1930  Convert(id,bb,nbb);
1931  MulLong(x1,na,aa,naa,x3,&nd);
1932  t = c; na = j = *nc; r = x1;
1933  NCOPY(r,t,j);
1934  MulLong(x2,nb,bb,nbb,c,nc);
1935  AddLong(x3,nd,c,*nc,x2,&nb);
1936  }
1937 }
1938 
1939 #endif
1940 
1941 /*
1942  #] GCD :
1943  #[ GcdLong : WORD GcdLong(a,na,b,nb,c,nc)
1944 
1945  Returns the Greatest Common Divider of a and b in c.
1946  If a and or b are zero an error message will be returned.
1947  The answer is always positive.
1948  In principle a and c can be the same.
1949 */
1950 
1951 #ifndef NEWTRICK
1952 /*
1953  #[ Old Routine :
1954 */
1955 
1956 WORD GcdLong(PHEAD UWORD *a, WORD na, UWORD *b, WORD nb, UWORD *c, WORD *nc)
1957 {
1958  GETBIDENTITY
1959  if ( !na || !nb ) {
1960  if ( !na && !nb ) {
1961  MLOCK(ErrorMessageLock);
1962  MesPrint("Cannot take gcd");
1963  MUNLOCK(ErrorMessageLock);
1964  return(-1);
1965  }
1966 
1967  if ( !na ) {
1968  *nc = abs(nb);
1969  NCOPY(c,b,*nc);
1970  *nc = abs(nb);
1971  return(0);
1972  }
1973 
1974  *nc = abs(na);
1975  NCOPY(c,a,*nc);
1976  *nc = abs(na);
1977  return(0);
1978  }
1979  if ( na < 0 ) na = -na;
1980  if ( nb < 0 ) nb = -nb;
1981  if ( na == 1 && nb == 1 ) {
1982 #ifdef EXTRAGCD2
1983  *c = (UWORD)GCD2((ULONG)*a,(ULONG)*b);
1984 #else
1985  UWORD x,y,z;
1986  x = *a;
1987  y = *b;
1988  do { z = x % y; x = y; } while ( ( y = z ) != 0 );
1989  *c = x;
1990 #endif
1991  *nc = 1;
1992  }
1993  else if ( na <= 2 && nb <= 2 ) {
1994  RLONG lx,ly,lz;
1995  if ( na == 2 ) { lx = (((RLONG)(a[1]))<<BITSINWORD) + *a; }
1996  else { lx = *a; }
1997  if ( nb == 2 ) { ly = (((RLONG)(b[1]))<<BITSINWORD) + *b; }
1998  else { ly = *b; }
1999 #ifdef LONGWORD
2000 #ifdef EXTRAGCD2
2001  lx = GCD2(lx,ly);
2002 #else
2003  do { lz = lx % ly; lx = ly; } while ( ( ly = lz ) != 0 );
2004 #endif
2005 #else
2006  if ( lx < ly ) { lz = lx; lx = ly; ly = lz; }
2007  do {
2008  lz = lx % ly; lx = ly;
2009  } while ( ( ly = lz ) != 0 && ( lx & AWORDMASK ) != 0 );
2010  if ( ly ) {
2011  do { *c = ((UWORD)lx)%((UWORD)ly); lx = ly; } while ( ( ly = *c ) != 0 );
2012  *c = (UWORD)lx;
2013  *nc = 1;
2014  }
2015  else
2016 #endif
2017  {
2018  *c++ = (UWORD)lx;
2019  if ( ( *c = (UWORD)(lx >> BITSINWORD) ) != 0 ) *nc = 2;
2020  else *nc = 1;
2021  }
2022  }
2023  else {
2024 #ifdef EXTRAGCD
2025  GCD(a,na,b,nb,c,nc);
2026 #else
2027 #ifdef NEWGCD
2028  UWORD *x3,*x1,*x2, *GLscrat7, *GLscrat8;
2029  WORD n1,n2,n3,n4;
2030  WORD i, j;
2031  x1 = c; x3 = a; n1 = i = na;
2032  NCOPY(x1,x3,i);
2033  GLscrat7 = NumberMalloc("GcdLong"); GLscrat8 = NumberMalloc("GcdLong");
2034  x2 = GLscrat8; x3 = b; n2 = i = nb;
2035  NCOPY(x2,x3,i);
2036  x1 = c; i = 0;
2037  while ( x1[0] == 0 ) { i += BITSINWORD; x1++; n1--; }
2038  while ( ( x1[0] & 1 ) == 0 ) { i++; SCHUIF(x1,n1) }
2039  x2 = GLscrat8; j = 0;
2040  while ( x2[0] == 0 ) { j += BITSINWORD; x2++; n2--; }
2041  while ( ( x2[0] & 1 ) == 0 ) { j++; SCHUIF(x2,n2) }
2042  if ( j > i ) j = i; /* powers of two in GCD */
2043  for(;;){
2044  if ( n1 > n2 ) {
2045 firstbig:
2046  SubPLon(x1,n1,x2,n2,x1,&n3);
2047  n1 = n3;
2048  if ( n1 == 0 ) {
2049  x1 = c;
2050  n1 = i = n2; NCOPY(x1,x2,i);
2051  break;
2052  }
2053  while ( ( x1[0] & 1 ) == 0 ) SCHUIF(x1,n1)
2054  if ( n1 == 1 ) {
2055  if ( DivLong(x2,n2,x1,n1,GLscrat7,&n3,x2,&n4) ) goto GcdErr;
2056  n2 = n4;
2057  if ( n2 == 0 ) {
2058  i = n1; x2 = c; NCOPY(x2,x1,i);
2059  break;
2060  }
2061 #ifdef EXTRAGCD2
2062  *c = (UWORD)GCD2((ULONG)x1[0],(ULONG)x2[0]);
2063 #else
2064  {
2065  UWORD x,y,z;
2066  x = x1[0];
2067  y = x2[0];
2068  do { z = x % y; x = y; } while ( ( y = z ) != 0 );
2069  *c = x;
2070  }
2071 #endif
2072  n1 = 1;
2073  break;
2074  }
2075  }
2076  else if ( n1 < n2 ) {
2077 lastbig:
2078  SubPLon(x2,n2,x1,n1,x2,&n3);
2079  n2 = n3;
2080  if ( n2 == 0 ) {
2081  i = n1; x2 = c; NCOPY(x2,x1,i);
2082  break;
2083  }
2084  while ( ( x2[0] & 1 ) == 0 ) SCHUIF(x2,n2)
2085  if ( n2 == 1 ) {
2086  if ( DivLong(x1,n1,x2,n2,GLscrat7,&n3,x1,&n4) ) goto GcdErr;
2087  n1 = n4;
2088  if ( n1 == 0 ) {
2089  x1 = c;
2090  n1 = i = n2; NCOPY(x1,x2,i);
2091  break;
2092  }
2093 #ifdef EXTRAGCD2
2094  *c = (UWORD)GCD2((ULONG)x2[0],(ULONG)x1[0]);
2095 #else
2096  {
2097  UWORD x,y,z;
2098  x = x2[0];
2099  y = x1[0];
2100  do { z = x % y; x = y; } while ( ( y = z ) != 0 );
2101  *c = x;
2102  }
2103 #endif
2104  n1 = 1;
2105  break;
2106  }
2107  }
2108  else {
2109  for ( i = n1-1; i >= 0; i-- ) {
2110  if ( x1[i] > x2[i] ) goto firstbig;
2111  else if ( x1[i] < x2[i] ) goto lastbig;
2112  }
2113  i = n1; x2 = c; NCOPY(x2,x1,i);
2114  break;
2115  }
2116  }
2117 /*
2118  Now the GCD is in c but still needs j powers of 2.
2119 */
2120  x1 = c;
2121  while ( j >= BITSINWORD ) {
2122  for ( i = n1; i > 0; i-- ) x1[i] = x1[i-1];
2123  x1[0] = 0; n1++;
2124  j -= BITSINWORD;
2125  }
2126  if ( j > 0 ) {
2127  ULONG a1,a2 = 0;
2128  for ( i = 0; i < n1; i++ ) {
2129  a1 = x1[i]; a1 <<= j;
2130  a2 += a1;
2131  x1[i] = a2;
2132  a2 >>= BITSINWORD;
2133  }
2134  if ( a2 != 0 ) {
2135  x1[n1++] = a2;
2136  }
2137  }
2138  *nc = n1;
2139  NumberFree(GLscrat7,"GcdLong"); NumberFree(GLscrat8,"GcdLong");
2140 #else
2141  UWORD *x1,*x2,*x3,*x4,*c1,*c2;
2142  WORD n1,n2,n3,n4,i;
2143  x1 = c; x3 = a; n1 = i = na;
2144  NCOPY(x1,x3,i);
2145  x1 = c; c1 = x2 = NumberMalloc("GcdLong"); x3 = NumberMalloc("GcdLong"); x4 = NumberMalloc("GcdLong");
2146  c2 = b; n2 = i = nb;
2147  NCOPY(c1,c2,i);
2148  for(;;){
2149  if ( DivLong(x1,n1,x2,n2,x4,&n4,x3,&n3) ) goto GcdErr;
2150  if ( !n3 ) { x1 = x2; n1 = n2; break; }
2151  if ( DivLong(x2,n2,x3,n3,x4,&n4,x1,&n1) ) goto GcdErr;
2152  if ( !n1 ) { x1 = x3; n1 = n3; break; }
2153  if ( DivLong(x3,n3,x1,n1,x4,&n4,x2,&n2) ) goto GcdErr;
2154  if ( !n2 ) {
2155  *nc = n1;
2156  NumberFree(x2,"GcdLong"); NumberFree(x3,"GcdLong"); NumberFree(x4,"GcdLong");
2157  return(0);
2158  }
2159  }
2160  *nc = i = n1;
2161  NCOPY(c,x1,i);
2162  NumberFree(x2,"GcdLong"); NumberFree(x3,"GcdLong"); NumberFree(x4,"GcdLong");
2163 #endif
2164 #endif
2165  }
2166  return(0);
2167 GcdErr:
2168  MLOCK(ErrorMessageLock);
2169  MesCall("GcdLong");
2170  MUNLOCK(ErrorMessageLock);
2171  SETERROR(-1)
2172 }
2173 /*
2174  #] Old Routine :
2175 */
2176 #else
2177 
2178 /*
2179  New routine for GcdLong that uses smart shortcuts.
2180  Algorithm by J. Vermaseren 15-nov-2006.
2181  It runs faster for very big numbers but only by a fixed factor.
2182  There is no improvement in the power behaviour.
2183  Improvement on the whole of hf9 (multiple zeta values at weight 9):
2184  Better than a factor 2 on a 32 bits architecture and 2.76 on a
2185  64 bits architecture.
2186  On hf10 (MZV's at weight 10), 64 bits architecture: factor 7.
2187 
2188  If we have two long numbers (na,nb > GCDMAX) we will work in a
2189  truncated way. At the moment of writing (15-nov-2006) it isn't
2190  clear whether this algorithm is an invention or a reinvention.
2191  A short search on the web didn't show anything.
2192 
2193  31-jul-2007:
2194  A better search shows that this is an adaptation of the Lehmer-Euclid
2195  algorithm, already described in Knuth. Here we can work without upper
2196  and lower limit because we are only interested in the GCD, not the
2197  extra numbers. Also it takes already some features of the double
2198  digit Lehmer-Euclid algorithm of Jebelean it seems.
2199 
2200  Maybe this can be programmed slightly better and we can get another
2201  few percent speed increase. Further improvements for the assymptotic
2202  case come from splitting the calculation as in Karatsuba and working
2203  with FFT divisions and multiplications etc. But this is when hundreds
2204  of words are involved at the least.
2205 
2206  Algorithm
2207 
2208  1: while ( na > nb || nb < GCDMAX ) {
2209  if ( nb == 0 ) { result in a }
2210  c = a % b;
2211  a = b;
2212  b = c;
2213  }
2214  2: Make the truncated values in which a and b are the combinations
2215  of the top two words of a and b. The whole numbers are aa and bb now.
2216  3: ma1 = 1; ma2 = 0; mb1 = 0; mb2 = 1;
2217  4: A = a; B = b; m = a/b; c = a - m*b;
2218  c = ma1*a+ma2*b-m*(mb1*a+mb2*b) = (ma1-m*mb1)*a+(ma2-m*mb2)*b
2219  mc1 = ma1-m*mb1; mc2 = ma2-m*mb2;
2220  5: a = b; ma1 = mb1; ma2 = mb2;
2221  b = c; mb1 = mc1; mb2 = mc2;
2222  6: if ( b != 0 && nb >= FULLMAX ) goto 4;
2223  7: Now construct the new quantities
2224  ma1*aa+ma2*bb and mb1*aa+mb2*bb
2225  8: goto 1;
2226 
2227  The essence of the above algorithm is that we do the divisions only
2228  on relatively short numbers. Also usually there are many steps 4&5
2229  for each step 7. This eliminates many operations.
2230  The termination at FULLMAX is that we make errors by not considering
2231  the tail of the number. If we run b down all the way, the errors combine
2232  in such a way that the new numbers may be of the same order as the old
2233  numbers. By stopping halfway we don't get the error beyond halfway
2234  either. Unfortunately this means that a >= FULLMAX and hence na > nb
2235  which means that next we will have a complete division. But just once.
2236  Running the steps 4-6 till a < FULLMAX runs already into problems.
2237  It may be necessary to experiment a bit to obtain the optimum value
2238  of GCDMAX.
2239 */
2240 
2241 WORD GcdLong(PHEAD UWORD *a, WORD na, UWORD *b, WORD nb, UWORD *c, WORD *nc)
2242 {
2243  GETBIDENTITY
2244  UWORD x,y,z;
2245  UWORD *x1,*x2,*x3,*x4,*x5,*d;
2246  UWORD *GLscrat6, *GLscrat7, *GLscrat8, *GLscrat9, *GLscrat10;
2247  WORD n1,n2,n3,n4,n5,i;
2248  RLONG lx,ly,lz;
2249  LONG ma1, ma2, mb1, mb2, mc1, mc2, m;
2250  if ( !na || !nb ) {
2251  if ( !na && !nb ) {
2252  MLOCK(ErrorMessageLock);
2253  MesPrint("Cannot take gcd");
2254  MUNLOCK(ErrorMessageLock);
2255  return(-1);
2256  }
2257 
2258  if ( !na ) {
2259  *nc = abs(nb);
2260  NCOPY(c,b,*nc);
2261  *nc = abs(nb);
2262  return(0);
2263  }
2264 
2265  *nc = abs(na);
2266  NCOPY(c,a,*nc);
2267  *nc = abs(na);
2268  return(0);
2269  }
2270  if ( na < 0 ) na = -na;
2271  if ( nb < 0 ) nb = -nb;
2272 /*
2273  #[ GMP stuff :
2274 */
2275 #ifdef WITHGMP
2276  if ( na > 3 && nb > 3 ) {
2277  int ii;
2278  mp_limb_t *upa, *upb, *upc, xx;
2279  UWORD *uw, *u1, *u2;
2280  unsigned int tcounta, tcountb, tcounta1, tcountb1;
2281  mp_size_t ana, anb, anc;
2282 
2283  u1 = uw = NumberMalloc("GcdLong");
2284  upa = (mp_limb_t *)u1;
2285  ana = na; tcounta1 = 0;
2286  while ( a[0] == 0 ) { a++; ana--; tcounta1++; }
2287  for ( ii = 0; ii < ana; ii++ ) { *uw++ = *a++; }
2288  if ( ( ana & 1 ) != 0 ) { *uw = 0; ana++; }
2289  ana >>= 1;
2290 
2291  u2 = uw = NumberMalloc("GcdLong");
2292  upb = (mp_limb_t *)u2;
2293  anb = nb; tcountb1 = 0;
2294  while ( b[0] == 0 ) { b++; anb--; tcountb1++; }
2295  for ( ii = 0; ii < anb; ii++ ) { *uw++ = *b++; }
2296  if ( ( anb & 1 ) != 0 ) { *uw = 0; anb++; }
2297  anb >>= 1;
2298 
2299  xx = upa[0]; tcounta = 0;
2300  while ( ( xx & 15 ) == 0 ) { tcounta += 4; xx >>= 4; }
2301  while ( ( xx & 1 ) == 0 ) { tcounta += 1; xx >>= 1; }
2302  xx = upb[0]; tcountb = 0;
2303  while ( ( xx & 15 ) == 0 ) { tcountb += 4; xx >>= 4; }
2304  while ( ( xx & 1 ) == 0 ) { tcountb += 1; xx >>= 1; }
2305 
2306  if ( tcounta ) {
2307  mpn_rshift(upa,upa,ana,tcounta);
2308  if ( upa[ana-1] == 0 ) ana--;
2309  }
2310  if ( tcountb ) {
2311  mpn_rshift(upb,upb,anb,tcountb);
2312  if ( upb[anb-1] == 0 ) anb--;
2313  }
2314 
2315  upc = (mp_limb_t *)(NumberMalloc("GcdLong"));
2316  if ( ( ana > anb ) || ( ( ana == anb ) && ( upa[ana-1] >= upb[ana-1] ) ) ) {
2317  anc = mpn_gcd(upc,upa,ana,upb,anb);
2318  }
2319  else {
2320  anc = mpn_gcd(upc,upb,anb,upa,ana);
2321  }
2322 
2323  tcounta = tcounta1*BITSINWORD + tcounta;
2324  tcountb = tcountb1*BITSINWORD + tcountb;
2325  if ( tcountb > tcounta ) tcountb = tcounta;
2326  tcounta = tcountb/BITSINWORD;
2327  tcountb = tcountb%BITSINWORD;
2328 
2329  if ( tcountb ) {
2330  xx = mpn_lshift(upc,upc,anc,tcountb);
2331  if ( xx ) { upc[anc] = xx; anc++; }
2332  }
2333 
2334  uw = (UWORD *)upc; anc *= 2;
2335  while ( uw[anc-1] == 0 ) anc--;
2336  for ( ii = 0; ii < (int)tcounta; ii++ ) *c++ = 0;
2337  for ( ii = 0; ii < anc; ii++ ) *c++ = *uw++;
2338  *nc = anc + tcounta;
2339  NumberFree(u1,"GcdLong"); NumberFree(u2,"GcdLong"); NumberFree((UWORD *)(upc),"GcdLong");
2340  return(0);
2341  }
2342 #endif
2343 /*
2344  #] GMP stuff :
2345 */
2346 /*
2347  #[ Easy cases :
2348 */
2349  if ( na == 1 && nb == 1 ) {
2350  x = *a;
2351  y = *b;
2352  do { z = x % y; x = y; } while ( ( y = z ) != 0 );
2353  *c = x;
2354  *nc = 1;
2355  return(0);
2356  }
2357  else if ( na <= 2 && nb <= 2 ) {
2358  if ( na == 2 ) { lx = (((RLONG)(a[1]))<<BITSINWORD) + *a; }
2359  else { lx = *a; }
2360  if ( nb == 2 ) { ly = (((RLONG)(b[1]))<<BITSINWORD) + *b; }
2361  else { ly = *b; }
2362  if ( lx < ly ) { lz = lx; lx = ly; ly = lz; }
2363 #if ( BITSINWORD == 16 )
2364  do {
2365  lz = lx % ly; lx = ly;
2366  } while ( ( ly = lz ) != 0 );
2367 #else
2368  do {
2369  lz = lx % ly; lx = ly;
2370  } while ( ( ly = lz ) != 0 && ( lx & AWORDMASK ) != 0 );
2371  if ( ly ) {
2372  x = (UWORD)lx; y = (UWORD)ly;
2373  do { *c = x % y; x = y; } while ( ( y = *c ) != 0 );
2374  *c = x;
2375  *nc = 1;
2376  }
2377  else
2378 #endif
2379  {
2380  *c++ = (UWORD)lx;
2381  if ( ( *c = (UWORD)(lx >> BITSINWORD) ) != 0 ) *nc = 2;
2382  else *nc = 1;
2383  }
2384  return(0);
2385  }
2386 /*
2387  #] Easy cases :
2388 */
2389  GLscrat6 = NumberMalloc("GcdLong"); GLscrat7 = NumberMalloc("GcdLong");
2390  GLscrat8 = NumberMalloc("GcdLong");
2391  GLscrat9 = NumberMalloc("GcdLong"); GLscrat10 = NumberMalloc("GcdLong");
2392 restart:;
2393 /*
2394  #[ Easy cases :
2395 */
2396  if ( na == 1 && nb == 1 ) {
2397  x = *a;
2398  y = *b;
2399  do { z = x % y; x = y; } while ( ( y = z ) != 0 );
2400  *c = x;
2401  *nc = 1;
2402  }
2403  else if ( na <= 2 && nb <= 2 ) {
2404  if ( na == 2 ) { lx = (((RLONG)(a[1]))<<BITSINWORD) + *a; }
2405  else { lx = *a; }
2406  if ( nb == 2 ) { ly = (((RLONG)(b[1]))<<BITSINWORD) + *b; }
2407  else { ly = *b; }
2408  if ( lx < ly ) { lz = lx; lx = ly; ly = lz; }
2409 #if ( BITSINWORD == 16 )
2410  do {
2411  lz = lx % ly; lx = ly;
2412  } while ( ( ly = lz ) != 0 );
2413 #else
2414  do {
2415  lz = lx % ly; lx = ly;
2416  } while ( ( ly = lz ) != 0 && ( lx & AWORDMASK ) != 0 );
2417  if ( ly ) {
2418  x = (UWORD)lx; y = (UWORD)ly;
2419  do { *c = x % y; x = y; } while ( ( y = *c ) != 0 );
2420  *c = x;
2421  *nc = 1;
2422  }
2423  else
2424 #endif
2425  {
2426  *c++ = (UWORD)lx;
2427  if ( ( *c = (UWORD)(lx >> BITSINWORD) ) != 0 ) *nc = 2;
2428  else *nc = 1;
2429  }
2430  }
2431 /*
2432  #] Easy cases :
2433  #[ Original code :
2434 */
2435  else if ( na < GCDMAX || nb < GCDMAX || na != nb ) {
2436  if ( na < nb ) {
2437  x2 = GLscrat8; x3 = a; n2 = i = na;
2438  NCOPY(x2,x3,i);
2439  x1 = c; x3 = b; n1 = i = nb;
2440  NCOPY(x1,x3,i);
2441  }
2442  else {
2443  x1 = c; x3 = a; n1 = i = na;
2444  NCOPY(x1,x3,i);
2445  x2 = GLscrat8; x3 = b; n2 = i = nb;
2446  NCOPY(x2,x3,i);
2447  }
2448  x1 = c; x2 = GLscrat8; x3 = GLscrat7; x4 = GLscrat6;
2449  for(;;){
2450  if ( DivLong(x1,n1,x2,n2,x4,&n4,x3,&n3) ) goto GcdErr;
2451  if ( !n3 ) { x1 = x2; n1 = n2; break; }
2452  if ( n2 <= 2 ) { a = x2; b = x3; na = n2; nb = n3; goto restart; }
2453  if ( n3 >= GCDMAX && n2 == n3 ) {
2454  a = GLscrat9; b = GLscrat10; na = n2; nb = n3;
2455  for ( i = 0; i < na; i++ ) a[i] = x2[i];
2456  for ( i = 0; i < nb; i++ ) b[i] = x3[i];
2457  goto newtrick;
2458  }
2459  if ( DivLong(x2,n2,x3,n3,x4,&n4,x1,&n1) ) goto GcdErr;
2460  if ( !n1 ) { x1 = x3; n1 = n3; break; }
2461  if ( n3 <= 2 ) { a = x3; b = x1; na = n3; nb = n1; goto restart; }
2462  if ( n1 >= GCDMAX && n1 == n3 ) {
2463  a = GLscrat9; b = GLscrat10; na = n3; nb = n1;
2464  for ( i = 0; i < na; i++ ) a[i] = x3[i];
2465  for ( i = 0; i < nb; i++ ) b[i] = x1[i];
2466  goto newtrick;
2467  }
2468  if ( DivLong(x3,n3,x1,n1,x4,&n4,x2,&n2) ) goto GcdErr;
2469  if ( !n2 ) { *nc = n1; goto normalend; }
2470  if ( n1 <= 2 ) { a = x1; b = x2; na = n1; nb = n2; goto restart; }
2471  if ( n2 >= GCDMAX && n2 == n1 ) {
2472  a = GLscrat9; b = GLscrat10; na = n1; nb = n2;
2473  for ( i = 0; i < na; i++ ) a[i] = x1[i];
2474  for ( i = 0; i < nb; i++ ) b[i] = x2[i];
2475  goto newtrick;
2476  }
2477  }
2478  *nc = i = n1;
2479  NCOPY(c,x1,i);
2480  }
2481 /*
2482  #] Original code :
2483  #[ New code :
2484 */
2485  else {
2486 /*
2487  This is the new algorithm starting at step 3.
2488 
2489  3: ma1 = 1; ma2 = 0; mb1 = 0; mb2 = 1;
2490  4: A = a; B = b; m = a/b; c = a - m*b;
2491  c = ma1*a+ma2*b-m*(mb1*a+mb2*b) = (ma1-m*mb1)*a+(ma2-m*mb2)*b
2492  mc1 = ma1-m*mb1; mc2 = ma2-m*mb2;
2493  5: a = b; ma1 = mb1; ma2 = mb2;
2494  b = c; mb1 = mc1; mb2 = mc2;
2495  6: if ( b != 0 ) goto 4;
2496 */
2497 newtrick:;
2498  ma1 = 1; ma2 = 0; mb1 = 0; mb2 = 1;
2499  lx = (((RLONG)(a[na-1]))<<BITSINWORD) + a[na-2];
2500  ly = (((RLONG)(b[nb-1]))<<BITSINWORD) + b[nb-2];
2501  if ( ly > lx ) { lz = lx; lx = ly; ly = lz; d = a; a = b; b = d; }
2502  do {
2503  m = lx/ly;
2504  mc1 = ma1-m*mb1; mc2 = ma2-m*mb2;
2505  ma1 = mb1; ma2 = mb2; mb1 = mc1; mb2 = mc2;
2506  lz = lx - m*ly; lx = ly; ly = lz;
2507  } while ( ly >= FULLMAX );
2508 /*
2509  Next the construction of the two new numbers
2510 
2511  7: Now construct the new quantities
2512  a = ma1*aa+ma2*bb and b = mb1*aa+mb2*bb
2513 */
2514  x1 = GLscrat6;
2515  x2 = GLscrat7;
2516  x3 = GLscrat8;
2517  x5 = GLscrat10;
2518  if ( ma1 < 0 ) {
2519  ma1 = -ma1;
2520  x1[0] = (UWORD)ma1;
2521  x1[1] = (UWORD)(ma1 >> BITSINWORD);
2522  if ( x1[1] ) n1 = -2;
2523  else n1 = -1;
2524  }
2525  else {
2526  x1[0] = (UWORD)ma1;
2527  x1[1] = (UWORD)(ma1 >> BITSINWORD);
2528  if ( x1[1] ) n1 = 2;
2529  else n1 = 1;
2530  }
2531  if ( MulLong(a,na,x1,n1,x2,&n2) ) goto GcdErr;
2532  if ( ma2 < 0 ) {
2533  ma2 = -ma2;
2534  x1[0] = (UWORD)ma2;
2535  x1[1] = (UWORD)(ma2 >> BITSINWORD);
2536  if ( x1[1] ) n1 = -2;
2537  else n1 = -1;
2538  }
2539  else {
2540  x1[0] = (UWORD)ma2;
2541  x1[1] = (UWORD)(ma2 >> BITSINWORD);
2542  if ( x1[1] ) n1 = 2;
2543  else n1 = 1;
2544  }
2545  if ( MulLong(b,nb,x1,n1,x3,&n3) ) goto GcdErr;
2546  if ( AddLong(x2,n2,x3,n3,c,&n4) ) goto GcdErr;
2547  if ( mb1 < 0 ) {
2548  mb1 = -mb1;
2549  x1[0] = (UWORD)mb1;
2550  x1[1] = (UWORD)(mb1 >> BITSINWORD);
2551  if ( x1[1] ) n1 = -2;
2552  else n1 = -1;
2553  }
2554  else {
2555  x1[0] = (UWORD)mb1;
2556  x1[1] = (UWORD)(mb1 >> BITSINWORD);
2557  if ( x1[1] ) n1 = 2;
2558  else n1 = 1;
2559  }
2560  if ( MulLong(a,na,x1,n1,x2,&n2) ) goto GcdErr;
2561  if ( mb2 < 0 ) {
2562  mb2 = -mb2;
2563  x1[0] = (UWORD)mb2;
2564  x1[1] = (UWORD)(mb2 >> BITSINWORD);
2565  if ( x1[1] ) n1 = -2;
2566  else n1 = -1;
2567  }
2568  else {
2569  x1[0] = (UWORD)mb2;
2570  x1[1] = (UWORD)(mb2 >> BITSINWORD);
2571  if ( x1[1] ) n1 = 2;
2572  else n1 = 1;
2573  }
2574  if ( MulLong(b,nb,x1,n1,x3,&n3) ) goto GcdErr;
2575  if ( AddLong(x2,n2,x3,n3,x5,&n5) ) goto GcdErr;
2576  a = c; na = n4; b = x5; nb = n5;
2577  if ( nb == 0 ) { *nc = n4; goto normalend; }
2578  x4 = GLscrat9;
2579  for ( i = 0; i < na; i++ ) x4[i] = a[i];
2580  a = x4;
2581  if ( na < 0 ) na = -na;
2582  if ( nb < 0 ) nb = -nb;
2583 /*
2584  The typical case now is that in a we have the last step to go
2585  to loose the leading word, while in b we have lost the leading word.
2586  We could go to DivLong now but we can also add an extra step that
2587  is less wasteful.
2588  In the case that the new leading word of b is extrememly short (like 1)
2589  we make a rather large error of course. In the worst case the whole
2590  will be intercepted by DivLong after all, but that is so rare that
2591  it shouldn't influence any timing in a measurable way.
2592 */
2593  if ( nb >= GCDMAX && na == nb+1 && b[nb-1] >= HALFMAX && b[nb-1] > a[na-1] ) {
2594  lx = (((RLONG)(a[na-1]))<<BITSINWORD) + a[na-2];
2595  x1[0] = lx/b[nb-1]; n1 = 1;
2596  MulLong(b,nb,x1,n1,x2,&n2);
2597  n2 = -n2;
2598  AddLong(a,na,x2,n2,x4,&n4);
2599  if ( n4 == 0 ) {
2600  *nc = nb;
2601  for ( i = 0; i < nb; i++ ) c[i] = b[i];
2602  goto normalend;
2603  }
2604  if ( n4 < 0 ) n4 = -n4;
2605  a = b; na = nb; b = x4; nb = n4;
2606  }
2607  goto restart;
2608 /*
2609  #] New code :
2610 */
2611  }
2612 normalend:
2613  NumberFree(GLscrat6,"GcdLong"); NumberFree(GLscrat7,"GcdLong"); NumberFree(GLscrat8,"GcdLong");
2614  NumberFree(GLscrat9,"GcdLong"); NumberFree(GLscrat10,"GcdLong");
2615  return(0);
2616 GcdErr:
2617  MLOCK(ErrorMessageLock);
2618  MesCall("GcdLong");
2619  MUNLOCK(ErrorMessageLock);
2620  NumberFree(GLscrat6,"GcdLong"); NumberFree(GLscrat7,"GcdLong"); NumberFree(GLscrat8,"GcdLong");
2621  NumberFree(GLscrat9,"GcdLong"); NumberFree(GLscrat10,"GcdLong");
2622  SETERROR(-1)
2623 }
2624 
2625 #endif
2626 
2627 /*
2628  #] GcdLong :
2629  #[ GetBinom : WORD GetBinom(a,na,i1,i2)
2630 */
2631 
2632 WORD GetBinom(UWORD *a, WORD *na, WORD i1, WORD i2)
2633 {
2634  GETIDENTITY
2635  WORD j, k, l;
2636  UWORD *GBscrat3, *GBscrat4;
2637  if ( i1-i2 < i2 ) i2 = i1-i2;
2638  if ( i2 == 0 ) { *a = 1; *na = 1; return(0); }
2639  if ( i2 > i1 ) { *a = 0; *na = 0; return(0); }
2640  *a = i1; *na = 1;
2641  GBscrat3 = NumberMalloc("GetBinom"); GBscrat4 = NumberMalloc("GetBinom");
2642  for ( j = 2; j <= i2; j++ ) {
2643  GBscrat3[0] = i1+1-j;
2644  if ( MulLong(a,*na,GBscrat3,(WORD)1,GBscrat4,&k) ) goto CalledFrom;
2645  GBscrat3[0] = j;
2646  if ( DivLong(GBscrat4,k,GBscrat3,(WORD)1,a,na,GBscrat3,&l) ) goto CalledFrom;
2647  }
2648  NumberFree(GBscrat3,"GetBinom"); NumberFree(GBscrat4,"GetBinom");
2649  return(0);
2650 CalledFrom:
2651  MLOCK(ErrorMessageLock);
2652  MesCall("GetBinom");
2653  MUNLOCK(ErrorMessageLock);
2654  NumberFree(GBscrat3,"GetBinom"); NumberFree(GBscrat4,"GetBinom");
2655  SETERROR(-1)
2656 }
2657 
2658 /*
2659  #] GetBinom :
2660  #[ LcmLong : WORD LcmLong(a,na,b,nb)
2661 
2662  Computes the LCM of the long numbers a and b and puts the result
2663  in c. c is allowed to be equal to a.
2664 */
2665 
2666 WORD LcmLong(PHEAD UWORD *a, WORD na, UWORD *b, WORD nb, UWORD *c, WORD *nc)
2667 {
2668  WORD error = 0;
2669  UWORD *d = NumberMalloc("LcmLong");
2670  UWORD *e = NumberMalloc("LcmLong");
2671  UWORD *f = NumberMalloc("LcmLong");
2672  WORD nd, ne, nf;
2673  GcdLong(BHEAD a, na, b, nb, d, &nd);
2674  DivLong(a,na,d,nd,e,&ne,f,&nf);
2675  if ( MulLong(b,nb,e,ne,c,nc) ) {
2676  MLOCK(ErrorMessageLock);
2677  MesCall("LcmLong");
2678  MUNLOCK(ErrorMessageLock);
2679  error = -1;
2680  }
2681  NumberFree(f,"LcmLong");
2682  NumberFree(e,"LcmLong");
2683  NumberFree(d,"LcmLong");
2684  return(error);
2685 }
2686 
2687 /*
2688  #] LcmLong :
2689  #[ TakeLongRoot: int TakeLongRoot(a,n,power)
2690 
2691  Takes the 'power'-root of the long number in a.
2692  If the root could be taken the return value is zero.
2693  If the root could not be taken, the return value is 1.
2694  The root will be in a if it could be taken, otherwise there will be garbage
2695  Algorithm: (assume b is guess of root, b' better guess)
2696  b' = (a-(power-1)*b^power)/(n*b^(power-1))
2697  Note: power should be positive!
2698 */
2699 
2700 int TakeLongRoot(UWORD *a, WORD *n, WORD power)
2701 {
2702  GETIDENTITY
2703  int numbits, guessbits, i, retval = 0;
2704  UWORD x, *b, *c, *d, *e;
2705  WORD na, nb, nc, nd, ne;
2706  if ( *n < 0 && ( power & 1 ) == 0 ) return(1);
2707  if ( power == 1 ) return(0);
2708  if ( *n < 0 ) { na = -*n; }
2709  else { na = *n; }
2710  if ( na == 1 ) {
2711 /* Special cases that are the most frequent */
2712  if ( a[0] == 1 ) return(0);
2713  if ( power < BITSINWORD && na == 1 && a[0] == (UWORD)(1<<power) ) {
2714  a[0] = 2; return(0);
2715  }
2716  if ( 2*power < BITSINWORD && na == 1 && a[0] == (UWORD)(1<<(2*power)) ) {
2717  a[0] = 4; return(0);
2718  }
2719  }
2720 /*
2721  Step 1: make a guess. We count the number of bits.
2722  numbits will be the 1+2log(a)
2723 */
2724  numbits = BITSINWORD*(na-1);
2725  x = a[na-1];
2726  while ( ( x >> 8 ) != 0 ) { numbits += 8; x >>= 8; }
2727  if ( ( x >> 4 ) != 0 ) { numbits += 4; x >>= 4; }
2728  if ( ( x >> 2 ) != 0 ) { numbits += 2; x >>= 2; }
2729  if ( ( x >> 1 ) != 0 ) numbits++;
2730  guessbits = numbits / power;
2731  if ( guessbits <= 0 ) return(1); /* root < 2 and 1 we did already */
2732  nb = guessbits/BITSINWORD;
2733 /*
2734  The recursion is:
2735  (b'-b) = (a/b^(power-1)-b)/n
2736  = (a/c-b)/n
2737  = (d-b)/n (remainder of a/c is e)
2738  = c/n (we reuse the scratch array c)
2739  Termination can be tricky. When a/c has no remainder and = b we have a root.
2740  When d = b but the remainder of a/c != 0, there is definitely no root.
2741 */
2742  b = NumberMalloc("TakeLongRoot"); c = NumberMalloc("TakeLongRoot");
2743  d = NumberMalloc("TakeLongRoot"); e = NumberMalloc("TakeLongRoot");
2744  for ( i = 0; i < nb; i++ ) { b[i] = 0; }
2745  b[nb] = 1 << (guessbits%BITSINWORD);
2746  nb++;
2747  for(;;) {
2748  nc = nb;
2749  for ( i = 0; i < nb; i++ ) c[i] = b[i];
2750  if ( RaisPow(BHEAD c,&nc,power-1) ) goto TLcall;
2751  if ( DivLong(a,na,c,nc,d,&nd,e,&ne) ) goto TLcall;
2752  nb = -nb;
2753  if ( AddLong(d,nd,b,nb,c,&nc) ) goto TLcall;
2754  nb = -nb;
2755  if ( nc == 0 ) {
2756  if ( ne == 0 ) break;
2757  retval = 1; break;
2758 /*
2759  else {
2760  NumberFree(b,"TakeLongRoot"); NumberFree(c,"TakeLongRoot");
2761  NumberFree(d,"TakeLongRoot"); NumberFree(e,"TakeLongRoot");
2762  return(1);
2763  }
2764 */
2765  }
2766  DivLong(c,nc,(UWORD *)(&power),1,d,&nd,e,&ne);
2767  if ( nd == 0 ) {
2768  retval = 1;
2769  break;
2770 /*
2771  NumberFree(b,"TakeLongRoot"); NumberFree(c,"TakeLongRoot");
2772  NumberFree(d,"TakeLongRoot"); NumberFree(e,"TakeLongRoot");
2773  return(1);
2774 */
2775 /*
2776  This code tries b+1 as a final possibility.
2777  We believe this is not needed
2778  UWORD one = 1;
2779  if ( AddLong(b,nb,&one,1,c,&nc) ) goto TLcall;
2780  if ( RaisPow(BHEAD c,&nc,power-1) ) goto TLcall;
2781  if ( DivLong(a,na,c,nc,d,&nd,e,&ne) ) goto TLcall;
2782  if ( ne != 0 ) return(1);
2783  nb = -nb;
2784  if ( SubLong(d,nd,b,nb,c,&nc) ) goto TLcall;
2785  nb = -nb;
2786  if ( nc != 0 ) {
2787  NumberFree(b,"TakeLongRoot"); NumberFree(c,"TakeLongRoot");
2788  NumberFree(d,"TakeLongRoot"); NumberFree(e,"TakeLongRoot");
2789  return(1);
2790  }
2791  break;
2792 */
2793  }
2794  if ( AddLong(b,nb,d,nd,b,&nb) ) goto TLcall;
2795  }
2796  for ( i = 0; i < nb; i++ ) a[i] = b[i];
2797  if ( *n < 0 ) *n = -nb;
2798  else *n = nb;
2799  NumberFree(b,"TakeLongRoot"); NumberFree(c,"TakeLongRoot");
2800  NumberFree(d,"TakeLongRoot"); NumberFree(e,"TakeLongRoot");
2801  return(retval);
2802 TLcall:
2803  MLOCK(ErrorMessageLock);
2804  MesCall("TakeLongRoot");
2805  MUNLOCK(ErrorMessageLock);
2806  NumberFree(b,"TakeLongRoot"); NumberFree(c,"TakeLongRoot");
2807  NumberFree(d,"TakeLongRoot"); NumberFree(e,"TakeLongRoot");
2808  Terminate(-1);
2809  return(-1);
2810 }
2811 
2812 /*
2813  #] TakeLongRoot:
2814  #[ MakeRational:
2815 
2816  Makes the integer a mod m into a traction b/c with |b|,|c| < sqrt(m)
2817  For the algorithm, see MakeLongRational.
2818 */
2819 
2820 int MakeRational(WORD a,WORD m, WORD *b, WORD *c)
2821 {
2822  LONG x1,x2,x3,x4,y1,y2;
2823  if ( a < 0 ) { a = a+m; }
2824  if ( a <= 1 ) {
2825  if ( a > m/2 ) a = a-m;
2826  *b = a; *c = 1; return(0);
2827  }
2828  x1 = m; x2 = a;
2829  if ( x2*x2 >= m ) {
2830  y1 = x1/x2; y2 = x1%x2; x3 = 1; x4 = -y1; x1 = x2; x2 = y2;
2831  while ( x2*x2 >= m ) {
2832  y1 = x1/x2; y2 = x1%x2; x1 = x2; x2 = y2; y2 = x3-y1*x4; x3 = x4; x4 = y2;
2833  }
2834  }
2835  else x4 = 1;
2836  if ( x2 == 0 ) { return(1); }
2837  if ( x2 > m/2 ) *b = x2-m;
2838  else *b = x2;
2839  if ( x4 > m/2 ) { *c = x4-m; *c = -*c; *b = -*b; }
2840  else if ( x4 <= -m/2 ) { x4 += m; *c = x4; }
2841  else if ( x4 < 0 ) { x4 = -x4; *c = x4; *b = -*b; }
2842  else *c = x4;
2843  return(0);
2844 }
2845 
2846 /*
2847  #] MakeRational:
2848  #[ MakeLongRational:
2849 
2850  Converts the long number a mod m into the fraction b
2851  One of the properties of b is that num,den < sqrt(m)
2852  The algorithm: Start with: m 0
2853  a 1
2854  Make now c=m%a, c1=m/a c c2=0-c1*1
2855  Make now d=a%c d1=a/c d d2=1-d1*c2
2856  Make now e=c%d e1=c/d e e2=1-e1*d2
2857  etc till in the first column we get a number < sqrt(m)
2858  We have then f,f2 and the fraction is f/f2.
2859  If at any moment we get a zero, m contained an unlucky prime.
2860 
2861  Note that this can be made a lot faster when we make the same
2862  improvements as in the GCD routine. That is something for later.
2863 #ifdef WITHMAKERATIONAL
2864 */
2865 
2866 #define COPYLONG(x1,nx1,x2,nx2) { int i; for(i=0;i<ABS(nx2);i++)x1[i]=x2[i];nx1=nx2; }
2867 
2868 int MakeLongRational(PHEAD UWORD *a, WORD na, UWORD *m, WORD nm, UWORD *b, WORD *nb)
2869 {
2870  UWORD *root = NumberMalloc("MakeRational");
2871  UWORD *x1 = NumberMalloc("MakeRational");
2872  UWORD *x2 = NumberMalloc("MakeRational");
2873  UWORD *x3 = NumberMalloc("MakeRational");
2874  UWORD *x4 = NumberMalloc("MakeRational");
2875  UWORD *y1 = NumberMalloc("MakeRational");
2876  UWORD *y2 = NumberMalloc("MakeRational");
2877  WORD nroot,nx1,nx2,nx3,nx4,ny1,ny2,retval = 0;
2878  WORD sign = 1;
2879 /*
2880  Step 1: Take the square root of m
2881 */
2882  COPYLONG(root,nroot,m,nm)
2883  TakeLongRoot(root,&nroot,2);
2884 /*
2885  Step 2: Set the start values
2886 */
2887  if ( na < 0 ) { na = -na; sign = -sign; }
2888  COPYLONG(x1,nx1,m,nm)
2889  COPYLONG(x2,nx2,a,na)
2890 /*
2891  x3[0] = 0, nx3 = 0;
2892  x4[0] = 1, nx4 = 1;
2893 */
2894 /*
2895  The start operation needs some special attention because of the zero.
2896 */
2897  if ( BigLong(x2,nx2,root,nroot) <= 0 ) {
2898  x4[0] = 1, nx4 = 1;
2899  goto gottheanswer;
2900  }
2901  DivLong(x1,nx1,x2,nx2,y1,&ny1,y2,&ny2);
2902  if ( ny2 == 0 ) { retval = 1; goto cleanup; }
2903  COPYLONG(x1,nx1,x2,nx2)
2904  COPYLONG(x2,nx2,y2,ny2)
2905  x3[0] = 1; nx3 = 1;
2906  COPYLONG(x4,nx4,y1,ny1)
2907  nx4 = -nx4;
2908 /*
2909  Now the loop.
2910 */
2911  while ( BigLong(x2,nx2,root,nroot) > 0 ) {
2912  DivLong(x1,nx1,x2,nx2,y1,&ny1,y2,&ny2);
2913  if ( ny2 == 0 ) { retval = 1; goto cleanup; }
2914  COPYLONG(x1,nx1,x2,nx2)
2915  COPYLONG(x2,nx2,y2,ny2)
2916  MulLong(y1,ny1,x4,nx4,y2,&ny2);
2917  ny2 = -ny2;
2918  AddLong(x3,nx3,y2,ny2,y1,&ny1);
2919  COPYLONG(x3,nx3,x4,nx4)
2920  COPYLONG(x4,nx4,y1,ny1)
2921  }
2922 /*
2923  Now we have the answer. It is x2/x4. It has to be packed into b.
2924 */
2925 gottheanswer:
2926  if ( nx4 < 0 ) { sign = -sign; nx4 = -nx4; }
2927  COPYLONG(b,*nb,x2,nx2)
2928  Pack(b,nb,x4,nx4);
2929  if ( sign < 0 ) *nb = -*nb;
2930 cleanup:
2931  NumberFree(y2,"MakeRational");
2932  NumberFree(y1,"MakeRational");
2933  NumberFree(x4,"MakeRational");
2934  NumberFree(x3,"MakeRational");
2935  NumberFree(x2,"MakeRational");
2936  NumberFree(x1,"MakeRational");
2937  NumberFree(root,"MakeRational");
2938  return(retval);
2939 }
2940 
2941 /*
2942 #endif
2943  #] MakeLongRational:
2944  #[ ChineseRemainder:
2945 */
2957 #ifdef WITHCHINESEREMAINDER
2958 
2959 int ChineseRemainder(PHEAD MODNUM *a1, MODNUM *a2, MODNUM *a)
2960 {
2961  UWORD *inv1 = NumberMalloc("ChineseRemainder");
2962  UWORD *inv2 = NumberMalloc("ChineseRemainder");
2963  UWORD *fac1 = NumberMalloc("ChineseRemainder");
2964  UWORD *fac2 = NumberMalloc("ChineseRemainder");
2965  UWORD two[1];
2966  WORD ninv1, ninv2, nfac1, nfac2;
2967  if ( a1->na < 0 ) {
2968  AddLong(a1->a,a1->na,a1->m,a1->nm,a1->a,&(a1->na));
2969  }
2970  if ( a2->na < 0 ) {
2971  AddLong(a2->a,a2->na,a2->m,a2->nm,a2->a,&(a2->na));
2972  }
2973  MulLong(a1->m,a1->nm,a2->m,a2->nm,a->m,&(a->nm));
2974 
2975  GetLongModInverses(BHEAD a1->m,a1->nm,a2->m,a2->nm,inv1,&ninv1,inv2,&ninv2);
2976  MulLong(inv1,ninv1,a1->m,a1->nm,fac1,&nfac1);
2977  MulLong(inv2,ninv2,a2->m,a2->nm,fac2,&nfac2);
2978 
2979  MulLong(fac1,nfac1,a2->a,a2->na,inv1,&ninv1);
2980  MulLong(fac2,nfac2,a1->a,a1->na,inv2,&ninv2);
2981  AddLong(inv1,ninv1,inv2,ninv2,a->a,&(a->na));
2982 
2983  two[0] = 2;
2984  MulLong(a->a,a->na,two,1,fac1,&nfac1);
2985  if ( BigLong(fac1,nfac1,a->m,a->nm) > 0 ) {
2986  a->nm = -a->nm;
2987  AddLong(a->a,a->na,a->m,a->nm,a->a,&(a->na));
2988  a->nm = -a->nm;
2989  }
2990  NumberFree(fac2,"ChineseRemainder");
2991  NumberFree(fac1,"ChineseRemainder");
2992  NumberFree(inv2,"ChineseRemainder");
2993  NumberFree(inv1,"ChineseRemainder");
2994  return(0);
2995 }
2996 
2997 #endif
2998 
2999 /*
3000  #] ChineseRemainder:
3001  #] RekenLong :
3002  #[ RekenTerms :
3003  #[ CompCoef : WORD CompCoef(term1,term2)
3004 
3005  Compares the coefficients of term1 and term2 by subtracting them.
3006  This does more work than needed but this routine is only called
3007  when sorting functions and function arguments.
3008  (and comparing values
3009 */
3010 /* #define 64SAVE */
3011 
3012 WORD CompCoef(WORD *term1, WORD *term2)
3013 {
3014  GETIDENTITY
3015  UWORD *c;
3016  WORD n1,n2,n3,*a;
3017  GETCOEF(term1,n1);
3018  GETCOEF(term2,n2);
3019  if ( term1[1] == 0 && n1 == 1 ) {
3020  if ( term2[1] == 0 && n2 == 1 ) return(0);
3021  if ( n2 < 0 ) return(1);
3022  return(-1);
3023  }
3024  else if ( term2[1] == 0 && n2 == 1 ) {
3025  if ( n1 < 0 ) return(-1);
3026  return(1);
3027  }
3028  if ( n1 > 0 ) {
3029  if ( n2 < 0 ) return(1);
3030  }
3031  else {
3032  if ( n2 > 0 ) return(-1);
3033  a = term1; term1 = term2; term2 = a;
3034  n3 = -n1; n1 = -n2; n2 = n3;
3035  }
3036  if ( term1[1] == 1 && term2[1] == 1 && n1 == 1 && n2 == 1 ) {
3037  if ( (UWORD)*term1 > (UWORD)*term2 ) return(1);
3038  else if ( (UWORD)*term1 < (UWORD)*term2 ) return(-1);
3039  else return(0);
3040  }
3041 
3042 /*
3043  The next call should get dedicated code, as AddRat does more than
3044  strictly needed. Also more attention should be given to overflow.
3045 */
3046  c = NumberMalloc("CompCoef");
3047  if ( AddRat(BHEAD (UWORD *)term1,n1,(UWORD *)term2,-n2,c,&n3) ) {
3048  MLOCK(ErrorMessageLock);
3049  MesCall("CompCoef");
3050  MUNLOCK(ErrorMessageLock);
3051  NumberFree(c,"CompCoef");
3052  SETERROR(-1)
3053  }
3054  NumberFree(c,"CompCoef");
3055  return(n3);
3056 }
3057 
3058 /*
3059  #] CompCoef :
3060  #[ Modulus : WORD Modulus(term)
3061 
3062  Routine takes the coefficient of term modulus b. The answer
3063  is in term again and the length of term is adjusted.
3064 
3065 */
3066 
3067 WORD Modulus(WORD *term)
3068 {
3069  WORD *t;
3070  WORD n1;
3071  t = term;
3072  GETCOEF(t,n1);
3073  if ( TakeModulus((UWORD *)t,&n1,AC.cmod,AC.ncmod,UNPACK) ) {
3074  MLOCK(ErrorMessageLock);
3075  MesCall("Modulus");
3076  MUNLOCK(ErrorMessageLock);
3077  SETERROR(-1)
3078  }
3079  if ( !n1 ) {
3080  *term = 0;
3081  return(0);
3082  }
3083  else if ( n1 > 0 ) {
3084  n1 <<= 1;
3085  t += n1; /* Note that n1 >= 0 */
3086  n1++;
3087  }
3088  else if ( n1 < 0 ) {
3089  n1 *= 2;
3090  t += -n1;
3091  n1--;
3092  }
3093  *t++ = n1;
3094  *term = WORDDIF(t,term);
3095  return(0);
3096 }
3097 
3098 /*
3099  #] Modulus :
3100  #[ TakeModulus : WORD TakeModulus(a,na,cmodvec,ncmod,par)
3101 
3102  Routine gets the rational number in a with reduced length na.
3103  It is called when AC.ncmod != 0 and the number in AC.cmod is the
3104  number wrt which we need the modulus.
3105  The result is returned in a and na again.
3106 
3107  If par == NOUNPACK we only do a single number, not a fraction.
3108  In addition we don't do fancy. We want a positive number and
3109  the input was supposed to be positive.
3110  We don't pack the result. The calling routine is responsible for that.
3111  This may not be a good idea. To be checked.
3112 */
3113 
3114 WORD TakeModulus(UWORD *a, WORD *na, UWORD *cmodvec, WORD ncmod, WORD par)
3115 {
3116  GETIDENTITY
3117  UWORD *c, *d, *e, *f, *g, *h;
3118  UWORD *x4,*x2;
3119  UWORD *x3,*x1,*x5,*x6,*x7,*x8;
3120  WORD y3,y1,y5,y6;
3121  WORD n1, i, y2, y4;
3122  WORD nh, tdenom, tnumer, nmod;
3123  LONG x;
3124  if ( ncmod == 0 ) return(0); /* No modulus operation */
3125  nmod = ABS(ncmod);
3126  n1 = *na;
3127  if ( ( par & UNPACK ) != 0 ) UnPack(a,n1,&tdenom,&tnumer);
3128  else { tnumer = n1; }
3129 /*
3130  We fish out the special case that the coefficient is short as well.
3131  There is no need to make lots of calls etc
3132 */
3133  if ( ( ( par & UNPACK ) == 0 ) && nmod == 1 && ( n1 == 1 || n1 == -1 ) ) {
3134  goto simplecase;
3135  }
3136  else if ( nmod == 1 && ( n1 == 1 || n1 == -1 ) ) {
3137  if ( a[1] != 1 ) {
3138  a[1] = a[1] % cmodvec[0];
3139  if ( a[1] == 0 ) {
3140  MesPrint("Division by zero in short modulus arithmetic");
3141  return(-1);
3142  }
3143  y1 = 0;
3144  if ( ( AC.modinverses != 0 ) && ( ( par & NOINVERSES ) == 0 ) ) {
3145  y1 = AC.modinverses[a[1]];
3146  }
3147  else {
3148  GetModInverses(a[1],cmodvec[0],&y1,&y2);
3149  }
3150  x = a[0];
3151  a[0] = (x*y1) % cmodvec[0];
3152  a[1] = 1;
3153  }
3154  else {
3155 simplecase:
3156  a[0] = a[0] % cmodvec[0];
3157  }
3158  if ( a[0] == 0 ) { *na = 0; return(0); }
3159  if ( ( AC.modmode & POSNEG ) != 0 ) {
3160  if ( a[0] > (UWORD)(cmodvec[0]/2) ) {
3161  a[0] = cmodvec[0] - a[0];
3162  *na = -*na;
3163  }
3164  }
3165  else if ( *na < 0 ) {
3166  *na = 1; a[0] = cmodvec[0] - a[0];
3167  }
3168  return(0);
3169  }
3170  c = NumberMalloc("TakeModulus"); d = NumberMalloc("TakeModulus"); e = NumberMalloc("TakeModulus");
3171  f = NumberMalloc("TakeModulus"); g = NumberMalloc("TakeModulus"); h = NumberMalloc("TakeModulus");
3172  n1 = ABS(n1);
3173  if ( DivLong(a,tnumer,(UWORD *)cmodvec,nmod,
3174  c,&nh,a,&tnumer) ) goto ModErr;
3175  if ( tnumer == 0 ) { *na = 0; goto normalreturn; }
3176  if ( ( par & UNPACK ) == 0 ) {
3177  if ( ( AC.modmode & POSNEG ) != 0 ) {
3178  NormalModulus(a,&tnumer);
3179  }
3180  else if ( tnumer < 0 ) {
3181  SubPLon((UWORD *)cmodvec,nmod,a,-tnumer,a,&tnumer);
3182  }
3183  *na = tnumer;
3184  goto normalreturn;
3185  }
3186  if ( tdenom == 1 && a[n1] == 1 ) {
3187  if ( ( AC.modmode & POSNEG ) != 0 ) {
3188  NormalModulus(a,&tnumer);
3189  }
3190  else if ( tnumer < 0 ) {
3191  SubPLon((UWORD *)cmodvec,nmod,a,-tnumer,a,&tnumer);
3192  }
3193  *na = tnumer;
3194  i = ABS(tnumer);
3195  a += i;
3196  *a++ = 1;
3197  while ( --i > 0 ) *a++ = 0;
3198  goto normalreturn;
3199  }
3200  if ( DivLong(a+n1,tdenom,(UWORD *)cmodvec,nmod,c,&nh,a+n1,&tdenom) ) goto ModErr;
3201  if ( !tdenom ) {
3202  MLOCK(ErrorMessageLock);
3203  MesPrint("Division by zero in modulus arithmetic");
3204  if ( AP.DebugFlag ) {
3205  AO.OutSkip = 3;
3206  FiniLine();
3207  i = *na;
3208  if ( i < 0 ) i = -i;
3209  while ( --i >= 0 ) { TalToLine((UWORD)(*a++)); TokenToLine((UBYTE *)" "); }
3210  i = *na;
3211  if ( i < 0 ) i = -i;
3212  while ( --i >= 0 ) { TalToLine((UWORD)(*a++)); TokenToLine((UBYTE *)" "); }
3213  TalToLine((UWORD)(*na));
3214  AO.OutSkip = 0;
3215  FiniLine();
3216  }
3217  MUNLOCK(ErrorMessageLock);
3218  NumberFree(c,"TakeModulus"); NumberFree(d,"TakeModulus"); NumberFree(e,"TakeModulus");
3219  NumberFree(f,"TakeModulus"); NumberFree(g,"TakeModulus"); NumberFree(h,"TakeModulus");
3220  return(-1);
3221  }
3222  if ( ( AC.modinverses != 0 ) && ( ( par & NOINVERSES ) == 0 )
3223  && ( tdenom == 1 || tdenom == -1 ) ) {
3224  *d = AC.modinverses[a[n1]]; y1 = 1; y2 = tdenom;
3225  if ( MulLong(a,tnumer,d,y1,c,&y3) ) goto ModErr;
3226  if ( DivLong(c,y3,(UWORD *)cmodvec,nmod,d,&y5,a,&tdenom) ) goto ModErr;
3227  if ( y2 < 0 ) tdenom = -tdenom;
3228  }
3229  else {
3230  x2 = (UWORD *)cmodvec; x1 = c; i = nmod; while ( --i >= 0 ) *x1++ = *x2++;
3231  x1 = c; x2 = a+n1; x3 = d; x4 = e; x5 = f; x6 = g;
3232  y1 = nmod; y2 = tdenom; y4 = 0; y5 = 1; *x5 = 1;
3233  for(;;) {
3234  if ( DivLong(x1,y1,x2,y2,h,&nh,x3,&y3) ) goto ModErr;
3235  if ( MulLong(x5,y5,h,nh,x6,&y6) ) goto ModErr;
3236  if ( AddLong(x4,y4,x6,-y6,x6,&y6) ) goto ModErr;
3237  if ( !y3 ) {
3238  if ( y2 != 1 || *x2 != 1 ) {
3239  MLOCK(ErrorMessageLock);
3240  MesPrint("Inverse in modulus arithmetic doesn't exist");
3241  MesPrint("Denominator and modulus are not relative prime");
3242  MUNLOCK(ErrorMessageLock);
3243  goto ModErr;
3244  }
3245  break;
3246  }
3247  x7 = x1; x1 = x2; y1 = y2; x2 = x3; y2 = y3; x3 = x7;
3248  x8 = x4; x4 = x5; y4 = y5; x5 = x6; y5 = y6; x6 = x8;
3249  }
3250  if ( y5 < 0 && AddLong((UWORD *)cmodvec,nmod,x5,y5,x5,&y5) ) goto ModErr;
3251  if ( MulLong(a,tnumer,x5,y5,c,&y3) ) goto ModErr;
3252  if ( DivLong(c,y3,(UWORD *)cmodvec,nmod,d,&y5,a,&tdenom) ) goto ModErr;
3253  }
3254  if ( !tdenom ) { *na = 0; goto normalreturn; }
3255  if ( ( ( AC.modmode & POSNEG ) != 0 ) && ( ( par & FROMFUNCTION ) == 0 ) ) {
3256  NormalModulus(a,&tdenom);
3257  }
3258  else if ( tdenom < 0 ) {
3259  SubPLon((UWORD *)cmodvec,nmod,a,-tdenom,a,&tdenom);
3260  }
3261  *na = tdenom;
3262  i = ABS(tdenom);
3263  a += i;
3264  *a++ = 1;
3265  while ( --i > 0 ) *a++ = 0;
3266 normalreturn:
3267  NumberFree(c,"TakeModulus"); NumberFree(d,"TakeModulus"); NumberFree(e,"TakeModulus");
3268  NumberFree(f,"TakeModulus"); NumberFree(g,"TakeModulus"); NumberFree(h,"TakeModulus");
3269  return(0);
3270 ModErr:
3271  MLOCK(ErrorMessageLock);
3272  MesCall("TakeModulus");
3273  MUNLOCK(ErrorMessageLock);
3274  NumberFree(c,"TakeModulus"); NumberFree(d,"TakeModulus"); NumberFree(e,"TakeModulus");
3275  NumberFree(f,"TakeModulus"); NumberFree(g,"TakeModulus"); NumberFree(h,"TakeModulus");
3276  SETERROR(-1)
3277 }
3278 
3279 /*
3280  #] TakeModulus :
3281  #[ TakeNormalModulus : WORD TakeNormalModulus(a,na,par)
3282 
3283  added by Jan [01-09-2010]
3284 */
3285 
3286 WORD TakeNormalModulus (UWORD *a, WORD *na, UWORD *c, WORD nc, WORD par)
3287 {
3288  WORD n;
3289  WORD nhalfc;
3290  UWORD *halfc;
3291 
3292  GETIDENTITY;
3293 
3294  /* determine c/2 by right shifting */
3295  halfc = NumberMalloc("TakeNormalModulus");
3296  nhalfc=nc;
3297  WCOPY(halfc,c,nc);
3298 
3299  for (n=0; n<nhalfc; n++) {
3300  halfc[n] >>= 1;
3301  if (n+1<nc) halfc[n] |= c[n+1] << (BITSINWORD-1);
3302  }
3303 
3304  if (halfc[nhalfc-1]==0)
3305  nhalfc--;
3306 
3307  /* takes care of the number never expanding, e.g., -1(mod 100) -> 99 -> -1 */
3308  if (BigLong(a,ABS(*na),halfc,nhalfc) > 0) {
3309 
3310  TakeModulus(a,na,c,nc,par);
3311 
3312  n = ABS(*na);
3313  if (BigLong(a,n,halfc,nhalfc) > 0) {
3314  SubPLon(c,nc,a,n,a,&n);
3315  *na = (*na > 0 ? -n : n);
3316  }
3317  }
3318 
3319  NumberFree(halfc,"TakeNormalModulus");
3320  return(0);
3321 }
3322 
3323 /*
3324  #] TakeNormalModulus :
3325  #[ MakeModTable : WORD MakeModTable()
3326 */
3327 
3328 WORD MakeModTable()
3329 {
3330  LONG size, i, j, n;
3331  n = ABS(AC.ncmod);
3332  if ( AC.modpowers ) {
3333  M_free(AC.modpowers,"AC.modpowers");
3334  AC.modpowers = NULL;
3335  }
3336  if ( n > 2 ) {
3337  MLOCK(ErrorMessageLock);
3338  MesPrint("&No memory for modulus generator power table");
3339  MUNLOCK(ErrorMessageLock);
3340  Terminate(-1);
3341  }
3342  if ( n == 0 ) return(0);
3343  size = (LONG)(*AC.cmod);
3344  if ( n == 2 ) size += (((LONG)AC.cmod[1])<<BITSINWORD);
3345  AC.modpowers = (UWORD *)Malloc1(size*n*sizeof(UWORD),"table for powers of modulus");
3346  if ( n == 1 ) {
3347  j = 1;
3348  for ( i = 0; i < size; i++ ) AC.modpowers[i] = 0;
3349  for ( i = 0; i < size; i++ ) {
3350  AC.modpowers[j] = (WORD)i;
3351  j *= *AC.powmod;
3352  j %= *AC.cmod;
3353  }
3354  for ( i = 2; i < size; i++ ) {
3355  if ( AC.modpowers[i] == 0 ) {
3356  MLOCK(ErrorMessageLock);
3357  MesPrint("&improper generator for this modulus");
3358  MUNLOCK(ErrorMessageLock);
3359  M_free(AC.modpowers,"AC.modpowers");
3360  return(-1);
3361  }
3362  }
3363  AC.modpowers[1] = 0;
3364  }
3365  else {
3366  GETIDENTITY
3367  WORD nScrat, n2;
3368  UWORD *MMscrat7 = NumberMalloc("MakeModTable"), *MMscratC = NumberMalloc("MakeModTable");
3369  *MMscratC = 1;
3370  nScrat = 1;
3371  j = size << 1;
3372  for ( i = 0; i < j; i+=2 ) { AC.modpowers[i] = 0; AC.modpowers[i+1] = 0; }
3373  for ( i = 0; i < size; i++ ) {
3374  j = *MMscratC + (((LONG)MMscratC[1])<<BITSINWORD);
3375  j <<= 1;
3376  AC.modpowers[j] = (WORD)(i & WORDMASK);
3377  AC.modpowers[j+1] = (WORD)(i >> BITSINWORD);
3378  MulLong((UWORD *)MMscratC,nScrat,(UWORD *)AC.powmod,
3379  AC.npowmod,(UWORD *)MMscrat7,&n2);
3380  TakeModulus(MMscrat7,&n2,AC.cmod,AC.ncmod,NOUNPACK);
3381  *MMscratC = *MMscrat7; MMscratC[1] = MMscrat7[1]; nScrat = n2;
3382  }
3383  NumberFree(MMscrat7,"MakeModTable"); NumberFree(MMscratC,"MakeModTable");
3384  j = size << 1;
3385  for ( i = 4; i < j; i+=2 ) {
3386  if ( AC.modpowers[i] == 0 && AC.modpowers[i+1] == 0 ) {
3387  MLOCK(ErrorMessageLock);
3388  MesPrint("&improper generator for this modulus");
3389  MUNLOCK(ErrorMessageLock);
3390  M_free(AC.modpowers,"AC.modpowers");
3391  return(-1);
3392  }
3393  }
3394  AC.modpowers[2] = AC.modpowers[3] = 0;
3395  }
3396  return(0);
3397 }
3398 
3399 /*
3400  #] MakeModTable :
3401  #] RekenTerms :
3402  #[ Functions :
3403  #[ Factorial : WORD Factorial(n,a,na)
3404 
3405  Starts with only the value of fac_(0).
3406  Builds up what is needed and remembers it for the next time.
3407 
3408  We have:
3409  AT.nfac: the number of the highest stored factorial
3410  AT.pfac: the array of locations in the array of stored factorials
3411  AT.factorials: the array with stored factorials
3412 */
3413 
3414 int Factorial(PHEAD WORD n, UWORD *a, WORD *na)
3415 {
3416  GETBIDENTITY
3417  UWORD *b, *c;
3418  WORD nc;
3419  int i, j;
3420  LONG ii;
3421  if ( n > AT.nfac ) {
3422  if ( AT.factorials == 0 ) {
3423  AT.nfac = 0; AT.mfac = 50; AT.sfact = 400;
3424  AT.pfac = (LONG *)Malloc1((AT.mfac+2)*sizeof(LONG),"factorials");
3425  AT.factorials = (UWORD *)Malloc1(AT.sfact*sizeof(UWORD),"factorials");
3426  AT.factorials[0] = 1; AT.pfac[0] = 0; AT.pfac[1] = 1;
3427  }
3428  b = a;
3429  c = AT.factorials+AT.pfac[AT.nfac];
3430  nc = i = AT.pfac[AT.nfac+1] - AT.pfac[AT.nfac];
3431  while ( --i >= 0 ) *b++ = *c++;
3432  for ( j = AT.nfac+1; j <= n; j++ ) {
3433  Product(a,&nc,j);
3434  if ( nc > AM.MaxTal ) {
3435  MLOCK(ErrorMessageLock);
3436  MesPrint("Overflow in factorial. MaxTal = %d",AM.MaxTal);
3437  MesPrint("Increase MaxTerm in %s",setupfilename);
3438  MUNLOCK(ErrorMessageLock);
3439  return(-1);
3440  }
3441  if ( j > AT.mfac ) { /* double the pfac buffer */
3442  LONG *p;
3443  p = (LONG *)Malloc1((AT.mfac*2+2)*sizeof(LONG),"factorials");
3444  i = AT.mfac;
3445  for ( i = AT.mfac+1; i >= 0; i-- ) p[i] = AT.pfac[i];
3446  M_free(AT.pfac,"factorial pointers"); AT.pfac = p; AT.mfac *= 2;
3447  }
3448  if ( AT.pfac[j] + nc >= AT.sfact ) { /* double the factorials buffer */
3449  UWORD *f;
3450  f = (UWORD *)Malloc1(AT.sfact*2*sizeof(UWORD),"factorials");
3451  ii = AT.sfact;
3452  c = AT.factorials; b = f;
3453  while ( --ii >= 0 ) *b++ = *c++;
3454  M_free(AT.factorials,"factorials");
3455  AT.factorials = f;
3456  AT.sfact *= 2;
3457  }
3458  b = a; c = AT.factorials + AT.pfac[j]; i = nc;
3459  while ( --i >= 0 ) *c++ = *b++;
3460  AT.pfac[j+1] = AT.pfac[j] + nc;
3461  }
3462  *na = nc;
3463  AT.nfac = n;
3464  }
3465  else if ( n == 0 ) {
3466  *a = 1; *na = 1;
3467  }
3468  else {
3469  *na = i = AT.pfac[n+1] - AT.pfac[n];
3470  b = AT.factorials + AT.pfac[n];
3471  while ( --i >= 0 ) *a++ = *b++;
3472  }
3473  return(0);
3474 }
3475 
3476 /*
3477  #] Factorial :
3478  #[ Bernoulli : WORD Bernoulli(n,a,na)
3479 
3480  Starts with only the value of bernoulli_(0).
3481  Builds up what is needed and remembers it for the next time.
3482  b_0 = 1
3483  (n+1)*b_n = -b_{n-1}-sum_(i,1,n-1,b_i*b_{n-i})
3484  The n-1 playes only a role for b_2.
3485  We have hard coded b_0,b_1,b_2 and b_odd. After that:
3486  (2n+1)*b_2n = -sum_(i,1,n-1,b_2i*b_{2n-2i})
3487 
3488  We have:
3489  AT.nBer: the number of the highest stored Bernoulli number
3490  AT.pBer: the array of locations in the array of stored Bernoulli numbers
3491  AT.bernoullis: the array with stored Bernoulli numbers
3492 */
3493 
3494 int Bernoulli(WORD n, UWORD *a, WORD *na)
3495 {
3496  GETIDENTITY
3497  UWORD *b, *c, *scrib, *ntop, *ntop1;
3498  WORD i, i1, i2, nhalf, nqua, nscrib, nntop, nntop1, *oldworkpointer;
3499  UWORD twee = 2, twonplus1;
3500  int j;
3501  LONG ii;
3502  if ( n <= 1 ) {
3503  if ( n == 0 ) { a[0] = a[1] = 1; *na = 3; }
3504  else if ( n == 1 ) { a[0] = 1; a[1] = 2; *na = 3; }
3505  return(0);
3506  }
3507  if ( ( n & 1 ) != 0 ) { a[0] = a[1] = 0; *na = 0; return(0); }
3508  nhalf = n/2;
3509  if ( nhalf > AT.nBer ) {
3510  oldworkpointer = AT.WorkPointer;
3511  if ( AT.bernoullis == 0 ) {
3512  AT.nBer = 1; AT.mBer = 50; AT.sBer = 400;
3513  AT.pBer = (LONG *)Malloc1((AT.mBer+2)*sizeof(LONG),"bernoullis");
3514  AT.bernoullis = (UWORD *)Malloc1(AT.sBer*sizeof(UWORD),"bernoullis");
3515  AT.pBer[1] = 0; AT.pBer[2] = 3;
3516  AT.bernoullis[0] = 3; AT.bernoullis[1] = 1; AT.bernoullis[2] = 12;
3517  if ( nhalf == 1 ) {
3518  a[0] = 1; a[1] = 12; *na = 3; return(0);
3519  }
3520  }
3521  while ( nhalf > AT.mBer ) {
3522  LONG *p;
3523  p = (LONG *)Malloc1((AT.mBer*2+1)*sizeof(LONG),"bernoullis");
3524  i = AT.mBer;
3525  for ( i = AT.mBer; i >= 0; i-- ) p[i] = AT.pBer[i];
3526  M_free(AT.pBer,"factorial pointers"); AT.pBer = p; AT.mBer *= 2;
3527  }
3528  for ( n = AT.nBer+1; n <= nhalf; n++ ) {
3529  scrib = (UWORD *)(AT.WorkPointer);
3530  nqua = n/2;
3531  if ( ( n & 1 ) == 1 ) {
3532  nscrib = 0; ntop = scrib;
3533  }
3534  else {
3535  b = AT.bernoullis + AT.pBer[nqua];
3536  nscrib = *b++;
3537  i = (WORD)(REDLENG(nscrib));
3538  MulRat(BHEAD b,i,b,i,scrib,&nscrib);
3539  ntop = scrib + 2*nscrib;
3540  nqua--;
3541  }
3542  for ( j = 1; j <= nqua; j++ ) {
3543  b = AT.bernoullis + AT.pBer[j];
3544  c = AT.bernoullis + AT.pBer[n-j];
3545  i1 = (WORD)(*b); i2 = (WORD)(*c);
3546  i1 = REDLENG(i1);
3547  i2 = REDLENG(i2);
3548  MulRat(BHEAD b+1,i1,c+1,i2,ntop,&nntop);
3549  Mully(BHEAD ntop,&nntop,&twee,1);
3550  if ( nscrib ) {
3551  i = (WORD)nntop; if ( i < 0 ) i = -i;
3552  ntop1 = ntop + 2*i;
3553  AddRat(BHEAD ntop,nntop,scrib,nscrib,ntop1,&nntop1);
3554  }
3555  else {
3556  ntop1 = ntop; nntop1 = nntop;
3557  }
3558  nscrib = i1 = (WORD)nntop1;
3559  if ( i1 < 0 ) i1 = - i1;
3560  i1 = 2*i1;
3561  for ( i = 0; i < i1; i++ ) scrib[i] = ntop1[i];
3562  ntop = scrib + i1;
3563  }
3564  twonplus1 = 2*n+1;
3565  Divvy(BHEAD scrib,&nscrib,&twonplus1,-1);
3566  i1 = INCLENG(nscrib);
3567  i2 = i1; if ( i2 < 0 ) i2 = -i2;
3568  i = (WORD)(AT.bernoullis[AT.pBer[n-1]]);
3569  if ( i < 0 ) i = -i;
3570  AT.pBer[n] = AT.pBer[n-1]+i;
3571  if ( AT.pBer[n] + i2 >= AT.sBer ) {
3572  UWORD *f;
3573  f = (UWORD *)Malloc1(AT.sBer*2*sizeof(UWORD),"bernoullis");
3574  ii = AT.sBer;
3575  c = AT.bernoullis; b = f;
3576  while ( --ii >= 0 ) *b++ = *c++;
3577  M_free(AT.bernoullis,"bernoullis");
3578  AT.bernoullis = f;
3579  AT.sBer *= 2;
3580  }
3581  c = AT.bernoullis + AT.pBer[n]; b = scrib;
3582  *c++ = i1;
3583  for ( i = 1; i < i2; i++ ) *c++ = *b++;
3584  }
3585  AT.nBer = nhalf;
3586  AT.WorkPointer = oldworkpointer;
3587  }
3588  b = AT.bernoullis + AT.pBer[nhalf];
3589  *na = i = (WORD)(*b++);
3590  if ( i < 0 ) i = -i;
3591  i--;
3592  while ( --i >= 0 ) *a++ = *b++;
3593  return(0);
3594 }
3595 
3596 /*
3597  #] Bernoulli :
3598  #[ NextPrime :
3599 */
3611 #if ( BITSINWORD == 32 )
3612 
3613 void StartPrimeList(PHEAD0)
3614 {
3615  int i, j;
3616  AR.PrimeList[AR.numinprimelist++] = 3;
3617  for ( i = 5; i < 46340; i += 2 ) {
3618  for ( j = 0; j < AR.numinprimelist && AR.PrimeList[j]*AR.PrimeList[j] <= i; j++ ) {
3619  if ( i % AR.PrimeList[j] == 0 ) goto nexti;
3620  }
3621  AR.PrimeList[AR.numinprimelist++] = i;
3622 nexti:;
3623  }
3624  AR.notfirstprime = 1;
3625 }
3626 
3627 #endif
3628 
3629 WORD NextPrime(PHEAD WORD num)
3630 {
3631  int i, j;
3632  WORD *newpl;
3633  LONG newsize, x;
3634 #if ( BITSINWORD == 32 )
3635  if ( AR.notfirstprime == 0 ) StartPrimeList(BHEAD0);
3636 #endif
3637  if ( num > AT.inprimelist ) {
3638  while ( AT.inprimelist < num ) {
3639  if ( num >= AT.sizeprimelist ) {
3640  if ( AT.sizeprimelist == 0 ) newsize = 32;
3641  else newsize = 2*AT.sizeprimelist;
3642  while ( num >= newsize ) newsize = newsize*2;
3643  newpl = (WORD *)Malloc1(newsize*sizeof(WORD),"NextPrime");
3644  for ( i = 0; i < AT.sizeprimelist; i++ ) {
3645  newpl[i] = AT.primelist[i];
3646  }
3647  if ( AT.sizeprimelist > 0 ) {
3648  M_free(AT.primelist,"NextPrime");
3649  }
3650  AT.sizeprimelist = newsize;
3651  AT.primelist = newpl;
3652  }
3653  if ( AT.inprimelist < 0 ) { i = MAXPOSITIVE; }
3654  else { i = AT.primelist[AT.inprimelist]; }
3655  while ( i > MAXPOWER ) {
3656  i -= 2; x = i;
3657 #if ( BITSINWORD == 32 )
3658  for ( j = 0; j < AR.numinprimelist && AR.PrimeList[j]*(LONG)(AR.PrimeList[j]) <= x; j++ ) {
3659  if ( x % AR.PrimeList[j] == 0 ) goto nexti;
3660  }
3661 #else
3662  for ( j = 3; j*((LONG)j) <= x; j += 2 ) {
3663  if ( x % j == 0 ) goto nexti;
3664  }
3665 #endif
3666  AT.inprimelist++;
3667  AT.primelist[AT.inprimelist] = i;
3668  break;
3669 nexti:;
3670  }
3671  if ( i < MAXPOWER ) {
3672  MLOCK(ErrorMessageLock);
3673  MesPrint("There are not enough short prime numbers for this calculation");
3674  MesPrint("Try to use a computer with a %d-bits architecture",
3675  (int)(BITSINWORD*4));
3676  MUNLOCK(ErrorMessageLock);
3677  Terminate(-1);
3678  }
3679  }
3680  }
3681  return(AT.primelist[num]);
3682 }
3683 
3684 /*
3685  #] NextPrime :
3686  #[ wranf :
3687 
3688  A random number generator that generates random WORDs with a very
3689  long sequence. It is based on the Knuth generator.
3690 
3691  We take some care that each thread can run its own, but each
3692  uses its own startup. Hence the seed includes the identity of
3693  the thread.
3694 
3695  For NPAIR1, NPAIR2 we can use any pair from the table on page 28.
3696  Candidates are 24,55 (the example on the pages 171,172)
3697  or (33,97) or (38,89)
3698  These values are defined in fsizes.h and used in startup.c and threads.c
3699 */
3700 
3701 #define WARMUP 6
3702 
3703 static void wranfnew(PHEAD0)
3704 {
3705  int i;
3706  LONG j;
3707  for ( i = 0; i < AR.wranfnpair1; i++ ) {
3708  j = AR.wranfia[i] - AR.wranfia[i+(AR.wranfnpair2-AR.wranfnpair1)];
3709  if ( j < 0 ) j += (LONG)1 << (2*BITSINWORD-2);
3710  AR.wranfia[i] = j;
3711  }
3712  for ( i = AR.wranfnpair1; i < AR.wranfnpair2; i++ ) {
3713  j = AR.wranfia[i] - AR.wranfia[i-AR.wranfnpair1];
3714  if ( j < 0 ) j += (LONG)1 << (2*BITSINWORD-2);
3715  AR.wranfia[i] = j;
3716  }
3717 }
3718 
3719 void iniwranf(PHEAD0)
3720 {
3721  int imax = AR.wranfnpair2-1;
3722  ULONG i, ii, seed = AR.wranfseed;
3723  LONG j, k;
3724  ULONG offset = 12345;
3725 #ifdef PARALLELCODE
3726  int id;
3727 #if defined(WITHPTHREADS)
3728  id = AT.identity;
3729 #elif defined(WITHMPI)
3730  id = PF.me;
3731 #endif
3732  seed += id;
3733  i = id + 1;
3734  if ( i > 1 ) {
3735  ULONG pow, accu;
3736  pow = offset; accu = 1;
3737  while ( i ) {
3738  if ( ( i & 1 ) != 0 ) accu *= pow;
3739  i >>= 1; pow = pow*pow;
3740  }
3741  offset = accu;
3742  }
3743 #endif
3744  if ( seed < ((LONG)1<<(BITSINWORD-1)) ) {
3745  j = ( (seed+31459L) << (BITSINWORD-2))+offset;
3746  }
3747  else if ( seed < ((LONG)1<<(BITSINWORD+10-1)) ) {
3748  j = ( (seed+31459L) << (BITSINWORD-10-2))+offset;
3749  }
3750  else {
3751  j = ( (seed+31459L) << 1)+offset;
3752  }
3753  if ( ( seed & 1 ) == 1 ) seed++;
3754  j += seed;
3755  AR.wranfia[imax] = j;
3756  k = 1;
3757  for ( i = 0; i <= (ULONG)(imax); i++ ) {
3758  ii = (AR.wranfnpair1*i)%AR.wranfnpair2;
3759  AR.wranfia[ii] = k;
3760  k = j - k;
3761  if ( k < 0 ) k += (LONG)1 << (2*BITSINWORD-2);
3762  j = AR.wranfia[ii];
3763  }
3764  for ( i = 0; i < WARMUP; i++ ) wranfnew(BHEAD0);
3765  AR.wranfcall = 0;
3766 }
3767 
3768 UWORD wranf(PHEAD0)
3769 {
3770  UWORD wval;
3771  if ( AR.wranfia == 0 ) {
3772  AR.wranfia = (ULONG *)Malloc1(AR.wranfnpair2*sizeof(ULONG),"wranf");
3773  iniwranf(BHEAD0);
3774  }
3775  if ( AR.wranfcall >= AR.wranfnpair2) {
3776  wranfnew(BHEAD0);
3777  AR.wranfcall = 0;
3778  }
3779  wval = (UWORD)(AR.wranfia[AR.wranfcall++]>>(BITSINWORD-1));
3780  return(wval);
3781 }
3782 
3783 /*
3784  Returns a random UWORD in the range (0,...,imax-1)
3785 */
3786 
3787 UWORD iranf(PHEAD UWORD imax)
3788 {
3789  UWORD i;
3790  ULONG x = (LONG)1 << BITSINWORD, xmax = x - x%imax;
3791  while ( ( i = wranf(BHEAD0) ) >= xmax ) {}
3792  return(i%imax);
3793 }
3794 
3795 /*
3796  #] wranf :
3797  #[ PreRandom :
3798 
3799  The random number generator of the preprocessor.
3800  This one is completely different from the execution time generator
3801  random_(number). In the preprocessor we generate a floating point
3802  number in a string according to a distribution.
3803  Currently allowed are:
3804  RANDOM_(log,min,max)
3805  RANDOM_(lin,min,max)
3806  The return value is a string with the floating point number.
3807 */
3808 
3809 UBYTE *PreRandom(UBYTE *s)
3810 {
3811  GETIDENTITY
3812  UBYTE *mode,*mins = 0,*maxs = 0, *outval;
3813  float num;
3814  double minval, maxval, value = 0;
3815  int linlog = -1;
3816  mode = s;
3817  while ( FG.cTable[*s] <= 1 ) s++;
3818  if ( *s == ',' ) { *s = 0; s++; }
3819  mins = s;
3820  while ( *s && *s != ',' ) s++;
3821  if ( *s == ',' ) { *s = 0; s++; }
3822  maxs = s;
3823  while ( *s && *s != ',' ) s++;
3824  if ( *s || *maxs == 0 || *mins == 0 ) {
3825  MesPrint("@Illegal arguments in macro RANDOM_");
3826  Terminate(-1);
3827  }
3828  if ( StrICmp(mode,(UBYTE *)"lin") == 0 ) {
3829  linlog = 0;
3830  }
3831  else if ( StrICmp(mode,(UBYTE *)"log") == 0 ) {
3832  linlog = 1;
3833  }
3834  else {
3835  MesPrint("@Illegal mode argument in macro RANDOM_");
3836  Terminate(-1);
3837  }
3838 
3839  sscanf((char *)mins,"%f",&num); minval = num;
3840  sscanf((char *)maxs,"%f",&num); maxval = num;
3841 
3842  /*
3843  * Note on ParFORM: we should use the same random number on all the
3844  * processes in the complication phase. The random number is generated
3845  * on the master and broadcast to the other processes.
3846  */
3847  {
3848  UWORD x;
3849  double xx;
3850 #ifdef WITHMPI
3851  if ( PF.me == MASTER ) {
3852  x = wranf(BHEAD0);
3853  }
3854  x = (UWORD)PF_BroadcastNumber((LONG)x);
3855 #else
3856  x = wranf(BHEAD0);
3857 #endif
3858  xx = x/pow(2.0,(double)(BITSINWORD-1));
3859  if ( linlog == 0 ) {
3860  value = minval + (maxval-minval)*xx;
3861  }
3862  else if ( linlog == 1 ) {
3863  value = minval * pow(maxval/minval,xx);
3864  }
3865  }
3866 
3867  outval = (UBYTE *)Malloc1(64,"PreRandom");
3868  if ( ABS(value) < 0.00001 || ABS(value) > 1000000. ) {
3869  sprintf((char *)outval,"%e",value);
3870  }
3871  else if ( ABS(value) < 0.0001 ) { sprintf((char *)outval,"%10f",value); }
3872  else if ( ABS(value) < 0.001 ) { sprintf((char *)outval,"%9f",value); }
3873  else if ( ABS(value) < 0.01 ) { sprintf((char *)outval,"%8f",value); }
3874  else if ( ABS(value) < 0.1 ) { sprintf((char *)outval,"%7f",value); }
3875  else if ( ABS(value) < 1. ) { sprintf((char *)outval,"%6f",value); }
3876  else if ( ABS(value) < 10. ) { sprintf((char *)outval,"%5f",value); }
3877  else if ( ABS(value) < 100. ) { sprintf((char *)outval,"%4f",value); }
3878  else if ( ABS(value) < 1000. ) { sprintf((char *)outval,"%3f",value); }
3879  else if ( ABS(value) < 10000. ) { sprintf((char *)outval,"%2f",value); }
3880  else { sprintf((char *)outval,"%1f",value); }
3881  return(outval);
3882 }
3883 
3884 /*
3885  #] PreRandom :
3886  #] Functions :
3887 */
#define PHEAD
Definition: ftypes.h:56
int GetModInverses(WORD m1, WORD m2, WORD *im1, WORD *im2)
Definition: reken.c:1443
WORD CompCoef(WORD *term1, WORD *term2)
Definition: reken.c:3012
int NormalModulus(UWORD *a, WORD *na)
Definition: reken.c:1370
LONG PF_BroadcastNumber(LONG x)
Definition: parallel.c:2096
WORD NextPrime(PHEAD WORD num)
Definition: reken.c:3629
int MakeInverses()
Definition: reken.c:1407
VOID RaisPowCached(PHEAD WORD x, WORD n, UWORD **c, WORD *nc)
Definition: reken.c:1286