source: git/libpolys/coeffs/modulop.h @ f081b9

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