source: git/libpolys/coeffs/rintegers2.cc @ 9e2cf88

spielwiese
Last change on this file since 9e2cf88 was 9e2cf88, checked in by Hans Schoenemann <hannes@…>, 5 years ago
split rintegers.cc
  • Property mode set to 100644
File size: 14.0 KB
Line 
1/****************************************
2*  Computer Algebra System SINGULAR     *
3****************************************/
4/*
5* ABSTRACT: numbers (integers)
6*/
7
8
9#ifdef HAVE_RINGS
10#if SI_INTEGER_VARIANT == 2
11
12#include "coeffs/si_gmp.h"
13
14/*
15 * Multiply two numbers
16 */
17static number nrzMult (number a, number b, const coeffs)
18{
19  mpz_ptr erg = (mpz_ptr) omAllocBin(gmp_nrz_bin);
20  mpz_init(erg);
21  mpz_mul(erg, (mpz_ptr) a, (mpz_ptr) b);
22  return (number) erg;
23}
24
25/*
26 * Give the smallest non unit k, such that a * x = k = b * y has a solution
27 */
28static number nrzLcm (number a,number b,const coeffs)
29{
30  mpz_ptr erg = (mpz_ptr) omAllocBin(gmp_nrz_bin);
31  mpz_init(erg);
32  mpz_lcm(erg, (mpz_ptr) a, (mpz_ptr) b);
33  return (number) erg;
34}
35
36/*
37 * Give the largest non unit k, such that a = x * k, b = y * k has
38 * a solution.
39 */
40static number nrzGcd (number a,number b,const coeffs)
41{
42  mpz_ptr erg = (mpz_ptr) omAllocBin(gmp_nrz_bin);
43  mpz_init(erg);
44  mpz_gcd(erg, (mpz_ptr) a, (mpz_ptr) b);
45  return (number) erg;
46}
47
48/*
49 * Give the largest non unit k, such that a = x * k, b = y * k has
50 * a solution and r, s, s.t. k = s*a + t*b
51 */
52static number nrzExtGcd (number a, number b, number *s, number *t, const coeffs)
53{
54  mpz_ptr erg = (mpz_ptr) omAllocBin(gmp_nrz_bin);
55  mpz_ptr bs = (mpz_ptr) omAllocBin(gmp_nrz_bin);
56  mpz_ptr bt = (mpz_ptr) omAllocBin(gmp_nrz_bin);
57  mpz_init(erg);
58  mpz_init(bs);
59  mpz_init(bt);
60  mpz_gcdext(erg, bs, bt, (mpz_ptr) a, (mpz_ptr) b);
61  *s = (number) bs;
62  *t = (number) bt;
63  return (number) erg;
64}
65
66static number nrzXExtGcd (number a, number b, number *s, number *t, number *u, number *v, const coeffs )
67{
68  mpz_ptr erg = (mpz_ptr) omAllocBin(gmp_nrz_bin);
69  mpz_ptr bs = (mpz_ptr) omAllocBin(gmp_nrz_bin);
70  mpz_ptr bt = (mpz_ptr) omAllocBin(gmp_nrz_bin);
71  mpz_init(erg);
72  mpz_init(bs);
73  mpz_init(bt);
74
75  mpz_gcdext(erg, bs, bt, (mpz_ptr)a, (mpz_ptr)b);
76
77  mpz_ptr bu = (mpz_ptr) omAllocBin(gmp_nrz_bin);
78  mpz_ptr bv = (mpz_ptr) omAllocBin(gmp_nrz_bin);
79
80  mpz_init_set(bu, (mpz_ptr) b);
81  mpz_init_set(bv, (mpz_ptr) a);
82
83  assume(mpz_cmp_si(erg, 0));
84
85  mpz_div(bu, bu, erg);
86  mpz_div(bv, bv, erg);
87
88  mpz_mul_si(bu, bu, -1);
89  *u = (number) bu;
90  *v = (number) bv;
91
92  *s = (number) bs;
93  *t = (number) bt;
94  return (number) erg;
95}
96
97static void nrzPower (number a, int i, number * result, const coeffs)
98{
99  mpz_ptr erg = (mpz_ptr) omAllocBin(gmp_nrz_bin);
100  mpz_init(erg);
101  mpz_pow_ui(erg, (mpz_ptr) a, i);
102  *result = (number) erg;
103}
104
105/*
106 * create a number from int
107 */
108number nrzInit (long i, const coeffs)
109{
110  mpz_ptr erg = (mpz_ptr) omAllocBin(gmp_nrz_bin);
111  mpz_init_set_si(erg, i);
112  return (number) erg;
113}
114
115void nrzDelete(number *a, const coeffs)
116{
117  if (*a != NULL)
118  {
119    mpz_clear((mpz_ptr) *a);
120    omFreeBin((ADDRESS) *a, gmp_nrz_bin);
121    *a = NULL;
122  }
123}
124
125static number nrzCopy(number a, const coeffs)
126{
127  if (a==NULL) return NULL;
128  mpz_ptr erg = (mpz_ptr) omAllocBin(gmp_nrz_bin);
129  mpz_init_set(erg, (mpz_ptr) a);
130  return (number) erg;
131}
132
133#if 0
134number nrzCopyMap(number a, const coeffs /*src*/, const coeffs dst)
135{
136  return nrzCopy(a,dst);
137}
138#endif
139
140int nrzSize(number a, const coeffs)
141{
142  mpz_ptr p=(mpz_ptr)a;
143  int s=p->_mp_alloc;
144  if (s==1) s=(mpz_cmp_ui(p,0)!=0);
145  return s;
146
147}
148
149/*
150 * convert a number to int
151 */
152static long nrzInt(number &n, const coeffs)
153{
154  return mpz_get_si( (mpz_ptr)n);
155}
156
157static number nrzAdd (number a, number b, const coeffs)
158{
159  mpz_ptr erg = (mpz_ptr) omAllocBin(gmp_nrz_bin);
160  mpz_init(erg);
161  mpz_add(erg, (mpz_ptr) a, (mpz_ptr) b);
162  return (number) erg;
163}
164
165static number nrzSub (number a, number b, const coeffs)
166{
167  mpz_ptr erg = (mpz_ptr) omAllocBin(gmp_nrz_bin);
168  mpz_init(erg);
169  mpz_sub(erg, (mpz_ptr) a, (mpz_ptr) b);
170  return (number) erg;
171}
172
173static number nrzGetUnit (number, const coeffs r)
174{
175  return nrzInit(1, r);
176}
177
178static BOOLEAN nrzIsUnit (number a, const coeffs)
179{
180  return 0 == mpz_cmpabs_ui((mpz_ptr) a, 1);
181}
182
183static BOOLEAN nrzIsZero (number  a, const coeffs)
184{
185  return 0 == mpz_cmpabs_ui((mpz_ptr) a, 0);
186}
187
188static BOOLEAN nrzIsOne (number a, const coeffs)
189{
190  return (0 == mpz_cmp_ui((mpz_ptr) a, 1));
191}
192
193static BOOLEAN nrzIsMOne (number a, const coeffs)
194{
195  return (0 == mpz_cmp_si((mpz_ptr) a, -1));
196}
197
198static BOOLEAN nrzEqual (number a,number b, const coeffs)
199{
200  return 0 == mpz_cmp((mpz_ptr) a, (mpz_ptr) b);
201}
202
203static BOOLEAN nrzGreater (number a,number b, const coeffs)
204{
205  return 0 < mpz_cmp((mpz_ptr) a, (mpz_ptr) b);
206}
207
208static BOOLEAN nrzGreaterZero (number k, const coeffs)
209{
210  return 0 < mpz_sgn1((mpz_ptr) k);
211}
212
213static BOOLEAN nrzDivBy (number a,number b, const coeffs)
214{
215  return mpz_divisible_p((mpz_ptr) a, (mpz_ptr) b) != 0;
216}
217
218static int nrzDivComp(number a, number b, const coeffs r)
219{
220  if (nrzDivBy(a, b, r))
221  {
222    if (nrzDivBy(b, a, r)) return 2;
223    return -1;
224  }
225  if (nrzDivBy(b, a, r)) return 1;
226  return 0;
227}
228
229static number nrzDiv (number a,number b, const coeffs)
230{
231  mpz_ptr erg = (mpz_ptr) omAllocBin(gmp_nrz_bin);
232  mpz_init(erg);
233  mpz_ptr r = (mpz_ptr) omAllocBin(gmp_nrz_bin);
234  mpz_init(r);
235  mpz_tdiv_qr(erg, r, (mpz_ptr) a, (mpz_ptr) b);
236  //if (!nrzIsZero((number) r, R))
237  //{
238  //  WerrorS("Division by non divisible element.");
239  //  WerrorS("Result is without remainder.");
240  //}
241  mpz_clear(r);
242  omFreeBin(r, gmp_nrz_bin);
243  return (number) erg;
244}
245
246static number nrzExactDiv (number a,number b, const coeffs)
247{
248  mpz_ptr erg = (mpz_ptr) omAllocBin(gmp_nrz_bin);
249  mpz_init(erg);
250  mpz_tdiv_q(erg, (mpz_ptr) a, (mpz_ptr) b);
251  return (number) erg;
252}
253
254static number nrzEucNorm (number a, const coeffs )
255{
256  mpz_ptr abs = (mpz_ptr) omAllocBin(gmp_nrz_bin);
257  mpz_init(abs);
258  mpz_abs(abs, (mpz_ptr)a);
259
260  return (number) abs;
261}
262
263static number nrzSmallestQuotRem (number a, number b, number * r, const coeffs )
264{
265  mpz_ptr qq = (mpz_ptr) omAllocBin(gmp_nrz_bin);
266  mpz_init(qq);
267  mpz_ptr rr = (mpz_ptr) omAllocBin(gmp_nrz_bin);
268  mpz_init(rr);
269  int gsign = mpz_sgn((mpz_ptr) b);
270  mpz_t gg, ghalf;
271  mpz_init(gg);
272  mpz_init(ghalf);
273  mpz_abs(gg, (mpz_ptr) b);
274  mpz_fdiv_qr(qq, rr, (mpz_ptr) a, gg);
275  mpz_tdiv_q_2exp(ghalf, gg, 1);
276  if (mpz_cmp(rr, ghalf) > 0)  // r > ghalf
277    {
278      mpz_sub(rr, rr, gg);
279      mpz_add_ui(qq, qq, 1);
280    }
281  if (gsign < 0) mpz_neg(qq, qq);
282
283  mpz_clear(gg);
284  mpz_clear(ghalf);
285  if (r==NULL)
286  {
287    mpz_clear(rr);
288    omFreeBin(rr,gmp_nrz_bin);
289  }
290  else
291  {
292    *r=(number)rr;
293  }
294  return (number) qq;
295}
296
297static number nrzQuotRem (number a, number b, number * r, const coeffs )
298{
299  mpz_ptr qq = (mpz_ptr) omAllocBin(gmp_nrz_bin);
300  mpz_init(qq);
301  mpz_ptr rr = (mpz_ptr) omAllocBin(gmp_nrz_bin);
302  mpz_init(rr);
303  mpz_tdiv_qr(qq, rr, (mpz_ptr) a, (mpz_ptr) b);
304  if (r==NULL)
305  {
306    mpz_clear(rr);
307    omFreeBin(rr,gmp_nrz_bin);
308  }
309  else
310  {
311    *r=(number)rr;
312  }
313  return (number) qq;
314}
315
316static number nrzIntMod (number a,number b, const coeffs)
317{
318  mpz_ptr erg = (mpz_ptr) omAllocBin(gmp_nrz_bin);
319  mpz_init(erg);
320  mpz_ptr r = (mpz_ptr) omAllocBin(gmp_nrz_bin);
321  mpz_init(r);
322  mpz_tdiv_qr(erg, r, (mpz_ptr) a, (mpz_ptr) b);
323  mpz_clear(erg);
324  omFreeBin(erg, gmp_nrz_bin);
325  return (number) r;
326}
327
328static number  nrzInvers (number c, const coeffs r)
329{
330  if (!nrzIsUnit((number) c, r))
331  {
332    WerrorS("Non invertible element.");
333    return (number)NULL;
334  }
335  return nrzCopy(c,r);
336}
337
338static number nrzNeg (number c, const coeffs)
339{
340// nNeg inplace !!!
341  mpz_mul_si((mpz_ptr) c, (mpz_ptr) c, -1);
342  return c;
343}
344
345static number nrzMapMachineInt(number from, const coeffs /*src*/, const coeffs /*dst*/)
346{
347  mpz_ptr erg = (mpz_ptr) omAllocBin(gmp_nrz_bin);
348  mpz_init_set_ui(erg, (unsigned long) from);
349  return (number) erg;
350}
351
352static number nrzMapZp(number from, const coeffs /*src*/, const coeffs /*dst*/)
353{
354  mpz_ptr erg = (mpz_ptr) omAllocBin(gmp_nrz_bin);
355  mpz_init_set_si(erg, (long) from);
356  return (number) erg;
357}
358
359static number nrzMapQ(number from, const coeffs src, const coeffs /*dst*/)
360{
361  mpz_ptr erg = (mpz_ptr) omAllocBin(gmp_nrz_bin);
362  mpz_init(erg);
363  nlGMP(from, erg, src); // FIXME? TODO? // extern void   nlGMP(number &i, number n, const coeffs r); // to be replaced with n_MPZ(erg, from, src); // ?
364  return (number) erg;
365}
366
367static nMapFunc nrzSetMap(const coeffs src, const coeffs /*dst*/)
368{
369  /* dst = currRing */
370  /* dst = nrn */
371  if ((src->rep==n_rep_gmp)
372  && (nCoeff_is_Z(src) || nCoeff_is_Zn(src) || nCoeff_is_Ring_PtoM(src)))
373  {
374    return ndCopyMap; //nrzCopyMap;
375  }
376  if ((src->rep==n_rep_gap_gmp) /*&& nCoeff_is_Z(src)*/)
377  {
378    return ndCopyMap; //nrzCopyMap;
379  }
380  if (nCoeff_is_Ring_2toM(src))
381  {
382    return nrzMapMachineInt;
383  }
384  if (nCoeff_is_Zp(src))
385  {
386    return nrzMapZp;
387  }
388  if (getCoeffType(src)==n_Q /*nCoeff_is_Q(src) or coeffs_BIGINT*/)
389  {
390    return nrzMapQ;
391  }
392  return NULL;      // default
393}
394
395/*
396 * set the exponent (allocate and init tables) (TODO)
397 */
398
399void nrzSetExp(int, coeffs)
400{
401}
402
403void nrzInitExp(int, coeffs)
404{
405}
406
407#ifdef LDEBUG
408static BOOLEAN nrzDBTest (number, const char *, const int, const coeffs)
409{
410  return TRUE;//TODO
411}
412#endif
413
414void nrzWrite (number a, const coeffs)
415{
416  char *s,*z;
417  if (a==NULL)
418  {
419    StringAppendS("o");
420  }
421  else
422  {
423    int l=mpz_sizeinbase((mpz_ptr) a, 10) + 2;
424    s=(char*)omAlloc(l);
425    z=mpz_get_str(s,10,(mpz_ptr) a);
426    StringAppendS(z);
427    omFreeSize((ADDRESS)s,l);
428  }
429}
430
431/*2
432* extracts a long integer from s, returns the rest    (COPY FROM longrat0.cc)
433*/
434static const char * nlEatLongC(char *s, mpz_ptr i)
435{
436  const char * start=s;
437
438  if (*s<'0' || *s>'9')
439  {
440    mpz_set_ui(i,1);
441    return s;
442  }
443  while (*s >= '0' && *s <= '9') s++;
444  if (*s=='\0')
445  {
446    mpz_set_str(i,start,10);
447  }
448  else
449  {
450    char c=*s;
451    *s='\0';
452    mpz_set_str(i,start,10);
453    *s=c;
454  }
455  return s;
456}
457
458
459static CanonicalForm nrzConvSingNFactoryN(number n, BOOLEAN setChar, const coeffs /*r*/)
460{
461  if (setChar) setCharacteristic( 0 );
462
463  CanonicalForm term;
464  mpz_t num;
465  mpz_init_set(num, *((mpz_t*)n));
466  term = make_cf(num);
467  return term;
468}
469
470static number nrzConvFactoryNSingN(const CanonicalForm n, const coeffs r)
471{
472  if (n.isImm())
473    return nrzInit(n.intval(),r);
474  else
475  {
476    mpz_ptr m = (mpz_ptr) omAllocBin(gmp_nrz_bin);
477    gmp_numerator(n,m);
478    if (!n.den().isOne())
479    {
480      WarnS("denominator is not 1 in factory");
481    }
482    return (number) m;
483  }
484}
485
486static const char * nrzRead (const char *s, number *a, const coeffs)
487{
488  mpz_ptr z = (mpz_ptr) omAllocBin(gmp_nrz_bin);
489  {
490    mpz_init(z);
491    s = nlEatLongC((char *) s, z);
492  }
493  *a = (number) z;
494  return s;
495}
496
497static coeffs nrzQuot1(number c, const coeffs r)
498{
499    long ch = r->cfInt(c, r);
500    mpz_t dummy;
501    mpz_init_set_ui(dummy, ch);
502    ZnmInfo info;
503    info.base = dummy;
504    info.exp = (unsigned long) 1;
505    coeffs rr = nInitChar(n_Zn, (void*)&info);
506    mpz_clear(dummy);
507    return(rr);
508}
509
510static number nrzInitMPZ(mpz_t m, const coeffs)
511{
512  mpz_ptr z = (mpz_ptr) omAllocBin(gmp_nrz_bin);
513  mpz_init_set(z, m);
514  return (number)z;
515}
516
517static number nrzFarey(number r, number N, const coeffs R)
518{
519  number a0 = nrzCopy(N, R);
520  number b0 = nrzInit(0, R);
521  number a1 = nrzCopy(r, R);
522  number b1 = nrzInit(1, R);
523  number two = nrzInit(2, R);
524#if 0
525  PrintS("Farey start with ");
526  n_Print(r, R);
527  PrintS(" mod ");
528  n_Print(N, R);
529  PrintLn();
530#endif
531  while (1)
532  {
533    number as = nrzMult(a1, a1, R);
534    n_InpMult(as, two, R);
535    if (nrzGreater(N, as, R))
536    {
537      nrzDelete(&as, R);
538      break;
539    }
540    nrzDelete(&as, R);
541    number q = nrzDiv(a0, a1, R);
542    number t = nrzMult(a1, q, R),
543           s = nrzSub(a0, t, R);
544    nrzDelete(&a0, R);
545    a0 = a1;
546    a1 = s;
547    nrzDelete(&t, R);
548
549    t = nrzMult(b1, q, R);
550    s = nrzSub(b0, t, R);
551    nrzDelete(&b0, R);
552    b0 = b1;
553    b1 = s;
554    nrzDelete(&t, R);
555    nrzDelete(&q, R);
556  }
557  number as = nrzMult(b1, b1, R);
558  n_InpMult(as, two, R);
559  nrzDelete(&two, R);
560  if (nrzGreater(as, N, R))
561  {
562    nrzDelete(&a0, R);
563    nrzDelete(&a1, R);
564    nrzDelete(&b0, R);
565    nrzDelete(&b1, R);
566    nrzDelete(&as, R);
567    return NULL;
568  }
569  nrzDelete(&as, R);
570  nrzDelete(&a0, R);
571  nrzDelete(&b0, R);
572
573  number a, b, ab;
574  coeffs Q = nInitChar(n_Q, 0);
575  nMapFunc f = n_SetMap(R, Q);
576  a = f(a1, R, Q);
577  b = f(b1, R, Q);
578  ab = n_Div(a, b, Q);
579  n_Delete(&a, Q);
580  n_Delete(&b, Q);
581  nKillChar(Q);
582
583  nrzDelete(&a1, R);
584  nrzDelete(&b1, R);
585  return ab;
586}
587
588void nrzWriteFd(number n, const ssiInfo* d, const coeffs)
589{
590  mpz_out_str (d->f_write,SSI_BASE, (mpz_ptr)n);
591  fputc(' ',d->f_write);
592}
593
594number nrzReadFd(const ssiInfo *d, const coeffs)
595{
596  mpz_ptr erg = (mpz_ptr) omAllocBin(gmp_nrz_bin);
597  mpz_init(erg);
598  s_readmpz_base(d->f_read,erg,SSI_BASE);
599  return (number)erg;
600}
601
602BOOLEAN nrzInitChar(coeffs r,  void *)
603{
604  assume( getCoeffType(r) == n_Z );
605
606  r->is_field=FALSE;
607  r->is_domain=TRUE;
608  r->rep=n_rep_gmp;
609
610  //r->nCoeffIsEqual = ndCoeffIsEqual;
611  r->cfCoeffString = nrzCoeffString;
612  r->cfCoeffName = nrzCoeffName;
613  r->cfCoeffWrite = nrzCoeffWrite;
614  //r->cfKillChar = ndKillChar;
615  r->cfMult  = nrzMult;
616  r->cfSub   = nrzSub;
617  r->cfAdd   = nrzAdd;
618  r->cfDiv   = nrzDiv;
619  r->cfIntMod= nrzIntMod;
620  r->cfExactDiv= nrzExactDiv;
621  r->cfInit = nrzInit;
622  r->cfInitMPZ = nrzInitMPZ;
623  r->cfSize  = nrzSize;
624  r->cfInt  = nrzInt;
625  r->cfDivComp = nrzDivComp;
626  r->cfIsUnit = nrzIsUnit;
627  r->cfGetUnit = nrzGetUnit;
628  r->cfExtGcd = nrzExtGcd;
629  r->cfXExtGcd = nrzXExtGcd;
630  r->cfDivBy = nrzDivBy;
631  r->cfEucNorm = nrzEucNorm;
632  r->cfQuotRem = nrzSmallestQuotRem;
633  r->cfInpNeg   = nrzNeg;
634  r->cfInvers= nrzInvers;
635  r->cfCopy  = nrzCopy;
636  r->cfWriteLong = nrzWrite;
637  r->cfRead = nrzRead;
638  r->cfGreater = nrzGreater;
639  r->cfEqual = nrzEqual;
640  r->cfIsZero = nrzIsZero;
641  r->cfIsOne = nrzIsOne;
642  r->cfIsMOne = nrzIsMOne;
643  r->cfGreaterZero = nrzGreaterZero;
644  r->cfPower = nrzPower;
645  r->cfGcd  = nrzGcd;
646  r->cfLcm  = nrzLcm;
647  r->cfDelete= nrzDelete;
648  r->cfSetMap = nrzSetMap;
649  r->cfQuot1 = nrzQuot1;
650  r->convSingNFactoryN=nrzConvSingNFactoryN;
651  r->convFactoryNSingN=nrzConvFactoryNSingN;
652  r->cfChineseRemainder=nlChineseRemainderSym;
653  r->cfFarey=nrzFarey;
654  r->cfWriteFd=nrzWriteFd;
655  r->cfReadFd=nrzReadFd;
656  // debug stuff
657
658#ifdef LDEBUG
659  r->cfDBTest=nrzDBTest;
660#endif
661
662  r->ch = 0;
663  r->has_simple_Alloc=FALSE;
664  r->has_simple_Inverse=FALSE;
665  return FALSE;
666}
667#endif
668#endif
Note: See TracBrowser for help on using the repository browser.