source: git/libpolys/coeffs/modulop.h

spielwiese
Last change on this file was 00c93d, checked in by Hans Schoenemann <hannes@…>, 2 years ago
adapt (better) to diffent cpu types
  • Property mode set to 100644
File size: 5.9 KB
RevLine 
[35aab3]1#ifndef MODULOP_H
2#define MODULOP_H
3/****************************************
4*  Computer Algebra System SINGULAR     *
5****************************************/
6/*
[a1299c]7* ABSTRACT: numbers modulo p (<=32749)
[35aab3]8*/
[f5f5c9]9#include "misc/auxiliary.h"
[00c93d]10#include "factory/factory.h"
11
12#ifdef HAVE_NTL
13  #include <NTL/config.h>
14  #ifdef NTL_AVOID_BRANCHING
15  #undef HAVE_GENERIC_ADD
16  #endif
17#endif
[35aab3]18
[2a4a23]19
[35aab3]20// define if a*b is with mod instead of tables
[2a4a23]21//#define HAVE_GENERIC_MULT
22// define if 1/b is from  tables
23//#define HAVE_INVTABLE
[00c93d]24// define if an if should be used(NTL does not define NTL_AVOID_BRANCHING)
[35aab3]25//#define HAVE_GENERIC_ADD
[c2bd74]26
[2a4a23]27//#undef HAVE_GENERIC_ADD
28//#undef HAVE_GENERIC_MULT
29//#undef HAVE_INVTABLE
30
31//#define HAVE_GENERIC_ADD 1
32//#define HAVE_GENERIC_MULT 1
33//#define HAVE_INVTABLE 1
34
[a1299c]35// enable large primes (32749 < p < 2^31-)
[c2bd74]36#define NV_OPS
[a1299c]37#define NV_MAX_PRIME 32749
[097353]38#define FACTORY_MAX_PRIME 536870909
[c2bd74]39
[2206753]40struct n_Procs_s; typedef struct  n_Procs_s  *coeffs;
41struct snumber; typedef struct snumber *   number;
[35aab3]42
[1cce47]43BOOLEAN npInitChar(coeffs r, void* p);
[8c484e]44
[e77676]45// inline number npMultM(number a, number b, int npPrimeM)
46// // return (a*b)%n
47// {
48//    double ab;
49//    long q, res;
[601105]50//
[e77676]51//    ab = ((double) ((int)a)) * ((double) ((int)b));
52//    q  = (long) (ab/((double) npPrimeM));  // q could be off by (+/-) 1
53//    res = (long) (ab - ((double) q)*((double) npPrimeM));
54//    res += (res >> 31) & npPrimeM;
55//    res -= npPrimeM;
56//    res += (res >> 31) & npPrimeM;
57//    return (number)res;
58// }
[2a4a23]59#ifdef HAVE_GENERIC_MULT
[7d90aa]60static inline number npMultM(number a, number b, const coeffs r)
[35aab3]61{
[601105]62  return (number)
[e77676]63    ((((unsigned long) a)*((unsigned long) b)) % ((unsigned long) r->ch));
[35aab3]64}
[bb865e]65static inline void npInpMultM(number &a, number b, const coeffs r)
66{
67  a=(number)
68    ((((unsigned long) a)*((unsigned long) b)) % ((unsigned long) r->ch));
69}
[35aab3]70#else
[7d90aa]71static inline number npMultM(number a, number b, const coeffs r)
[35aab3]72{
[7d90aa]73  long x = (long)r->npLogTable[(long)a]+ r->npLogTable[(long)b];
[55b747]74  #ifdef HAVE_GENERIC_ADD
[98a94b]75    if (x>=r->npPminus1M) x-=r->npPminus1M;
[55b747]76  #else
77    x-=r->npPminus1M;
78    #if SIZEOF_LONG == 8
[98a94b]79      x += (x >> 63) & r->npPminus1M;
[55b747]80    #else
[98a94b]81      x += (x >> 31) & r->npPminus1M;
[55b747]82    #endif
83  #endif
84  return (number)(long)r->npExpTable[x];
[35aab3]85}
[bb865e]86static inline void npInpMultM(number &a, number b, const coeffs r)
87{
88  long x = (long)r->npLogTable[(long)a]+ r->npLogTable[(long)b];
[55b747]89  #ifdef HAVE_GENERIC_ADD
[98a94b]90    if (x>=r->npPminus1M) x-=r->npPminus1M;
[55b747]91  #else
92    x-=r->npPminus1M;
93    #if SIZEOF_LONG == 8
[98a94b]94      x += (x >> 63) & r->npPminus1M;
[55b747]95    #else
[98a94b]96      x += (x >> 31) & r->npPminus1M;
[55b747]97    #endif
98  #endif
99  a=(number)(long)r->npExpTable[x];
[bb865e]100}
[35aab3]101#endif
102
103#if 0
104inline number npAddAsm(number a, number b, int m)
105{
106  number r;
107    asm ("addl %2, %1; cmpl %3, %1; jb 0f; subl %3, %1; 0:"
108         : "=&r" (r)
109         : "%0" (a), "g" (b), "g" (m)
110         : "cc");
111  return r;
112}
113inline number npSubAsm(number a, number b, int m)
114{
115  number r;
116  asm ("subl %2, %1; jnc 0f; addl %3, %1; 0:"
117        : "=&r" (r)
118        : "%0" (a), "g" (b), "g" (m)
119        : "cc");
120  return r;
121}
122#endif
123#ifdef HAVE_GENERIC_ADD
[7d90aa]124static inline number npAddM(number a, number b, const coeffs r)
[35aab3]125{
[f0af17]126  unsigned long R = (unsigned long)a + (unsigned long)b;
[755296]127  return (number)(R >= (unsigned long)r->ch ? R - (unsigned long)r->ch : R);
[35aab3]128}
[bb865e]129static inline void npInpAddM(number &a, number b, const coeffs r)
130{
131  unsigned long R = (unsigned long)a + (unsigned long)b;
[755296]132  a=(number)(R >= (unsigned long)r->ch ? R - (unsigned long)r->ch : R);
[bb865e]133}
[7d90aa]134static inline number npSubM(number a, number b, const coeffs r)
[35aab3]135{
[6c56a8]136  return (number)((long)a<(long)b ?
[e77676]137                       r->ch-(long)b+(long)a : (long)a-(long)b);
[35aab3]138}
139#else
[7d90aa]140static inline number npAddM(number a, number b, const coeffs r)
[35aab3]141{
[755296]142   unsigned long res = ((unsigned long)a + (unsigned long)b);
[e77676]143   res -= r->ch;
[6c56a8]144#if SIZEOF_LONG == 8
[f0af17]145   res += ((long)res >> 63) & r->ch;
[6c56a8]146#else
[f0af17]147   res += ((long)res >> 31) & r->ch;
[6c56a8]148#endif
[35aab3]149   return (number)res;
150}
[bb865e]151static inline void npInpAddM(number &a, number b, const coeffs r)
152{
[755296]153   unsigned long res = ((unsigned long)a + (unsigned long)b);
[bb865e]154   res -= r->ch;
155#if SIZEOF_LONG == 8
156   res += ((long)res >> 63) & r->ch;
157#else
158   res += ((long)res >> 31) & r->ch;
159#endif
160   a=(number)res;
161}
[7d90aa]162static inline number npSubM(number a, number b, const coeffs r)
[35aab3]163{
[6c56a8]164   long res = ((long)a - (long)b);
165#if SIZEOF_LONG == 8
[e77676]166   res += (res >> 63) & r->ch;
[6c56a8]167#else
[e77676]168   res += (res >> 31) & r->ch;
[6c56a8]169#endif
[35aab3]170   return (number)res;
171}
[40971d]172#endif
[e8cd08]173
[40971d]174static inline number npNegM(number a, const coeffs r)
[e8cd08]175{
[40971d]176  return (number)((long)(r->ch)-(long)(a));
[e8cd08]177}
178
[6e969c]179static inline BOOLEAN npIsOne (number a, const coeffs)
180{
181  return 1 == (long)a;
182}
183
184static inline long npInvMod(long a, const coeffs R)
185{
[755296]186   long s;
[6e969c]187
[755296]188   long  u, v, u0, u1, u2, q, r;
[6e969c]189
190   assume(a>0);
191   u1=1; u2=0;
192   u = a; v = R->ch;
193
[647b762]194   do
[6e969c]195   {
196      q = u / v;
197      //r = u % v;
198      r = u - q*v;
199      u = v;
200      v = r;
201      u0 = u2;
202      u2 = u1 - q*u2;
203      u1 = u0;
[647b762]204   } while (v != 0);
[6e969c]205
206   assume(u==1);
207   s = u1;
208#ifdef HAVE_GENERIC_ADD
209   if (s < 0)
210      return s + R->ch;
211   else
212      return s;
213#else
214  #if SIZEOF_LONG == 8
215  s += (s >> 63) & R->ch;
216  #else
217  s += (s >> 31) & R->ch;
218  #endif
219  return s;
220#endif
221}
222
223static inline number npInversM (number c, const coeffs r)
224{
225  n_Test(c, r);
226#ifndef HAVE_GENERIC_MULT
227  #ifndef HAVE_INVTABLE
228  number d = (number)(long)r->npExpTable[r->npPminus1M - r->npLogTable[(long)c]];
229  #else
230  long inv=(long)r->npInvTable[(long)c];
231  if (inv==0)
232  {
233    inv = (long)r->npExpTable[r->npPminus1M - r->npLogTable[(long)c]];
234    r->npInvTable[(long)c]=inv;
235  }
236  number d = (number)inv;
237  #endif
238#else
239  #ifdef HAVE_INVTABLE
240  long inv=(long)r->npInvTable[(long)c];
241  if (inv==0)
242  {
243    inv=npInvMod((long)c,r);
244    r->npInvTable[(long)c]=inv;
245  }
246  #else
247  long  inv=npInvMod((long)c,r);
248  #endif
249  number d = (number)inv;
250#endif
251  n_Test(d, r);
252  return d;
253}
[35aab3]254
255
[75f460]256// The folloing is reused inside gnumpc.cc, gnumpfl.cc and longrat.cc
[777f8b]257long    npInt         (number &n, const coeffs r);
[35aab3]258
[e77676]259#define npEqualM(A,B,r)  ((A)==(B))
[35aab3]260
261#endif
Note: See TracBrowser for help on using the repository browser.