1 /****************************************
2 * Computer Algebra System SINGULAR *
3 ****************************************/
4 /*
5 * ABSTRACT: numbers modulo p (<=32749)
6 */
7
8 #include "misc/auxiliary.h"
9
10 #include "factory/factory.h"
11
12 #include "misc/mylimits.h"
13 #include "misc/sirandom.h"
14
15 #include "reporter/reporter.h"
16
17 #include "coeffs/coeffs.h"
18 #include "coeffs/numbers.h"
19 #include "coeffs/mpr_complex.h"
20
21 #include "coeffs/longrat.h"
22 #include "coeffs/modulop.h"
23
24 #include <string.h>
25
26 BOOLEAN npGreaterZero (number k, const coeffs r);
27 long npInt (number &n, const coeffs r);
28 void npPower (number a, int i, number * result,const coeffs r);
29 BOOLEAN npIsMOne (number a,const coeffs r);
30 number npDiv (number a, number b,const coeffs r);
31 number npNeg (number c,const coeffs r);
32 number npInvers (number c,const coeffs r);
33 BOOLEAN npGreater (number a, number b,const coeffs r);
34 BOOLEAN npEqual (number a, number b,const coeffs r);
35 void npWrite (number a, const coeffs r);
36 const char * npRead (const char *s, number *a,const coeffs r);
37 void nvInpMult(number &a, number b, const coeffs r);
38
39 #ifdef LDEBUG
40 BOOLEAN npDBTest (number a, const char *f, const int l, const coeffs r);
41 #endif
42
43 nMapFunc npSetMap(const coeffs src, const coeffs dst);
44
45 #include "coeffs/modulop_inl.h" // npMult, npInit
46
47 #ifdef NV_OPS
48 number nvDiv (number a, number b, const coeffs r);
49 number nvInvers (number c, const coeffs r);
50 //void nvPower (number a, int i, number * result, const coeffs r);
51 #endif
52
npGreaterZero(number k,const coeffs r)53 BOOLEAN npGreaterZero (number k, const coeffs r)
54 {
55 n_Test(k, r);
56
57 int h = (int)((long) k);
58 return ((int)h !=0) && (h <= (r->ch>>1));
59 }
60
61 //unsigned long npMultMod(unsigned long a, unsigned long b, int npPrimeM)
62 //{
63 // unsigned long c = a*b;
64 // c = c % npPrimeM;
65 // assume(c == (unsigned long) npMultM((number) a, (number) b, npPrimeM));
66 // return c;
67 //}
68
69
npInpMult(number & a,number b,const coeffs r)70 void npInpMult (number &a,number b, const coeffs r)
71 {
72 n_Test(a, r);
73 n_Test(b, r);
74
75 if (((long)a == 0) || ((long)b == 0))
76 a=(number)0;
77 else
78 a = npMultM(a,b, r);
79 n_Test(a, r);
80 }
81
82 /*2
83 * convert a number to an int in (-p/2 .. p/2]
84 */
npInt(number & n,const coeffs r)85 long npInt(number &n, const coeffs r)
86 {
87 n_Test(n, r);
88
89 if ((long)n > (((long)r->ch) >>1)) return ((long)n -((long)r->ch));
90 else return ((long)n);
91 }
92
npIsMOne(number a,const coeffs r)93 BOOLEAN npIsMOne (number a, const coeffs r)
94 {
95 n_Test(a, r);
96
97 return ((r->npPminus1M == (long)a) &&(1L!=(long)a))/*for char 2*/;
98 }
99
npDiv(number a,number b,const coeffs r)100 number npDiv (number a,number b, const coeffs r)
101 {
102 n_Test(a, r);
103 n_Test(b, r);
104
105 if ((long)b==0L)
106 {
107 WerrorS(nDivBy0);
108 return (number)0L;
109 }
110 if ((long)a==0) return (number)0L;
111
112 number d;
113 #ifndef HAVE_GENERIC_MULT
114 int s = r->npLogTable[(long)a] - r->npLogTable[(long)b];
115 #ifdef HAVE_GENERIC_ADD
116 if (s < 0)
117 s += r->npPminus1M;
118 #else
119 #if SIZEOF_LONG == 8
120 s += ((long)s >> 63) & r->npPminus1M;
121 #else
122 s += ((long)s >> 31) & r->npPminus1M;
123 #endif
124 #endif
125 d = (number)(long)r->npExpTable[s];
126 #else
127 number inv=npInversM(b,r);
128 d = npMultM(a,inv,r);
129 #endif
130
131 n_Test(d, r);
132 return d;
133
134 }
npInvers(number c,const coeffs r)135 number npInvers (number c, const coeffs r)
136 {
137 n_Test(c, r);
138
139 if ((long)c==0L)
140 {
141 WerrorS("1/0");
142 return (number)0L;
143 }
144 number d = npInversM(c,r);
145
146 n_Test(d, r);
147 return d;
148 }
149
npNeg(number c,const coeffs r)150 number npNeg (number c, const coeffs r)
151 {
152 n_Test(c, r);
153
154 if ((long)c==0L) return c;
155
156 #if 0
157 number d = npNegM(c,r);
158 n_Test(d, r);
159 return d;
160 #else
161 c = npNegM(c,r);
162 n_Test(c, r);
163 return c;
164 #endif
165 }
166
npGreater(number a,number b,const coeffs r)167 BOOLEAN npGreater (number a,number b, const coeffs r)
168 {
169 n_Test(a, r);
170 n_Test(b, r);
171
172 //return (long)a != (long)b;
173 return ((long)a) > ((long)b);
174 }
175
npEqual(number a,number b,const coeffs r)176 BOOLEAN npEqual (number a,number b, const coeffs r)
177 {
178 n_Test(a, r);
179 n_Test(b, r);
180
181 // return (long)a == (long)b;
182
183 return npEqualM(a,b,r);
184 }
185
npWrite(number a,const coeffs r)186 void npWrite (number a, const coeffs r)
187 {
188 n_Test(a, r);
189
190 if ((long)a>(((long)r->ch) >>1)) StringAppend("-%d",(int)(((long)r->ch)-((long)a)));
191 else StringAppend("%d",(int)((long)a));
192 }
193
194 #if 0
195 void npPower (number a, int i, number * result, const coeffs r)
196 {
197 n_Test(a, r);
198
199 if (i==0)
200 {
201 //npInit(1,result);
202 *(long *)result = 1;
203 }
204 else if (i==1)
205 {
206 *result = a;
207 }
208 else
209 {
210 npPower(a,i-1,result,r);
211 *result = npMultM(a,*result,r);
212 }
213 }
214 #endif
215
npEati(const char * s,int * i,const coeffs r)216 static inline const char* npEati(const char *s, int *i, const coeffs r)
217 {
218 return nEati((char *)s,i,(int)r->ch);
219 }
220
npRead(const char * s,number * a,const coeffs r)221 const char * npRead (const char *s, number *a, const coeffs r)
222 {
223 int z;
224 int n=1;
225
226 s = npEati(s, &z, r);
227 if ((*s) == '/')
228 {
229 s++;
230 s = npEati(s, &n, r);
231 }
232 if (n == 1)
233 *a = (number)(long)z;
234 else
235 {
236 if ((z==0)&&(n==0))
237 {
238 WerrorS(nDivBy0);
239 *a=(number)0L;
240 }
241 else
242 {
243 #ifdef NV_OPS
244 if (r->ch>NV_MAX_PRIME)
245 *a = nvDiv((number)(long)z,(number)(long)n,r);
246 else
247 #endif
248 *a = npDiv((number)(long)z,(number)(long)n,r);
249 }
250 }
251 n_Test(*a, r);
252 return s;
253 }
254
255 /*2
256 * set the charcteristic (allocate and init tables)
257 */
258
npKillChar(coeffs r)259 void npKillChar(coeffs r)
260 {
261 #ifdef HAVE_INVTABLE
262 if (r->npInvTable!=NULL)
263 {
264 omFreeSize( (void *)r->npInvTable, r->ch*sizeof(unsigned short) );
265 r->npInvTable=NULL;
266 }
267 #endif
268 #ifndef HAVE_GENERIC_MULT
269 if (r->npExpTable!=NULL)
270 {
271 omFreeSize( (void *)r->npExpTable, r->ch*sizeof(unsigned short) );
272 omFreeSize( (void *)r->npLogTable, r->ch*sizeof(unsigned short) );
273 r->npExpTable=NULL; r->npLogTable=NULL;
274 }
275 #endif
276 }
277
npCoeffsEqual(const coeffs r,n_coeffType n,void * parameter)278 static BOOLEAN npCoeffsEqual(const coeffs r, n_coeffType n, void * parameter)
279 {
280 /* test, if r is an instance of nInitCoeffs(n,parameter) */
281 return (n==n_Zp) && (r->ch==(int)(long)parameter);
282 }
npConvSingNFactoryN(number n,BOOLEAN setChar,const coeffs r)283 CanonicalForm npConvSingNFactoryN( number n, BOOLEAN setChar, const coeffs r )
284 {
285 if (setChar) setCharacteristic( r->ch );
286 return CanonicalForm(npInt( n,r ));
287 }
288
npConvFactoryNSingN(const CanonicalForm n,const coeffs r)289 number npConvFactoryNSingN( const CanonicalForm n, const coeffs r)
290 {
291 if (n.isImm())
292 {
293 return npInit(n.intval(),r);
294 }
295 else
296 {
297 assume(0);
298 return NULL;
299 }
300 }
301
npCoeffName(const coeffs cf)302 static char* npCoeffName(const coeffs cf)
303 {
304 STATIC_VAR char npCoeffName_buf[15];
305 snprintf(npCoeffName_buf,14,"ZZ/%d",cf->ch);
306 return npCoeffName_buf;
307 }
308
npWriteFd(number n,const ssiInfo * d,const coeffs)309 static void npWriteFd(number n, const ssiInfo* d, const coeffs)
310 {
311 fprintf(d->f_write,"%d ",(int)(long)n);
312 }
313
npReadFd(const ssiInfo * d,const coeffs)314 static number npReadFd(const ssiInfo *d, const coeffs)
315 {
316 // read int
317 int dd;
318 dd=s_readint(d->f_read);
319 return (number)(long)dd;
320 }
321
npRandom(siRandProc p,number,number,const coeffs cf)322 static number npRandom(siRandProc p, number, number, const coeffs cf)
323 {
324 return npInit(p(),cf);
325 }
326
327
328 #ifndef HAVE_GENERIC_MULT
npPar(int,coeffs r)329 static number npPar(int, coeffs r)
330 {
331 return (number)(long)r->npExpTable[1];
332 }
333 #endif
334
npInitMPZ(mpz_t m,const coeffs r)335 static number npInitMPZ(mpz_t m, const coeffs r)
336 {
337 return (number)mpz_fdiv_ui(m, r->ch);
338 }
339
npInitChar(coeffs r,void * p)340 BOOLEAN npInitChar(coeffs r, void* p)
341 {
342 assume( getCoeffType(r) == n_Zp );
343 const int c = (int) (long) p;
344
345 assume( c > 0 );
346
347 int i, w;
348
349 r->is_field=TRUE;
350 r->is_domain=TRUE;
351 r->rep=n_rep_int;
352
353 r->ch = c;
354 r->npPminus1M = c /*r->ch*/ - 1;
355
356 //r->cfInitChar=npInitChar;
357 r->cfKillChar=npKillChar;
358 r->nCoeffIsEqual=npCoeffsEqual;
359 r->cfCoeffName=npCoeffName;
360
361 r->cfMult = npMult;
362 r->cfInpMult = npInpMult;
363 r->cfSub = npSubM;
364 r->cfAdd = npAddM;
365 r->cfInpAdd = npInpAddM;
366 r->cfDiv = npDiv;
367 r->cfInit = npInit;
368 //r->cfSize = ndSize;
369 r->cfInt = npInt;
370 r->cfInitMPZ = npInitMPZ;
371 #ifdef HAVE_RINGS
372 //r->cfDivComp = NULL; // only for ring stuff
373 //r->cfIsUnit = NULL; // only for ring stuff
374 //r->cfGetUnit = NULL; // only for ring stuff
375 //r->cfExtGcd = NULL; // only for ring stuff
376 // r->cfDivBy = NULL; // only for ring stuff
377 #endif
378 r->cfInpNeg = npNeg;
379 r->cfInvers= npInvers;
380 //r->cfCopy = ndCopy;
381 //r->cfRePart = ndCopy;
382 //r->cfImPart = ndReturn0;
383 r->cfWriteLong = npWrite;
384 r->cfRead = npRead;
385 //r->cfNormalize=ndNormalize;
386 r->cfGreater = npGreater;
387 r->cfEqual = npEqual;
388 r->cfIsZero = npIsZero;
389 r->cfIsOne = npIsOne;
390 r->cfIsMOne = npIsMOne;
391 r->cfGreaterZero = npGreaterZero;
392 //r->cfPower = npPower;
393 //r->cfGetDenom = ndGetDenom;
394 //r->cfGetNumerator = ndGetNumerator;
395 //r->cfGcd = ndGcd;
396 //r->cfLcm = ndGcd;
397 //r->cfDelete= ndDelete;
398 r->cfSetMap = npSetMap;
399 //r->cfName = ndName;
400 //r->cfInpMult=ndInpMult;
401 r->convSingNFactoryN=npConvSingNFactoryN;
402 r->convFactoryNSingN=npConvFactoryNSingN;
403 r->cfRandom=npRandom;
404 #ifdef LDEBUG
405 // debug stuff
406 r->cfDBTest=npDBTest;
407 #endif
408
409 // io via ssi
410 r->cfWriteFd=npWriteFd;
411 r->cfReadFd=npReadFd;
412
413 // the variables:
414 r->type = n_Zp;
415 r->has_simple_Alloc=TRUE;
416 r->has_simple_Inverse=TRUE;
417
418 // the tables
419 #ifdef NV_OPS
420 if (r->ch <=NV_MAX_PRIME)
421 #endif
422 {
423 #ifdef HAVE_INVTABLE
424 r->npInvTable=(unsigned short*)omAlloc0( r->ch*sizeof(unsigned short) );
425 #endif
426 #ifndef HAVE_GENERIC_MULT
427 r->cfParameter=npPar; /* Singular.jl */
428 r->npExpTable=(unsigned short *)omAlloc0( r->ch*sizeof(unsigned short) );
429 r->npLogTable=(unsigned short *)omAlloc0( r->ch*sizeof(unsigned short) );
430 r->npExpTable[0] = 1;
431 r->npLogTable[0] = 0;
432 if (r->ch > 2)
433 {
434 w = 1;
435 loop
436 {
437 r->npLogTable[1] = 0;
438 w++;
439 i = 0;
440 loop
441 {
442 i++;
443 r->npExpTable[i] =(int)(((long)w * (long)r->npExpTable[i-1]) % r->ch);
444 r->npLogTable[r->npExpTable[i]] = i;
445 if /*(i == r->ch - 1 ) ||*/ (/*(*/ r->npExpTable[i] == 1 /*)*/)
446 break;
447 }
448 if (i == r->ch - 1)
449 break;
450 }
451 }
452 else
453 {
454 r->npExpTable[1] = 1;
455 r->npLogTable[1] = 0;
456 }
457 #endif
458 }
459 #ifdef NV_OPS
460 else /*if (c>NV_MAX_PRIME)*/
461 {
462 r->cfMult = nvMult;
463 r->cfDiv = nvDiv;
464 r->cfExactDiv = nvDiv;
465 r->cfInvers = nvInvers;
466 r->cfInpMult = nvInpMult;
467 //r->cfPower= nvPower;
468 //if (c>FACTORY_MAX_PRIME) // factory will catch this error
469 //{
470 // r->convSingNFactoryN=ndConvSingNFactoryN;
471 //}
472 }
473 #endif
474 return FALSE;
475 }
476
477 #ifdef LDEBUG
npDBTest(number a,const char * f,const int l,const coeffs r)478 BOOLEAN npDBTest (number a, const char *f, const int l, const coeffs r)
479 {
480 if (((long)a<0L) || ((long)a>(long)r->ch))
481 {
482 Print("wrong mod p number %ld at %s,%d\n",(long)a,f,l);
483 return FALSE;
484 }
485 return TRUE;
486 }
487 #endif
488
npMapP(number from,const coeffs src,const coeffs dst_r)489 static number npMapP(number from, const coeffs src, const coeffs dst_r)
490 {
491 long i = (long)from;
492 if (i>src->ch/2)
493 {
494 i-=src->ch;
495 while (i < 0) i+=dst_r->ch;
496 }
497 i%=dst_r->ch;
498 return (number)i;
499 }
500
npMapLongR(number from,const coeffs,const coeffs dst_r)501 static number npMapLongR(number from, const coeffs /*src*/, const coeffs dst_r)
502 {
503 gmp_float *ff=(gmp_float*)from;
504 mpf_t *f=ff->_mpfp();
505 number res;
506 mpz_ptr dest,ndest;
507 int size,i;
508 int e,al,bl;
509 long iz;
510 mp_ptr qp,dd,nn;
511
512 size = (*f)[0]._mp_size;
513 if (size == 0)
514 return npInit(0,dst_r);
515 if(size<0)
516 size = -size;
517
518 qp = (*f)[0]._mp_d;
519 while(qp[0]==0)
520 {
521 qp++;
522 size--;
523 }
524
525 if(dst_r->ch>2)
526 e=(*f)[0]._mp_exp-size;
527 else
528 e=0;
529 res = ALLOC_RNUMBER();
530 #if defined(LDEBUG)
531 res->debug=123456;
532 #endif
533 dest = res->z;
534
535 long in=0;
536 if (e<0)
537 {
538 al = dest->_mp_size = size;
539 if (al<2) al = 2;
540 dd = (mp_ptr)omAlloc(sizeof(mp_limb_t)*al);
541 for (i=0;i<size;i++) dd[i] = qp[i];
542 bl = 1-e;
543 nn = (mp_ptr)omAlloc(sizeof(mp_limb_t)*bl);
544 nn[bl-1] = 1;
545 for (i=bl-2;i>=0;i--) nn[i] = 0;
546 ndest = res->n;
547 ndest->_mp_d = nn;
548 ndest->_mp_alloc = ndest->_mp_size = bl;
549 res->s = 0;
550 in=mpz_fdiv_ui(ndest,dst_r->ch);
551 mpz_clear(ndest);
552 }
553 else
554 {
555 al = dest->_mp_size = size+e;
556 if (al<2) al = 2;
557 dd = (mp_ptr)omAlloc(sizeof(mp_limb_t)*al);
558 for (i=0;i<size;i++) dd[i+e] = qp[i];
559 for (i=0;i<e;i++) dd[i] = 0;
560 res->s = 3;
561 }
562
563 dest->_mp_d = dd;
564 dest->_mp_alloc = al;
565 iz=mpz_fdiv_ui(dest,dst_r->ch);
566 mpz_clear(dest);
567 if(res->s==0)
568 iz=(long)npDiv((number)iz,(number)in,dst_r);
569 FREE_RNUMBER(res); // Q!?
570 return (number)iz;
571 }
572
573 #ifdef HAVE_RINGS
574 /*2
575 * convert from a GMP integer
576 */
npMapGMP(number from,const coeffs,const coeffs dst)577 static number npMapGMP(number from, const coeffs /*src*/, const coeffs dst)
578 {
579 return (number)mpz_fdiv_ui((mpz_ptr) from, dst->ch);
580 }
581
npMapZ(number from,const coeffs src,const coeffs dst)582 static number npMapZ(number from, const coeffs src, const coeffs dst)
583 {
584 if (SR_HDL(from) & SR_INT)
585 {
586 long f_i=SR_TO_INT(from);
587 return npInit(f_i,dst);
588 }
589 return npMapGMP(from,src,dst);
590 }
591
592 /*2
593 * convert from an machine long
594 */
npMapMachineInt(number from,const coeffs,const coeffs dst)595 static number npMapMachineInt(number from, const coeffs /*src*/,const coeffs dst)
596 {
597 long i = (long) (((unsigned long) from) % dst->ch);
598 return (number) i;
599 }
600 #endif
601
npMapCanonicalForm(number a,const coeffs,const coeffs dst)602 static number npMapCanonicalForm (number a, const coeffs /*src*/, const coeffs dst)
603 {
604 setCharacteristic (dst ->ch);
605 CanonicalForm f= CanonicalForm ((InternalCF*)(a));
606 return (number) (f.intval());
607 }
608
npSetMap(const coeffs src,const coeffs)609 nMapFunc npSetMap(const coeffs src, const coeffs)
610 {
611 #ifdef HAVE_RINGS
612 if ((src->rep==n_rep_int) && nCoeff_is_Ring_2toM(src))
613 {
614 return npMapMachineInt;
615 }
616 if (src->rep==n_rep_gmp) //nCoeff_is_Z(src) || nCoeff_is_Ring_PtoM(src) || nCoeff_is_Zn(src))
617 {
618 return npMapGMP;
619 }
620 if (src->rep==n_rep_gap_gmp) //nCoeff_is_Z(src)
621 {
622 return npMapZ;
623 }
624 #endif
625 if (src->rep==n_rep_gap_rat) /* Q, Z */
626 {
627 return nlModP; // npMap0; // FIXME? TODO? // extern number nlModP(number q, const coeffs Q, const coeffs Zp); // Map q \in QQ \to Zp // FIXME!
628 }
629 if ((src->rep==n_rep_int) && nCoeff_is_Zp(src) )
630 {
631 return npMapP;
632 }
633 if ((src->rep==n_rep_gmp_float) && nCoeff_is_long_R(src))
634 {
635 return npMapLongR;
636 }
637 if (nCoeff_is_CF (src))
638 {
639 return npMapCanonicalForm;
640 }
641 return NULL; /* default */
642 }
643
644 // -----------------------------------------------------------
645 // operation for very large primes (32749< p < 2^31-1)
646 // ----------------------------------------------------------
647 #ifdef NV_OPS
nvInpMult(number & a,number b,const coeffs r)648 void nvInpMult(number &a, number b, const coeffs r)
649 {
650 number n=nvMult(a,b,r);
651 a=n;
652 }
653
nvInversM(number c,const coeffs r)654 static inline number nvInversM (number c, const coeffs r)
655 {
656 long inv=npInvMod((long)c,r);
657 return (number)inv;
658 }
659
nvDiv(number a,number b,const coeffs r)660 number nvDiv (number a,number b, const coeffs r)
661 {
662 if ((long)a==0L)
663 return (number)0L;
664 else if ((long)b==0L)
665 {
666 WerrorS(nDivBy0);
667 return (number)0L;
668 }
669 else
670 {
671 number inv=nvInversM(b,r);
672 return nvMult(a,inv,r);
673 }
674 }
nvInvers(number c,const coeffs r)675 number nvInvers (number c, const coeffs r)
676 {
677 if ((long)c==0L)
678 {
679 WerrorS(nDivBy0);
680 return (number)0L;
681 }
682 return nvInversM(c,r);
683 }
684 #if 0
685 void nvPower (number a, int i, number * result, const coeffs r)
686 {
687 if (i==0)
688 {
689 //npInit(1,result);
690 *(long *)result = 1;
691 }
692 else if (i==1)
693 {
694 *result = a;
695 }
696 else
697 {
698 nvPower(a,i-1,result,r);
699 *result = nvMult(a,*result,r);
700 }
701 }
702 #endif
703 #endif
704
705