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