My Project
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 #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
18 
19 
20 // define if a*b is with mod instead of tables
21 //#define HAVE_GENERIC_MULT
22 // define if 1/b is from tables
23 //#define HAVE_INVTABLE
24 // define if an if should be used(NTL does not define NTL_AVOID_BRANCHING)
25 //#define HAVE_GENERIC_ADD
26 
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 
35 // enable large primes (32749 < p < 2^31-)
36 #define NV_OPS
37 #define NV_MAX_PRIME 32749
38 #define FACTORY_MAX_PRIME 536870909
39 
40 struct n_Procs_s; typedef struct n_Procs_s *coeffs;
41 struct snumber; typedef struct snumber * number;
42 
43 BOOLEAN npInitChar(coeffs r, void* p);
44 
45 // inline number npMultM(number a, number b, int npPrimeM)
46 // // return (a*b)%n
47 // {
48 // double ab;
49 // long q, res;
50 //
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 // }
59 #ifdef HAVE_GENERIC_MULT
60 static inline number npMultM(number a, number b, const coeffs r)
61 {
62  return (number)
63  ((((unsigned long) a)*((unsigned long) b)) % ((unsigned long) r->ch));
64 }
65 static 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 }
70 #else
71 static inline number npMultM(number a, number b, const coeffs r)
72 {
73  long x = (long)r->npLogTable[(long)a]+ r->npLogTable[(long)b];
74  #ifdef HAVE_GENERIC_ADD
75  if (x>=r->npPminus1M) x-=r->npPminus1M;
76  #else
77  x-=r->npPminus1M;
78  #if SIZEOF_LONG == 8
79  x += (x >> 63) & r->npPminus1M;
80  #else
81  x += (x >> 31) & r->npPminus1M;
82  #endif
83  #endif
84  return (number)(long)r->npExpTable[x];
85 }
86 static inline void npInpMultM(number &a, number b, const coeffs r)
87 {
88  long x = (long)r->npLogTable[(long)a]+ r->npLogTable[(long)b];
89  #ifdef HAVE_GENERIC_ADD
90  if (x>=r->npPminus1M) x-=r->npPminus1M;
91  #else
92  x-=r->npPminus1M;
93  #if SIZEOF_LONG == 8
94  x += (x >> 63) & r->npPminus1M;
95  #else
96  x += (x >> 31) & r->npPminus1M;
97  #endif
98  #endif
99  a=(number)(long)r->npExpTable[x];
100 }
101 #endif
102 
103 #if 0
104 inline 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 }
113 inline 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
124 static inline number npAddM(number a, number b, const coeffs r)
125 {
126  unsigned long R = (unsigned long)a + (unsigned long)b;
127  return (number)(R >= (unsigned long)r->ch ? R - (unsigned long)r->ch : R);
128 }
129 static inline void npInpAddM(number &a, number b, const coeffs r)
130 {
131  unsigned long R = (unsigned long)a + (unsigned long)b;
132  a=(number)(R >= (unsigned long)r->ch ? R - (unsigned long)r->ch : R);
133 }
134 static inline number npSubM(number a, number b, const coeffs r)
135 {
136  return (number)((long)a<(long)b ?
137  r->ch-(long)b+(long)a : (long)a-(long)b);
138 }
139 #else
140 static inline number npAddM(number a, number b, const coeffs r)
141 {
142  unsigned long res = ((unsigned long)a + (unsigned long)b);
143  res -= r->ch;
144 #if SIZEOF_LONG == 8
145  res += ((long)res >> 63) & r->ch;
146 #else
147  res += ((long)res >> 31) & r->ch;
148 #endif
149  return (number)res;
150 }
151 static inline void npInpAddM(number &a, number b, const coeffs r)
152 {
153  unsigned long res = ((unsigned long)a + (unsigned long)b);
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 }
162 static inline number npSubM(number a, number b, const coeffs r)
163 {
164  long res = ((long)a - (long)b);
165 #if SIZEOF_LONG == 8
166  res += (res >> 63) & r->ch;
167 #else
168  res += (res >> 31) & r->ch;
169 #endif
170  return (number)res;
171 }
172 #endif
173 
174 static inline number npNegM(number a, const coeffs r)
175 {
176  return (number)((long)(r->ch)-(long)(a));
177 }
178 
179 static inline BOOLEAN npIsOne (number a, const coeffs)
180 {
181  return 1 == (long)a;
182 }
183 
184 static inline long npInvMod(long a, const coeffs R)
185 {
186  long s;
187 
188  long u, v, u0, u1, u2, q, r;
189 
190  assume(a>0);
191  u1=1; u2=0;
192  u = a; v = R->ch;
193 
194  do
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;
204  } while (v != 0);
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 
223 static 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 }
254 
255 
256 // The folloing is reused inside gnumpc.cc, gnumpfl.cc and longrat.cc
257 long npInt (number &n, const coeffs r);
258 
259 #define npEqualM(A,B,r) ((A)==(B))
260 
261 #endif
All the auxiliary stuff.
int BOOLEAN
Definition: auxiliary.h:87
int m
Definition: cfEzgcd.cc:128
Variable x
Definition: cfModGcd.cc:4082
int p
Definition: cfModGcd.cc:4078
CanonicalForm b
Definition: cfModGcd.cc:4103
#define n_Test(a, r)
BOOLEAN n_Test(number a, const coeffs r)
Definition: coeffs.h:709
const CanonicalForm int s
Definition: facAbsFact.cc:51
CanonicalForm res
Definition: facAbsFact.cc:60
const Variable & v
< [in] a sqrfree bivariate poly
Definition: facBivar.h:39
mpz_t n
Definition: longrat.h:51
'SR_INT' is the type of those integers small enough to fit into 29 bits.
Definition: longrat.h:49
#define assume(x)
Definition: mod2.h:389
static BOOLEAN npIsOne(number a, const coeffs)
Definition: modulop.h:179
static number npAddM(number a, number b, const coeffs r)
Definition: modulop.h:124
static number npMultM(number a, number b, const coeffs r)
Definition: modulop.h:71
BOOLEAN npInitChar(coeffs r, void *p)
Definition: modulop.cc:338
static number npNegM(number a, const coeffs r)
Definition: modulop.h:174
static void npInpMultM(number &a, number b, const coeffs r)
Definition: modulop.h:86
static long npInvMod(long a, const coeffs R)
Definition: modulop.h:184
static number npInversM(number c, const coeffs r)
Definition: modulop.h:223
long npInt(number &n, const coeffs r)
Definition: modulop.cc:83
static void npInpAddM(number &a, number b, const coeffs r)
Definition: modulop.h:129
static number npSubM(number a, number b, const coeffs r)
Definition: modulop.h:134
The main handler for Singular numbers which are suitable for Singular polynomials.
#define R
Definition: sirandom.c:27