source: git/libpolys/coeffs/modulop.cc @ d8b0997

jengelh-datetimespielwiese
Last change on this file since d8b0997 was d8b0997, checked in by Hans Schoenemann <hannes@…>, 10 years ago
fix: npInit for large ints
  • Property mode set to 100644
File size: 15.7 KB
Line 
1/****************************************
2*  Computer Algebra System SINGULAR     *
3****************************************/
4/*
5* ABSTRACT: numbers modulo p (<=32003)
6*/
7
8#include "config.h"
9#include <misc/auxiliary.h>
10
11#ifdef HAVE_FACTORY
12#include <factory/factory.h>
13#endif
14
15#include <string.h>
16#include <omalloc/omalloc.h>
17#include <coeffs/coeffs.h>
18#include <reporter/reporter.h>
19#include <coeffs/numbers.h>
20#include <coeffs/longrat.h>
21#include <coeffs/mpr_complex.h>
22#include <misc/mylimits.h>
23#include <coeffs/modulop.h>
24
25// int npGen=0;
26
27/// Our Type!
28static const n_coeffType ID = n_Zp;
29
30#ifdef NV_OPS
31#pragma GCC diagnostic ignored "-Wlong-long"
32static inline number nvMultM(number a, number b, const coeffs r)
33{
34  assume( getCoeffType(r) == ID );
35 
36#if SIZEOF_LONG == 4
37#define ULONG64 (unsigned long long)(unsigned long)
38#else
39#define ULONG64 (unsigned long)
40#endif
41  return (number) 
42      (unsigned long)((ULONG64 a)*(ULONG64 b) % (ULONG64 r->ch));
43}
44number  nvMult        (number a, number b, const coeffs r);
45number  nvDiv         (number a, number b, const coeffs r);
46number  nvInvers      (number c, const coeffs r);
47void    nvPower       (number a, int i, number * result, const coeffs r);
48#endif
49
50
51
52
53BOOLEAN npGreaterZero (number k, const coeffs r)
54{
55  assume( n_Test(k, r) );
56 
57  int h = (int)((long) k);
58  return ((int)h !=0) && (h <= (r->ch>>1));
59}
60
61//unsigned long npMultMod(unsigned long a, unsigned long b, int npPrimeM)
62//{
63//  unsigned long c = a*b;
64//  c = c % npPrimeM;
65//  assume(c == (unsigned long) npMultM((number) a, (number) b, npPrimeM));
66//  return c;
67//}
68
69number npMult (number a,number b, const coeffs r)
70{
71  assume( n_Test(a, r) );
72  assume( n_Test(b, r) );
73
74  if (((long)a == 0) || ((long)b == 0))
75    return (number)0;
76  number c = npMultM(a,b, r);
77  assume( n_Test(c, r) );
78  return c; 
79}
80
81/*2
82* create a number from int
83*/
84number npInit (long i, const coeffs r)
85{
86  long ii=i % (long)r->ch;
87  if (ii <  0L)                         ii += (long)r->ch;
88
89  number c = (number)ii;
90  assume( n_Test(c, r) );
91  return c;
92 
93}
94
95
96/*2
97 * convert a number to an int in (-p/2 .. p/2]
98 */
99int npInt(number &n, const coeffs r)
100{
101  assume( n_Test(n, r) );
102 
103  if ((long)n > (((long)r->ch) >>1)) return (int)((long)n -((long)r->ch));
104  else                               return (int)((long)n);
105}
106
107number npAdd (number a, number b, const coeffs r)
108{
109  assume( n_Test(a, r) );
110  assume( n_Test(b, r) );
111 
112  number c = npAddM(a,b, r);
113
114  assume( n_Test(c, r) );
115 
116  return c;
117}
118
119number npSub (number a, number b, const coeffs r)
120{
121  assume( n_Test(a, r) );
122  assume( n_Test(b, r) );
123 
124  number c = npSubM(a,b,r);
125
126  assume( n_Test(c, r) );
127
128  return c;
129}
130
131BOOLEAN npIsZero (number  a, const coeffs r)
132{
133  assume( n_Test(a, r) );
134 
135  return 0 == (long)a;
136}
137
138BOOLEAN npIsOne (number a, const coeffs r)
139{
140  assume( n_Test(a, r) );
141 
142  return 1 == (long)a;
143}
144
145BOOLEAN npIsMOne (number a, const coeffs r)
146{
147  assume( n_Test(a, r) );
148 
149  return ((r->npPminus1M == (long)a)&&((long)1!=(long)a));
150}
151
152#ifdef HAVE_DIV_MOD
153
154#ifdef USE_NTL_XGCD
155
156//ifdef HAVE_NTL // in ntl.a
157//extern void XGCD(long& d, long& s, long& t, long a, long b);
158#include <NTL/ZZ.h>
159#ifdef NTL_CLIENT
160NTL_CLIENT
161#endif
162
163#endif
164
165long InvMod(long a, const coeffs R)
166{
167   long d, s, t;
168
169#ifdef USE_NTL_XGCD
170   XGCD(d, s, t, a, R->ch);
171   assume (d == 1);
172#else
173   long  u, v, u0, v0, u1, v1, u2, v2, q, r;
174
175   assume(a>0);
176   u1=1; u2=0;
177   u = a; v = R->ch;
178
179   while (v != 0)
180   {
181      q = u / v;
182      r = u % v;
183      u = v;
184      v = r;
185      u0 = u2;
186      u2 = u1 - q*u2;
187      u1 = u0;
188   }
189
190   assume(u==1);
191   s = u1;
192#endif
193   if (s < 0)
194      return s + R->ch;
195   else
196      return s;
197}
198#endif
199
200inline number npInversM (number c, const coeffs r)
201{
202  assume( n_Test(c, r) );
203#ifndef HAVE_DIV_MOD
204  number d = (number)(long)r->npExpTable[r->npPminus1M - r->npLogTable[(long)c]];
205#else
206  long inv=(long)r->npInvTable[(long)c];
207  if (inv==0)
208  {
209    inv=InvMod((long)c,r);
210    r->npInvTable[(long)c]=inv;
211  }
212  number d = (number)inv;
213#endif
214  assume( n_Test(d, r) );
215  return d;
216
217}
218
219number npDiv (number a,number b, const coeffs r)
220{
221  assume( n_Test(a, r) );
222  assume( n_Test(b, r) );
223
224//#ifdef NV_OPS
225//  if (r->ch>NV_MAX_PRIME)
226//    return nvDiv(a,b);
227//#endif
228  if ((long)a==0)
229    return (number)0;
230  number d;
231 
232#ifndef HAVE_DIV_MOD
233  if ((long)b==0)
234  {
235    WerrorS(nDivBy0);
236    return (number)0;
237  }
238 
239  int s = r->npLogTable[(long)a] - r->npLogTable[(long)b];
240  if (s < 0)
241    s += r->npPminus1M;
242  d = (number)(long)r->npExpTable[s];
243#else
244  number inv=npInversM(b,r);
245  d = npMultM(a,inv,r);
246#endif
247
248  assume( n_Test(d, r) );
249  return d;
250 
251}
252number  npInvers (number c, const coeffs r)
253{
254  assume( n_Test(c, r) );
255 
256  if ((long)c==0)
257  {
258    WerrorS("1/0");
259    return (number)0;
260  }
261  number d = npInversM(c,r);
262 
263  assume( n_Test(d, r) );
264  return d;
265 
266}
267
268number npNeg (number c, const coeffs r)
269{
270  assume( n_Test(c, r) );
271 
272  if ((long)c==0) return c;
273
274#if 0 
275  number d = npNegM(c,r); 
276  assume( n_Test(d, r) );
277  return d;
278#else
279  c = npNegM(c,r); 
280  assume( n_Test(c, r) );
281  return c;
282#endif 
283}
284
285BOOLEAN npGreater (number a,number b, const coeffs r)
286{
287  assume( n_Test(a, r) );
288  assume( n_Test(b, r) );
289
290  //return (long)a != (long)b;
291  return (long)a > (long)b;
292}
293
294BOOLEAN npEqual (number a,number b, const coeffs r)
295{
296  assume( n_Test(a, r) );
297  assume( n_Test(b, r) );
298 
299//  return (long)a == (long)b;
300 
301  return npEqualM(a,b,r);
302}
303
304void npWrite (number &a, const coeffs r)
305{
306  assume( n_Test(a, r) );
307 
308  if ((long)a>(((long)r->ch) >>1)) StringAppend("-%d",(int)(((long)r->ch)-((long)a)));
309  else                             StringAppend("%d",(int)((long)a));
310}
311
312void npPower (number a, int i, number * result, const coeffs r)
313{
314  assume( n_Test(a, r) );
315
316  if (i==0)
317  {
318    //npInit(1,result);
319    *(long *)result = 1;
320  }
321  else if (i==1)
322  {
323    *result = a;
324  }
325  else
326  {
327    npPower(a,i-1,result,r);
328    *result = npMultM(a,*result,r);
329  }
330}
331
332static const char* npEati(const char *s, int *i, const coeffs r)
333{
334
335  if (((*s) >= '0') && ((*s) <= '9'))
336  {
337    unsigned long ii=0L;
338    do
339    {
340      ii *= 10;
341      ii += *s++ - '0';
342      if (ii >= (MAX_INT_VAL / 10)) ii = ii % r->ch;
343    }
344    while (((*s) >= '0') && ((*s) <= '9'));
345    if (ii >= r->ch) ii = ii % r->ch;
346    *i=(int)ii;
347  }
348  else (*i) = 1;
349  return s;
350}
351
352const char * npRead (const char *s, number *a, const coeffs r)
353{
354  int z;
355  int n=1;
356
357  s = npEati(s, &z, r);
358  if ((*s) == '/')
359  {
360    s++;
361    s = npEati(s, &n, r);
362  }
363  if (n == 1)
364    *a = (number)z;
365  else
366  {
367    if ((z==0)&&(n==0)) WerrorS(nDivBy0);
368    else
369    {
370#ifdef NV_OPS
371      if (r->ch>NV_MAX_PRIME)
372        *a = nvDiv((number)z,(number)n,r);
373      else
374#endif
375        *a = npDiv((number)z,(number)n,r);
376    }
377  }
378  assume( n_Test(*a, r) );
379  return s;
380}
381
382/*2
383* set the charcteristic (allocate and init tables)
384*/
385
386void npKillChar(coeffs r)
387{
388  #ifdef HAVE_DIV_MOD
389  if (r->npInvTable!=NULL)
390  omFreeSize( (void *)r->npInvTable, r->ch*sizeof(unsigned short) );
391  r->npInvTable=NULL;
392  #else
393  if (r->npExpTable!=NULL)
394  {
395    omFreeSize( (void *)r->npExpTable, r->ch*sizeof(unsigned short) );
396    omFreeSize( (void *)r->npLogTable, r->ch*sizeof(unsigned short) );
397    r->npExpTable=NULL; r->npLogTable=NULL;
398  }
399  #endif
400}
401
402static BOOLEAN npCoeffsEqual(const coeffs r, n_coeffType n, void * parameter)
403{
404  /* test, if r is an instance of nInitCoeffs(n,parameter) */
405  return (n==n_Zp) && (r->ch==(int)(long)parameter);
406}
407#ifdef HAVE_FACTORY
408CanonicalForm npConvSingNFactoryN( number n, BOOLEAN setChar, const coeffs r )
409{
410  if (setChar) setCharacteristic( r->ch );
411  CanonicalForm term(npInt( n,r ));
412  return term;
413}
414
415number npConvFactoryNSingN( const CanonicalForm n, const coeffs r)
416{
417  if (n.isImm())
418  {
419    return npInit(n.intval(),r);
420  }
421  else
422  {
423    assume(0);
424    return NULL;
425  }
426}
427#endif
428
429
430BOOLEAN npInitChar(coeffs r, void* p)
431{
432  assume( getCoeffType(r) == ID );
433  const int c = (int) (long) p;
434
435  assume( c > 0 );
436 
437  int i, w;
438
439  r->ch = c;
440  r->npPminus1M = c /*r->ch*/ - 1;
441
442  //r->cfInitChar=npInitChar;
443  r->cfKillChar=npKillChar;
444  r->nCoeffIsEqual=npCoeffsEqual;
445
446  r->cfMult  = npMult;
447  r->cfSub   = npSub;
448  r->cfAdd   = npAdd;
449  r->cfDiv   = npDiv;
450  r->cfIntDiv= npDiv;
451  //r->cfIntMod= ndIntMod;
452  r->cfExactDiv= npDiv;
453  r->cfInit = npInit;
454  //r->cfSize  = ndSize;
455  r->cfInt  = npInt;
456  #ifdef HAVE_RINGS
457  //r->cfDivComp = NULL; // only for ring stuff
458  //r->cfIsUnit = NULL; // only for ring stuff
459  //r->cfGetUnit = NULL; // only for ring stuff
460  //r->cfExtGcd = NULL; // only for ring stuff
461  // r->cfDivBy = NULL; // only for ring stuff
462  #endif
463  r->cfNeg   = npNeg;
464  r->cfInvers= npInvers;
465  //r->cfCopy  = ndCopy;
466  //r->cfRePart = ndCopy;
467  //r->cfImPart = ndReturn0;
468  r->cfWriteLong = npWrite;
469  r->cfRead = npRead;
470  //r->cfNormalize=ndNormalize;
471  r->cfGreater = npGreater;
472  r->cfEqual = npEqual;
473  r->cfIsZero = npIsZero;
474  r->cfIsOne = npIsOne;
475  r->cfIsMOne = npIsMOne;
476  r->cfGreaterZero = npGreaterZero;
477  r->cfPower = npPower;
478  r->cfGetDenom = ndGetDenom;
479  r->cfGetNumerator = ndGetNumerator;
480  //r->cfGcd  = ndGcd;
481  //r->cfLcm  = ndGcd;
482  //r->cfDelete= ndDelete;
483  r->cfSetMap = npSetMap;
484  //r->cfName = ndName;
485  r->cfInpMult=ndInpMult;
486  r->cfInit_bigint= nlModP; // npMap0;
487#ifdef NV_OPS
488  if (c>NV_MAX_PRIME)
489  {
490    r->cfMult  = nvMult;
491    r->cfDiv   = nvDiv;
492    r->cfExactDiv= nvDiv;
493    r->cfInvers= nvInvers;
494    r->cfPower= nvPower;
495  }
496#endif
497  r->cfCoeffWrite=npCoeffWrite;
498#ifdef LDEBUG
499  // debug stuff
500  r->cfDBTest=npDBTest;
501#endif
502
503#ifdef HAVE_FACTORY
504  r->convSingNFactoryN=npConvSingNFactoryN;
505  r->convFactoryNSingN=npConvFactoryNSingN;
506#endif
507 
508  // the variables:
509  r->nNULL = (number)0;
510  r->type = n_Zp;
511  r->ch = c;
512  r->has_simple_Alloc=TRUE;
513  r->has_simple_Inverse=TRUE;
514
515  // the tables
516#ifdef NV_OPS
517  if (r->ch <=NV_MAX_PRIME)
518#endif
519  {
520#if !defined(HAVE_DIV_MOD) || !defined(HAVE_MULT_MOD)
521    r->npExpTable=(unsigned short *)omAlloc( r->ch*sizeof(unsigned short) );
522    r->npLogTable=(unsigned short *)omAlloc( r->ch*sizeof(unsigned short) );
523    r->npExpTable[0] = 1;
524    r->npLogTable[0] = 0;
525    if (r->ch > 2)
526    {
527      w = 1;
528      loop
529      {
530        r->npLogTable[1] = 0;
531        w++;
532        i = 0;
533        loop
534        {
535          i++;
536          r->npExpTable[i] =(int)(((long)w * (long)r->npExpTable[i-1])
537                               % r->ch);
538          r->npLogTable[r->npExpTable[i]] = i;
539          if (/*(i == r->ch - 1 ) ||*/ (r->npExpTable[i] == 1))
540            break;
541        }
542        if (i == r->ch - 1)
543          break;
544      }
545    }
546    else
547    {
548      r->npExpTable[1] = 1;
549      r->npLogTable[1] = 0;
550    }
551#endif
552#ifdef HAVE_DIV_MOD
553    r->npInvTable=(unsigned short*)omAlloc0( r->ch*sizeof(unsigned short) );
554#endif
555  }
556  return FALSE;
557}
558
559#ifdef LDEBUG
560BOOLEAN npDBTest (number a, const char *f, const int l, const coeffs r)
561{
562  if (((long)a<0) || ((long)a>r->ch))
563  {
564    Print("wrong mod p number %ld at %s,%d\n",(long)a,f,l);
565    return FALSE;
566  }
567  return TRUE;
568}
569#endif
570
571number npMapP(number from, const coeffs src, const coeffs dst_r)
572{
573  long i = (long)from;
574  if (i>src->ch/2)
575  {
576    i-=src->ch;
577    while (i < 0) i+=dst_r->ch;
578  }
579  i%=dst_r->ch;
580  return (number)i;
581}
582
583static number npMapLongR(number from, const coeffs /*src*/, const coeffs dst_r)
584{
585  gmp_float *ff=(gmp_float*)from;
586  mpf_t *f=ff->_mpfp();
587  number res;
588  mpz_ptr dest,ndest;
589  int size,i;
590  int e,al,bl;
591  long iz;
592  mp_ptr qp,dd,nn;
593
594  size = (*f)[0]._mp_size;
595  if (size == 0)
596    return npInit(0,dst_r);
597  if(size<0)
598    size = -size;
599
600  qp = (*f)[0]._mp_d;
601  while(qp[0]==0)
602  {
603    qp++;
604    size--;
605  }
606
607  if(dst_r->ch>2)
608    e=(*f)[0]._mp_exp-size;
609  else
610    e=0;
611  res = ALLOC_RNUMBER();
612#if defined(LDEBUG)
613  res->debug=123456;
614#endif
615  dest = res->z;
616
617  int in=0;
618  if (e<0)
619  {
620    al = dest->_mp_size = size;
621    if (al<2) al = 2;
622    dd = (mp_ptr)omAlloc(sizeof(mp_limb_t)*al);
623    for (i=0;i<size;i++) dd[i] = qp[i];
624    bl = 1-e;
625    nn = (mp_ptr)omAlloc(sizeof(mp_limb_t)*bl);
626    nn[bl-1] = 1;
627    for (i=bl-2;i>=0;i--) nn[i] = 0;
628    ndest = res->n;
629    ndest->_mp_d = nn;
630    ndest->_mp_alloc = ndest->_mp_size = bl;
631    res->s = 0;
632    in=mpz_fdiv_ui(ndest,dst_r->ch);
633    mpz_clear(ndest);
634  }
635  else
636  {
637    al = dest->_mp_size = size+e;
638    if (al<2) al = 2;
639    dd = (mp_ptr)omAlloc(sizeof(mp_limb_t)*al);
640    for (i=0;i<size;i++) dd[i+e] = qp[i];
641    for (i=0;i<e;i++) dd[i] = 0;
642    res->s = 3;
643  }
644
645  dest->_mp_d = dd;
646  dest->_mp_alloc = al;
647  iz=mpz_fdiv_ui(dest,dst_r->ch);
648  mpz_clear(dest);
649  if(res->s==0)
650    iz=(long)npDiv((number)iz,(number)in,dst_r);
651  FREE_RNUMBER(res); // Q!?
652  return (number)iz;
653}
654
655#ifdef HAVE_RINGS
656/*2
657* convert from a GMP integer
658*/
659number npMapGMP(number from, const coeffs /*src*/, const coeffs dst)
660{
661  int_number erg = (int_number) omAlloc(sizeof(mpz_t)); // evtl. spaeter mit bin
662  mpz_init(erg);
663
664  mpz_mod_ui(erg, (int_number) from, dst->ch);
665  number r = (number) mpz_get_si(erg);
666
667  mpz_clear(erg);
668  omFree((void *) erg);
669  return (number) r;
670}
671
672/*2
673* convert from an machine long
674*/
675number npMapMachineInt(number from, const coeffs /*src*/,const coeffs dst)
676{
677  long i = (long) (((unsigned long) from) % dst->ch);
678  return (number) i;
679}
680#endif
681
682#ifdef HAVE_FACTORY
683number npMapCanonicalForm (number a, const coeffs /*src*/, const coeffs dst)
684{
685  setCharacteristic (dst ->ch);
686  CanonicalForm f= CanonicalForm ((InternalCF*)(a));
687  return (number) (f.intval());
688}
689#endif
690
691nMapFunc npSetMap(const coeffs src, const coeffs dst)
692{
693#ifdef HAVE_RINGS
694  if (nCoeff_is_Ring_2toM(src))
695  {
696    return npMapMachineInt;
697  }
698  if (nCoeff_is_Ring_Z(src) || nCoeff_is_Ring_PtoM(src) || nCoeff_is_Ring_ModN(src))
699  {
700    return npMapGMP;
701  }
702#endif
703  if (nCoeff_is_Q(src))
704  {
705    return nlModP; // npMap0;
706  }
707  if ( nCoeff_is_Zp(src) )
708  {
709    if (n_GetChar(src) == n_GetChar(dst))
710    {
711      return ndCopyMap;
712    }
713    else
714    {
715      return npMapP;
716    }
717  }
718  if (nCoeff_is_long_R(src))
719  {
720    return npMapLongR;
721  }
722#ifdef HAVE_FACTORY
723  if (nCoeff_is_CF (src))
724  {
725    return npMapCanonicalForm;
726  }
727#endif
728  return NULL;      /* default */
729}
730
731// -----------------------------------------------------------
732//  operation for very large primes (32003< p < 2^31-1)
733// ----------------------------------------------------------
734#ifdef NV_OPS
735
736number nvMult (number a,number b, const coeffs r)
737{
738  //if (((long)a == 0) || ((long)b == 0))
739  //  return (number)0;
740  //else
741    return nvMultM(a,b,r);
742}
743
744void   nvInpMult(number &a, number b, const coeffs r)
745{
746  number n=nvMultM(a,b,r);
747  a=n;
748}
749
750
751inline long nvInvMod(long a, const coeffs R)
752{
753#ifdef HAVE_DIV_MOD
754  return InvMod(a, R);
755#else
756/// TODO: use "long InvMod(long a, const coeffs R)"?!
757 
758   long  s;
759
760   long  u, u0, u1, u2, q, r; // v0, v1, v2,
761
762   u1=1; // v1=0;
763   u2=0; // v2=1;
764   u = a;
765
766   long v = R->ch;
767
768   while (v != 0)
769   {
770      q = u / v;
771      r = u % v;
772      u = v;
773      v = r;
774      u0 = u2;
775//      v0 = v2;
776      u2 = u1 - q*u2;
777//      v2 = v1 - q*v2;
778      u1 = u0;
779//      v1 = v0;
780   }
781
782   s = u1;
783   //t = v1;
784   if (s < 0)
785      return s + R->ch;
786   else
787     return s;
788#endif
789}
790
791inline number nvInversM (number c, const coeffs r)
792{
793  long inv=nvInvMod((long)c,r);
794  return (number)inv;
795}
796
797number nvDiv (number a,number b, const coeffs r)
798{
799  if ((long)a==0)
800    return (number)0;
801  else if ((long)b==0)
802  {
803    WerrorS(nDivBy0);
804    return (number)0;
805  }
806  else
807  {
808    number inv=nvInversM(b,r);
809    return nvMultM(a,inv,r);
810  }
811}
812number  nvInvers (number c, const coeffs r)
813{
814  if ((long)c==0)
815  {
816    WerrorS(nDivBy0);
817    return (number)0;
818  }
819  return nvInversM(c,r);
820}
821void nvPower (number a, int i, number * result, const coeffs r)
822{
823  if (i==0)
824  {
825    //npInit(1,result);
826    *(long *)result = 1;
827  }
828  else if (i==1)
829  {
830    *result = a;
831  }
832  else
833  {
834    nvPower(a,i-1,result,r);
835    *result = nvMultM(a,*result,r);
836  }
837}
838#endif
839
840void    npCoeffWrite  (const coeffs r, BOOLEAN /*details*/)
841{
842  Print("//   characteristic : %d\n",r->ch);
843}
844
Note: See TracBrowser for help on using the repository browser.