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