source: git/libpolys/coeffs/modulop.cc @ 85bcd6

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