source: git/libpolys/coeffs/rintegers2.cc @ 96ce32

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