My Project
Loading...
Searching...
No Matches
longrat.cc
Go to the documentation of this file.
1/****************************************
2* Computer Algebra System SINGULAR *
3****************************************/
4/*
5* ABSTRACT: computation with long rational numbers (Hubert Grassmann)
6*/
7
8#include "misc/auxiliary.h"
9
10#include "factory/factory.h"
11
12#include "misc/sirandom.h"
13#include "misc/prime.h"
14#include "reporter/reporter.h"
15
16#include "coeffs/coeffs.h"
17#include "coeffs/numbers.h"
18#include "coeffs/rmodulon.h" // ZnmInfo
19#include "coeffs/longrat.h"
20#include "coeffs/shortfl.h"
21#include "coeffs/modulop.h"
22#include "coeffs/mpr_complex.h"
23
24#include <string.h>
25#include <float.h>
26
27// allow inlining only from p_Numbers.h and if ! LDEBUG
28#if defined(DO_LINLINE) && defined(P_NUMBERS_H) && !defined(LDEBUG)
29#define LINLINE static FORCE_INLINE
30#else
31#define LINLINE
32#undef DO_LINLINE
33#endif // DO_LINLINE
34
36LINLINE number nlInit(long i, const coeffs r);
41LINLINE void nlDelete(number *a, const coeffs r);
46LINLINE void nlInpAdd(number &a, number b, const coeffs r);
47LINLINE void nlInpMult(number &a, number b, const coeffs r);
48
49number nlRInit (long i);
50
51
52// number nlInitMPZ(mpz_t m, const coeffs r);
53// void nlMPZ(mpz_t m, number &n, const coeffs r);
54
55void nlNormalize(number &x, const coeffs r);
56
57number nlGcd(number a, number b, const coeffs r);
59number nlNormalizeHelper(number a, number b, const coeffs r); /*special routine !*/
61BOOLEAN nlIsMOne(number a, const coeffs r);
62long nlInt(number &n, const coeffs r);
64
66number nlInvers(number a, const coeffs r);
67number nlDiv(number a, number b, const coeffs r);
69number nlIntDiv(number a, number b, const coeffs r);
70number nlIntMod(number a, number b, const coeffs r);
71void nlPower(number x, int exp, number *lu, const coeffs r);
72const char * nlRead (const char *s, number *a, const coeffs r);
73void nlWrite(number a, const coeffs r);
74
76
77#ifdef LDEBUG
78BOOLEAN nlDBTest(number a, const char *f, const int l);
79#endif
80
81nMapFunc nlSetMap(const coeffs src, const coeffs dst);
82
83// in-place operations
84void nlInpIntDiv(number &a, number b, const coeffs r);
85
86#ifdef LDEBUG
87#define nlTest(a, r) nlDBTest(a,__FILE__,__LINE__, r)
88BOOLEAN nlDBTest(number a, const char *f,int l, const coeffs r);
89#else
90#define nlTest(a, r) do {} while (0)
91#endif
92
93
94// 64 bit version:
95//#if SIZEOF_LONG == 8
96#if 0
97#define MAX_NUM_SIZE 60
98#define POW_2_28 (1L<<60)
99#define POW_2_28_32 (1L<<28)
100#define LONG long
101#else
102#define MAX_NUM_SIZE 28
103#define POW_2_28 (1L<<28)
104#define POW_2_28_32 (1L<<28)
105#define LONG int
106#endif
107
108
109static inline number nlShort3(number x) // assume x->s==3
110{
111 assume(x->s==3);
112 if (mpz_sgn1(x->z)==0)
113 {
114 mpz_clear(x->z);
116 return INT_TO_SR(0);
117 }
118 if (mpz_size1(x->z)<=MP_SMALL)
119 {
120 LONG ui=mpz_get_si(x->z);
121 if ((((ui<<3)>>3)==ui)
122 && (mpz_cmp_si(x->z,(long)ui)==0))
123 {
124 mpz_clear(x->z);
126 return INT_TO_SR(ui);
127 }
128 }
129 return x;
130}
131
132#ifndef LONGRAT_CC
133#define LONGRAT_CC
134
135#ifndef BYTES_PER_MP_LIMB
136#define BYTES_PER_MP_LIMB sizeof(mp_limb_t)
137#endif
138
139//#define SR_HDL(A) ((long)(A))
140/*#define SR_INT 1L*/
141/*#define INT_TO_SR(INT) ((number) (((long)INT << 2) + SR_INT))*/
142// #define SR_TO_INT(SR) (((long)SR) >> 2)
143
144#define MP_SMALL 1
145//#define mpz_isNeg(A) (mpz_sgn1(A)<0)
146#define mpz_isNeg(A) ((A)->_mp_size<0)
147#define mpz_limb_size(A) ((A)->_mp_size)
148#define mpz_limb_d(A) ((A)->_mp_d)
149
150void _nlDelete_NoImm(number *a);
151
152/***************************************************************
153 *
154 * Routines which are never inlined by p_Numbers.h
155 *
156 *******************************************************************/
157#ifndef P_NUMBERS_H
158
160{
161 return nlShort3(x);
162}
163
165{
166 number z = ALLOC_RNUMBER();
167 z->s = 3;
168 #ifdef LDEBUG
169 z->debug=123456;
170 #endif
171 mpz_init_set(z->z, m);
172 z=nlShort3(z);
173 return z;
174}
175
176#if (__GNU_MP_VERSION*10+__GNU_MP_VERSION_MINOR < 31)
177void mpz_mul_si (mpz_ptr r, mpz_srcptr s, long int si)
178{
179 if (si>=0)
180 mpz_mul_ui(r,s,si);
181 else
182 {
183 mpz_mul_ui(r,s,-si);
184 mpz_neg(r,r);
185 }
186}
187#endif
188
189static number nlMapP(number from, const coeffs src, const coeffs dst)
190{
191 assume( getCoeffType(src) == n_Zp );
192
193 number to = nlInit(npInt(from,src), dst); // FIXME? TODO? // extern long npInt (number &n, const coeffs r);
194
195 return to;
196}
197
198static number nlMapLongR(number from, const coeffs src, const coeffs dst);
199static number nlMapR(number from, const coeffs src, const coeffs dst);
200
201
202/*2
203* convert from a GMP integer
204*/
205static inline number nlMapGMP(number from, const coeffs /*src*/, const coeffs dst)
206{
207 return nlInitMPZ((mpz_ptr)from,dst);
208}
209
210number nlMapZ(number from, const coeffs /*src*/, const coeffs dst)
211{
212 if (SR_HDL(from) & SR_INT)
213 {
214 return from;
215 }
216 return nlInitMPZ((mpz_ptr)from,dst);
217}
218
219/*2
220* convert from an machine long
221*/
222number nlMapMachineInt(number from, const coeffs /*src*/, const coeffs /*dst*/)
223{
225#if defined(LDEBUG)
226 z->debug=123456;
227#endif
228 mpz_init_set_ui(z->z,(unsigned long) from);
229 z->s = 3;
230 z=nlShort3(z);
231 return z;
232}
233
234#ifdef LDEBUG
235BOOLEAN nlDBTest(number a, const char *f,const int l, const coeffs /*r*/)
236{
237 if (a==NULL)
238 {
239 Print("!!longrat: NULL in %s:%d\n",f,l);
240 return FALSE;
241 }
242 //if ((int)a==1) Print("!! 0x1 as number ? %s %d\n",f,l);
243 if ((((long)a)&3L)==3L)
244 {
245 Print(" !!longrat:ptr(3) in %s:%d\n",f,l);
246 return FALSE;
247 }
248 if ((((long)a)&3L)==1L)
249 {
250 if (((((LONG)(long)a)<<1)>>1)!=((LONG)(long)a))
251 {
252 Print(" !!longrat:arith:%lx in %s:%d\n",(long)a, f,l);
253 return FALSE;
254 }
255 return TRUE;
256 }
257 /* TODO: If next line is active, then computations in algebraic field
258 extensions over Q will throw a lot of assume violations although
259 everything is computed correctly and no seg fault appears.
260 Maybe the test is not appropriate in this case. */
261 omCheckIf(omCheckAddrSize(a,sizeof(*a)), return FALSE);
262 if (a->debug!=123456)
263 {
264 Print("!!longrat:debug:%d in %s:%d\n",a->debug,f,l);
265 a->debug=123456;
266 return FALSE;
267 }
268 if ((a->s<0)||(a->s>4))
269 {
270 Print("!!longrat:s=%d in %s:%d\n",a->s,f,l);
271 return FALSE;
272 }
273 /* TODO: If next line is active, then computations in algebraic field
274 extensions over Q will throw a lot of assume violations although
275 everything is computed correctly and no seg fault appears.
276 Maybe the test is not appropriate in this case. */
277 //omCheckAddrSize(a->z[0]._mp_d,a->z[0]._mp_alloc*BYTES_PER_MP_LIMB);
278 if (a->z[0]._mp_alloc==0)
279 Print("!!longrat:z->alloc=0 in %s:%d\n",f,l);
280
281 if (a->s<2)
282 {
283 if ((a->n[0]._mp_d[0]==0)&&(a->n[0]._mp_alloc<=1))
284 {
285 Print("!!longrat: n==0 in %s:%d\n",f,l);
286 return FALSE;
287 }
288 /* TODO: If next line is active, then computations in algebraic field
289 extensions over Q will throw a lot of assume violations although
290 everything is computed correctly and no seg fault appears.
291 Maybe the test is not appropriate in this case. */
292 //omCheckIf(omCheckAddrSize(a->n[0]._mp_d,a->n[0]._mp_alloc*BYTES_PER_MP_LIMB), return FALSE);
293 if (a->z[0]._mp_alloc==0)
294 Print("!!longrat:n->alloc=0 in %s:%d\n",f,l);
295 if ((mpz_size1(a->n) ==1) && (mpz_cmp_si(a->n,1L)==0))
296 {
297 Print("!!longrat:integer as rational in %s:%d\n",f,l);
298 mpz_clear(a->n); a->s=3;
299 return FALSE;
300 }
301 else if (mpz_isNeg(a->n))
302 {
303 Print("!!longrat:div. by negative in %s:%d\n",f,l);
304 mpz_neg(a->z,a->z);
305 mpz_neg(a->n,a->n);
306 return FALSE;
307 }
308 return TRUE;
309 }
310 //if (a->s==2)
311 //{
312 // Print("!!longrat:s=2 in %s:%d\n",f,l);
313 // return FALSE;
314 //}
315 if (mpz_size1(a->z)>MP_SMALL) return TRUE;
316 LONG ui=(LONG)mpz_get_si(a->z);
317 if ((((ui<<3)>>3)==ui)
318 && (mpz_cmp_si(a->z,(long)ui)==0))
319 {
320 Print("!!longrat:im int %d in %s:%d\n",ui,f,l);
321 return FALSE;
322 }
323 return TRUE;
324}
325#endif
326
328{
329 if (setChar) setCharacteristic( 0 );
330
332 if ( SR_HDL(n) & SR_INT )
333 {
334 long nn=SR_TO_INT(n);
335 term = nn;
336 }
337 else
338 {
339 if ( n->s == 3 )
340 {
341 mpz_t dummy;
342 long lz=mpz_get_si(n->z);
343 if (mpz_cmp_si(n->z,lz)==0) term=lz;
344 else
345 {
346 mpz_init_set( dummy,n->z );
347 term = make_cf( dummy );
348 }
349 }
350 else
351 {
352 // assume s==0 or s==1
353 mpz_t num, den;
355 mpz_init_set( num, n->z );
356 mpz_init_set( den, n->n );
357 term = make_cf( num, den, ( n->s != 1 ));
358 }
359 }
360 return term;
361}
362
363number nlRInit (long i);
364
366{
367 if (f.isImm())
368 {
369 return nlInit(f.intval(),r);
370 }
371 else
372 {
373 number z = ALLOC_RNUMBER();
374#if defined(LDEBUG)
375 z->debug=123456;
376#endif
377 gmp_numerator( f, z->z );
378 if ( f.den().isOne() )
379 {
380 z->s = 3;
381 z=nlShort3(z);
382 }
383 else
384 {
385 gmp_denominator( f, z->n );
386 z->s = 1;
387 }
388 return z;
389 }
390}
391
392static number nlMapR(number from, const coeffs src, const coeffs dst)
393{
394 assume( getCoeffType(src) == n_R );
395
396 double f=nrFloat(from); // FIXME? TODO? // extern float nrFloat(number n);
397 if (f==0.0) return INT_TO_SR(0);
398 int f_sign=1;
399 if (f<0.0)
400 {
401 f_sign=-1;
402 f=-f;
403 }
404 int i=0;
405 mpz_t h1;
407 while((FLT_RADIX*f) < DBL_MAX && i<DBL_MANT_DIG)
408 {
409 f*=FLT_RADIX;
411 i++;
412 }
413 number re=nlRInit(1);
414 mpz_set_d(re->z,f);
415 memcpy(&(re->n),&h1,sizeof(h1));
416 re->s=0; /* not normalized */
417 if(f_sign==-1) re=nlNeg(re,dst);
419 return re;
420}
421
422static number nlMapR_BI(number from, const coeffs src, const coeffs dst)
423{
424 assume( getCoeffType(src) == n_R );
425
426 double f=nrFloat(from); // FIXME? TODO? // extern float nrFloat(number n);
427 if (f==0.0) return INT_TO_SR(0);
428 long l=long(f);
429 return nlInit(l,dst);
430}
431
432static number nlMapLongR(number from, const coeffs src, const coeffs dst)
433{
434 assume( getCoeffType(src) == n_long_R );
435
436 gmp_float *ff=(gmp_float*)from;
437 mpf_t *f=ff->_mpfp();
438 number res;
439 mpz_ptr dest,ndest;
440 int size, i,negative;
441 int e,al,bl;
442 mp_ptr qp,dd,nn;
443
444 size = (*f)[0]._mp_size;
445 if (size == 0)
446 return INT_TO_SR(0);
447 if(size<0)
448 {
449 negative = 1;
450 size = -size;
451 }
452 else
453 negative = 0;
454
455 qp = (*f)[0]._mp_d;
456 while(qp[0]==0)
457 {
458 qp++;
459 size--;
460 }
461
462 e=(*f)[0]._mp_exp-size;
463 res = ALLOC_RNUMBER();
464#if defined(LDEBUG)
465 res->debug=123456;
466#endif
467 dest = res->z;
468
469 void* (*allocfunc) (size_t);
471 if (e<0)
472 {
473 al = dest->_mp_size = size;
474 if (al<2) al = 2;
475 dd = (mp_ptr)allocfunc(sizeof(mp_limb_t)*al);
476 for (i=0;i<size;i++) dd[i] = qp[i];
477 bl = 1-e;
478 nn = (mp_ptr)allocfunc(sizeof(mp_limb_t)*bl);
479 memset(nn,0,sizeof(mp_limb_t)*bl);
480 nn[bl-1] = 1;
481 ndest = res->n;
482 ndest->_mp_d = nn;
483 ndest->_mp_alloc = ndest->_mp_size = bl;
484 res->s = 0;
485 }
486 else
487 {
488 al = dest->_mp_size = size+e;
489 if (al<2) al = 2;
490 dd = (mp_ptr)allocfunc(sizeof(mp_limb_t)*al);
491 memset(dd,0,sizeof(mp_limb_t)*al);
492 for (i=0;i<size;i++) dd[i+e] = qp[i];
493 for (i=0;i<e;i++) dd[i] = 0;
494 res->s = 3;
495 }
496
497 dest->_mp_d = dd;
498 dest->_mp_alloc = al;
499 if (negative) mpz_neg(dest,dest);
500
501 if (res->s==0)
503 else if (mpz_size1(res->z)<=MP_SMALL)
504 {
505 // res is new, res->ref is 1
507 }
508 nlTest(res, dst);
509 return res;
510}
511
512static number nlMapLongR_BI(number from, const coeffs src, const coeffs dst)
513{
514 assume( getCoeffType(src) == n_long_R );
515
516 gmp_float *ff=(gmp_float*)from;
517 if (mpf_fits_slong_p(ff->t))
518 {
519 long l=mpf_get_si(ff->t);
520 return nlInit(l,dst);
521 }
522 char *out=floatToStr(*(gmp_float*)from, src->float_len);
523 char *p=strchr(out,'.');
524 *p='\0';
525 number res;
526 res = ALLOC_RNUMBER();
527#if defined(LDEBUG)
528 res->debug=123456;
529#endif
530 res->s=3;
531 mpz_init(res->z);
532 if (out[0]=='-')
533 {
534 mpz_set_str(res->z,out+1,10);
535 res=nlNeg(res,dst);
536 }
537 else
538 {
539 mpz_set_str(res->z,out,10);
540 }
541 omFree( (void *)out );
542 return res;
543}
544
545static number nlMapC(number from, const coeffs src, const coeffs dst)
546{
547 assume( getCoeffType(src) == n_long_C );
548 if ( ! ((gmp_complex*)from)->imag().isZero() )
549 return INT_TO_SR(0);
550
551 if (dst->is_field==FALSE) /* ->ZZ */
552 {
553 char *s=floatToStr(((gmp_complex*)from)->real(),src->float_len);
554 mpz_t z;
555 mpz_init(z);
556 char *ss=nEatLong(s,z);
557 if (*ss=='\0')
558 {
559 omFree(s);
560 number n=nlInitMPZ(z,dst);
561 mpz_clear(z);
562 return n;
563 }
564 omFree(s);
565 mpz_clear(z);
566 WarnS("conversion problem in CC -> ZZ mapping");
567 return INT_TO_SR(0);
568 }
569
570 mpf_t *f = ((gmp_complex*)from)->real()._mpfp();
571
572 number res;
573 mpz_ptr dest,ndest;
574 int size, i,negative;
575 int e,al,bl;
576 mp_ptr qp,dd,nn;
577
578 size = (*f)[0]._mp_size;
579 if (size == 0)
580 return INT_TO_SR(0);
581 if(size<0)
582 {
583 negative = 1;
584 size = -size;
585 }
586 else
587 negative = 0;
588
589 qp = (*f)[0]._mp_d;
590 while(qp[0]==0)
591 {
592 qp++;
593 size--;
594 }
595
596 e=(*f)[0]._mp_exp-size;
597 res = ALLOC_RNUMBER();
598#if defined(LDEBUG)
599 res->debug=123456;
600#endif
601 dest = res->z;
602
603 void* (*allocfunc) (size_t);
605 if (e<0)
606 {
607 al = dest->_mp_size = size;
608 if (al<2) al = 2;
609 dd = (mp_ptr)allocfunc(sizeof(mp_limb_t)*al);
610 for (i=0;i<size;i++) dd[i] = qp[i];
611 bl = 1-e;
612 nn = (mp_ptr)allocfunc(sizeof(mp_limb_t)*bl);
613 memset(nn,0,sizeof(mp_limb_t)*bl);
614 nn[bl-1] = 1;
615 ndest = res->n;
616 ndest->_mp_d = nn;
617 ndest->_mp_alloc = ndest->_mp_size = bl;
618 res->s = 0;
619 }
620 else
621 {
622 al = dest->_mp_size = size+e;
623 if (al<2) al = 2;
624 dd = (mp_ptr)allocfunc(sizeof(mp_limb_t)*al);
625 memset(dd,0,sizeof(mp_limb_t)*al);
626 for (i=0;i<size;i++) dd[i+e] = qp[i];
627 for (i=0;i<e;i++) dd[i] = 0;
628 res->s = 3;
629 }
630
631 dest->_mp_d = dd;
632 dest->_mp_alloc = al;
633 if (negative) mpz_neg(dest,dest);
634
635 if (res->s==0)
637 else if (mpz_size1(res->z)<=MP_SMALL)
638 {
639 // res is new, res->ref is 1
641 }
642 nlTest(res, dst);
643 return res;
644}
645
646//static number nlMapLongR(number from)
647//{
648// gmp_float *ff=(gmp_float*)from;
649// const mpf_t *f=ff->mpfp();
650// int f_size=ABS((*f)[0]._mp_size);
651// if (f_size==0)
652// return nlInit(0);
653// int f_sign=1;
654// number work=ngcCopy(from);
655// if (!ngcGreaterZero(work))
656// {
657// f_sign=-1;
658// work=ngcNeg(work);
659// }
660// int i=0;
661// mpz_t h1;
662// mpz_init_set_ui(h1,1);
663// while((FLT_RADIX*f) < DBL_MAX && i<DBL_MANT_DIG)
664// {
665// f*=FLT_RADIX;
666// mpz_mul_ui(h1,h1,FLT_RADIX);
667// i++;
668// }
669// number r=nlRInit(1);
670// mpz_set_d(&(r->z),f);
671// memcpy(&(r->n),&h1,sizeof(h1));
672// r->s=0; /* not normalized */
673// nlNormalize(r);
674// return r;
675//
676//
677// number r=nlRInit(1);
678// int f_shift=f_size+(*f)[0]._mp_exp;
679// if ( f_shift > 0)
680// {
681// r->s=0;
682// mpz_init(&r->n);
683// mpz_setbit(&r->n,f_shift*BYTES_PER_MP_LIMB*8);
684// mpz_setbit(&r->z,f_size*BYTES_PER_MP_LIMB*8-1);
685// // now r->z has enough space
686// memcpy(mpz_limb_d(&r->z),((*f)[0]._mp_d),f_size*BYTES_PER_MP_LIMB);
687// nlNormalize(r);
688// }
689// else
690// {
691// r->s=3;
692// if (f_shift==0)
693// {
694// mpz_setbit(&r->z,f_size*BYTES_PER_MP_LIMB*8-1);
695// // now r->z has enough space
696// memcpy(mpz_limb_d(&r->z),((*f)[0]._mp_d),f_size*BYTES_PER_MP_LIMB);
697// }
698// else /* f_shift < 0 */
699// {
700// mpz_setbit(&r->z,(f_size-f_shift)*BYTES_PER_MP_LIMB*8-1);
701// // now r->z has enough space
702// memcpy(mpz_limb_d(&r->z)-f_shift,((*f)[0]._mp_d),
703// f_size*BYTES_PER_MP_LIMB);
704// }
705// }
706// if ((*f)[0]._mp_size<0);
707// r=nlNeg(r);
708// return r;
709//}
710
711int nlSize(number a, const coeffs)
712{
713 if (a==INT_TO_SR(0))
714 return 0; /* rational 0*/
715 if (SR_HDL(a) & SR_INT)
716 return 1; /* immediate int */
717 int s=a->z[0]._mp_alloc;
718// while ((s>0) &&(a->z._mp_d[s]==0L)) s--;
719//#if SIZEOF_LONG == 8
720// if (a->z._mp_d[s] < (unsigned long)0x100000000L) s=s*2-1;
721// else s *=2;
722//#endif
723// s++;
724 if (a->s<2)
725 {
726 int d=a->n[0]._mp_alloc;
727// while ((d>0) && (a->n._mp_d[d]==0L)) d--;
728//#if SIZEOF_LONG == 8
729// if (a->n._mp_d[d] < (unsigned long)0x100000000L) d=d*2-1;
730// else d *=2;
731//#endif
732 s+=d;
733 }
734 return s;
735}
736
737/*2
738* convert number to int
739*/
740long nlInt(number &i, const coeffs r)
741{
742 nlTest(i, r);
743 nlNormalize(i,r);
744 if (SR_HDL(i) & SR_INT)
745 {
746 return SR_TO_INT(i);
747 }
748 if (i->s==3)
749 {
750 if(mpz_size1(i->z)>MP_SMALL) return 0;
751 long ul=mpz_get_si(i->z);
752 if (mpz_cmp_si(i->z,ul)!=0) return 0;
753 return ul;
754 }
755 mpz_t tmp;
756 long ul;
757 mpz_init(tmp);
758 mpz_tdiv_q(tmp,i->z,i->n);
759 if(mpz_size1(tmp)>MP_SMALL) ul=0;
760 else
761 {
763 if (mpz_cmp_si(tmp,ul)!=0) ul=0;
764 }
765 mpz_clear(tmp);
766 return ul;
767}
768
769/*2
770* convert number to bigint
771*/
773{
774 nlTest(i, r);
775 nlNormalize(i,r);
776 if (SR_HDL(i) & SR_INT) return (i);
777 if (i->s==3)
778 {
779 return nlCopy(i,r);
780 }
781 number tmp=nlRInit(1);
782 mpz_tdiv_q(tmp->z,i->z,i->n);
784 return tmp;
785}
786
787/*
788* 1/a
789*/
791{
792 nlTest(a, r);
793 number n;
794 if (SR_HDL(a) & SR_INT)
795 {
796 if ((a==INT_TO_SR(1L)) || (a==INT_TO_SR(-1L)))
797 {
798 return a;
799 }
800 if (nlIsZero(a,r))
801 {
803 return INT_TO_SR(0);
804 }
805 n=ALLOC_RNUMBER();
806#if defined(LDEBUG)
807 n->debug=123456;
808#endif
809 n->s=1;
810 if (((long)a)>0L)
811 {
812 mpz_init_set_ui(n->z,1L);
813 mpz_init_set_si(n->n,(long)SR_TO_INT(a));
814 }
815 else
816 {
817 mpz_init_set_si(n->z,-1L);
818 mpz_init_set_si(n->n,(long)-SR_TO_INT(a));
819 }
820 nlTest(n, r);
821 return n;
822 }
823 n=ALLOC_RNUMBER();
824#if defined(LDEBUG)
825 n->debug=123456;
826#endif
827 {
828 mpz_init_set(n->n,a->z);
829 switch (a->s)
830 {
831 case 0:
832 case 1:
833 n->s=a->s;
834 mpz_init_set(n->z,a->n);
835 if (mpz_isNeg(n->n)) /* && n->s<2*/
836 {
837 mpz_neg(n->z,n->z);
838 mpz_neg(n->n,n->n);
839 }
840 if (mpz_cmp_ui(n->n,1L)==0)
841 {
842 mpz_clear(n->n);
843 n->s=3;
844 n=nlShort3(n);
845 }
846 break;
847 case 3:
848 // i.e. |a| > 2^...
849 n->s=1;
850 if (mpz_isNeg(n->n)) /* && n->s<2*/
851 {
852 mpz_neg(n->n,n->n);
853 mpz_init_set_si(n->z,-1L);
854 }
855 else
856 {
857 mpz_init_set_ui(n->z,1L);
858 }
859 break;
860 }
861 }
862 nlTest(n, r);
863 return n;
864}
865
866
867/*2
868* u := a / b in Z, if b | a (else undefined)
869*/
871{
872 if (b==INT_TO_SR(0))
873 {
875 return INT_TO_SR(0);
876 }
877 number u;
878 if (SR_HDL(a) & SR_HDL(b) & SR_INT)
879 {
880 /* the small int -(1<<28) divided by -1 is the large int (1<<28) */
881 if ((a==INT_TO_SR(-(POW_2_28)))&&(b==INT_TO_SR(-1L)))
882 {
883 return nlRInit(POW_2_28);
884 }
885 long aa=SR_TO_INT(a);
886 long bb=SR_TO_INT(b);
887 return INT_TO_SR(aa/bb);
888 }
889 number aa=NULL;
890 number bb=NULL;
891 if (SR_HDL(a) & SR_INT)
892 {
893 aa=nlRInit(SR_TO_INT(a));
894 a=aa;
895 }
896 if (SR_HDL(b) & SR_INT)
897 {
899 b=bb;
900 }
901 u=ALLOC_RNUMBER();
902#if defined(LDEBUG)
903 u->debug=123456;
904#endif
905 mpz_init(u->z);
906 /* u=a/b */
907 u->s = 3;
908 assume(a->s==3);
909 assume(b->s==3);
910 mpz_divexact(u->z,a->z,b->z);
911 if (aa!=NULL)
912 {
913 mpz_clear(aa->z);
914#if defined(LDEBUG)
915 aa->debug=654324;
916#endif
917 FREE_RNUMBER(aa); // omFreeBin((void *)aa, rnumber_bin);
918 }
919 if (bb!=NULL)
920 {
921 mpz_clear(bb->z);
922#if defined(LDEBUG)
923 bb->debug=654324;
924#endif
925 FREE_RNUMBER(bb); // omFreeBin((void *)bb, rnumber_bin);
926 }
927 u=nlShort3(u);
928 nlTest(u, r);
929 return u;
930}
931
932/*2
933* u := a / b in Z
934*/
936{
937 if (b==INT_TO_SR(0))
938 {
940 return INT_TO_SR(0);
941 }
942 number u;
943 if (SR_HDL(a) & SR_HDL(b) & SR_INT)
944 {
945 /* the small int -(1<<28) divided by -1 is the large int (1<<28) */
946 if ((a==INT_TO_SR(-(POW_2_28)))&&(b==INT_TO_SR(-1L)))
947 {
948 return nlRInit(POW_2_28);
949 }
950 LONG aa=SR_TO_INT(a);
952 LONG rr=aa%bb;
953 if (rr<0) rr+=ABS(bb);
954 LONG cc=(aa-rr)/bb;
955 return INT_TO_SR(cc);
956 }
957 number aa=NULL;
958 if (SR_HDL(a) & SR_INT)
959 {
960 /* the small int -(1<<28) divided by 2^28 is 1 */
961 if (a==INT_TO_SR(-(POW_2_28)))
962 {
963 if(mpz_cmp_si(b->z,(POW_2_28))==0)
964 {
965 return INT_TO_SR(-1);
966 }
967 }
968 aa=nlRInit(SR_TO_INT(a));
969 a=aa;
970 }
971 number bb=NULL;
972 if (SR_HDL(b) & SR_INT)
973 {
975 b=bb;
976 }
977 u=ALLOC_RNUMBER();
978#if defined(LDEBUG)
979 u->debug=123456;
980#endif
981 assume(a->s==3);
982 assume(b->s==3);
983 /* u=u/b */
984 mpz_t rr;
985 mpz_init(rr);
986 mpz_mod(rr,a->z,b->z);
987 u->s = 3;
988 mpz_init(u->z);
989 mpz_sub(u->z,a->z,rr);
990 mpz_clear(rr);
991 mpz_divexact(u->z,u->z,b->z);
992 if (aa!=NULL)
993 {
994 mpz_clear(aa->z);
995#if defined(LDEBUG)
996 aa->debug=654324;
997#endif
999 }
1000 if (bb!=NULL)
1001 {
1002 mpz_clear(bb->z);
1003#if defined(LDEBUG)
1004 bb->debug=654324;
1005#endif
1007 }
1008 u=nlShort3(u);
1009 nlTest(u,r);
1010 return u;
1011}
1012
1013/*2
1014* u := a mod b in Z, u>=0
1015*/
1017{
1018 if (b==INT_TO_SR(0))
1019 {
1021 return INT_TO_SR(0);
1022 }
1023 if (a==INT_TO_SR(0))
1024 return INT_TO_SR(0);
1025 number u;
1026 if (SR_HDL(a) & SR_HDL(b) & SR_INT)
1027 {
1028 LONG aa=SR_TO_INT(a);
1029 LONG bb=SR_TO_INT(b);
1030 LONG c=aa % bb;
1031 if (c<0) c+=ABS(bb);
1032 return INT_TO_SR(c);
1033 }
1034 if (SR_HDL(a) & SR_INT)
1035 {
1036 LONG ai=SR_TO_INT(a);
1037 mpz_t aa;
1039 u=ALLOC_RNUMBER();
1040#if defined(LDEBUG)
1041 u->debug=123456;
1042#endif
1043 u->s = 3;
1044 mpz_init(u->z);
1045 mpz_mod(u->z, aa, b->z);
1046 mpz_clear(aa);
1047 u=nlShort3(u);
1048 nlTest(u,r);
1049 return u;
1050 }
1051 number bb=NULL;
1052 if (SR_HDL(b) & SR_INT)
1053 {
1055 b=bb;
1056 }
1057 u=ALLOC_RNUMBER();
1058#if defined(LDEBUG)
1059 u->debug=123456;
1060#endif
1061 mpz_init(u->z);
1062 u->s = 3;
1063 mpz_mod(u->z, a->z, b->z);
1064 if (bb!=NULL)
1065 {
1066 mpz_clear(bb->z);
1067#if defined(LDEBUG)
1068 bb->debug=654324;
1069#endif
1071 }
1072 u=nlShort3(u);
1073 nlTest(u,r);
1074 return u;
1075}
1076
1078{
1079 if (SR_HDL(a) & SR_HDL(b) & SR_INT)
1080 {
1081 return ((SR_TO_INT(a) % SR_TO_INT(b))==0);
1082 }
1083 if (SR_HDL(b) & SR_INT)
1084 {
1085 return (mpz_divisible_ui_p(a->z,SR_TO_INT(b))!=0);
1086 }
1087 if (SR_HDL(a) & SR_INT) return FALSE;
1088 return mpz_divisible_p(a->z, b->z) != 0;
1089}
1090
1092{
1093 if (nlDivBy(a, b, r))
1094 {
1095 if (nlDivBy(b, a, r)) return 2;
1096 return -1;
1097 }
1098 if (nlDivBy(b, a, r)) return 1;
1099 return 0;
1100}
1101
1103{
1104 if (nlGreaterZero(n,cf)) return INT_TO_SR(1);
1105 else return INT_TO_SR(-1);
1106}
1107
1109{
1110 long ch = r->cfInt(c, r);
1111 int p=IsPrime(ch);
1112 coeffs rr=NULL;
1113 if (((long)p)==ch)
1114 {
1115 rr = nInitChar(n_Zp,(void*)ch);
1116 }
1117 else
1118 {
1119 mpz_t dummy;
1121 ZnmInfo info;
1122 info.base = dummy;
1123 info.exp = (unsigned long) 1;
1124 rr = nInitChar(n_Zn, (void*)&info);
1126 }
1127 return(rr);
1128}
1129
1130
1132{
1133 return ((SR_HDL(a) & SR_INT) && (ABS(SR_TO_INT(a))==1));
1134}
1135
1136
1137/*2
1138* u := a / b
1139*/
1141{
1142 if (nlIsZero(b,r))
1143 {
1145 return INT_TO_SR(0);
1146 }
1147 number u;
1148// ---------- short / short ------------------------------------
1149 if (SR_HDL(a) & SR_HDL(b) & SR_INT)
1150 {
1151 LONG i=SR_TO_INT(a);
1152 LONG j=SR_TO_INT(b);
1153 if (j==1L) return a;
1154 if ((i==-POW_2_28) && (j== -1L))
1155 {
1156 return nlRInit(POW_2_28);
1157 }
1158 LONG r=i%j;
1159 if (r==0)
1160 {
1161 return INT_TO_SR(i/j);
1162 }
1163 u=ALLOC_RNUMBER();
1164 u->s=0;
1165 #if defined(LDEBUG)
1166 u->debug=123456;
1167 #endif
1168 mpz_init_set_si(u->z,(long)i);
1169 mpz_init_set_si(u->n,(long)j);
1170 }
1171 else
1172 {
1173 u=ALLOC_RNUMBER();
1174 u->s=0;
1175 #if defined(LDEBUG)
1176 u->debug=123456;
1177 #endif
1178 mpz_init(u->z);
1179// ---------- short / long ------------------------------------
1180 if (SR_HDL(a) & SR_INT)
1181 {
1182 // short a / (z/n) -> (a*n)/z
1183 if (b->s<2)
1184 {
1185 mpz_mul_si(u->z,b->n,SR_TO_INT(a));
1186 }
1187 else
1188 // short a / long z -> a/z
1189 {
1190 mpz_set_si(u->z,SR_TO_INT(a));
1191 }
1192 if (mpz_cmp(u->z,b->z)==0)
1193 {
1194 mpz_clear(u->z);
1195 FREE_RNUMBER(u);
1196 return INT_TO_SR(1);
1197 }
1198 mpz_init_set(u->n,b->z);
1199 }
1200// ---------- long / short ------------------------------------
1201 else if (SR_HDL(b) & SR_INT)
1202 {
1203 mpz_set(u->z,a->z);
1204 // (z/n) / b -> z/(n*b)
1205 if (a->s<2)
1206 {
1207 mpz_init_set(u->n,a->n);
1208 if (((long)b)>0L)
1209 mpz_mul_ui(u->n,u->n,SR_TO_INT(b));
1210 else
1211 {
1212 mpz_mul_ui(u->n,u->n,-SR_TO_INT(b));
1213 mpz_neg(u->z,u->z);
1214 }
1215 }
1216 else
1217 // long z / short b -> z/b
1218 {
1219 //mpz_set(u->z,a->z);
1221 }
1222 }
1223// ---------- long / long ------------------------------------
1224 else
1225 {
1226 mpz_set(u->z,a->z);
1227 mpz_init_set(u->n,b->z);
1228 if (a->s<2) mpz_mul(u->n,u->n,a->n);
1229 if (b->s<2) mpz_mul(u->z,u->z,b->n);
1230 }
1231 }
1232 if (mpz_isNeg(u->n))
1233 {
1234 mpz_neg(u->z,u->z);
1235 mpz_neg(u->n,u->n);
1236 }
1237 if (mpz_cmp_si(u->n,1L)==0)
1238 {
1239 mpz_clear(u->n);
1240 u->s=3;
1241 u=nlShort3(u);
1242 }
1243 nlTest(u, r);
1244 return u;
1245}
1246
1247/*2
1248* u:= x ^ exp
1249*/
1250void nlPower (number x,int exp,number * u, const coeffs r)
1251{
1252 *u = INT_TO_SR(0); // 0^e, e!=0
1253 if (exp==0)
1254 *u= INT_TO_SR(1);
1255 else if (!nlIsZero(x,r))
1256 {
1257 nlTest(x, r);
1258 number aa=NULL;
1259 if (SR_HDL(x) & SR_INT)
1260 {
1262 x=aa;
1263 }
1264 else if (x->s==0)
1265 nlNormalize(x,r);
1266 *u=ALLOC_RNUMBER();
1267#if defined(LDEBUG)
1268 (*u)->debug=123456;
1269#endif
1270 mpz_init((*u)->z);
1271 mpz_pow_ui((*u)->z,x->z,(unsigned long)exp);
1272 if (x->s<2)
1273 {
1274 if (mpz_cmp_si(x->n,1L)==0)
1275 {
1276 x->s=3;
1277 mpz_clear(x->n);
1278 }
1279 else
1280 {
1281 mpz_init((*u)->n);
1282 mpz_pow_ui((*u)->n,x->n,(unsigned long)exp);
1283 }
1284 }
1285 (*u)->s = x->s;
1286 if ((*u)->s==3) *u=nlShort3(*u);
1287 if (aa!=NULL)
1288 {
1289 mpz_clear(aa->z);
1291 }
1292 }
1293#ifdef LDEBUG
1294 if (exp<0) Print("nlPower: neg. exp. %d\n",exp);
1295 nlTest(*u, r);
1296#endif
1297}
1298
1299
1300/*2
1301* za >= 0 ?
1302*/
1304{
1305 nlTest(a, r);
1306 if (SR_HDL(a) & SR_INT) return SR_HDL(a)>1L /* represents number(0) */;
1307 return (!mpz_isNeg(a->z));
1308}
1309
1310/*2
1311* a > b ?
1312*/
1314{
1315 nlTest(a, r);
1316 nlTest(b, r);
1317 number re;
1318 BOOLEAN rr;
1319 re=nlSub(a,b,r);
1320 rr=(!nlIsZero(re,r)) && (nlGreaterZero(re,r));
1321 nlDelete(&re,r);
1322 return rr;
1323}
1324
1325/*2
1326* a == -1 ?
1327*/
1329{
1330#ifdef LDEBUG
1331 if (a==NULL) return FALSE;
1332 nlTest(a, r);
1333#endif
1334 return (a==INT_TO_SR(-1L));
1335}
1336
1337/*2
1338* result =gcd(a,b)
1339*/
1341{
1342 number result;
1343 nlTest(a, r);
1344 nlTest(b, r);
1345 //nlNormalize(a);
1346 //nlNormalize(b);
1347 if ((a==INT_TO_SR(1L))||(a==INT_TO_SR(-1L))
1348 || (b==INT_TO_SR(1L))||(b==INT_TO_SR(-1L)))
1349 return INT_TO_SR(1L);
1350 if (a==INT_TO_SR(0)) /* gcd(0,b) ->b */
1351 return nlCopy(b,r);
1352 if (b==INT_TO_SR(0)) /* gcd(a,0) -> a */
1353 return nlCopy(a,r);
1354 if (SR_HDL(a) & SR_HDL(b) & SR_INT)
1355 {
1356 long i=SR_TO_INT(a);
1357 long j=SR_TO_INT(b);
1358 long l;
1359 i=ABS(i);
1360 j=ABS(j);
1361 do
1362 {
1363 l=i%j;
1364 i=j;
1365 j=l;
1366 } while (l!=0L);
1367 if (i==POW_2_28)
1369 else
1371 nlTest(result,r);
1372 return result;
1373 }
1374 if (((!(SR_HDL(a) & SR_INT))&&(a->s<2))
1375 || ((!(SR_HDL(b) & SR_INT))&&(b->s<2))) return INT_TO_SR(1);
1376 if (SR_HDL(a) & SR_INT)
1377 {
1378 LONG aa=ABS(SR_TO_INT(a));
1379 unsigned long t=mpz_gcd_ui(NULL,b->z,(long)aa);
1380 if (t==POW_2_28)
1382 else
1383 result=INT_TO_SR(t);
1384 }
1385 else
1386 if (SR_HDL(b) & SR_INT)
1387 {
1388 LONG bb=ABS(SR_TO_INT(b));
1389 unsigned long t=mpz_gcd_ui(NULL,a->z,(long)bb);
1390 if (t==POW_2_28)
1392 else
1393 result=INT_TO_SR(t);
1394 }
1395 else
1396 {
1398 result->s = 3;
1399 #ifdef LDEBUG
1400 result->debug=123456;
1401 #endif
1402 mpz_init(result->z);
1403 mpz_gcd(result->z,a->z,b->z);
1405 }
1406 nlTest(result, r);
1407 return result;
1408}
1409
1410static int int_extgcd(int a, int b, int * u, int* x, int * v, int* y)
1411{
1412 int q, r;
1413 if (a==0)
1414 {
1415 *u = 0;
1416 *v = 1;
1417 *x = -1;
1418 *y = 0;
1419 return b;
1420 }
1421 if (b==0)
1422 {
1423 *u = 1;
1424 *v = 0;
1425 *x = 0;
1426 *y = 1;
1427 return a;
1428 }
1429 *u=1;
1430 *v=0;
1431 *x=0;
1432 *y=1;
1433 do
1434 {
1435 q = a/b;
1436 r = a%b;
1437 assume (q*b+r == a);
1438 a = b;
1439 b = r;
1440
1441 r = -(*v)*q+(*u);
1442 (*u) =(*v);
1443 (*v) = r;
1444
1445 r = -(*y)*q+(*x);
1446 (*x) = (*y);
1447 (*y) = r;
1448 } while (b);
1449
1450 return a;
1451}
1452
1453//number nlGcd_dummy(number a, number b, const coeffs r)
1454//{
1455// extern char my_yylinebuf[80];
1456// Print("nlGcd in >>%s<<\n",my_yylinebuf);
1457// return nlGcd(a,b,r);;
1458//}
1459
1460number nlShort1(number x) // assume x->s==0/1
1461{
1462 assume(x->s<2);
1463 if (mpz_sgn1(x->z)==0)
1464 {
1466 return INT_TO_SR(0);
1467 }
1468 if (x->s<2)
1469 {
1470 if (mpz_cmp(x->z,x->n)==0)
1471 {
1473 return INT_TO_SR(1);
1474 }
1475 }
1476 return x;
1477}
1478/*2
1479* simplify x
1480*/
1481void nlNormalize (number &x, const coeffs r)
1482{
1483 if ((SR_HDL(x) & SR_INT) ||(x==NULL))
1484 return;
1485 if (x->s==3)
1486 {
1488 nlTest(x,r);
1489 return;
1490 }
1491 else if (x->s==0)
1492 {
1493 if (mpz_cmp_si(x->n,1L)==0)
1494 {
1495 mpz_clear(x->n);
1496 x->s=3;
1497 x=nlShort3(x);
1498 }
1499 else
1500 {
1501 mpz_t gcd;
1502 mpz_init(gcd);
1503 mpz_gcd(gcd,x->z,x->n);
1504 x->s=1;
1505 if (mpz_cmp_si(gcd,1L)!=0)
1506 {
1507 mpz_divexact(x->z,x->z,gcd);
1508 mpz_divexact(x->n,x->n,gcd);
1509 if (mpz_cmp_si(x->n,1L)==0)
1510 {
1511 mpz_clear(x->n);
1512 x->s=3;
1514 }
1515 }
1516 mpz_clear(gcd);
1517 }
1518 }
1519 nlTest(x, r);
1520}
1521
1522/*2
1523* returns in result->z the lcm(a->z,b->n)
1524*/
1526{
1527 number result;
1528 nlTest(a, r);
1529 nlTest(b, r);
1530 if ((SR_HDL(b) & SR_INT)
1531 || (b->s==3))
1532 {
1533 // b is 1/(b->n) => b->n is 1 => result is a
1534 return nlCopy(a,r);
1535 }
1537#if defined(LDEBUG)
1538 result->debug=123456;
1539#endif
1540 result->s=3;
1541 mpz_t gcd;
1542 mpz_init(gcd);
1543 mpz_init(result->z);
1544 if (SR_HDL(a) & SR_INT)
1545 mpz_gcd_ui(gcd,b->n,ABS(SR_TO_INT(a)));
1546 else
1547 mpz_gcd(gcd,a->z,b->n);
1548 if (mpz_cmp_si(gcd,1L)!=0)
1549 {
1550 mpz_t bt;
1551 mpz_init(bt);
1552 mpz_divexact(bt,b->n,gcd);
1553 if (SR_HDL(a) & SR_INT)
1555 else
1556 mpz_mul(result->z,bt,a->z);
1557 mpz_clear(bt);
1558 }
1559 else
1560 if (SR_HDL(a) & SR_INT)
1561 mpz_mul_si(result->z,b->n,SR_TO_INT(a));
1562 else
1563 mpz_mul(result->z,b->n,a->z);
1564 mpz_clear(gcd);
1566 nlTest(result, r);
1567 return result;
1568}
1569
1570// Map q \in QQ or ZZ \to Zp or an extension of it
1571// src = Q or Z, dst = Zp (or an extension of Zp)
1572number nlModP(number q, const coeffs /*Q*/, const coeffs Zp)
1573{
1574 const int p = n_GetChar(Zp);
1575 assume( p > 0 );
1576
1577 const long P = p;
1578 assume( P > 0 );
1579
1580 // embedded long within q => only long numerator has to be converted
1581 // to int (modulo char.)
1582 if (SR_HDL(q) & SR_INT)
1583 {
1584 long i = SR_TO_INT(q);
1585 return n_Init( i, Zp );
1586 }
1587
1588 const unsigned long PP = p;
1589
1590 // numerator modulo char. should fit into int
1591 number z = n_Init( static_cast<long>(mpz_fdiv_ui(q->z, PP)), Zp );
1592
1593 // denominator != 1?
1594 if (q->s!=3)
1595 {
1596 // denominator modulo char. should fit into int
1597 number n = n_Init( static_cast<long>(mpz_fdiv_ui(q->n, PP)), Zp );
1598
1599 number res = n_Div( z, n, Zp );
1600
1601 n_Delete(&z, Zp);
1602 n_Delete(&n, Zp);
1603
1604 return res;
1605 }
1606
1607 return z;
1608}
1609
1610/*2
1611* convert number i (from Q) to GMP and warn if denom != 1
1612*/
1613void nlGMP(number &i, mpz_t n, const coeffs r)
1614{
1615 // Hier brauche ich einfach die GMP Zahl
1616 nlTest(i, r);
1617 nlNormalize(i, r);
1618 if (SR_HDL(i) & SR_INT)
1619 {
1620 mpz_set_si(n, SR_TO_INT(i));
1621 return;
1622 }
1623 if (i->s!=3)
1624 {
1625 WarnS("Omitted denominator during coefficient mapping !");
1626 }
1627 mpz_set(n, i->z);
1628}
1629
1630/*2
1631* access to denominator, other 1 for integers
1632*/
1634{
1635 if (!(SR_HDL(n) & SR_INT))
1636 {
1637 if (n->s==0)
1638 {
1639 nlNormalize(n,r);
1640 }
1641 if (!(SR_HDL(n) & SR_INT))
1642 {
1643 if (n->s!=3)
1644 {
1646 u->s=3;
1647#if defined(LDEBUG)
1648 u->debug=123456;
1649#endif
1650 mpz_init_set(u->z,n->n);
1651 u=nlShort3_noinline(u);
1652 return u;
1653 }
1654 }
1655 }
1656 return INT_TO_SR(1);
1657}
1658
1659/*2
1660* access to Nominator, nlCopy(n) for integers
1661*/
1663{
1664 if (!(SR_HDL(n) & SR_INT))
1665 {
1666 if (n->s==0)
1667 {
1668 nlNormalize(n,r);
1669 }
1670 if (!(SR_HDL(n) & SR_INT))
1671 {
1673#if defined(LDEBUG)
1674 u->debug=123456;
1675#endif
1676 u->s=3;
1677 mpz_init_set(u->z,n->z);
1678 if (n->s!=3)
1679 {
1680 u=nlShort3_noinline(u);
1681 }
1682 return u;
1683 }
1684 }
1685 return n; // imm. int
1686}
1687
1688/***************************************************************
1689 *
1690 * routines which are needed by Inline(d) routines
1691 *
1692 *******************************************************************/
1694{
1695 assume(! (SR_HDL(a) & SR_HDL(b) & SR_INT));
1696// long - short
1697 BOOLEAN bo;
1698 if (SR_HDL(b) & SR_INT)
1699 {
1700 if (a->s!=0) return FALSE;
1701 number n=b; b=a; a=n;
1702 }
1703// short - long
1704 if (SR_HDL(a) & SR_INT)
1705 {
1706 if (b->s!=0)
1707 return FALSE;
1708 if ((((long)a) > 0L) && (mpz_isNeg(b->z)))
1709 return FALSE;
1710 if ((((long)a) < 0L) && (!mpz_isNeg(b->z)))
1711 return FALSE;
1712 mpz_t bb;
1713 mpz_init(bb);
1714 mpz_mul_si(bb,b->n,(long)SR_TO_INT(a));
1715 bo=(mpz_cmp(bb,b->z)==0);
1716 mpz_clear(bb);
1717 return bo;
1718 }
1719// long - long
1720 if (((a->s==1) && (b->s==3))
1721 || ((b->s==1) && (a->s==3)))
1722 return FALSE;
1723 if (mpz_isNeg(a->z)&&(!mpz_isNeg(b->z)))
1724 return FALSE;
1725 if (mpz_isNeg(b->z)&&(!mpz_isNeg(a->z)))
1726 return FALSE;
1727 mpz_t aa;
1728 mpz_t bb;
1729 mpz_init_set(aa,a->z);
1730 mpz_init_set(bb,b->z);
1731 if (a->s<2) mpz_mul(bb,bb,a->n);
1732 if (b->s<2) mpz_mul(aa,aa,b->n);
1733 bo=(mpz_cmp(aa,bb)==0);
1734 mpz_clear(aa);
1735 mpz_clear(bb);
1736 return bo;
1737}
1738
1739// copy not immediate number a
1741{
1742 assume(!(SR_HDL(a) & SR_INT));
1743 //nlTest(a, r);
1745#if defined(LDEBUG)
1746 b->debug=123456;
1747#endif
1748 switch (a->s)
1749 {
1750 case 0:
1751 case 1:
1752 mpz_init_set(b->n,a->n); /*no break*/
1753 case 3:
1754 mpz_init_set(b->z,a->z);
1755 break;
1756 }
1757 b->s = a->s;
1758 return b;
1759}
1760
1762{
1763 {
1764 switch ((*a)->s)
1765 {
1766 case 0:
1767 case 1:
1768 mpz_clear((*a)->n); /*no break*/
1769 case 3:
1770 mpz_clear((*a)->z);
1771 }
1772 #ifdef LDEBUG
1773 memset(*a,0,sizeof(**a));
1774 #endif
1775 FREE_RNUMBER(*a); // omFreeBin((void *) *a, rnumber_bin);
1776 }
1777}
1778
1780{
1781 mpz_neg(a->z,a->z);
1782 if (a->s==3)
1783 {
1784 a=nlShort3(a);
1785 }
1786 return a;
1787}
1788
1789// conditio to use nlNormalize_Gcd in intermediate computations:
1790#define GCD_NORM_COND(OLD,NEW) (mpz_size1(NEW->z)>mpz_size1(OLD->z))
1791
1793{
1794 mpz_t gcd;
1795 mpz_init(gcd);
1796 mpz_gcd(gcd,x->z,x->n);
1797 x->s=1;
1798 if (mpz_cmp_si(gcd,1L)!=0)
1799 {
1800 mpz_divexact(x->z,x->z,gcd);
1801 mpz_divexact(x->n,x->n,gcd);
1802 if (mpz_cmp_si(x->n,1L)==0)
1803 {
1804 mpz_clear(x->n);
1805 x->s=3;
1807 }
1808 }
1809 mpz_clear(gcd);
1810}
1811
1813{
1815#if defined(LDEBUG)
1816 u->debug=123456;
1817#endif
1818 mpz_init(u->z);
1819 if (SR_HDL(b) & SR_INT)
1820 {
1821 number x=a;
1822 a=b;
1823 b=x;
1824 }
1825 if (SR_HDL(a) & SR_INT)
1826 {
1827 switch (b->s)
1828 {
1829 case 0:
1830 case 1:/* a:short, b:1 */
1831 {
1832 mpz_t x;
1833 mpz_init(x);
1834 mpz_mul_si(x,b->n,SR_TO_INT(a));
1835 mpz_add(u->z,b->z,x);
1836 mpz_clear(x);
1837 if (mpz_sgn1(u->z)==0)
1838 {
1839 mpz_clear(u->z);
1840 FREE_RNUMBER(u);
1841 return INT_TO_SR(0);
1842 }
1843 if (mpz_cmp(u->z,b->n)==0)
1844 {
1845 mpz_clear(u->z);
1846 FREE_RNUMBER(u);
1847 return INT_TO_SR(1);
1848 }
1849 mpz_init_set(u->n,b->n);
1850 u->s = 0;
1851 if (GCD_NORM_COND(b,u)) { nlNormalize_Gcd(u); }
1852 break;
1853 }
1854 case 3:
1855 {
1856 if (((long)a)>0L)
1857 mpz_add_ui(u->z,b->z,SR_TO_INT(a));
1858 else
1859 mpz_sub_ui(u->z,b->z,-SR_TO_INT(a));
1860 u->s = 3;
1861 u=nlShort3(u);
1862 break;
1863 }
1864 }
1865 }
1866 else
1867 {
1868 switch (a->s)
1869 {
1870 case 0:
1871 case 1:
1872 {
1873 switch(b->s)
1874 {
1875 case 0:
1876 case 1:
1877 {
1878 mpz_t x;
1879 mpz_init(x);
1880
1881 mpz_mul(x,b->z,a->n);
1882 mpz_mul(u->z,a->z,b->n);
1883 mpz_add(u->z,u->z,x);
1884 mpz_clear(x);
1885
1886 if (mpz_sgn1(u->z)==0)
1887 {
1888 mpz_clear(u->z);
1889 FREE_RNUMBER(u);
1890 return INT_TO_SR(0);
1891 }
1892 mpz_init(u->n);
1893 mpz_mul(u->n,a->n,b->n);
1894 if (mpz_cmp(u->z,u->n)==0)
1895 {
1896 mpz_clear(u->z);
1897 mpz_clear(u->n);
1898 FREE_RNUMBER(u);
1899 return INT_TO_SR(1);
1900 }
1901 u->s = 0;
1902 if (GCD_NORM_COND(b,u)) { nlNormalize_Gcd(u); }
1903 break;
1904 }
1905 case 3: /* a:1 b:3 */
1906 {
1907 mpz_mul(u->z,b->z,a->n);
1908 mpz_add(u->z,u->z,a->z);
1909 if (mpz_sgn1(u->z)==0)
1910 {
1911 mpz_clear(u->z);
1912 FREE_RNUMBER(u);
1913 return INT_TO_SR(0);
1914 }
1915 if (mpz_cmp(u->z,a->n)==0)
1916 {
1917 mpz_clear(u->z);
1918 FREE_RNUMBER(u);
1919 return INT_TO_SR(1);
1920 }
1921 mpz_init_set(u->n,a->n);
1922 u->s = 0;
1923 if (GCD_NORM_COND(a,u)) { nlNormalize_Gcd(u); }
1924 break;
1925 }
1926 } /*switch (b->s) */
1927 break;
1928 }
1929 case 3:
1930 {
1931 switch(b->s)
1932 {
1933 case 0:
1934 case 1:/* a:3, b:1 */
1935 {
1936 mpz_mul(u->z,a->z,b->n);
1937 mpz_add(u->z,u->z,b->z);
1938 if (mpz_sgn1(u->z)==0)
1939 {
1940 mpz_clear(u->z);
1941 FREE_RNUMBER(u);
1942 return INT_TO_SR(0);
1943 }
1944 if (mpz_cmp(u->z,b->n)==0)
1945 {
1946 mpz_clear(u->z);
1947 FREE_RNUMBER(u);
1948 return INT_TO_SR(1);
1949 }
1950 mpz_init_set(u->n,b->n);
1951 u->s = 0;
1952 if (GCD_NORM_COND(b,u)) { nlNormalize_Gcd(u); }
1953 break;
1954 }
1955 case 3:
1956 {
1957 mpz_add(u->z,a->z,b->z);
1958 u->s = 3;
1959 u=nlShort3(u);
1960 break;
1961 }
1962 }
1963 break;
1964 }
1965 }
1966 }
1967 return u;
1968}
1969
1971{
1972 if (SR_HDL(b) & SR_INT)
1973 {
1974 switch (a->s)
1975 {
1976 case 0:
1977 case 1:/* b:short, a:1 */
1978 {
1979 mpz_t x;
1980 mpz_init(x);
1981 mpz_mul_si(x,a->n,SR_TO_INT(b));
1982 mpz_add(a->z,a->z,x);
1983 mpz_clear(x);
1984 nlNormalize_Gcd(a);
1985 break;
1986 }
1987 case 3:
1988 {
1989 if (((long)b)>0L)
1990 mpz_add_ui(a->z,a->z,SR_TO_INT(b));
1991 else
1992 mpz_sub_ui(a->z,a->z,-SR_TO_INT(b));
1993 a->s = 3;
1994 a=nlShort3_noinline(a);
1995 break;
1996 }
1997 }
1998 return;
1999 }
2000 else if (SR_HDL(a) & SR_INT)
2001 {
2003 #if defined(LDEBUG)
2004 u->debug=123456;
2005 #endif
2006 mpz_init(u->z);
2007 switch (b->s)
2008 {
2009 case 0:
2010 case 1:/* a:short, b:1 */
2011 {
2012 mpz_t x;
2013 mpz_init(x);
2014
2015 mpz_mul_si(x,b->n,SR_TO_INT(a));
2016 mpz_add(u->z,b->z,x);
2017 mpz_clear(x);
2018 // result cannot be 0, if coeffs are normalized
2019 mpz_init_set(u->n,b->n);
2020 u->s=0;
2021 if (GCD_NORM_COND(b,u)) { nlNormalize_Gcd(u); }
2022 else { u=nlShort1(u); }
2023 break;
2024 }
2025 case 3:
2026 {
2027 if (((long)a)>0L)
2028 mpz_add_ui(u->z,b->z,SR_TO_INT(a));
2029 else
2030 mpz_sub_ui(u->z,b->z,-SR_TO_INT(a));
2031 // result cannot be 0, if coeffs are normalized
2032 u->s = 3;
2033 u=nlShort3_noinline(u);
2034 break;
2035 }
2036 }
2037 a=u;
2038 }
2039 else
2040 {
2041 switch (a->s)
2042 {
2043 case 0:
2044 case 1:
2045 {
2046 switch(b->s)
2047 {
2048 case 0:
2049 case 1: /* a:1 b:1 */
2050 {
2051 mpz_t x;
2052 mpz_t y;
2053 mpz_init(x);
2054 mpz_init(y);
2055 mpz_mul(x,b->z,a->n);
2056 mpz_mul(y,a->z,b->n);
2057 mpz_add(a->z,x,y);
2058 mpz_clear(x);
2059 mpz_clear(y);
2060 mpz_mul(a->n,a->n,b->n);
2061 a->s=0;
2062 if (GCD_NORM_COND(b,a)) { nlNormalize_Gcd(a); }
2063 else { a=nlShort1(a);}
2064 break;
2065 }
2066 case 3: /* a:1 b:3 */
2067 {
2068 mpz_t x;
2069 mpz_init(x);
2070 mpz_mul(x,b->z,a->n);
2071 mpz_add(a->z,a->z,x);
2072 mpz_clear(x);
2073 a->s=0;
2074 if (GCD_NORM_COND(b,a)) { nlNormalize_Gcd(a); }
2075 else { a=nlShort1(a);}
2076 break;
2077 }
2078 } /*switch (b->s) */
2079 break;
2080 }
2081 case 3:
2082 {
2083 switch(b->s)
2084 {
2085 case 0:
2086 case 1:/* a:3, b:1 */
2087 {
2088 mpz_t x;
2089 mpz_init(x);
2090 mpz_mul(x,a->z,b->n);
2091 mpz_add(a->z,b->z,x);
2092 mpz_clear(x);
2093 mpz_init_set(a->n,b->n);
2094 a->s=0;
2095 if (GCD_NORM_COND(b,a)) { nlNormalize_Gcd(a); }
2096 else { a=nlShort1(a);}
2097 break;
2098 }
2099 case 3:
2100 {
2101 mpz_add(a->z,a->z,b->z);
2102 a->s = 3;
2103 a=nlShort3_noinline(a);
2104 break;
2105 }
2106 }
2107 break;
2108 }
2109 }
2110 }
2111}
2112
2114{
2116#if defined(LDEBUG)
2117 u->debug=123456;
2118#endif
2119 mpz_init(u->z);
2120 if (SR_HDL(a) & SR_INT)
2121 {
2122 switch (b->s)
2123 {
2124 case 0:
2125 case 1:/* a:short, b:1 */
2126 {
2127 mpz_t x;
2128 mpz_init(x);
2129 mpz_mul_si(x,b->n,SR_TO_INT(a));
2130 mpz_sub(u->z,x,b->z);
2131 mpz_clear(x);
2132 if (mpz_sgn1(u->z)==0)
2133 {
2134 mpz_clear(u->z);
2135 FREE_RNUMBER(u);
2136 return INT_TO_SR(0);
2137 }
2138 if (mpz_cmp(u->z,b->n)==0)
2139 {
2140 mpz_clear(u->z);
2141 FREE_RNUMBER(u);
2142 return INT_TO_SR(1);
2143 }
2144 mpz_init_set(u->n,b->n);
2145 u->s=0;
2146 if (GCD_NORM_COND(b,u)) { nlNormalize_Gcd(u); }
2147 break;
2148 }
2149 case 3:
2150 {
2151 if (((long)a)>0L)
2152 {
2153 mpz_sub_ui(u->z,b->z,SR_TO_INT(a));
2154 mpz_neg(u->z,u->z);
2155 }
2156 else
2157 {
2158 mpz_add_ui(u->z,b->z,-SR_TO_INT(a));
2159 mpz_neg(u->z,u->z);
2160 }
2161 u->s = 3;
2162 u=nlShort3(u);
2163 break;
2164 }
2165 }
2166 }
2167 else if (SR_HDL(b) & SR_INT)
2168 {
2169 switch (a->s)
2170 {
2171 case 0:
2172 case 1:/* b:short, a:1 */
2173 {
2174 mpz_t x;
2175 mpz_init(x);
2176 mpz_mul_si(x,a->n,SR_TO_INT(b));
2177 mpz_sub(u->z,a->z,x);
2178 mpz_clear(x);
2179 if (mpz_sgn1(u->z)==0)
2180 {
2181 mpz_clear(u->z);
2182 FREE_RNUMBER(u);
2183 return INT_TO_SR(0);
2184 }
2185 if (mpz_cmp(u->z,a->n)==0)
2186 {
2187 mpz_clear(u->z);
2188 FREE_RNUMBER(u);
2189 return INT_TO_SR(1);
2190 }
2191 mpz_init_set(u->n,a->n);
2192 u->s=0;
2193 if (GCD_NORM_COND(a,u)) { nlNormalize_Gcd(u); }
2194 break;
2195 }
2196 case 3:
2197 {
2198 if (((long)b)>0L)
2199 {
2200 mpz_sub_ui(u->z,a->z,SR_TO_INT(b));
2201 }
2202 else
2203 {
2204 mpz_add_ui(u->z,a->z,-SR_TO_INT(b));
2205 }
2206 u->s = 3;
2207 u=nlShort3(u);
2208 break;
2209 }
2210 }
2211 }
2212 else
2213 {
2214 switch (a->s)
2215 {
2216 case 0:
2217 case 1:
2218 {
2219 switch(b->s)
2220 {
2221 case 0:
2222 case 1:
2223 {
2224 mpz_t x;
2225 mpz_t y;
2226 mpz_init(x);
2227 mpz_init(y);
2228 mpz_mul(x,b->z,a->n);
2229 mpz_mul(y,a->z,b->n);
2230 mpz_sub(u->z,y,x);
2231 mpz_clear(x);
2232 mpz_clear(y);
2233 if (mpz_sgn1(u->z)==0)
2234 {
2235 mpz_clear(u->z);
2236 FREE_RNUMBER(u);
2237 return INT_TO_SR(0);
2238 }
2239 mpz_init(u->n);
2240 mpz_mul(u->n,a->n,b->n);
2241 if (mpz_cmp(u->z,u->n)==0)
2242 {
2243 mpz_clear(u->z);
2244 mpz_clear(u->n);
2245 FREE_RNUMBER(u);
2246 return INT_TO_SR(1);
2247 }
2248 u->s=0;
2249 if (GCD_NORM_COND(a,u)) { nlNormalize_Gcd(u); }
2250 break;
2251 }
2252 case 3: /* a:1, b:3 */
2253 {
2254 mpz_t x;
2255 mpz_init(x);
2256 mpz_mul(x,b->z,a->n);
2257 mpz_sub(u->z,a->z,x);
2258 mpz_clear(x);
2259 if (mpz_sgn1(u->z)==0)
2260 {
2261 mpz_clear(u->z);
2262 FREE_RNUMBER(u);
2263 return INT_TO_SR(0);
2264 }
2265 if (mpz_cmp(u->z,a->n)==0)
2266 {
2267 mpz_clear(u->z);
2268 FREE_RNUMBER(u);
2269 return INT_TO_SR(1);
2270 }
2271 mpz_init_set(u->n,a->n);
2272 u->s=0;
2273 if (GCD_NORM_COND(a,u)) { nlNormalize_Gcd(u); }
2274 break;
2275 }
2276 }
2277 break;
2278 }
2279 case 3:
2280 {
2281 switch(b->s)
2282 {
2283 case 0:
2284 case 1: /* a:3, b:1 */
2285 {
2286 mpz_t x;
2287 mpz_init(x);
2288 mpz_mul(x,a->z,b->n);
2289 mpz_sub(u->z,x,b->z);
2290 mpz_clear(x);
2291 if (mpz_sgn1(u->z)==0)
2292 {
2293 mpz_clear(u->z);
2294 FREE_RNUMBER(u);
2295 return INT_TO_SR(0);
2296 }
2297 if (mpz_cmp(u->z,b->n)==0)
2298 {
2299 mpz_clear(u->z);
2300 FREE_RNUMBER(u);
2301 return INT_TO_SR(1);
2302 }
2303 mpz_init_set(u->n,b->n);
2304 u->s=0;
2305 if (GCD_NORM_COND(b,u)) { nlNormalize_Gcd(u); }
2306 break;
2307 }
2308 case 3: /* a:3 , b:3 */
2309 {
2310 mpz_sub(u->z,a->z,b->z);
2311 u->s = 3;
2312 u=nlShort3(u);
2313 break;
2314 }
2315 }
2316 break;
2317 }
2318 }
2319 }
2320 return u;
2321}
2322
2323// a and b are intermediate, but a*b not
2325{
2327#if defined(LDEBUG)
2328 u->debug=123456;
2329#endif
2330 u->s=3;
2331 mpz_init_set_si(u->z,SR_TO_INT(a));
2332 mpz_mul_si(u->z,u->z,SR_TO_INT(b));
2333 return u;
2334}
2335
2336// a or b are not immediate
2338{
2339 assume(! (SR_HDL(a) & SR_HDL(b) & SR_INT));
2341#if defined(LDEBUG)
2342 u->debug=123456;
2343#endif
2344 mpz_init(u->z);
2345 if (SR_HDL(b) & SR_INT)
2346 {
2347 number x=a;
2348 a=b;
2349 b=x;
2350 }
2351 if (SR_HDL(a) & SR_INT)
2352 {
2353 u->s=b->s;
2354 if (u->s==1) u->s=0;
2355 if (((long)a)>0L)
2356 {
2357 mpz_mul_ui(u->z,b->z,(unsigned long)SR_TO_INT(a));
2358 }
2359 else
2360 {
2361 if (a==INT_TO_SR(-1))
2362 {
2363 mpz_set(u->z,b->z);
2364 mpz_neg(u->z,u->z);
2365 u->s=b->s;
2366 }
2367 else
2368 {
2369 mpz_mul_ui(u->z,b->z,(unsigned long)-SR_TO_INT(a));
2370 mpz_neg(u->z,u->z);
2371 }
2372 }
2373 if (u->s<2)
2374 {
2375 if (mpz_cmp(u->z,b->n)==0)
2376 {
2377 mpz_clear(u->z);
2378 FREE_RNUMBER(u);
2379 return INT_TO_SR(1);
2380 }
2381 mpz_init_set(u->n,b->n);
2382 if (GCD_NORM_COND(b,u)) { nlNormalize_Gcd(u); }
2383 }
2384 else //u->s==3
2385 {
2386 u=nlShort3(u);
2387 }
2388 }
2389 else
2390 {
2391 mpz_mul(u->z,a->z,b->z);
2392 u->s = 0;
2393 if(a->s==3)
2394 {
2395 if(b->s==3)
2396 {
2397 u->s = 3;
2398 }
2399 else
2400 {
2401 if (mpz_cmp(u->z,b->n)==0)
2402 {
2403 mpz_clear(u->z);
2404 FREE_RNUMBER(u);
2405 return INT_TO_SR(1);
2406 }
2407 mpz_init_set(u->n,b->n);
2408 if (GCD_NORM_COND(b,u)) { nlNormalize_Gcd(u); }
2409 }
2410 }
2411 else
2412 {
2413 if(b->s==3)
2414 {
2415 if (mpz_cmp(u->z,a->n)==0)
2416 {
2417 mpz_clear(u->z);
2418 FREE_RNUMBER(u);
2419 return INT_TO_SR(1);
2420 }
2421 mpz_init_set(u->n,a->n);
2422 if (GCD_NORM_COND(a,u)) { nlNormalize_Gcd(u); }
2423 }
2424 else
2425 {
2426 mpz_init(u->n);
2427 mpz_mul(u->n,a->n,b->n);
2428 if (mpz_cmp(u->z,u->n)==0)
2429 {
2430 mpz_clear(u->z);
2431 mpz_clear(u->n);
2432 FREE_RNUMBER(u);
2433 return INT_TO_SR(1);
2434 }
2435 if (GCD_NORM_COND(a,u)) { nlNormalize_Gcd(u); }
2436 }
2437 }
2438 }
2439 return u;
2440}
2441
2442/*2
2443* copy a to b for mapping
2444*/
2445number nlCopyMap(number a, const coeffs /*src*/, const coeffs /*dst*/)
2446{
2447 if ((SR_HDL(a) & SR_INT)||(a==NULL))
2448 {
2449 return a;
2450 }
2451 return _nlCopy_NoImm(a);
2452}
2453
2455{
2456 if ((SR_HDL(a) & SR_INT)||(a==NULL))
2457 {
2458 return a;
2459 }
2460 if (a->s==3) return _nlCopy_NoImm(a);
2461 number a0=a;
2462 BOOLEAN a1=FALSE;
2463 if (a->s==0) { a0=_nlCopy_NoImm(a); a1=TRUE; }
2465 number b2=nlGetDenom(a0,src);
2467 nlDelete(&b1,src);
2468 nlDelete(&b2,src);
2469 if (a1) _nlDelete_NoImm(&a0);
2470 return b;
2471}
2472
2474{
2475 if (src->rep==n_rep_gap_rat) /*Q, coeffs_BIGINT */
2476 {
2477 if ((src->is_field==dst->is_field) /* Q->Q, Z->Z*/
2478 || (src->is_field==FALSE)) /* Z->Q */
2479 return nlCopyMap;
2480 return nlMapQtoZ; /* Q->Z */
2481 }
2482 if ((src->rep==n_rep_int) && nCoeff_is_Zp(src))
2483 {
2484 return nlMapP;
2485 }
2486 if ((src->rep==n_rep_float) && nCoeff_is_R(src))
2487 {
2488 if (dst->is_field) /* R -> Q */
2489 return nlMapR;
2490 else
2491 return nlMapR_BI; /* R -> bigint */
2492 }
2493 if ((src->rep==n_rep_gmp_float) && nCoeff_is_long_R(src))
2494 {
2495 if (dst->is_field)
2496 return nlMapLongR; /* long R -> Q */
2497 else
2498 return nlMapLongR_BI;
2499 }
2500 if (nCoeff_is_long_C(src))
2501 {
2502 return nlMapC; /* C -> Q */
2503 }
2504 if (src->rep==n_rep_gmp) // nCoeff_is_Z(src) || nCoeff_is_Ring_PtoM(src) || nCoeff_is_Zn(src))
2505 {
2506 return nlMapGMP;
2507 }
2508 if (src->rep==n_rep_gap_gmp)
2509 {
2510 return nlMapZ;
2511 }
2512 if ((src->rep==n_rep_int) && nCoeff_is_Ring_2toM(src))
2513 {
2514 return nlMapMachineInt;
2515 }
2516 return NULL;
2517}
2518/*2
2519* z := i
2520*/
2522{
2524#if defined(LDEBUG)
2525 z->debug=123456;
2526#endif
2527 mpz_init_set_si(z->z,i);
2528 z->s = 3;
2529 return z;
2530}
2531
2532/*2
2533* z := i/j
2534*/
2535number nlInit2 (int i, int j, const coeffs r)
2536{
2538#if defined(LDEBUG)
2539 z->debug=123456;
2540#endif
2541 mpz_init_set_si(z->z,(long)i);
2542 mpz_init_set_si(z->n,(long)j);
2543 z->s = 0;
2544 nlNormalize(z,r);
2545 return z;
2546}
2547
2549{
2551#if defined(LDEBUG)
2552 z->debug=123456;
2553#endif
2554 mpz_init_set(z->z,i);
2555 mpz_init_set(z->n,j);
2556 z->s = 0;
2557 nlNormalize(z,r);
2558 return z;
2559}
2560
2561#else // DO_LINLINE
2562
2563// declare immediate routines
2564number nlRInit (long i);
2573
2574#endif
2575
2576/***************************************************************
2577 *
2578 * Routines which might be inlined by p_Numbers.h
2579 *
2580 *******************************************************************/
2581#if defined(DO_LINLINE) || !defined(P_NUMBERS_H)
2582
2583// routines which are always inlined/static
2584
2585/*2
2586* a = b ?
2587*/
2589{
2590 nlTest(a, r);
2591 nlTest(b, r);
2592// short - short
2593 if (SR_HDL(a) & SR_HDL(b) & SR_INT) return a==b;
2594 return _nlEqual_aNoImm_OR_bNoImm(a, b);
2595}
2596
2598{
2599 number n;
2600 #if MAX_NUM_SIZE == 60
2601 if (((i << 3) >> 3) == i) n=INT_TO_SR(i);
2602 else n=nlRInit(i);
2603 #else
2604 LONG ii=(LONG)i;
2605 if ( ((((long)ii)==i) && ((ii << 3) >> 3) == ii )) n=INT_TO_SR(ii);
2606 else n=nlRInit(i);
2607 #endif
2608 nlTest(n, r);
2609 return n;
2610}
2611
2612/*2
2613* a == 1 ?
2614*/
2616{
2617#ifdef LDEBUG
2618 if (a==NULL) return FALSE;
2619 nlTest(a, r);
2620#endif
2621 return (a==INT_TO_SR(1));
2622}
2623
2625{
2626 #if 0
2627 if (a==INT_TO_SR(0)) return TRUE;
2628 if ((SR_HDL(a) & SR_INT)||(a==NULL)) return FALSE;
2629 if (mpz_cmp_si(a->z,0L)==0)
2630 {
2631 printf("gmp-0 in nlIsZero\n");
2632 dErrorBreak();
2633 return TRUE;
2634 }
2635 return FALSE;
2636 #else
2637 return (a==NULL)|| (a==INT_TO_SR(0));
2638 #endif
2639}
2640
2641/*2
2642* copy a to b
2643*/
2645{
2646 if (SR_HDL(a) & SR_INT)
2647 {
2648 return a;
2649 }
2650 return _nlCopy_NoImm(a);
2651}
2652
2653
2654/*2
2655* delete a
2656*/
2657LINLINE void nlDelete (number * a, const coeffs r)
2658{
2659 if (*a!=NULL)
2660 {
2661 nlTest(*a, r);
2662 if ((SR_HDL(*a) & SR_INT)==0)
2663 {
2664 _nlDelete_NoImm(a);
2665 }
2666 *a=NULL;
2667 }
2668}
2669
2670/*2
2671* za:= - za
2672*/
2674{
2675 nlTest(a, R);
2676 if(SR_HDL(a) &SR_INT)
2677 {
2678 LONG r=SR_TO_INT(a);
2679 if (r==(-(POW_2_28))) a=nlRInit(POW_2_28);
2680 else a=INT_TO_SR(-r);
2681 return a;
2682 }
2683 a = _nlNeg_NoImm(a);
2684 nlTest(a, R);
2685 return a;
2686
2687}
2688
2689/*2
2690* u:= a + b
2691*/
2693{
2694 if (SR_HDL(a) & SR_HDL(b) & SR_INT)
2695 {
2696 LONG r=SR_HDL(a)+SR_HDL(b)-1L;
2697 if ( ((r << 1) >> 1) == r )
2698 return (number)(long)r;
2699 else
2700 return nlRInit(SR_TO_INT(r));
2701 }
2703 nlTest(u, R);
2704 return u;
2705}
2706
2709
2711{
2712 // a=a+b
2713 if (SR_HDL(a) & SR_HDL(b) & SR_INT)
2714 {
2715 LONG r=SR_HDL(a)+SR_HDL(b)-1L;
2716 if ( ((r << 1) >> 1) == r )
2717 a=(number)(long)r;
2718 else
2719 a=nlRInit(SR_TO_INT(r));
2720 }
2721 else
2722 {
2724 nlTest(a,r);
2725 }
2726}
2727
2729{
2730 nlTest(a, R);
2731 nlTest(b, R);
2732 if (a==INT_TO_SR(0)) return INT_TO_SR(0);
2733 if (b==INT_TO_SR(0)) return INT_TO_SR(0);
2734 if (SR_HDL(a) & SR_HDL(b) & SR_INT)
2735 {
2736 LONG r=(LONG)((unsigned LONG)(SR_HDL(a)-1L))*((unsigned LONG)(SR_HDL(b)>>1));
2737 if ((r/(SR_HDL(b)>>1))==(SR_HDL(a)-1L))
2738 {
2739 number u=((number) ((r>>1)+SR_INT));
2740 if (((((LONG)SR_HDL(u))<<1)>>1)==SR_HDL(u)) return (u);
2741 return nlRInit(SR_HDL(u)>>2);
2742 }
2744 nlTest(u, R);
2745 return u;
2746
2747 }
2749 nlTest(u, R);
2750 return u;
2751
2752}
2753
2754
2755/*2
2756* u:= a - b
2757*/
2759{
2760 if (SR_HDL(a) & SR_HDL(b) & SR_INT)
2761 {
2762 LONG r=SR_HDL(a)-SR_HDL(b)+1;
2763 if ( ((r << 1) >> 1) == r )
2764 {
2765 return (number)(long)r;
2766 }
2767 else
2768 return nlRInit(SR_TO_INT(r));
2769 }
2771 nlTest(u, r);
2772 return u;
2773
2774}
2775
2777{
2778 number aa=a;
2779 if (((SR_HDL(b)|SR_HDL(aa))&SR_INT))
2780 {
2781 number n=nlMult(aa,b,r);
2782 nlDelete(&a,r);
2783 a=n;
2784 }
2785 else
2786 {
2787 mpz_mul(aa->z,a->z,b->z);
2788 if (aa->s==3)
2789 {
2790 if(b->s!=3)
2791 {
2792 mpz_init_set(a->n,b->n);
2793 a->s=0;
2794 }
2795 }
2796 else
2797 {
2798 if(b->s!=3)
2799 {
2800 mpz_mul(a->n,a->n,b->n);
2801 }
2802 a->s=0;
2803 }
2804 }
2805}
2806#endif // DO_LINLINE
2807
2808#ifndef P_NUMBERS_H
2809
2810void nlMPZ(mpz_t m, number &n, const coeffs r)
2811{
2812 nlTest(n, r);
2813 nlNormalize(n, r);
2814 if (SR_HDL(n) & SR_INT) mpz_init_set_si(m, SR_TO_INT(n)); /* n fits in an int */
2815 else mpz_init_set(m, (mpz_ptr)n->z);
2816}
2817
2818
2820{
2821 if (SR_HDL(a) & SR_HDL(b) & SR_INT)
2822 {
2823 int uu, vv, x, y;
2824 int g = int_extgcd(SR_TO_INT(a), SR_TO_INT(b), &uu, &vv, &x, &y);
2825 *s = INT_TO_SR(uu);
2826 *t = INT_TO_SR(vv);
2827 *u = INT_TO_SR(x);
2828 *v = INT_TO_SR(y);
2829 return INT_TO_SR(g);
2830 }
2831 else
2832 {
2833 mpz_t aa, bb;
2834 if (SR_HDL(a) & SR_INT)
2835 {
2837 }
2838 else
2839 {
2840 mpz_init_set(aa, a->z);
2841 }
2842 if (SR_HDL(b) & SR_INT)
2843 {
2845 }
2846 else
2847 {
2848 mpz_init_set(bb, b->z);
2849 }
2851 mpz_init(erg);
2852 mpz_init(bs);
2853 mpz_init(bt);
2854
2855 mpz_gcdext(erg, bs, bt, aa, bb);
2856
2857 mpz_div(aa, aa, erg);
2858 *u=nlInitMPZ(bb,r);
2859 *u=nlNeg(*u,r);
2860 *v=nlInitMPZ(aa,r);
2861
2862 mpz_clear(aa);
2863 mpz_clear(bb);
2864
2865 *s = nlInitMPZ(bs,r);
2866 *t = nlInitMPZ(bt,r);
2867 return nlInitMPZ(erg,r);
2868 }
2869}
2870
2872{
2873 assume(SR_TO_INT(b)!=0);
2874 if (SR_HDL(a) & SR_HDL(b) & SR_INT)
2875 {
2876 if (r!=NULL)
2877 *r = INT_TO_SR(SR_TO_INT(a) % SR_TO_INT(b));
2878 return INT_TO_SR(SR_TO_INT(a)/SR_TO_INT(b));
2879 }
2880 else if (SR_HDL(a) & SR_INT)
2881 {
2882 // -2^xx / 2^xx
2883 if ((a==INT_TO_SR(-(POW_2_28)))&&(b==INT_TO_SR(-1L)))
2884 {
2885 if (r!=NULL) *r=INT_TO_SR(0);
2886 return nlRInit(POW_2_28);
2887 }
2888 //a is small, b is not, so q=0, r=a
2889 if (r!=NULL)
2890 *r = a;
2891 return INT_TO_SR(0);
2892 }
2893 else if (SR_HDL(b) & SR_INT)
2894 {
2895 unsigned long rr;
2896 mpz_t qq;
2897 mpz_init(qq);
2898 mpz_t rrr;
2899 mpz_init(rrr);
2900 rr = mpz_divmod_ui(qq, rrr, a->z, (unsigned long)ABS(SR_TO_INT(b)));
2901 mpz_clear(rrr);
2902
2903 if (r!=NULL)
2904 *r = INT_TO_SR(rr);
2905 if (SR_TO_INT(b)<0)
2906 {
2907 mpz_neg(qq, qq);
2908 }
2909 return nlInitMPZ(qq,R);
2910 }
2911 mpz_t qq,rr;
2912 mpz_init(qq);
2913 mpz_init(rr);
2914 mpz_divmod(qq, rr, a->z, b->z);
2915 if (r!=NULL)
2916 *r = nlInitMPZ(rr,R);
2917 else
2918 {
2919 mpz_clear(rr);
2920 }
2921 return nlInitMPZ(qq,R);
2922}
2923
2924void nlInpGcd(number &a, number b, const coeffs r)
2925{
2926 if ((SR_HDL(b)|SR_HDL(a))&SR_INT)
2927 {
2928 number n=nlGcd(a,b,r);
2929 nlDelete(&a,r);
2930 a=n;
2931 }
2932 else
2933 {
2934 mpz_gcd(a->z,a->z,b->z);
2935 a=nlShort3_noinline(a);
2936 }
2937}
2938
2940{
2941 if ((SR_HDL(b)|SR_HDL(a))&SR_INT)
2942 {
2943 number n=nlIntDiv(a,b, r);
2944 nlDelete(&a,r);
2945 a=n;
2946 }
2947 else
2948 {
2949 mpz_t rr;
2950 mpz_init(rr);
2951 mpz_mod(rr,a->z,b->z);
2952 mpz_sub(a->z,a->z,rr);
2953 mpz_clear(rr);
2954 mpz_divexact(a->z,a->z,b->z);
2955 a=nlShort3_noinline(a);
2956 }
2957}
2958
2960{
2961 mpz_t A,B,C,D,E,N,P,tmp;
2963 else mpz_init_set(P,nP->z);
2964 const mp_bitcnt_t bits=2*(mpz_size1(P)+1)*GMP_LIMB_BITS;
2965 mpz_init2(N,bits);
2967 else mpz_set(N,nN->z);
2968 assume(!mpz_isNeg(P));
2969 if (mpz_isNeg(N)) mpz_add(N,N,P);
2970 mpz_init2(A,bits); mpz_set_ui(A,0L);
2971 mpz_init2(B,bits); mpz_set_ui(B,1L);
2972 mpz_init2(C,bits); mpz_set_ui(C,0L);
2973 mpz_init2(D,bits);
2974 mpz_init2(E,bits); mpz_set(E,P);
2976 number z=INT_TO_SR(0);
2977 while(mpz_sgn1(N)!=0)
2978 {
2979 mpz_mul(tmp,N,N);
2980 mpz_add(tmp,tmp,tmp);
2981 if (mpz_cmp(tmp,P)<0)
2982 {
2983 if (mpz_isNeg(B))
2984 {
2985 mpz_neg(B,B);
2986 mpz_neg(N,N);
2987 }
2988 // check for gcd(N,B)==1
2989 mpz_gcd(tmp,N,B);
2990 if (mpz_cmp_ui(tmp,1)==0)
2991 {
2992 // return N/B
2993 z=ALLOC_RNUMBER();
2994 #ifdef LDEBUG
2995 z->debug=123456;
2996 #endif
2997 memcpy(z->z,N,sizeof(mpz_t));
2998 memcpy(z->n,B,sizeof(mpz_t));
2999 z->s = 0;
3000 nlNormalize(z,r);
3001 }
3002 else
3003 {
3004 // return nN (the input) instead of "fail"
3005 z=nlCopy(nN,r);
3006 mpz_clear(B);
3007 mpz_clear(N);
3008 }
3009 break;
3010 }
3011 //mpz_mod(D,E,N);
3012 //mpz_div(tmp,E,N);
3013 mpz_divmod(tmp,D,E,N);
3014 mpz_mul(tmp,tmp,B);
3015 mpz_sub(C,A,tmp);
3016 mpz_set(E,N);
3017 mpz_set(N,D);
3018 mpz_set(A,B);
3019 mpz_set(B,C);
3020 }
3021 mpz_clear(tmp);
3022 mpz_clear(A);
3023 mpz_clear(C);
3024 mpz_clear(D);
3025 mpz_clear(E);
3026 mpz_clear(P);
3027 return z;
3028}
3029
3031{
3032 mpz_ptr aa,bb;
3033 *s=ALLOC_RNUMBER();
3034 mpz_init((*s)->z); (*s)->s=3;
3035 (*t)=ALLOC_RNUMBER();
3036 mpz_init((*t)->z); (*t)->s=3;
3038 mpz_init(g->z); g->s=3;
3039 #ifdef LDEBUG
3040 g->debug=123456;
3041 (*s)->debug=123456;
3042 (*t)->debug=123456;
3043 #endif
3044 if (SR_HDL(a) & SR_INT)
3045 {
3046 aa=(mpz_ptr)omAlloc(sizeof(mpz_t));
3048 }
3049 else
3050 {
3051 aa=a->z;
3052 }
3053 if (SR_HDL(b) & SR_INT)
3054 {
3055 bb=(mpz_ptr)omAlloc(sizeof(mpz_t));
3057 }
3058 else
3059 {
3060 bb=b->z;
3061 }
3062 mpz_gcdext(g->z,(*s)->z,(*t)->z,aa,bb);
3063 g=nlShort3(g);
3064 (*s)=nlShort3((*s));
3065 (*t)=nlShort3((*t));
3066 if (SR_HDL(a) & SR_INT)
3067 {
3068 mpz_clear(aa);
3069 omFreeSize(aa, sizeof(mpz_t));
3070 }
3071 if (SR_HDL(b) & SR_INT)
3072 {
3073 mpz_clear(bb);
3074 omFreeSize(bb, sizeof(mpz_t));
3075 }
3076 return g;
3077}
3078
3079//void nlCoeffWrite (const coeffs r, BOOLEAN /*details*/)
3080//{
3081// if (r->is_field) PrintS("QQ");
3082// else PrintS("ZZ");
3083//}
3084
3087// elements in the array are x[0..(rl-1)], q[0..(rl-1)]
3088{
3089 setCharacteristic( 0 ); // only in char 0
3091 CFArray X(rl), Q(rl);
3092 int i;
3093 for(i=rl-1;i>=0;i--)
3094 {
3095 X[i]=CF->convSingNFactoryN(x[i],FALSE,CF); // may be larger MAX_INT
3096 Q[i]=CF->convSingNFactoryN(q[i],FALSE,CF); // may be larger MAX_INT
3097 }
3099 if (n_SwitchChinRem)
3101 else
3103 number n=CF->convFactoryNSingN(xnew,CF);
3104 if (sym)
3105 {
3106 number p=CF->convFactoryNSingN(qnew,CF);
3107 number p2;
3108 if (getCoeffType(CF) == n_Q) p2=nlIntDiv(p,nlInit(2, CF),CF);
3109 else p2=CF->cfDiv(p,CF->cfInit(2, CF),CF);
3110 if (CF->cfGreater(n,p2,CF))
3111 {
3112 number n2=CF->cfSub(n,p,CF);
3113 CF->cfDelete(&n,CF);
3114 n=n2;
3115 }
3116 CF->cfDelete(&p2,CF);
3117 CF->cfDelete(&p,CF);
3118 }
3119 CF->cfNormalize(n,CF);
3120 return n;
3121}
3122#if 0
3123number nlChineseRemainder(number *x, number *q,int rl, const coeffs C)
3124{
3125 CFArray inv(rl);
3126 return nlChineseRemainderSym(x,q,rl,TRUE,inv,C);
3127}
3128#endif
3129
3131{
3132 assume(cf != NULL);
3133
3135
3136 if( !numberCollectionEnumerator.MoveNext() ) // empty zero polynomial?
3137 {
3138 c = nlInit(1, cf);
3139 return;
3140 }
3141
3142 // all coeffs are given by integers!!!
3143
3144 // part 1, find a small candidate for gcd
3146 int s1,s;
3147 s=2147483647; // max. int
3148
3150
3151 int normalcount = 0;
3152 do
3153 {
3154 number& n = numberCollectionEnumerator.Current();
3155 nlNormalize(n, cf); ++normalcount;
3156 cand1 = n;
3157
3158 if (SR_HDL(cand1)&SR_INT) { cand=cand1; break; }
3159 assume(cand1->s==3); // all coeffs should be integers // ==0?!! after printing
3160 s1=mpz_size1(cand1->z);
3161 if (s>s1)
3162 {
3163 cand=cand1;
3164 s=s1;
3165 }
3166 } while (numberCollectionEnumerator.MoveNext() );
3167
3168// assume( nlGreaterZero(cand,cf) ); // cand may be a negative integer!
3169
3170 cand=nlCopy(cand,cf);
3171 // part 2: compute gcd(cand,all coeffs)
3172
3174
3175 while (numberCollectionEnumerator.MoveNext() )
3176 {
3177 number& n = numberCollectionEnumerator.Current();
3178
3179 if( (--normalcount) <= 0)
3180 nlNormalize(n, cf);
3181
3182 nlInpGcd(cand, n, cf);
3184
3185 if(nlIsOne(cand,cf))
3186 {
3187 c = cand;
3188
3189 if(!lc_is_pos)
3190 {
3191 // make the leading coeff positive
3192 c = nlNeg(c, cf);
3194
3195 while (numberCollectionEnumerator.MoveNext() )
3196 {
3198 nn = nlNeg(nn, cf);
3199 }
3200 }
3201 return;
3202 }
3203 }
3204
3205 // part3: all coeffs = all coeffs / cand
3206 if (!lc_is_pos)
3207 cand = nlNeg(cand,cf);
3208
3209 c = cand;
3211
3212 while (numberCollectionEnumerator.MoveNext() )
3213 {
3214 number& n = numberCollectionEnumerator.Current();
3215 number t=nlExactDiv(n, cand, cf); // simple integer exact division, no ratios to remain
3216 nlDelete(&n, cf);
3217 n = t;
3218 }
3219}
3220
3222{
3223 assume(cf != NULL);
3224
3226
3227 if( !numberCollectionEnumerator.MoveNext() ) // empty zero polynomial?
3228 {
3229 c = nlInit(1, cf);
3230// assume( n_GreaterZero(c, cf) );
3231 return;
3232 }
3233
3234 // all coeffs are given by integers after returning from this routine
3235
3236 // part 1, collect product of all denominators /gcds
3237 number cand;
3239#if defined(LDEBUG)
3240 cand->debug=123456;
3241#endif
3242 cand->s=3;
3243
3244 int s=0;
3245
3247
3248 do
3249 {
3251
3252 if (!(SR_HDL(cand1)&SR_INT))
3253 {
3255 if ((!(SR_HDL(cand1)&SR_INT)) // not a short int
3256 && (cand1->s==1)) // and is a normalised rational
3257 {
3258 if (s==0) // first denom, we meet
3259 {
3260 mpz_init_set(cand->z, cand1->n); // cand->z = cand1->n
3261 s=1;
3262 }
3263 else // we have already something
3264 {
3265 mpz_lcm(cand->z, cand->z, cand1->n);
3266 }
3267 }
3268 }
3269 }
3270 while (numberCollectionEnumerator.MoveNext() );
3271
3272
3273 if (s==0) // nothing to do, all coeffs are already integers
3274 {
3275// mpz_clear(tmp);
3277 if (lc_is_pos)
3278 c=nlInit(1,cf);
3279 else
3280 {
3281 // make the leading coeff positive
3282 c=nlInit(-1,cf);
3283
3284 // TODO: incorporate the following into the loop below?
3286 while (numberCollectionEnumerator.MoveNext() )
3287 {
3288 number& n = numberCollectionEnumerator.Current();
3289 n = nlNeg(n, cf);
3290 }
3291 }
3292// assume( n_GreaterZero(c, cf) );
3293 return;
3294 }
3295
3296 cand = nlShort3(cand);
3297
3298 // part2: all coeffs = all coeffs * cand
3299 // make the lead coeff positive
3301
3302 if (!lc_is_pos)
3303 cand = nlNeg(cand, cf);
3304
3305 c = cand;
3306
3307 while (numberCollectionEnumerator.MoveNext() )
3308 {
3309 number &n = numberCollectionEnumerator.Current();
3310 nlInpMult(n, cand, cf);
3311 }
3312
3313}
3314
3315char * nlCoeffName(const coeffs r)
3316{
3317 if (r->cfDiv==nlDiv) return (char*)"QQ";
3318 else return (char*)"ZZ";
3319}
3320
3321void nlWriteFd(number n, const ssiInfo* d, const coeffs)
3322{
3323 if(SR_HDL(n) & SR_INT)
3324 {
3325 #if SIZEOF_LONG == 4
3326 fprintf(d->f_write,"4 %ld ",SR_TO_INT(n));
3327 #else
3328 long nn=SR_TO_INT(n);
3329 if ((nn<POW_2_28_32)&&(nn>= -POW_2_28_32))
3330 {
3331 int nnn=(int)nn;
3332 fprintf(d->f_write,"4 %d ",nnn);
3333 }
3334 else
3335 {
3336 mpz_t tmp;
3338 fputs("8 ",d->f_write);
3340 fputc(' ',d->f_write);
3341 mpz_clear(tmp);
3342 }
3343 #endif
3344 }
3345 else if (n->s<2)
3346 {
3347 //gmp_fprintf(f,"%d %Zd %Zd ",n->s,n->z,n->n);
3348 fprintf(d->f_write,"%d ",n->s+5);
3349 mpz_out_str (d->f_write,SSI_BASE, n->z);
3350 fputc(' ',d->f_write);
3351 mpz_out_str (d->f_write,SSI_BASE, n->n);
3352 fputc(' ',d->f_write);
3353
3354 //if (d->f_debug!=NULL) gmp_fprintf(d->f_debug,"number: s=%d gmp/gmp \"%Zd %Zd\" ",n->s,n->z,n->n);
3355 }
3356 else /*n->s==3*/
3357 {
3358 //gmp_fprintf(d->f_write,"3 %Zd ",n->z);
3359 fputs("8 ",d->f_write);
3360 mpz_out_str (d->f_write,SSI_BASE, n->z);
3361 fputc(' ',d->f_write);
3362
3363 //if (d->f_debug!=NULL) gmp_fprintf(d->f_debug,"number: gmp \"%Zd\" ",n->z);
3364 }
3365}
3366
3368{
3369 int sub_type=-1;
3371 switch(sub_type)
3372 {
3373 case 0:
3374 case 1:
3375 {// read mpz_t, mpz_t
3376 number n=nlRInit(0);
3377 mpz_init(n->n);
3378 s_readmpz(d->f_read,n->z);
3379 s_readmpz(d->f_read,n->n);
3380 n->s=sub_type;
3381 return n;
3382 }
3383
3384 case 3:
3385 {// read mpz_t
3386 number n=nlRInit(0);
3387 s_readmpz(d->f_read,n->z);
3388 n->s=3; /*sub_type*/
3389 #if SIZEOF_LONG == 8
3390 n=nlShort3(n);
3391 #endif
3392 return n;
3393 }
3394 case 4:
3395 {
3397 //#if SIZEOF_LONG == 8
3398 return INT_TO_SR(dd);
3399 //#else
3400 //return nlInit(dd,NULL);
3401 //#endif
3402 }
3403 case 5:
3404 case 6:
3405 {// read raw mpz_t, mpz_t
3406 number n=nlRInit(0);
3407 mpz_init(n->n);
3408 s_readmpz_base (d->f_read,n->z, SSI_BASE);
3409 s_readmpz_base (d->f_read,n->n, SSI_BASE);
3410 n->s=sub_type-5;
3411 return n;
3412 }
3413 case 8:
3414 {// read raw mpz_t
3415 number n=nlRInit(0);
3416 s_readmpz_base (d->f_read,n->z, SSI_BASE);
3417 n->s=sub_type=3; /*subtype-5*/
3418 #if SIZEOF_LONG == 8
3419 n=nlShort3(n);
3420 #endif
3421 return n;
3422 }
3423
3424 default: Werror("error in reading number: invalid subtype %d",sub_type);
3425 return NULL;
3426 }
3427 return NULL;
3428}
3429
3431{
3432 /* test, if r is an instance of nInitCoeffs(n,parameter) */
3433 /* if parameter is not needed */
3434 if (n==r->type)
3435 {
3436 if ((p==NULL)&&(r->cfDiv==nlDiv)) return TRUE;
3437 if ((p!=NULL)&&(r->cfDiv!=nlDiv)) return TRUE;
3438 }
3439 return FALSE;
3440}
3441
3443{
3444 number g=nlGcd(a,b,r);
3445 number n1=nlMult(a,b,r);
3446 number n2=nlExactDiv(n1,g,r);
3447 nlDelete(&g,r);
3448 nlDelete(&n1,r);
3449 return n2;
3450}
3451
3453{
3454 number a=nlInit(p(),cf);
3455 if (v2!=NULL)
3456 {
3457 number b=nlInit(p(),cf);
3458 number c=nlDiv(a,b,cf);
3459 nlDelete(&b,cf);
3460 nlDelete(&a,cf);
3461 a=c;
3462 }
3463 return a;
3464}
3465
3467{
3468 r->is_domain=TRUE;
3469 r->rep=n_rep_gap_rat;
3470
3471 r->nCoeffIsEqual=nlCoeffIsEqual;
3472 //r->cfKillChar = ndKillChar; /* dummy */
3473 //r->cfCoeffString=nlCoeffString;
3474 r->cfCoeffName=nlCoeffName;
3475
3476 r->cfInitMPZ = nlInitMPZ;
3477 r->cfMPZ = nlMPZ;
3478
3479 r->cfMult = nlMult;
3480 r->cfSub = nlSub;
3481 r->cfAdd = nlAdd;
3482 r->cfExactDiv= nlExactDiv;
3483 if (p==NULL) /* Q */
3484 {
3485 r->is_field=TRUE;
3486 r->cfDiv = nlDiv;
3487 //r->cfGcd = ndGcd_dummy;
3488 r->cfSubringGcd = nlGcd;
3489 }
3490 else /* Z: coeffs_BIGINT */
3491 {
3492 r->is_field=FALSE;
3493 r->cfDiv = nlIntDiv;
3494 r->cfIntMod= nlIntMod;
3495 r->cfGcd = nlGcd;
3496 r->cfDivBy=nlDivBy;
3497 r->cfDivComp = nlDivComp;
3498 r->cfIsUnit = nlIsUnit;
3499 r->cfGetUnit = nlGetUnit;
3500 r->cfQuot1 = nlQuot1;
3501 r->cfLcm = nlLcm;
3502 r->cfXExtGcd=nlXExtGcd;
3503 r->cfQuotRem=nlQuotRem;
3504 }
3505 r->cfInit = nlInit;
3506 r->cfSize = nlSize;
3507 r->cfInt = nlInt;
3508
3509 r->cfChineseRemainder=nlChineseRemainderSym;
3510 r->cfFarey=nlFarey;
3511 r->cfInpNeg = nlNeg;
3512 r->cfInvers= nlInvers;
3513 r->cfCopy = nlCopy;
3514 r->cfRePart = nlCopy;
3515 //r->cfImPart = ndReturn0;
3516 r->cfWriteLong = nlWrite;
3517 r->cfRead = nlRead;
3518 r->cfNormalize=nlNormalize;
3519 r->cfGreater = nlGreater;
3520 r->cfEqual = nlEqual;
3521 r->cfIsZero = nlIsZero;
3522 r->cfIsOne = nlIsOne;
3523 r->cfIsMOne = nlIsMOne;
3524 r->cfGreaterZero = nlGreaterZero;
3525 r->cfPower = nlPower;
3526 r->cfGetDenom = nlGetDenom;
3527 r->cfGetNumerator = nlGetNumerator;
3528 r->cfExtGcd = nlExtGcd; // only for ring stuff and Z
3529 r->cfNormalizeHelper = nlNormalizeHelper;
3530 r->cfDelete= nlDelete;
3531 r->cfSetMap = nlSetMap;
3532 //r->cfName = ndName;
3533 r->cfInpMult=nlInpMult;
3534 r->cfInpAdd=nlInpAdd;
3535 //r->cfCoeffWrite=nlCoeffWrite;
3536
3537 r->cfClearContent = nlClearContent;
3538 r->cfClearDenominators = nlClearDenominators;
3539
3540#ifdef LDEBUG
3541 // debug stuff
3542 r->cfDBTest=nlDBTest;
3543#endif
3544 r->convSingNFactoryN=nlConvSingNFactoryN;
3545 r->convFactoryNSingN=nlConvFactoryNSingN;
3546
3547 r->cfRandom=nlRandom;
3548
3549 // io via ssi
3550 r->cfWriteFd=nlWriteFd;
3551 r->cfReadFd=nlReadFd;
3552
3553 //r->type = n_Q;
3554 r->ch = 0;
3555 r->has_simple_Alloc=FALSE;
3556 r->has_simple_Inverse=FALSE;
3557
3558 // variables for this type of coeffs:
3559 // (none)
3560 return FALSE;
3561}
3562#if 0
3564{
3565 if (((SR_HDL(b)&SR_HDL(a))&SR_INT)
3566 {
3567 int bi=SR_TO_INT(b);
3568 int ai=SR_TO_INT(a);
3569 int bb=ABS(bi);
3570 int c=ai%bb;
3571 if (c<0) c+=bb;
3572 return (INT_TO_SR(c));
3573 }
3574 number al;
3575 number bl;
3576 if (SR_HDL(a))&SR_INT)
3577 al=nlRInit(SR_TO_INT(a));
3578 else
3579 al=nlCopy(a);
3580 if (SR_HDL(b))&SR_INT)
3582 else
3583 bl=nlCopy(b);
3584 number r=nlRInit(0);
3585 mpz_mod(r->z,al->z,bl->z);
3586 nlDelete(&al);
3587 nlDelete(&bl);
3588 if (mpz_size1(&r->z)<=MP_SMALL)
3589 {
3590 LONG ui=(int)mpz_get_si(&r->z);
3591 if ((((ui<<3)>>3)==ui)
3592 && (mpz_cmp_si(x->z,(long)ui)==0))
3593 {
3594 mpz_clear(&r->z);
3595 FREE_RNUMBER(r); // omFreeBin((void *)r, rnumber_bin);
3596 r=INT_TO_SR(ui);
3597 }
3598 }
3599 return r;
3600}
3601#endif
3602#endif // not P_NUMBERS_H
3603#endif // LONGRAT_CC
All the auxiliary stuff.
#define SSI_BASE
Definition auxiliary.h:135
static int ABS(int v)
Definition auxiliary.h:112
int BOOLEAN
Definition auxiliary.h:87
#define TRUE
Definition auxiliary.h:100
#define FALSE
Definition auxiliary.h:96
void On(int sw)
switches
void Off(int sw)
switches
int size(const CanonicalForm &f, const Variable &v)
int size ( const CanonicalForm & f, const Variable & v )
Definition cf_ops.cc:600
void FACTORY_PUBLIC setCharacteristic(int c)
Definition cf_char.cc:28
CanonicalForm num(const CanonicalForm &f)
CanonicalForm den(const CanonicalForm &f)
const CanonicalForm CFMap CFMap & N
Definition cfEzgcd.cc:56
int l
Definition cfEzgcd.cc:100
int m
Definition cfEzgcd.cc:128
int i
Definition cfEzgcd.cc:132
const CanonicalForm const CanonicalForm const CanonicalForm const CanonicalForm & cand
Definition cfModGcd.cc:70
Variable x
Definition cfModGcd.cc:4090
int p
Definition cfModGcd.cc:4086
g
Definition cfModGcd.cc:4098
CanonicalForm cf
Definition cfModGcd.cc:4091
CanonicalForm b
Definition cfModGcd.cc:4111
void FACTORY_PUBLIC chineseRemainder(const CanonicalForm &x1, const CanonicalForm &q1, const CanonicalForm &x2, const CanonicalForm &q2, CanonicalForm &xnew, CanonicalForm &qnew)
void chineseRemainder ( const CanonicalForm & x1, const CanonicalForm & q1, const CanonicalForm & x2,...
Definition cf_chinese.cc:57
void FACTORY_PUBLIC chineseRemainderCached(const CanonicalForm &x1, const CanonicalForm &q1, const CanonicalForm &x2, const CanonicalForm &q2, CanonicalForm &xnew, CanonicalForm &qnew, CFArray &inv)
static const int SW_RATIONAL
set to 1 for computations over Q
Definition cf_defs.h:31
FILE * f
Definition checklibs.c:9
factory's main class
gmp_complex numbers based on
Coefficient rings, fields and other domains suitable for Singular polynomials.
static FORCE_INLINE BOOLEAN nCoeff_is_long_R(const coeffs r)
Definition coeffs.h:884
n_coeffType
Definition coeffs.h:27
@ n_R
single prescision (6,6) real numbers
Definition coeffs.h:31
@ n_Q
rational (GMP) numbers
Definition coeffs.h:30
@ n_Zn
only used if HAVE_RINGS is defined
Definition coeffs.h:44
@ n_long_R
real floating point (GMP) numbers
Definition coeffs.h:33
@ n_Zp
\F{p < 2^31}
Definition coeffs.h:29
@ n_long_C
complex floating point (GMP) numbers
Definition coeffs.h:41
static FORCE_INLINE number n_Div(number a, number b, const coeffs r)
return the quotient of 'a' and 'b', i.e., a/b; raises an error if 'b' is not invertible in r exceptio...
Definition coeffs.h:616
coeffs nInitChar(n_coeffType t, void *parameter)
one-time initialisations for new coeffs in case of an error return NULL
Definition numbers.cc:406
#define ALLOC_RNUMBER()
Definition coeffs.h:94
static FORCE_INLINE n_coeffType getCoeffType(const coeffs r)
Returns the type of coeffs domain.
Definition coeffs.h:429
static FORCE_INLINE int n_GetChar(const coeffs r)
Return the characteristic of the coeff. domain.
Definition coeffs.h:448
static FORCE_INLINE void n_Delete(number *p, const coeffs r)
delete 'p'
Definition coeffs.h:459
static FORCE_INLINE BOOLEAN nCoeff_is_Zp(const coeffs r)
Definition coeffs.h:793
static FORCE_INLINE number n_Init(long i, const coeffs r)
a number representing i in the given coeff field/ring r
Definition coeffs.h:539
static FORCE_INLINE BOOLEAN nCoeff_is_Ring_2toM(const coeffs r)
Definition coeffs.h:724
#define FREE_RNUMBER(x)
Definition coeffs.h:93
@ n_rep_gap_rat
(number), see longrat.h
Definition coeffs.h:118
@ n_rep_gap_gmp
(), see rinteger.h, new impl.
Definition coeffs.h:119
@ n_rep_float
(float), see shortfl.h
Definition coeffs.h:123
@ n_rep_int
(int), see modulop.h
Definition coeffs.h:117
@ n_rep_gmp_float
(gmp_float), see
Definition coeffs.h:124
@ n_rep_gmp
(mpz_ptr), see rmodulon,h
Definition coeffs.h:122
#define ALLOC0_RNUMBER()
Definition coeffs.h:95
number(* nMapFunc)(number a, const coeffs src, const coeffs dst)
maps "a", which lives in src, into dst
Definition coeffs.h:80
static FORCE_INLINE BOOLEAN nCoeff_is_R(const coeffs r)
Definition coeffs.h:829
static FORCE_INLINE BOOLEAN nCoeff_is_long_C(const coeffs r)
Definition coeffs.h:887
#define Print
Definition emacs.cc:80
#define WarnS
Definition emacs.cc:78
return result
const CanonicalForm int s
Definition facAbsFact.cc:51
const CanonicalForm int const CFList const Variable & y
Definition facAbsFact.cc:53
CanonicalForm res
Definition facAbsFact.cc:60
REvaluation E(1, terms.length(), IntRandom(25))
b *CanonicalForm B
Definition facBivar.cc:52
const Variable & v
< [in] a sqrfree bivariate poly
Definition facBivar.h:39
int j
Definition facHensel.cc:110
bool isZero(const CFArray &A)
checks if entries of A are zero
CanonicalForm FACTORY_PUBLIC make_cf(const mpz_ptr n)
Definition singext.cc:66
void FACTORY_PUBLIC gmp_numerator(const CanonicalForm &f, mpz_ptr result)
Definition singext.cc:20
void FACTORY_PUBLIC gmp_denominator(const CanonicalForm &f, mpz_ptr result)
Definition singext.cc:40
void WerrorS(const char *s)
Definition feFopen.cc:24
#define D(A)
Definition gentable.cc:128
#define VAR
Definition globaldefs.h:5
#define info
Definition libparse.cc:1256
static number nlMapP(number from, const coeffs src, const coeffs dst)
Definition longrat.cc:189
#define nlTest(a, r)
Definition longrat.cc:87
void nlWriteFd(number n, const ssiInfo *d, const coeffs)
Definition longrat.cc:3321
LINLINE void nlInpMult(number &a, number b, const coeffs r)
Definition longrat.cc:2776
LINLINE BOOLEAN nlEqual(number a, number b, const coeffs r)
Definition longrat.cc:2588
LINLINE number nlAdd(number la, number li, const coeffs r)
Definition longrat.cc:2692
number nlMapZ(number from, const coeffs, const coeffs dst)
Definition longrat.cc:210
long nlInt(number &n, const coeffs r)
Definition longrat.cc:740
static number nlLcm(number a, number b, const coeffs r)
Definition longrat.cc:3442
static number nlMapLongR_BI(number from, const coeffs src, const coeffs dst)
Definition longrat.cc:512
number nlInit2(int i, int j, const coeffs r)
create a rational i/j (implicitly) over Q NOTE: make sure to use correct Q in debug mode
Definition longrat.cc:2535
#define POW_2_28
Definition longrat.cc:103
LINLINE number nl_Copy(number a, const coeffs r)
number nlInit2gmp(mpz_t i, mpz_t j, const coeffs r)
create a rational i/j (implicitly) over Q NOTE: make sure to use correct Q in debug mode
Definition longrat.cc:2548
void _nlInpAdd_aNoImm_OR_bNoImm(number &a, number b)
Definition longrat.cc:1970
LINLINE number nlSub(number la, number li, const coeffs r)
Definition longrat.cc:2758
number nlIntMod(number a, number b, const coeffs r)
Definition longrat.cc:1016
number _nlCopy_NoImm(number a)
Definition longrat.cc:1740
number _nlSub_aNoImm_OR_bNoImm(number a, number b)
Definition longrat.cc:2113
LINLINE number nlCopy(number a, const coeffs r)
Definition longrat.cc:2644
LINLINE number nlNeg(number za, const coeffs r)
Definition longrat.cc:2673
number nlXExtGcd(number a, number b, number *s, number *t, number *u, number *v, const coeffs r)
Definition longrat.cc:2819
void nlPower(number x, int exp, number *lu, const coeffs r)
Definition longrat.cc:1250
number nlQuotRem(number a, number b, number *r, const coeffs R)
Definition longrat.cc:2871
number nlFarey(number nN, number nP, const coeffs CF)
Definition longrat.cc:2959
LINLINE BOOLEAN nlIsOne(number a, const coeffs r)
Definition longrat.cc:2615
#define mpz_isNeg(A)
Definition longrat.cc:146
static number nlMapC(number from, const coeffs src, const coeffs dst)
Definition longrat.cc:545
number nlNormalizeHelper(number a, number b, const coeffs r)
Definition longrat.cc:1525
LINLINE void nlDelete(number *a, const coeffs r)
Definition longrat.cc:2657
BOOLEAN nlGreaterZero(number za, const coeffs r)
Definition longrat.cc:1303
number _nlNeg_NoImm(number a)
Definition longrat.cc:1779
number nlModP(number q, const coeffs, const coeffs Zp)
Definition longrat.cc:1572
LINLINE void nlInpAdd(number &a, number b, const coeffs r)
Definition longrat.cc:2710
number nlExactDiv(number a, number b, const coeffs r)
Definition longrat.cc:870
void mpz_mul_si(mpz_ptr r, mpz_srcptr s, long int si)
Definition longrat.cc:177
VAR int n_SwitchChinRem
Definition longrat.cc:3085
const char * nlRead(const char *s, number *a, const coeffs r)
Definition longrat0.cc:31
void nlMPZ(mpz_t m, number &n, const coeffs r)
Definition longrat.cc:2810
number nlInvers(number a, const coeffs r)
Definition longrat.cc:790
BOOLEAN nlIsUnit(number a, const coeffs)
Definition longrat.cc:1131
void nlInpIntDiv(number &a, number b, const coeffs r)
Definition longrat.cc:2939
static void nlNormalize_Gcd(number &x)
Definition longrat.cc:1792
static number nlConvFactoryNSingN(const CanonicalForm f, const coeffs r)
Definition longrat.cc:365
number nlChineseRemainderSym(number *x, number *q, int rl, BOOLEAN sym, CFArray &inv_cache, const coeffs CF)
Definition longrat.cc:3086
int nlDivComp(number a, number b, const coeffs r)
Definition longrat.cc:1091
void _nlDelete_NoImm(number *a)
Definition longrat.cc:1761
#define LINLINE
Definition longrat.cc:31
char * nlCoeffName(const coeffs r)
Definition longrat.cc:3315
#define POW_2_28_32
Definition longrat.cc:104
BOOLEAN nlInitChar(coeffs r, void *p)
Definition longrat.cc:3466
number nlCopyMap(number a, const coeffs, const coeffs)
Definition longrat.cc:2445
number nlExtGcd(number a, number b, number *s, number *t, const coeffs)
Definition longrat.cc:3030
static number nlMapGMP(number from, const coeffs, const coeffs dst)
Definition longrat.cc:205
LINLINE number nlMult(number a, number b, const coeffs r)
Definition longrat.cc:2728
static number nlInitMPZ(mpz_t m, const coeffs)
Definition longrat.cc:164
number nlIntDiv(number a, number b, const coeffs r)
Definition longrat.cc:935
static void nlClearDenominators(ICoeffsEnumerator &numberCollectionEnumerator, number &c, const coeffs cf)
Definition longrat.cc:3221
static number nlMapLongR(number from, const coeffs src, const coeffs dst)
Definition longrat.cc:432
LINLINE BOOLEAN nlIsZero(number za, const coeffs r)
Definition longrat.cc:2624
number nlGetDenom(number &n, const coeffs r)
Definition longrat.cc:1633
number nlGcd(number a, number b, const coeffs r)
Definition longrat.cc:1340
number _nlMult_aImm_bImm_rNoImm(number a, number b)
Definition longrat.cc:2324
number nlReadFd(const ssiInfo *d, const coeffs)
Definition longrat.cc:3367
int nlSize(number a, const coeffs)
Definition longrat.cc:711
number nlMapMachineInt(number from, const coeffs, const coeffs)
Definition longrat.cc:222
nMapFunc nlSetMap(const coeffs src, const coeffs dst)
Definition longrat.cc:2473
number nlBigInt(number &n)
static number nlShort3(number x)
Definition longrat.cc:109
#define GCD_NORM_COND(OLD, NEW)
Definition longrat.cc:1790
BOOLEAN nlDBTest(number a, const char *f, const int l)
number nlDiv(number a, number b, const coeffs r)
Definition longrat.cc:1140
number nlRInit(long i)
Definition longrat.cc:2521
BOOLEAN nlIsMOne(number a, const coeffs r)
Definition longrat.cc:1328
static void nlClearContent(ICoeffsEnumerator &numberCollectionEnumerator, number &c, const coeffs cf)
Definition longrat.cc:3130
number _nlMult_aNoImm_OR_bNoImm(number a, number b)
Definition longrat.cc:2337
LINLINE number nlInit(long i, const coeffs r)
Definition longrat.cc:2597
number nlShort3_noinline(number x)
Definition longrat.cc:159
number nlGetNumerator(number &n, const coeffs r)
Definition longrat.cc:1662
number _nlAdd_aNoImm_OR_bNoImm(number a, number b)
Definition longrat.cc:1812
#define LONG
Definition longrat.cc:105
BOOLEAN nlCoeffIsEqual(const coeffs r, n_coeffType n, void *p)
Definition longrat.cc:3430
static CanonicalForm nlConvSingNFactoryN(number n, const BOOLEAN setChar, const coeffs)
Definition longrat.cc:327
static number nlMapR(number from, const coeffs src, const coeffs dst)
Definition longrat.cc:392
number nlGetUnit(number n, const coeffs cf)
Definition longrat.cc:1102
coeffs nlQuot1(number c, const coeffs r)
Definition longrat.cc:1108
BOOLEAN _nlEqual_aNoImm_OR_bNoImm(number a, number b)
Definition longrat.cc:1693
number nlShort1(number x)
Definition longrat.cc:1460
#define MP_SMALL
Definition longrat.cc:144
BOOLEAN nlGreater(number a, number b, const coeffs r)
Definition longrat.cc:1313
static number nlMapR_BI(number from, const coeffs src, const coeffs dst)
Definition longrat.cc:422
void nlGMP(number &i, mpz_t n, const coeffs r)
Definition longrat.cc:1613
void nlNormalize(number &x, const coeffs r)
Definition longrat.cc:1481
BOOLEAN nlDivBy(number a, number b, const coeffs)
Definition longrat.cc:1077
static int int_extgcd(int a, int b, int *u, int *x, int *v, int *y)
Definition longrat.cc:1410
void nlWrite(number a, const coeffs r)
Definition longrat0.cc:90
void nlInpGcd(number &a, number b, const coeffs r)
Definition longrat.cc:2924
static number nlRandom(siRandProc p, number v2, number, const coeffs cf)
Definition longrat.cc:3452
number nlMapQtoZ(number a, const coeffs src, const coeffs dst)
Definition longrat.cc:2454
#define SR_INT
Definition longrat.h:67
#define INT_TO_SR(INT)
Definition longrat.h:68
#define SR_TO_INT(SR)
Definition longrat.h:69
void dErrorBreak(void)
Definition dError.cc:140
#define assume(x)
Definition mod2.h:389
long npInt(number &n, const coeffs r)
Definition modulop.cc:83
char * floatToStr(const gmp_float &r, const unsigned int oprec)
gmp_float exp(const gmp_float &a)
The main handler for Singular numbers which are suitable for Singular polynomials.
char * nEatLong(char *s, mpz_ptr i)
extracts a long integer from s, returns the rest
Definition numbers.cc:706
const char *const nDivBy0
Definition numbers.h:89
#define omFreeSize(addr, size)
#define omAlloc(size)
#define omCheckIf(cond, test)
#define omCheckAddrSize(addr, size)
#define omFree(addr)
#define NULL
Definition omList.c:12
int IsPrime(int p)
Definition prime.cc:61
void Werror(const char *fmt,...)
Definition reporter.cc:189
void s_readmpz(s_buff F, mpz_t a)
Definition s_buff.cc:184
void s_readmpz_base(s_buff F, mpz_ptr a, int base)
Definition s_buff.cc:209
int s_readint(s_buff F)
Definition s_buff.cc:112
long s_readlong(s_buff F)
Definition s_buff.cc:140
s_buff f_read
Definition s_buff.h:22
FILE * f_write
Definition s_buff.h:23
SI_FLOAT nrFloat(number n)
Converts a n_R number into a float. Needed by Maps.
Definition shortfl.cc:48
#define mpz_size1(A)
Definition si_gmp.h:17
#define mpz_sgn1(A)
Definition si_gmp.h:18
#define R
Definition sirandom.c:27
#define A
Definition sirandom.c:24
#define Q
Definition sirandom.c:26
int(* siRandProc)(void)
Definition sirandom.h:9
#define SR_HDL(A)
Definition tgb.cc:35
int gcd(int a, int b)