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