modulop.h
Go to the documentation of this file.
1 #ifndef MODULOP_H
2 #define MODULOP_H
3 /****************************************
4 * Computer Algebra System SINGULAR *
5 ****************************************/
6 /*
7 * ABSTRACT: numbers modulo p (<=32749)
8 */
9 #include "misc/auxiliary.h"
10 
11 
12 // define if a*b is with mod instead of tables
13 //#define HAVE_GENERIC_MULT
14 // define if 1/b is from tables
15 //#define HAVE_INVTABLE
16 // define if an if should be used
17 //#define HAVE_GENERIC_ADD
18 
19 //#undef HAVE_GENERIC_ADD
20 //#undef HAVE_GENERIC_MULT
21 //#undef HAVE_INVTABLE
22 
23 //#define HAVE_GENERIC_ADD 1
24 //#define HAVE_GENERIC_MULT 1
25 //#define HAVE_INVTABLE 1
26 
27 // enable large primes (32749 < p < 2^31-)
28 #define NV_OPS
29 #define NV_MAX_PRIME 32749
30 #define FACTORY_MAX_PRIME 536870909
31 
32 #ifdef USE_NTL_XGCD
33 // in ntl.a
34 extern void XGCD(long& d, long& s, long& t, long a, long b);
35 #endif
36 
37 struct n_Procs_s; typedef struct n_Procs_s *coeffs;
38 struct snumber; typedef struct snumber * number;
39 
40 BOOLEAN npInitChar(coeffs r, void* p);
41 
42 // inline number npMultM(number a, number b, int npPrimeM)
43 // // return (a*b)%n
44 // {
45 // double ab;
46 // long q, res;
47 //
48 // ab = ((double) ((int)a)) * ((double) ((int)b));
49 // q = (long) (ab/((double) npPrimeM)); // q could be off by (+/-) 1
50 // res = (long) (ab - ((double) q)*((double) npPrimeM));
51 // res += (res >> 31) & npPrimeM;
52 // res -= npPrimeM;
53 // res += (res >> 31) & npPrimeM;
54 // return (number)res;
55 // }
56 #ifdef HAVE_GENERIC_MULT
57 static inline number npMultM(number a, number b, const coeffs r)
58 {
59  return (number)
60  ((((unsigned long) a)*((unsigned long) b)) % ((unsigned long) r->ch));
61 }
62 static inline void npInpMultM(number &a, number b, const coeffs r)
63 {
64  a=(number)
65  ((((unsigned long) a)*((unsigned long) b)) % ((unsigned long) r->ch));
66 }
67 #else
68 static inline number npMultM(number a, number b, const coeffs r)
69 {
70  long x = (long)r->npLogTable[(long)a]+ r->npLogTable[(long)b];
71  #ifdef HAVE_GENERIC_ADD
72  if (x>=r->npPminus1M) x-=r->npPminus1M;
73  #else
74  x-=r->npPminus1M;
75  #if SIZEOF_LONG == 8
76  x += (x >> 63) & r->npPminus1M;
77  #else
78  x += (x >> 31) & r->npPminus1M;
79  #endif
80  #endif
81  return (number)(long)r->npExpTable[x];
82 }
83 static inline void npInpMultM(number &a, number b, const coeffs r)
84 {
85  long x = (long)r->npLogTable[(long)a]+ r->npLogTable[(long)b];
86  #ifdef HAVE_GENERIC_ADD
87  if (x>=r->npPminus1M) x-=r->npPminus1M;
88  #else
89  x-=r->npPminus1M;
90  #if SIZEOF_LONG == 8
91  x += (x >> 63) & r->npPminus1M;
92  #else
93  x += (x >> 31) & r->npPminus1M;
94  #endif
95  #endif
96  a=(number)(long)r->npExpTable[x];
97 }
98 #endif
99 
100 #if 0
101 inline number npAddAsm(number a, number b, int m)
102 {
103  number r;
104  asm ("addl %2, %1; cmpl %3, %1; jb 0f; subl %3, %1; 0:"
105  : "=&r" (r)
106  : "%0" (a), "g" (b), "g" (m)
107  : "cc");
108  return r;
109 }
110 inline number npSubAsm(number a, number b, int m)
111 {
112  number r;
113  asm ("subl %2, %1; jnc 0f; addl %3, %1; 0:"
114  : "=&r" (r)
115  : "%0" (a), "g" (b), "g" (m)
116  : "cc");
117  return r;
118 }
119 #endif
120 #ifdef HAVE_GENERIC_ADD
121 static inline number npAddM(number a, number b, const coeffs r)
122 {
123  unsigned long R = (unsigned long)a + (unsigned long)b;
124  return (number)(R >= r->ch ? R - r->ch : R);
125 }
126 static inline void npInpAddM(number &a, number b, const coeffs r)
127 {
128  unsigned long R = (unsigned long)a + (unsigned long)b;
129  a=(number)(R >= r->ch ? R - r->ch : R);
130 }
131 static inline number npSubM(number a, number b, const coeffs r)
132 {
133  return (number)((long)a<(long)b ?
134  r->ch-(long)b+(long)a : (long)a-(long)b);
135 }
136 #else
137 static inline number npAddM(number a, number b, const coeffs r)
138 {
139  unsigned long res = (long)((unsigned long)a + (unsigned long)b);
140  res -= r->ch;
141 #if SIZEOF_LONG == 8
142  res += ((long)res >> 63) & r->ch;
143 #else
144  res += ((long)res >> 31) & r->ch;
145 #endif
146  return (number)res;
147 }
148 static inline void npInpAddM(number &a, number b, const coeffs r)
149 {
150  unsigned long res = (long)((unsigned long)a + (unsigned long)b);
151  res -= r->ch;
152 #if SIZEOF_LONG == 8
153  res += ((long)res >> 63) & r->ch;
154 #else
155  res += ((long)res >> 31) & r->ch;
156 #endif
157  a=(number)res;
158 }
159 static inline number npSubM(number a, number b, const coeffs r)
160 {
161  long res = ((long)a - (long)b);
162 #if SIZEOF_LONG == 8
163  res += (res >> 63) & r->ch;
164 #else
165  res += (res >> 31) & r->ch;
166 #endif
167  return (number)res;
168 }
169 #endif
170 
171 static inline number npNegM(number a, const coeffs r)
172 {
173  return (number)((long)(r->ch)-(long)(a));
174 }
175 
176 static inline BOOLEAN npIsZeroM (number a, const coeffs)
177 {
178  return 0 == (long)a;
179 }
180 static inline BOOLEAN npIsOne (number a, const coeffs)
181 {
182  return 1 == (long)a;
183 }
184 
185 static inline long npInvMod(long a, const coeffs R)
186 {
187  long s, t;
188 
189 #ifdef USE_NTL_XGCD
190  long d;
191  XGCD(d, s, t, a, R->ch);
192  assume (d == 1);
193 #else
194  long u, v, u0, v0, u1, u2, q, r;
195 
196  assume(a>0);
197  u1=1; u2=0;
198  u = a; v = R->ch;
199 
200  while (v != 0)
201  {
202  q = u / v;
203  //r = u % v;
204  r = u - q*v;
205  u = v;
206  v = r;
207  u0 = u2;
208  u2 = u1 - q*u2;
209  u1 = u0;
210  }
211 
212  assume(u==1);
213  s = u1;
214 #endif
215 #ifdef HAVE_GENERIC_ADD
216  if (s < 0)
217  return s + R->ch;
218  else
219  return s;
220 #else
221  #if SIZEOF_LONG == 8
222  s += (s >> 63) & R->ch;
223  #else
224  s += (s >> 31) & R->ch;
225  #endif
226  return s;
227 #endif
228 }
229 
230 static inline number npInversM (number c, const coeffs r)
231 {
232  n_Test(c, r);
233 #ifndef HAVE_GENERIC_MULT
234  #ifndef HAVE_INVTABLE
235  number d = (number)(long)r->npExpTable[r->npPminus1M - r->npLogTable[(long)c]];
236  #else
237  long inv=(long)r->npInvTable[(long)c];
238  if (inv==0)
239  {
240  inv = (long)r->npExpTable[r->npPminus1M - r->npLogTable[(long)c]];
241  r->npInvTable[(long)c]=inv;
242  }
243  number d = (number)inv;
244  #endif
245 #else
246  #ifdef HAVE_INVTABLE
247  long inv=(long)r->npInvTable[(long)c];
248  if (inv==0)
249  {
250  inv=npInvMod((long)c,r);
251  r->npInvTable[(long)c]=inv;
252  }
253  #else
254  long inv=npInvMod((long)c,r);
255  #endif
256  number d = (number)inv;
257 #endif
258  n_Test(d, r);
259  return d;
260 }
261 
262 
263 // The folloing is reused inside gnumpc.cc, gnumpfl.cc and longrat.cc
264 long npInt (number &n, const coeffs r);
265 
266 // The folloing is reused inside tgb*.cc
267 number npMult (number a, number b, const coeffs r);
268 // The following is currently used in OPAE.cc, OPAEQ.cc and OPAEp.cc for setting their SetMap...
269 nMapFunc npSetMap(const coeffs src, const coeffs dst); // FIXME! BUG?
270 
271 #define npEqualM(A,B,r) ((A)==(B))
272 
273 
274 #endif
const CanonicalForm int s
Definition: facAbsFact.cc:55
static number npInversM(number c, const coeffs r)
Definition: modulop.h:230
static number npMultM(number a, number b, const coeffs r)
Definition: modulop.h:68
&#39;SR_INT&#39; is the type of those integers small enough to fit into 29 bits.
Definition: longrat.h:47
number npMult(number a, number b, const coeffs r)
Definition: modulop.cc:86
BOOLEAN npInitChar(coeffs r, void *p)
Definition: modulop.cc:377
static number npNegM(number a, const coeffs r)
Definition: modulop.h:171
static number npSubM(number a, number b, const coeffs r)
Definition: modulop.h:131
CanonicalForm b
Definition: cfModGcd.cc:4044
mpz_t n
Definition: longrat.h:50
CanonicalForm res
Definition: facAbsFact.cc:64
#define assume(x)
Definition: mod2.h:390
The main handler for Singular numbers which are suitable for Singular polynomials.
static BOOLEAN npIsOne(number a, const coeffs)
Definition: modulop.h:180
number(* nMapFunc)(number a, const coeffs src, const coeffs dst)
maps "a", which lives in src, into dst
Definition: coeffs.h:73
#define n_Test(a, r)
BOOLEAN n_Test(number a, const coeffs r)
Definition: coeffs.h:738
static number npAddM(number a, number b, const coeffs r)
Definition: modulop.h:121
All the auxiliary stuff.
int m
Definition: cfEzgcd.cc:121
static long npInvMod(long a, const coeffs R)
Definition: modulop.h:185
static void npInpAddM(number &a, number b, const coeffs r)
Definition: modulop.h:126
const Variable & v
< [in] a sqrfree bivariate poly
Definition: facBivar.h:37
static BOOLEAN npIsZeroM(number a, const coeffs)
Definition: modulop.h:176
#define R
Definition: sirandom.c:27
Variable x
Definition: cfModGcd.cc:4023
long npInt(number &n, const coeffs r)
Definition: modulop.cc:127
nMapFunc npSetMap(const coeffs src, const coeffs dst)
Definition: modulop.cc:654
int p
Definition: cfModGcd.cc:4019
static void npInpMultM(number &a, number b, const coeffs r)
Definition: modulop.h:83
int BOOLEAN
Definition: auxiliary.h:87