source: git/libpolys/coeffs/modulop.cc @ 0461f0

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