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