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

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