- Timestamp:
- Nov 9, 2018, 11:53:46 AM (5 years ago)
- Branches:
- (u'fieker-DuVal', '117eb8c30fc9e991c4decca4832b1d19036c4c65')(u'spielwiese', '4bd32dfef92ec9f5ed8dceee82d14318ae147107')
- Children:
- 0f9f55fe5c567ee77fac4106ff7d41e33dd40d89
- Parents:
- 03ff5f33fb7a5e051b5bbb6dbc0866eddba7c169
- Location:
- libpolys/coeffs
- Files:
-
- 2 added
- 3 edited
-
Makefile.am (modified) (1 diff)
-
rintegers.cc (modified) (3 diffs)
-
rintegers.h (modified) (1 diff)
-
rintegers2.cc (added)
-
rintegers3.cc (added)
Legend:
- Unmodified
- Added
- Removed
-
libpolys/coeffs/Makefile.am
r03ff5f r9e2cf88 43 43 test_SOURCES = test.cc 44 44 test_LDADD = libcoeffs.la $(libcoeffs_la_LIBADD) 45 46 EXTRA_DIST = rintegers2.cc rinteger3.cc -
libpolys/coeffs/rintegers.cc
r03ff5f r9e2cf88 8 8 #include "misc/auxiliary.h" 9 9 10 #ifdef HAVE_RINGS 10 11 #include "factory/factory.h" 11 12 … … 13 14 #include "reporter/reporter.h" 14 15 15 #include "coeffs/coeffs.h" 16 #include "coeffs/rintegers.h" 17 #include "coeffs/mpr_complex.h" 18 #include "coeffs/rmodulon.h" 16 19 #include "coeffs/longrat.h" 17 20 #include "coeffs/numbers.h" 18 21 19 #include "coeffs/si_gmp.h"20 21 #include "coeffs/mpr_complex.h"22 #include "coeffs/rintegers.h"23 #include "coeffs/rmodulon.h"24 #include "coeffs/longrat.h"25 26 22 #include <string.h> 27 23 28 #ifdef HAVE_RINGS29 30 24 omBin gmp_nrz_bin = omGetSpecBin(sizeof(mpz_t)); 31 32 #if SI_INTEGER_VARIANT == 233 /*34 * Multiply two numbers35 */36 static number nrzMult (number a, number b, const coeffs)37 {38 mpz_ptr erg = (mpz_ptr) omAllocBin(gmp_nrz_bin);39 mpz_init(erg);40 mpz_mul(erg, (mpz_ptr) a, (mpz_ptr) b);41 return (number) erg;42 }43 44 /*45 * Give the smallest non unit k, such that a * x = k = b * y has a solution46 */47 static number nrzLcm (number a,number b,const coeffs)48 {49 mpz_ptr erg = (mpz_ptr) omAllocBin(gmp_nrz_bin);50 mpz_init(erg);51 mpz_lcm(erg, (mpz_ptr) a, (mpz_ptr) b);52 return (number) erg;53 }54 55 /*56 * Give the largest non unit k, such that a = x * k, b = y * k has57 * a solution.58 */59 static number nrzGcd (number a,number b,const coeffs)60 {61 mpz_ptr erg = (mpz_ptr) omAllocBin(gmp_nrz_bin);62 mpz_init(erg);63 mpz_gcd(erg, (mpz_ptr) a, (mpz_ptr) b);64 return (number) erg;65 }66 67 /*68 * Give the largest non unit k, such that a = x * k, b = y * k has69 * a solution and r, s, s.t. k = s*a + t*b70 */71 static number nrzExtGcd (number a, number b, number *s, number *t, const coeffs)72 {73 mpz_ptr erg = (mpz_ptr) omAllocBin(gmp_nrz_bin);74 mpz_ptr bs = (mpz_ptr) omAllocBin(gmp_nrz_bin);75 mpz_ptr bt = (mpz_ptr) omAllocBin(gmp_nrz_bin);76 mpz_init(erg);77 mpz_init(bs);78 mpz_init(bt);79 mpz_gcdext(erg, bs, bt, (mpz_ptr) a, (mpz_ptr) b);80 *s = (number) bs;81 *t = (number) bt;82 return (number) erg;83 }84 85 static number nrzXExtGcd (number a, number b, number *s, number *t, number *u, number *v, const coeffs )86 {87 mpz_ptr erg = (mpz_ptr) omAllocBin(gmp_nrz_bin);88 mpz_ptr bs = (mpz_ptr) omAllocBin(gmp_nrz_bin);89 mpz_ptr bt = (mpz_ptr) omAllocBin(gmp_nrz_bin);90 mpz_init(erg);91 mpz_init(bs);92 mpz_init(bt);93 94 mpz_gcdext(erg, bs, bt, (mpz_ptr)a, (mpz_ptr)b);95 96 mpz_ptr bu = (mpz_ptr) omAllocBin(gmp_nrz_bin);97 mpz_ptr bv = (mpz_ptr) omAllocBin(gmp_nrz_bin);98 99 mpz_init_set(bu, (mpz_ptr) b);100 mpz_init_set(bv, (mpz_ptr) a);101 102 assume(mpz_cmp_si(erg, 0));103 104 mpz_div(bu, bu, erg);105 mpz_div(bv, bv, erg);106 107 mpz_mul_si(bu, bu, -1);108 *u = (number) bu;109 *v = (number) bv;110 111 *s = (number) bs;112 *t = (number) bt;113 return (number) erg;114 }115 116 static void nrzPower (number a, int i, number * result, const coeffs)117 {118 mpz_ptr erg = (mpz_ptr) omAllocBin(gmp_nrz_bin);119 mpz_init(erg);120 mpz_pow_ui(erg, (mpz_ptr) a, i);121 *result = (number) erg;122 }123 124 /*125 * create a number from int126 */127 number nrzInit (long i, const coeffs)128 {129 mpz_ptr erg = (mpz_ptr) omAllocBin(gmp_nrz_bin);130 mpz_init_set_si(erg, i);131 return (number) erg;132 }133 134 void nrzDelete(number *a, const coeffs)135 {136 if (*a != NULL)137 {138 mpz_clear((mpz_ptr) *a);139 omFreeBin((ADDRESS) *a, gmp_nrz_bin);140 *a = NULL;141 }142 }143 144 static number nrzCopy(number a, const coeffs)145 {146 if (a==NULL) return NULL;147 mpz_ptr erg = (mpz_ptr) omAllocBin(gmp_nrz_bin);148 mpz_init_set(erg, (mpz_ptr) a);149 return (number) erg;150 }151 152 #if 0153 number nrzCopyMap(number a, const coeffs /*src*/, const coeffs dst)154 {155 return nrzCopy(a,dst);156 }157 #endif158 159 int nrzSize(number a, const coeffs)160 {161 mpz_ptr p=(mpz_ptr)a;162 int s=p->_mp_alloc;163 if (s==1) s=(mpz_cmp_ui(p,0)!=0);164 return s;165 166 }167 168 /*169 * convert a number to int170 */171 static long nrzInt(number &n, const coeffs)172 {173 return mpz_get_si( (mpz_ptr)n);174 }175 176 static number nrzAdd (number a, number b, const coeffs)177 {178 mpz_ptr erg = (mpz_ptr) omAllocBin(gmp_nrz_bin);179 mpz_init(erg);180 mpz_add(erg, (mpz_ptr) a, (mpz_ptr) b);181 return (number) erg;182 }183 184 static number nrzSub (number a, number b, const coeffs)185 {186 mpz_ptr erg = (mpz_ptr) omAllocBin(gmp_nrz_bin);187 mpz_init(erg);188 mpz_sub(erg, (mpz_ptr) a, (mpz_ptr) b);189 return (number) erg;190 }191 192 static number nrzGetUnit (number, const coeffs r)193 {194 return nrzInit(1, r);195 }196 197 static BOOLEAN nrzIsUnit (number a, const coeffs)198 {199 return 0 == mpz_cmpabs_ui((mpz_ptr) a, 1);200 }201 202 static BOOLEAN nrzIsZero (number a, const coeffs)203 {204 return 0 == mpz_cmpabs_ui((mpz_ptr) a, 0);205 }206 207 static BOOLEAN nrzIsOne (number a, const coeffs)208 {209 return (0 == mpz_cmp_ui((mpz_ptr) a, 1));210 }211 212 static BOOLEAN nrzIsMOne (number a, const coeffs)213 {214 return (0 == mpz_cmp_si((mpz_ptr) a, -1));215 }216 217 static BOOLEAN nrzEqual (number a,number b, const coeffs)218 {219 return 0 == mpz_cmp((mpz_ptr) a, (mpz_ptr) b);220 }221 222 static BOOLEAN nrzGreater (number a,number b, const coeffs)223 {224 return 0 < mpz_cmp((mpz_ptr) a, (mpz_ptr) b);225 }226 227 static BOOLEAN nrzGreaterZero (number k, const coeffs)228 {229 return 0 < mpz_sgn1((mpz_ptr) k);230 }231 232 static BOOLEAN nrzDivBy (number a,number b, const coeffs)233 {234 return mpz_divisible_p((mpz_ptr) a, (mpz_ptr) b) != 0;235 }236 237 static int nrzDivComp(number a, number b, const coeffs r)238 {239 if (nrzDivBy(a, b, r))240 {241 if (nrzDivBy(b, a, r)) return 2;242 return -1;243 }244 if (nrzDivBy(b, a, r)) return 1;245 return 0;246 }247 248 static number nrzDiv (number a,number b, const coeffs)249 {250 mpz_ptr erg = (mpz_ptr) omAllocBin(gmp_nrz_bin);251 mpz_init(erg);252 mpz_ptr r = (mpz_ptr) omAllocBin(gmp_nrz_bin);253 mpz_init(r);254 mpz_tdiv_qr(erg, r, (mpz_ptr) a, (mpz_ptr) b);255 //if (!nrzIsZero((number) r, R))256 //{257 // WerrorS("Division by non divisible element.");258 // WerrorS("Result is without remainder.");259 //}260 mpz_clear(r);261 omFreeBin(r, gmp_nrz_bin);262 return (number) erg;263 }264 265 static number nrzExactDiv (number a,number b, const coeffs)266 {267 mpz_ptr erg = (mpz_ptr) omAllocBin(gmp_nrz_bin);268 mpz_init(erg);269 mpz_tdiv_q(erg, (mpz_ptr) a, (mpz_ptr) b);270 return (number) erg;271 }272 273 static number nrzEucNorm (number a, const coeffs )274 {275 mpz_ptr abs = (mpz_ptr) omAllocBin(gmp_nrz_bin);276 mpz_init(abs);277 mpz_abs(abs, (mpz_ptr)a);278 279 return (number) abs;280 }281 282 static number nrzSmallestQuotRem (number a, number b, number * r, const coeffs )283 {284 mpz_ptr qq = (mpz_ptr) omAllocBin(gmp_nrz_bin);285 mpz_init(qq);286 mpz_ptr rr = (mpz_ptr) omAllocBin(gmp_nrz_bin);287 mpz_init(rr);288 int gsign = mpz_sgn((mpz_ptr) b);289 mpz_t gg, ghalf;290 mpz_init(gg);291 mpz_init(ghalf);292 mpz_abs(gg, (mpz_ptr) b);293 mpz_fdiv_qr(qq, rr, (mpz_ptr) a, gg);294 mpz_tdiv_q_2exp(ghalf, gg, 1);295 if (mpz_cmp(rr, ghalf) > 0) // r > ghalf296 {297 mpz_sub(rr, rr, gg);298 mpz_add_ui(qq, qq, 1);299 }300 if (gsign < 0) mpz_neg(qq, qq);301 302 mpz_clear(gg);303 mpz_clear(ghalf);304 if (r==NULL)305 {306 mpz_clear(rr);307 omFreeBin(rr,gmp_nrz_bin);308 }309 else310 {311 *r=(number)rr;312 }313 return (number) qq;314 }315 316 static number nrzQuotRem (number a, number b, number * r, const coeffs )317 {318 mpz_ptr qq = (mpz_ptr) omAllocBin(gmp_nrz_bin);319 mpz_init(qq);320 mpz_ptr rr = (mpz_ptr) omAllocBin(gmp_nrz_bin);321 mpz_init(rr);322 mpz_tdiv_qr(qq, rr, (mpz_ptr) a, (mpz_ptr) b);323 if (r==NULL)324 {325 mpz_clear(rr);326 omFreeBin(rr,gmp_nrz_bin);327 }328 else329 {330 *r=(number)rr;331 }332 return (number) qq;333 }334 335 static number nrzIntMod (number a,number b, const coeffs)336 {337 mpz_ptr erg = (mpz_ptr) omAllocBin(gmp_nrz_bin);338 mpz_init(erg);339 mpz_ptr r = (mpz_ptr) omAllocBin(gmp_nrz_bin);340 mpz_init(r);341 mpz_tdiv_qr(erg, r, (mpz_ptr) a, (mpz_ptr) b);342 mpz_clear(erg);343 omFreeBin(erg, gmp_nrz_bin);344 return (number) r;345 }346 347 static number nrzInvers (number c, const coeffs r)348 {349 if (!nrzIsUnit((number) c, r))350 {351 WerrorS("Non invertible element.");352 return (number)NULL;353 }354 return nrzCopy(c,r);355 }356 357 static number nrzNeg (number c, const coeffs)358 {359 // nNeg inplace !!!360 mpz_mul_si((mpz_ptr) c, (mpz_ptr) c, -1);361 return c;362 }363 364 static number nrzMapMachineInt(number from, const coeffs /*src*/, const coeffs /*dst*/)365 {366 mpz_ptr erg = (mpz_ptr) omAllocBin(gmp_nrz_bin);367 mpz_init_set_ui(erg, (unsigned long) from);368 return (number) erg;369 }370 371 static number nrzMapZp(number from, const coeffs /*src*/, const coeffs /*dst*/)372 {373 mpz_ptr erg = (mpz_ptr) omAllocBin(gmp_nrz_bin);374 mpz_init_set_si(erg, (long) from);375 return (number) erg;376 }377 378 static number nrzMapQ(number from, const coeffs src, const coeffs /*dst*/)379 {380 mpz_ptr erg = (mpz_ptr) omAllocBin(gmp_nrz_bin);381 mpz_init(erg);382 nlGMP(from, erg, src); // FIXME? TODO? // extern void nlGMP(number &i, number n, const coeffs r); // to be replaced with n_MPZ(erg, from, src); // ?383 return (number) erg;384 }385 386 static nMapFunc nrzSetMap(const coeffs src, const coeffs /*dst*/)387 {388 /* dst = currRing */389 /* dst = nrn */390 if ((src->rep==n_rep_gmp)391 && (nCoeff_is_Z(src) || nCoeff_is_Zn(src) || nCoeff_is_Ring_PtoM(src)))392 {393 return ndCopyMap; //nrzCopyMap;394 }395 if ((src->rep==n_rep_gap_gmp) /*&& nCoeff_is_Z(src)*/)396 {397 return ndCopyMap; //nrzCopyMap;398 }399 if (nCoeff_is_Ring_2toM(src))400 {401 return nrzMapMachineInt;402 }403 if (nCoeff_is_Zp(src))404 {405 return nrzMapZp;406 }407 if (getCoeffType(src)==n_Q /*nCoeff_is_Q(src) or coeffs_BIGINT*/)408 {409 return nrzMapQ;410 }411 return NULL; // default412 }413 414 /*415 * set the exponent (allocate and init tables) (TODO)416 */417 418 void nrzSetExp(int, coeffs)419 {420 }421 422 void nrzInitExp(int, coeffs)423 {424 }425 426 #ifdef LDEBUG427 static BOOLEAN nrzDBTest (number, const char *, const int, const coeffs)428 {429 return TRUE;//TODO430 }431 #endif432 433 void nrzWrite (number a, const coeffs)434 {435 char *s,*z;436 if (a==NULL)437 {438 StringAppendS("o");439 }440 else441 {442 int l=mpz_sizeinbase((mpz_ptr) a, 10) + 2;443 s=(char*)omAlloc(l);444 z=mpz_get_str(s,10,(mpz_ptr) a);445 StringAppendS(z);446 omFreeSize((ADDRESS)s,l);447 }448 }449 450 /*2451 * extracts a long integer from s, returns the rest (COPY FROM longrat0.cc)452 */453 static const char * nlEatLongC(char *s, mpz_ptr i)454 {455 const char * start=s;456 457 if (*s<'0' || *s>'9')458 {459 mpz_set_ui(i,1);460 return s;461 }462 while (*s >= '0' && *s <= '9') s++;463 if (*s=='\0')464 {465 mpz_set_str(i,start,10);466 }467 else468 {469 char c=*s;470 *s='\0';471 mpz_set_str(i,start,10);472 *s=c;473 }474 return s;475 }476 477 478 static CanonicalForm nrzConvSingNFactoryN(number n, BOOLEAN setChar, const coeffs /*r*/)479 {480 if (setChar) setCharacteristic( 0 );481 482 CanonicalForm term;483 mpz_t num;484 mpz_init_set(num, *((mpz_t*)n));485 term = make_cf(num);486 return term;487 }488 489 static number nrzConvFactoryNSingN(const CanonicalForm n, const coeffs r)490 {491 if (n.isImm())492 return nrzInit(n.intval(),r);493 else494 {495 mpz_ptr m = (mpz_ptr) omAllocBin(gmp_nrz_bin);496 gmp_numerator(n,m);497 if (!n.den().isOne())498 {499 WarnS("denominator is not 1 in factory");500 }501 return (number) m;502 }503 }504 505 static const char * nrzRead (const char *s, number *a, const coeffs)506 {507 mpz_ptr z = (mpz_ptr) omAllocBin(gmp_nrz_bin);508 {509 mpz_init(z);510 s = nlEatLongC((char *) s, z);511 }512 *a = (number) z;513 return s;514 }515 25 516 26 static void nrzCoeffWrite (const coeffs, BOOLEAN /*details*/) … … 529 39 } 530 40 531 static coeffs nrzQuot1(number c, const coeffs r) 532 { 533 long ch = r->cfInt(c, r); 534 mpz_t dummy; 535 mpz_init_set_ui(dummy, ch); 536 ZnmInfo info; 537 info.base = dummy; 538 info.exp = (unsigned long) 1; 539 coeffs rr = nInitChar(n_Zn, (void*)&info); 540 mpz_clear(dummy); 541 return(rr); 542 } 543 544 static number nrzInitMPZ(mpz_t m, const coeffs) 545 { 546 mpz_ptr z = (mpz_ptr) omAllocBin(gmp_nrz_bin); 547 mpz_init_set(z, m); 548 return (number)z; 549 } 550 551 static number nrzFarey(number r, number N, const coeffs R) 552 { 553 number a0 = nrzCopy(N, R); 554 number b0 = nrzInit(0, R); 555 number a1 = nrzCopy(r, R); 556 number b1 = nrzInit(1, R); 557 number two = nrzInit(2, R); 558 #if 0 559 PrintS("Farey start with "); 560 n_Print(r, R); 561 PrintS(" mod "); 562 n_Print(N, R); 563 PrintLn(); 564 #endif 565 while (1) 566 { 567 number as = nrzMult(a1, a1, R); 568 n_InpMult(as, two, R); 569 if (nrzGreater(N, as, R)) 570 { 571 nrzDelete(&as, R); 572 break; 573 } 574 nrzDelete(&as, R); 575 number q = nrzDiv(a0, a1, R); 576 number t = nrzMult(a1, q, R), 577 s = nrzSub(a0, t, R); 578 nrzDelete(&a0, R); 579 a0 = a1; 580 a1 = s; 581 nrzDelete(&t, R); 582 583 t = nrzMult(b1, q, R); 584 s = nrzSub(b0, t, R); 585 nrzDelete(&b0, R); 586 b0 = b1; 587 b1 = s; 588 nrzDelete(&t, R); 589 nrzDelete(&q, R); 590 } 591 number as = nrzMult(b1, b1, R); 592 n_InpMult(as, two, R); 593 nrzDelete(&two, R); 594 if (nrzGreater(as, N, R)) 595 { 596 nrzDelete(&a0, R); 597 nrzDelete(&a1, R); 598 nrzDelete(&b0, R); 599 nrzDelete(&b1, R); 600 nrzDelete(&as, R); 601 return NULL; 602 } 603 nrzDelete(&as, R); 604 nrzDelete(&a0, R); 605 nrzDelete(&b0, R); 606 607 number a, b, ab; 608 coeffs Q = nInitChar(n_Q, 0); 609 nMapFunc f = n_SetMap(R, Q); 610 a = f(a1, R, Q); 611 b = f(b1, R, Q); 612 ab = n_Div(a, b, Q); 613 n_Delete(&a, Q); 614 n_Delete(&b, Q); 615 nKillChar(Q); 616 617 nrzDelete(&a1, R); 618 nrzDelete(&b1, R); 619 return ab; 620 } 621 622 void nrzWriteFd(number n, const ssiInfo* d, const coeffs) 623 { 624 mpz_out_str (d->f_write,SSI_BASE, (mpz_ptr)n); 625 fputc(' ',d->f_write); 626 } 627 628 number nrzReadFd(const ssiInfo *d, const coeffs) 629 { 630 mpz_ptr erg = (mpz_ptr) omAllocBin(gmp_nrz_bin); 631 mpz_init(erg); 632 s_readmpz_base(d->f_read,erg,SSI_BASE); 633 return (number)erg; 634 } 635 636 BOOLEAN nrzInitChar(coeffs r, void *) 637 { 638 assume( getCoeffType(r) == n_Z ); 639 640 r->is_field=FALSE; 641 r->is_domain=TRUE; 642 r->rep=n_rep_gmp; 643 644 //r->nCoeffIsEqual = ndCoeffIsEqual; 645 r->cfCoeffString = nrzCoeffString; 646 r->cfCoeffName = nrzCoeffName; 647 r->cfCoeffWrite = nrzCoeffWrite; 648 //r->cfKillChar = ndKillChar; 649 r->cfMult = nrzMult; 650 r->cfSub = nrzSub; 651 r->cfAdd = nrzAdd; 652 r->cfDiv = nrzDiv; 653 r->cfIntMod= nrzIntMod; 654 r->cfExactDiv= nrzExactDiv; 655 r->cfInit = nrzInit; 656 r->cfInitMPZ = nrzInitMPZ; 657 r->cfSize = nrzSize; 658 r->cfInt = nrzInt; 659 r->cfDivComp = nrzDivComp; 660 r->cfIsUnit = nrzIsUnit; 661 r->cfGetUnit = nrzGetUnit; 662 r->cfExtGcd = nrzExtGcd; 663 r->cfXExtGcd = nrzXExtGcd; 664 r->cfDivBy = nrzDivBy; 665 r->cfEucNorm = nrzEucNorm; 666 r->cfQuotRem = nrzSmallestQuotRem; 667 r->cfInpNeg = nrzNeg; 668 r->cfInvers= nrzInvers; 669 r->cfCopy = nrzCopy; 670 r->cfWriteLong = nrzWrite; 671 r->cfRead = nrzRead; 672 r->cfGreater = nrzGreater; 673 r->cfEqual = nrzEqual; 674 r->cfIsZero = nrzIsZero; 675 r->cfIsOne = nrzIsOne; 676 r->cfIsMOne = nrzIsMOne; 677 r->cfGreaterZero = nrzGreaterZero; 678 r->cfPower = nrzPower; 679 r->cfGcd = nrzGcd; 680 r->cfLcm = nrzLcm; 681 r->cfDelete= nrzDelete; 682 r->cfSetMap = nrzSetMap; 683 r->cfQuot1 = nrzQuot1; 684 r->convSingNFactoryN=nrzConvSingNFactoryN; 685 r->convFactoryNSingN=nrzConvFactoryNSingN; 686 r->cfChineseRemainder=nlChineseRemainderSym; 687 r->cfFarey=nrzFarey; 688 r->cfWriteFd=nrzWriteFd; 689 r->cfReadFd=nrzReadFd; 690 // debug stuff 691 692 #ifdef LDEBUG 693 r->cfDBTest=nrzDBTest; 694 #endif 695 696 r->ch = 0; 697 r->has_simple_Alloc=FALSE; 698 r->has_simple_Inverse=FALSE; 699 return FALSE; 700 } 701 41 #if SI_INTEGER_VARIANT == 2 42 #include "rintegers2.cc" 702 43 #elif SI_INTEGER_VARIANT == 3 703 704 //make sure that a small number is an immediate integer 705 //bascially coped from longrat.cc nlShort3 706 //TODO: is there any point in checking 0 first??? 707 //TODO: it is not clear that this works in 32/64 bit everywhere. 708 // too many hacks. 709 #ifdef LDEBUG 710 #define nrzTest(A) nrzDBTest(A,__FILE__,__LINE__,NULL) 711 BOOLEAN nrzDBTest (number x, const char *f, const int l, const coeffs); 712 #else 713 #define nrzTest(A) 714 #endif 715 716 #undef CF_DEBUG 717 static inline number nrz_short(number x) 718 { 719 #if CF_DEBUG 720 StringAppendS("short("); 721 nrzWrite(x, NULL); 722 #endif 723 if (mpz_sgn1((mpz_ptr) x)==0) 724 { 725 mpz_clear((mpz_ptr)x); 726 omFreeBin(x, gmp_nrz_bin); 727 #if CF_DEBUG 728 StringAppendS(")=0"); 729 #endif 730 return INT_TO_SR(0); 731 } 732 if (mpz_size1((mpz_ptr)x)<=MP_SMALL) 733 { 734 long ui=mpz_get_si((mpz_ptr)x); 735 if ((((ui<<3)>>3)==ui) 736 && (mpz_cmp_si((mpz_ptr)x,ui)==0)) 737 { 738 mpz_clear((mpz_ptr)x); 739 omFreeBin(x, gmp_nrz_bin); 740 #if CF_DEBUG 741 StringAppendS(")=imm"); 742 #endif 743 return INT_TO_SR(ui); 744 } 745 } 746 #if CF_DEBUG 747 StringAppendS(")"); 748 #endif 749 return x; 750 } 751 752 753 static int nrzSize(number a, const coeffs) 754 { 755 if (a==INT_TO_SR(0)) return 0; 756 if (n_Z_IS_SMALL(a)) return 1; 757 return ((mpz_ptr)a)->_mp_alloc; 758 } 759 760 761 /* 762 * Multiply two numbers 763 * check for 0, 1, -1 maybe 764 */ 765 #if CF_DEBUG 766 number _nrzMult(number, number, const coeffs); 767 number nrzMult(number a, number b, const coeffs R) 768 { 769 StringSetS("Mult: "); 770 nrzWrite(a, R); 771 StringAppendS(" by "); 772 nrzWrite(b, R); 773 number c = _nrzMult(a, b, R); 774 StringAppendS(" is "); 775 nrzWrite(c, R); 776 char * s = StringEndS(); 777 Print("%s\n", s); 778 omFree(s); 779 return c; 780 } 781 number _nrzMult (number a, number b, const coeffs R) 782 #else 783 number nrzMult (number a, number b, const coeffs R) 784 #endif 785 { 786 if (n_Z_IS_SMALL(a) && n_Z_IS_SMALL(b)) 787 { 788 //from longrat.cc 789 if (SR_TO_INT(a)==0) 790 return a; 791 if (SR_TO_INT(b)==0) 792 return b; 793 long r=(long)((unsigned long)(SR_HDL(a)-1L))*((unsigned long)(SR_HDL(b)>>1)); 794 if ((r/(SR_HDL(b)>>1))==(SR_HDL(a)-1L)) 795 { 796 number u=((number) ((r>>1)+SR_INT)); 797 // if (((((long)SR_HDL(u))<<1)>>1)==SR_HDL(u)) return (u); 798 return nrzInit(SR_HDL(u)>>2, R); 799 } 800 mpz_ptr erg = (mpz_ptr) omAllocBin(gmp_nrz_bin); 801 mpz_init(erg); 802 mpz_set_si(erg, SR_TO_INT(a)); 803 mpz_mul_si(erg, erg, SR_TO_INT(b)); 804 nrzTest((number)erg); 805 return (number) erg; 806 } 807 else if (n_Z_IS_SMALL(a)) 808 { 809 if (SR_TO_INT(a)==0) 810 return a; 811 mpz_ptr erg = (mpz_ptr) omAllocBin(gmp_nrz_bin); 812 mpz_init_set(erg, (mpz_ptr) b); 813 mpz_mul_si(erg, erg, SR_TO_INT(a)); 814 nrzTest((number)erg); 815 return (number) erg; 816 } 817 else if (n_Z_IS_SMALL(b)) 818 { 819 if (SR_TO_INT(b)==0) 820 return b; 821 mpz_ptr erg = (mpz_ptr) omAllocBin(gmp_nrz_bin); 822 mpz_init_set(erg, (mpz_ptr) a); 823 mpz_mul_si(erg, erg, SR_TO_INT(b)); 824 nrzTest((number)erg); 825 return (number) erg; 826 } 827 else 828 { 829 mpz_ptr erg = (mpz_ptr) omAllocBin(gmp_nrz_bin); 830 mpz_init(erg); 831 mpz_mul(erg, (mpz_ptr) a, (mpz_ptr) b); 832 nrzTest((number)erg); 833 return (number) erg; 834 } 835 } 836 837 838 static long int_gcd(long a, long b) 839 { 840 long r; 841 a = ABS(a); 842 b = ABS(b); 843 if (!a) return b; 844 if (!b) return a; 845 do 846 { 847 r = a % b; 848 a = b; 849 b = r; 850 } while (b); 851 return ABS(a); // % in c doeas not imply a signn 852 // it would be unlikely to see a negative here 853 // but who knows 854 } 855 856 /* 857 * Give the smallest non unit k, such that a * x = k = b * y has a solution 858 */ 859 static number nrzLcm (number a, number b, const coeffs R) 860 { 861 #ifdef CF_DEBUG 862 PrintS("nrzLcm\n"); 863 #endif 864 mpz_ptr erg; 865 if (n_Z_IS_SMALL(a) && n_Z_IS_SMALL(b)) 866 { 867 long g = int_gcd(SR_TO_INT(a), SR_TO_INT(b)); 868 return nrzMult(a, INT_TO_SR(SR_TO_INT(b)/g), R); 869 } 870 else 871 { 872 erg = (mpz_ptr) omAllocBin(gmp_nrz_bin); 873 if (n_Z_IS_SMALL(a)) 874 { 875 mpz_init_set(erg, (mpz_ptr) b); 876 unsigned long g = mpz_gcd_ui(NULL, erg, (unsigned long) ABS(SR_TO_INT(a))); 877 mpz_mul_si(erg, erg, SR_TO_INT(a)/g); 878 } 879 else if (n_Z_IS_SMALL(b)) 880 { 881 mpz_init_set(erg, (mpz_ptr) a); 882 unsigned long g = mpz_gcd_ui(NULL, erg, (unsigned long) ABS(SR_TO_INT(b))); 883 mpz_mul_si(erg, erg, SR_TO_INT(b)/g); 884 } 885 else 886 { 887 mpz_init(erg); 888 mpz_lcm(erg, (mpz_ptr) a, (mpz_ptr) b); 889 } 890 } 891 return (number) erg; 892 } 893 894 static number nrzCopy(number a, const coeffs) 895 { 896 if (n_Z_IS_SMALL(a)) return a; 897 mpz_ptr erg = (mpz_ptr) omAllocBin(gmp_nrz_bin); 898 mpz_init_set(erg, (mpz_ptr) a); 899 return (number) erg; 900 } 901 902 /* 903 * Give the largest non unit k, such that a = x * k, b = y * k has 904 * a solution. 905 */ 906 static number nrzGcd (number a,number b,const coeffs R) 907 { 908 if (n_Z_IS_SMALL(a) && n_Z_IS_SMALL(b)) 909 { 910 long g = int_gcd(SR_TO_INT(a), SR_TO_INT(b)); 911 return INT_TO_SR(g); 912 } 913 else if (n_Z_IS_SMALL(a)) 914 { 915 if (a==INT_TO_SR(0)) 916 return nrzCopy(b, R); 917 unsigned long g = mpz_gcd_ui(NULL, (mpz_ptr)b, (unsigned long) ABS(SR_TO_INT(a))); 918 return INT_TO_SR( g); 919 } 920 else if (n_Z_IS_SMALL(b)) 921 { 922 if (b==INT_TO_SR(0)) 923 return nrzCopy(a, R); 924 unsigned long g = mpz_gcd_ui(NULL, (mpz_ptr)a, (unsigned long) ABS(SR_TO_INT(b))); 925 return INT_TO_SR(g); 926 } 927 else 928 { 929 mpz_ptr erg = (mpz_ptr) omAllocBin(gmp_nrz_bin); 930 mpz_init(erg); 931 mpz_gcd(erg, (mpz_ptr) a, (mpz_ptr) b); 932 return (number) erg; 933 } 934 } 935 936 /* 937 * Give the largest non unit k, such that a = x * k, b = y * k has 938 * a solution and r, s, s.t. k = s*a + t*b 939 */ 940 static long int_extgcd(long a, long b, long * u, long* x, long * v, long* y) 941 { 942 long q, r; 943 if (!a) 944 { 945 *u = 0; 946 *v = 1; 947 *x = -1; 948 *y = 0; 949 return b; 950 } 951 if (!b) 952 { 953 *u = 1; 954 *v = 0; 955 *x = 0; 956 *y = 1; 957 return a; 958 } 959 *u=1; 960 *v=0; 961 *x=0; 962 *y=1; 963 do 964 { 965 q = a/b; 966 r = a%b; 967 assume (q*b+r == a); 968 a = b; 969 b = r; 970 971 r = -(*v)*q+(*u); 972 (*u) =(*v); 973 (*v) = r; 974 975 r = -(*y)*q+(*x); 976 (*x) = (*y); 977 (*y) = r; 978 } while (b); 979 980 return a; 981 } 982 983 static number nrzExtGcd (number a, number b, number *s, number *t, const coeffs) 984 { 985 if (n_Z_IS_SMALL(a) && n_Z_IS_SMALL(b)) 986 { 987 long u, v, x, y; 988 long g = int_extgcd(SR_TO_INT(a), SR_TO_INT(b), &u, &v, &x, &y); 989 *s = INT_TO_SR(u); 990 *t = INT_TO_SR(v); 991 return INT_TO_SR(g); 992 } 993 else 994 { 995 mpz_t aa, bb; 996 if (n_Z_IS_SMALL(a)) 997 { 998 mpz_init_set_si(aa, SR_TO_INT(a)); 999 } 1000 else 1001 { 1002 mpz_init_set(aa, (mpz_ptr) a); 1003 } 1004 if (n_Z_IS_SMALL(b)) 1005 { 1006 mpz_init_set_si(bb, SR_TO_INT(b)); 1007 } 1008 else 1009 { 1010 mpz_init_set(bb, (mpz_ptr) b); 1011 } 1012 mpz_ptr erg = (mpz_ptr) omAllocBin(gmp_nrz_bin); 1013 mpz_ptr bs = (mpz_ptr) omAllocBin(gmp_nrz_bin); 1014 mpz_ptr bt = (mpz_ptr) omAllocBin(gmp_nrz_bin); 1015 mpz_init(erg); 1016 mpz_init(bs); 1017 mpz_init(bt); 1018 mpz_gcdext(erg, bs, bt, aa, bb); 1019 *s = nrz_short((number) bs); 1020 *t = nrz_short((number) bt); 1021 mpz_clear(aa); 1022 mpz_clear(bb); 1023 return nrz_short((number) erg); 1024 } 1025 } 1026 #if CF_DEBUG 1027 static number _nrzXExtGcd(number, number, number *, number *, number *, number *, const coeffs); 1028 static number nrzXExtGcd(number a, number b, number *x, number * y, number * u, number * v, const coeffs R) 1029 { 1030 char * s; 1031 StringSetS("XExtGcd: "); 1032 nrzWrite(a, R); 1033 StringAppendS(" by "); 1034 nrzWrite(b, R); 1035 number c = _nrzXExtGcd(a, b, x, y, u, v, R); 1036 StringAppendS(" is "); 1037 nrzWrite(c, R); 1038 StringAppendS("[["); 1039 nrzWrite(*x, R); 1040 StringAppendS(", "); 1041 nrzWrite(*y, R); 1042 StringAppendS("], "); 1043 nrzWrite(*u, R); 1044 StringAppendS(", "); 1045 nrzWrite(*v, R); 1046 s=StringEndS(); 1047 Print("%s]]\n", s); 1048 omFree(s); 1049 return c; 1050 } 1051 static number _nrzXExtGcd (number a, number b, number *s, number *t, number *u, number *v, const coeffs ) 1052 #else 1053 static number nrzXExtGcd (number a, number b, number *s, number *t, number *u, number *v, const coeffs ) 1054 #endif 1055 { 1056 if (n_Z_IS_SMALL(a) && n_Z_IS_SMALL(b)) 1057 { 1058 long uu, vv, x, y; 1059 long g = int_extgcd(SR_TO_INT(a), SR_TO_INT(b), &uu, &vv, &x, &y); 1060 *s = INT_TO_SR(uu); 1061 *t = INT_TO_SR(vv); 1062 *u = INT_TO_SR(x); 1063 *v = INT_TO_SR(y); 1064 return INT_TO_SR(g); 1065 } 1066 else 1067 { 1068 mpz_t aa, bb; 1069 if (n_Z_IS_SMALL(a)) 1070 { 1071 mpz_init_set_si(aa, SR_TO_INT(a)); 1072 } 1073 else 1074 { 1075 mpz_init_set(aa, (mpz_ptr) a); 1076 } 1077 if (n_Z_IS_SMALL(b)) 1078 { 1079 mpz_init_set_si(bb, SR_TO_INT(b)); 1080 } 1081 else 1082 { 1083 mpz_init_set(bb, (mpz_ptr) b); 1084 } 1085 mpz_ptr erg = (mpz_ptr) omAllocBin(gmp_nrz_bin); 1086 mpz_ptr bs = (mpz_ptr) omAllocBin(gmp_nrz_bin); 1087 mpz_ptr bt = (mpz_ptr) omAllocBin(gmp_nrz_bin); 1088 mpz_init(erg); 1089 mpz_init(bs); 1090 mpz_init(bt); 1091 1092 mpz_gcdext(erg, bs, bt, aa, bb); 1093 1094 mpz_ptr bu = (mpz_ptr) omAllocBin(gmp_nrz_bin); 1095 mpz_ptr bv = (mpz_ptr) omAllocBin(gmp_nrz_bin); 1096 1097 mpz_init_set(bu, (mpz_ptr) bb); 1098 mpz_init_set(bv, (mpz_ptr) aa); 1099 1100 mpz_clear(aa); 1101 mpz_clear(bb); 1102 assume(mpz_cmp_si(erg, 0)); 1103 1104 mpz_div(bu, bu, erg); 1105 mpz_div(bv, bv, erg); 1106 1107 mpz_mul_si(bu, bu, -1); 1108 *u = nrz_short((number) bu); 1109 *v = nrz_short((number) bv); 1110 1111 *s = nrz_short((number) bs); 1112 *t = nrz_short((number) bt); 1113 return nrz_short((number) erg); 1114 } 1115 } 1116 #if CF_DEBUG 1117 static number _nrzQuotRem(number, number, number *, const coeffs); 1118 static number nrzQuotRem(number a, number b, number * r, const coeffs R) 1119 { 1120 StringSetS("QuotRem: "); 1121 nrzWrite(a, R); 1122 StringAppendS(" by "); 1123 nrzWrite(b, R); 1124 number c = _nrzQuotRem(a, b, r, R); 1125 StringAppendS(" is "); 1126 nrzWrite(c, R); 1127 if (r) { 1128 StringAppendS("+R("); 1129 nrzWrite(*r, R); 1130 StringAppendS(")"); 1131 } 1132 char * s = StringEndS(); 1133 Print("%s\n", s); 1134 omFree(s); 1135 return c; 1136 } 1137 static number _nrzQuotRem (number a, number b, number * r, const coeffs ) 1138 #else 1139 static number nrzQuotRem (number a, number b, number * r, const coeffs ) 1140 #endif 1141 { 1142 assume(SR_TO_INT(b)); 1143 if (n_Z_IS_SMALL(a) && n_Z_IS_SMALL(b)) 1144 { 1145 if (r) 1146 *r = INT_TO_SR(SR_TO_INT(a) % SR_TO_INT(b)); 1147 return INT_TO_SR(SR_TO_INT(a)/SR_TO_INT(b)); 1148 } 1149 else if (n_Z_IS_SMALL(a)) 1150 { 1151 //a is small, b is not, so q=0, r=a 1152 if (r) 1153 *r = a; 1154 return INT_TO_SR(0); 1155 } 1156 else if (n_Z_IS_SMALL(b)) 1157 { 1158 unsigned long rr; 1159 mpz_ptr qq = (mpz_ptr) omAllocBin(gmp_nrz_bin); 1160 mpz_init(qq); 1161 mpz_t rrr; 1162 mpz_init(rrr); 1163 rr = mpz_divmod_ui(qq, rrr, (mpz_ptr) a, (unsigned long)ABS(SR_TO_INT(b))); 1164 mpz_clear(rrr); 1165 1166 if (r) 1167 *r = INT_TO_SR(rr); 1168 if (SR_TO_INT(b)<0) 1169 { 1170 mpz_mul_si(qq, qq, -1); 1171 } 1172 return nrz_short((number)qq); 1173 } 1174 mpz_ptr qq = (mpz_ptr) omAllocBin(gmp_nrz_bin), 1175 rr = (mpz_ptr) omAllocBin(gmp_nrz_bin); 1176 mpz_init(qq); 1177 mpz_init(rr); 1178 mpz_divmod(qq, rr, (mpz_ptr)a, (mpz_ptr)b); 1179 if (r) 1180 *r = (number) rr; 1181 else 1182 { 1183 mpz_clear(rr); 1184 } 1185 nrzTest((number)qq); 1186 return (number) qq; 1187 } 1188 1189 static void nrzPower (number a, int i, number * result, const coeffs) 1190 { 1191 mpz_ptr erg = (mpz_ptr) omAllocBin(gmp_nrz_bin); 1192 mpz_init(erg); 1193 mpz_t aa; 1194 if (n_Z_IS_SMALL(a)) 1195 mpz_init_set_si(aa, SR_TO_INT(a)); 1196 else 1197 mpz_init_set(aa, (mpz_ptr) a); 1198 mpz_pow_ui(erg, aa, i); 1199 *result = nrz_short((number) erg); 1200 } 1201 1202 /* 1203 * create a number from int 1204 * TODO: do not create an mpz if not necessary 1205 */ 1206 number nrzInit (long i, const coeffs) 1207 { 1208 mpz_ptr erg = (mpz_ptr) omAllocBin(gmp_nrz_bin); 1209 mpz_init_set_si(erg, i); 1210 return nrz_short((number) erg); 1211 } 1212 1213 static number nrzInitMPZ(mpz_t m, const coeffs) 1214 { 1215 mpz_ptr erg = (mpz_ptr) omAllocBin(gmp_nrz_bin); 1216 mpz_init_set(erg, m); 1217 return nrz_short((number) erg); 1218 } 1219 1220 1221 static void nrzDelete(number *a, const coeffs) 1222 { 1223 if (*a == NULL) return; 1224 if (n_Z_IS_SMALL(*a)==0) 1225 { 1226 mpz_clear((mpz_ptr) *a); 1227 omFreeBin((ADDRESS) *a, gmp_nrz_bin); 1228 } 1229 *a = NULL; 1230 } 1231 1232 /* 1233 * convert a number to int 1234 */ 1235 static long nrzInt(number &n, const coeffs) 1236 { 1237 if (n_Z_IS_SMALL(n)) return SR_TO_INT(n); 1238 return mpz_get_si( (mpz_ptr)n); 1239 } 1240 #if CF_DEBUG 1241 static number _nrzAdd(number, number, const coeffs); 1242 static number nrzAdd(number a, number b, const coeffs R) 1243 { 1244 StringSetS("Add: "); 1245 nrzWrite(a, R); 1246 StringAppendS(" to "); 1247 nrzWrite(b, R); 1248 number c = _nrzAdd(a, b, R); 1249 StringAppendS(" is "); 1250 nrzWrite(c, R); 1251 char * s = StringEndS(); 1252 Print("%s\n", s); 1253 omFree(s); 1254 return c; 1255 } 1256 static number _nrzAdd (number a, number b, const coeffs ) 1257 #else 1258 static number nrzAdd (number a, number b, const coeffs ) 1259 #endif 1260 { 1261 if (n_Z_IS_SMALL(a) && n_Z_IS_SMALL(b)) 1262 { 1263 long c = SR_TO_INT(a) + SR_TO_INT(b); 1264 if (INT_IS_SMALL(c)) 1265 return INT_TO_SR(c); 1266 mpz_ptr erg = (mpz_ptr) omAllocBin(gmp_nrz_bin); 1267 mpz_init_set_si(erg, c); 1268 1269 nrzTest((number)erg); 1270 return (number) erg; 1271 } 1272 else if (n_Z_IS_SMALL(a)) 1273 { 1274 mpz_ptr erg = (mpz_ptr) omAllocBin(gmp_nrz_bin); 1275 mpz_init(erg); 1276 if (SR_TO_INT(a)>0) 1277 mpz_add_ui(erg, (mpz_ptr) b, (unsigned long)SR_TO_INT(a)); 1278 else 1279 mpz_sub_ui(erg, (mpz_ptr) b, (unsigned long)-(SR_TO_INT(a))); 1280 return nrz_short((number) erg); 1281 } 1282 else if (n_Z_IS_SMALL(b)) 1283 { 1284 mpz_ptr erg = (mpz_ptr) omAllocBin(gmp_nrz_bin); 1285 mpz_init(erg); 1286 if (SR_TO_INT(b)>0) 1287 mpz_add_ui(erg, (mpz_ptr) a, (unsigned long)SR_TO_INT(b)); 1288 else 1289 mpz_sub_ui(erg, (mpz_ptr) a, (unsigned long)-(SR_TO_INT(b))); 1290 return nrz_short((number) erg); 1291 } 1292 else 1293 { 1294 mpz_ptr erg = (mpz_ptr) omAllocBin(gmp_nrz_bin); 1295 mpz_init(erg); 1296 mpz_add(erg, (mpz_ptr) a, (mpz_ptr) b); 1297 return nrz_short((number) erg); 1298 } 1299 } 1300 1301 static number nrzSub (number a, number b, const coeffs ) 1302 { 1303 if (n_Z_IS_SMALL(a) && n_Z_IS_SMALL(b)) 1304 { 1305 long c = SR_TO_INT(a) - SR_TO_INT(b); 1306 if (INT_IS_SMALL(c)) 1307 return INT_TO_SR(c); 1308 mpz_ptr erg = (mpz_ptr) omAllocBin(gmp_nrz_bin); 1309 mpz_init_set_si(erg, c); 1310 nrzTest((number)erg); 1311 return (number) erg; 1312 } 1313 else if (n_Z_IS_SMALL(a)) 1314 { 1315 mpz_ptr erg = (mpz_ptr) omAllocBin(gmp_nrz_bin); 1316 mpz_init(erg); 1317 1318 if (SR_TO_INT(a)>0) 1319 mpz_ui_sub(erg, (unsigned long)SR_TO_INT(a), (mpz_ptr) b); 1320 else 1321 { 1322 mpz_add_ui(erg, (mpz_ptr) b, (unsigned long)-SR_TO_INT(a)); 1323 mpz_neg(erg, erg); 1324 } 1325 return nrz_short((number) erg); 1326 } 1327 else if (n_Z_IS_SMALL(b)) 1328 { 1329 mpz_ptr erg = (mpz_ptr) omAllocBin(gmp_nrz_bin); 1330 mpz_init(erg); 1331 if (SR_TO_INT(b)>0) 1332 mpz_sub_ui(erg, (mpz_ptr) a, (unsigned long)SR_TO_INT(b)); 1333 else 1334 mpz_add_ui(erg, (mpz_ptr) a, (unsigned long)-SR_TO_INT(b)); 1335 return nrz_short((number) erg); 1336 } 1337 else 1338 { 1339 mpz_ptr erg = (mpz_ptr) omAllocBin(gmp_nrz_bin); 1340 mpz_init(erg); 1341 mpz_sub(erg, (mpz_ptr) a, (mpz_ptr) b); 1342 return nrz_short((number) erg); 1343 } 1344 } 1345 1346 static BOOLEAN nrzGreater (number a,number b, const coeffs) 1347 { 1348 if (n_Z_IS_SMALL(a) && n_Z_IS_SMALL(b)) 1349 return ((long)a)>((long)b); 1350 else if (n_Z_IS_SMALL(a)) 1351 return 0 > mpz_cmp_si((mpz_ptr)b,SR_TO_INT(a)); 1352 else if (n_Z_IS_SMALL(b)) 1353 return 0 < mpz_cmp_si((mpz_ptr)a,SR_TO_INT(b)); 1354 return 0 < mpz_cmp((mpz_ptr) a, (mpz_ptr) b); 1355 } 1356 1357 static BOOLEAN nrzGreaterZero (number k, const coeffs C) 1358 { 1359 return nrzGreater(k, INT_TO_SR(0), C); 1360 } 1361 1362 static number nrzGetUnit (number n, const coeffs r) 1363 { 1364 if (nrzGreaterZero(n, r)) 1365 return INT_TO_SR(1); 1366 /*else*/ 1367 return INT_TO_SR(-1); 1368 } 1369 1370 static number nrzAnn(number n, const coeffs) 1371 { 1372 if (SR_TO_INT(n)) // in Z: the annihilator of !=0 is 0 1373 return INT_TO_SR(0); 1374 else 1375 return INT_TO_SR(1); 1376 } 1377 1378 static BOOLEAN nrzIsUnit (number a, const coeffs) 1379 { 1380 return ABS(SR_TO_INT(a))==1; 1381 } 1382 1383 static BOOLEAN nrzIsZero (number a, const coeffs) 1384 { 1385 return (a==INT_TO_SR(0)); 1386 } 1387 1388 static BOOLEAN nrzIsOne (number a, const coeffs) 1389 { 1390 return a==INT_TO_SR(1); 1391 } 1392 1393 static BOOLEAN nrzIsMOne (number a, const coeffs) 1394 { 1395 return a==INT_TO_SR(-1); 1396 } 1397 1398 static BOOLEAN nrzEqual (number a,number b, const coeffs) 1399 { 1400 if (n_Z_IS_SMALL(a) && n_Z_IS_SMALL(b)) 1401 return a==b; 1402 else if (n_Z_IS_SMALL(a) || n_Z_IS_SMALL(b)) 1403 return FALSE; 1404 else 1405 return 0 == mpz_cmp((mpz_ptr) a, (mpz_ptr) b); 1406 } 1407 1408 static BOOLEAN nrzDivBy (number a,number b, const coeffs) 1409 { 1410 if (n_Z_IS_SMALL(a) && n_Z_IS_SMALL(b)) 1411 { 1412 return SR_TO_INT(a) %SR_TO_INT(b) ==0; 1413 } 1414 else if (n_Z_IS_SMALL(a)) 1415 { 1416 return a==INT_TO_SR(0); 1417 } 1418 else if (n_Z_IS_SMALL(b)) 1419 { 1420 return mpz_divisible_ui_p((mpz_ptr)a, (unsigned long)ABS(SR_TO_INT(b))) != 0; 1421 } 1422 else 1423 return mpz_divisible_p((mpz_ptr) a, (mpz_ptr) b) != 0; 1424 } 1425 1426 static int nrzDivComp(number a, number b, const coeffs r) 1427 { 1428 if (nrzDivBy(a, b, r)) 1429 { 1430 if (nrzDivBy(b, a, r)) return 2; 1431 return -1; 1432 } 1433 if (nrzDivBy(b, a, r)) return 1; 1434 return 0; 1435 } 1436 1437 static number nrzDiv (number a,number b, const coeffs) 1438 { 1439 assume(SR_TO_INT(b)); 1440 if (n_Z_IS_SMALL(a) && n_Z_IS_SMALL(b)) 1441 { 1442 //if (SR_TO_INT(a) % SR_TO_INT(b)) 1443 //{ 1444 // WerrorS("1:Division by non divisible element."); 1445 // WerrorS("Result is without remainder."); 1446 //} 1447 return INT_TO_SR(SR_TO_INT(a)/SR_TO_INT(b)); 1448 } 1449 else if (n_Z_IS_SMALL(a)) 1450 { 1451 //if (SR_TO_INT(a)) 1452 //{ 1453 // WerrorS("2:Division by non divisible element."); 1454 // WerrorS("Result is without remainder."); 1455 //} 1456 return INT_TO_SR(0); 1457 } 1458 else if (n_Z_IS_SMALL(b)) 1459 { 1460 mpz_ptr erg = (mpz_ptr) omAllocBin(gmp_nrz_bin); 1461 mpz_t r; 1462 mpz_init(r); 1463 mpz_init(erg); 1464 if (mpz_divmod_ui(erg, r, (mpz_ptr) a, (unsigned long)ABS(SR_TO_INT(b)))) { 1465 // WerrorS("3:Division by non divisible element."); 1466 // WerrorS("Result is without remainder."); 1467 } 1468 mpz_clear(r); 1469 if (SR_TO_INT(b)<0) 1470 mpz_neg(erg, erg); 1471 return nrz_short((number) erg); 1472 } 1473 mpz_ptr erg = (mpz_ptr) omAllocBin(gmp_nrz_bin); 1474 mpz_init(erg); 1475 mpz_t r; 1476 mpz_init(r); 1477 mpz_tdiv_qr(erg, r, (mpz_ptr) a, (mpz_ptr) b); 1478 #if CF_DEBUG 1479 StringSetS("division of"); 1480 nrzWrite(a, R); 1481 StringAppendS(" by "); 1482 nrzWrite(b, R); 1483 StringAppendS(" is "); 1484 number du; 1485 nrzWrite(du = (number)erg, R); 1486 StringAppendS(" rest "); 1487 nrzWrite(du = (number)r, R); 1488 char * s = StringEndS(); 1489 Print("%s\n", s); 1490 omFree(s); 1491 #endif 1492 1493 if (mpz_sgn1(r)!=0) 1494 { 1495 //WerrorS("4:Division by non divisible element."); 1496 //WerrorS("Result is without remainder."); 1497 } 1498 mpz_clear(r); 1499 return nrz_short((number) erg); 1500 } 1501 1502 static number nrzExactDiv (number a,number b, const coeffs) 1503 { 1504 assume(SR_TO_INT(b)); 1505 mpz_t aa, bb; 1506 if (n_Z_IS_SMALL(a)) 1507 mpz_init_set_si(aa, SR_TO_INT(a)); 1508 else 1509 mpz_init_set(aa, (mpz_ptr) a); 1510 if (n_Z_IS_SMALL(b)) 1511 mpz_init_set_si(bb, SR_TO_INT(b)); 1512 else 1513 mpz_init_set(bb, (mpz_ptr) b); 1514 mpz_ptr erg = (mpz_ptr) omAllocBin(gmp_nrz_bin); 1515 mpz_init(erg); 1516 mpz_tdiv_q(erg, (mpz_ptr) aa, (mpz_ptr) bb); 1517 mpz_clear(aa); 1518 mpz_clear(bb); 1519 nrzTest((number)erg); 1520 return (number) erg; 1521 } 1522 1523 static number nrzIntMod (number a,number b, const coeffs) 1524 { 1525 mpz_t aa, bb; 1526 assume(SR_TO_INT(b)); 1527 if (n_Z_IS_SMALL(a)) 1528 mpz_init_set_si(aa, SR_TO_INT(a)); 1529 else 1530 mpz_init_set(aa, (mpz_ptr) a); 1531 if (n_Z_IS_SMALL(b)) 1532 mpz_init_set_si(bb, SR_TO_INT(b)); 1533 else 1534 mpz_init_set(bb, (mpz_ptr) b); 1535 1536 mpz_t erg; 1537 mpz_init(erg); 1538 mpz_ptr r = (mpz_ptr) omAllocBin(gmp_nrz_bin); 1539 mpz_init(r); 1540 mpz_tdiv_qr(erg, r, (mpz_ptr) aa, (mpz_ptr) bb); 1541 mpz_clear(erg); 1542 mpz_clear(aa); 1543 mpz_clear(bb); 1544 1545 return nrz_short((number) r); 1546 } 1547 1548 static number nrzInvers (number c, const coeffs r) 1549 { 1550 if (!nrzIsUnit((number) c, r)) 1551 { 1552 WerrorS("Non invertible element."); 1553 return (number)NULL; 1554 } 1555 return c; // has to be 1 or -1.... 1556 } 1557 1558 static number nrzNeg (number c, const coeffs) 1559 { 1560 // nNeg inplace !!! 1561 if (n_Z_IS_SMALL(c)) 1562 return INT_TO_SR(-SR_TO_INT(c)); 1563 mpz_mul_si((mpz_ptr) c, (mpz_ptr) c, -1); 1564 return c; 1565 } 1566 1567 static number nrzFarey(number r, number N, const coeffs R) 1568 { 1569 number a0 = nrzCopy(N, R); 1570 number b0 = nrzInit(0, R); 1571 number a1 = nrzCopy(r, R); 1572 number b1 = nrzInit(1, R); 1573 number two = nrzInit(2, R); 1574 #if 0 1575 PrintS("Farey start with "); 1576 n_Print(r, R); 1577 PrintS(" mod "); 1578 n_Print(N, R); 1579 PrintLn(); 1580 #endif 1581 while (1) 1582 { 1583 number as = nrzMult(a1, a1, R); 1584 n_InpMult(as, two, R); 1585 if (nrzGreater(N, as, R)) 1586 { 1587 nrzDelete(&as, R); 1588 break; 1589 } 1590 nrzDelete(&as, R); 1591 number q = nrzDiv(a0, a1, R); 1592 number t = nrzMult(a1, q, R), 1593 s = nrzSub(a0, t, R); 1594 nrzDelete(&a0, R); 1595 a0 = a1; 1596 a1 = s; 1597 nrzDelete(&t, R); 1598 1599 t = nrzMult(b1, q, R); 1600 s = nrzSub(b0, t, R); 1601 nrzDelete(&b0, R); 1602 b0 = b1; 1603 b1 = s; 1604 nrzDelete(&t, R); 1605 nrzDelete(&q, R); 1606 } 1607 number as = nrzMult(b1, b1, R); 1608 n_InpMult(as, two, R); 1609 nrzDelete(&two, R); 1610 if (nrzGreater(as, N, R)) 1611 { 1612 nrzDelete(&a0, R); 1613 nrzDelete(&a1, R); 1614 nrzDelete(&b0, R); 1615 nrzDelete(&b1, R); 1616 nrzDelete(&as, R); 1617 return NULL; 1618 } 1619 nrzDelete(&as, R); 1620 nrzDelete(&a0, R); 1621 nrzDelete(&b0, R); 1622 1623 number a, b, ab; 1624 coeffs Q = nInitChar(n_Q, 0); 1625 nMapFunc f = n_SetMap(R, Q); 1626 a = f(a1, R, Q); 1627 b = f(b1, R, Q); 1628 ab = n_Div(a, b, Q); 1629 n_Delete(&a, Q); 1630 n_Delete(&b, Q); 1631 nKillChar(Q); 1632 1633 nrzDelete(&a1, R); 1634 nrzDelete(&b1, R); 1635 return ab; 1636 } 1637 1638 static number nrzMapMachineInt(number from, const coeffs /*src*/, const coeffs /*dst*/) 1639 { 1640 mpz_ptr erg = (mpz_ptr) omAllocBin(gmp_nrz_bin); 1641 mpz_init_set_ui(erg, (unsigned long) from); 1642 return nrz_short((number) erg); 1643 } 1644 1645 static number nrzMapZp(number from, const coeffs /*src*/, const coeffs /*dst*/) 1646 { 1647 mpz_ptr erg = (mpz_ptr) omAllocBin(gmp_nrz_bin); 1648 mpz_init_set_si(erg, (long) from); 1649 return nrz_short((number) erg); 1650 } 1651 1652 static number nrzModNMap(number from, const coeffs /* src */, const coeffs /*dst*/) 1653 { 1654 mpz_ptr erg = (mpz_ptr) omAllocBin(gmp_nrz_bin); 1655 mpz_init_set(erg, (mpz_ptr) from); 1656 return nrz_short((number) erg); 1657 } 1658 1659 static number nrzMapQ(number from, const coeffs /* src */, const coeffs dst) 1660 { 1661 if (SR_HDL(from) & SR_INT) 1662 return nrzInit(SR_TO_INT(from),dst); 1663 if (from->s!=3) 1664 { 1665 WerrorS("rational in map to integer"); 1666 return NULL; 1667 } 1668 mpz_ptr erg = (mpz_ptr) omAllocBin(gmp_nrz_bin); 1669 mpz_init_set(erg, from->z); 1670 return nrz_short((number) erg); 1671 } 1672 1673 static nMapFunc nrzSetMap(const coeffs src, const coeffs /*dst*/) 1674 { 1675 /* dst = rintegers */ 1676 if (src->rep==n_rep_gmp) //nCoeff_is_Zn(src) || nCoeff_is_Ring_PtoM(src)) 1677 return nrzModNMap; 1678 1679 if ((src->rep==n_rep_gap_gmp) && nCoeff_is_Z(src)) 1680 { 1681 return ndCopyMap; //nrzCopyMap; 1682 } 1683 if (src->rep==n_rep_gap_rat) /*&& nCoeff_is_Z(src)) Q, bigint*/ 1684 { 1685 return nrzMapQ; 1686 } 1687 if ((src->rep==n_rep_int) && nCoeff_is_Ring_2toM(src)) 1688 { 1689 return nrzMapMachineInt; 1690 } 1691 if ((src->rep==n_rep_int) && nCoeff_is_Zp(src)) 1692 { 1693 return nrzMapZp; 1694 } 1695 return NULL; // default 1696 } 1697 1698 1699 /* 1700 * set the exponent (allocate and init tables) (TODO) 1701 */ 1702 1703 void nrzSetExp(int, coeffs) 1704 { 1705 } 1706 1707 void nrzInitExp(int, coeffs) 1708 { 1709 } 1710 1711 #ifdef LDEBUG 1712 BOOLEAN nrzDBTest (number x, const char *f, const int l, const coeffs) 1713 { 1714 if (SR_HDL(x) & SR_INT) return TRUE; 1715 if (mpz_sgn1((mpz_ptr) x)==0) 1716 { 1717 Print("gmp-0 %s:%d\n",f,l); 1718 return FALSE; 1719 } 1720 if (mpz_size1((mpz_ptr)x)<=MP_SMALL) 1721 { 1722 long ui=mpz_get_si((mpz_ptr)x); 1723 if ((((ui<<3)>>3)==ui) 1724 && (mpz_cmp_si((mpz_ptr)x,ui)==0)) 1725 { 1726 Print("gmp-small %s:%d\n",f,l); 1727 return FALSE; 1728 } 1729 } 1730 return TRUE; 1731 } 1732 #endif 1733 1734 void nrzWrite (number a, const coeffs) 1735 { 1736 char *s,*z; 1737 if (a==NULL) 1738 { 1739 StringAppendS("o"); 1740 } 1741 else 1742 { 1743 if (n_Z_IS_SMALL(a)) 1744 { 1745 StringAppend("%d", SR_TO_INT(a)); 1746 } 1747 else 1748 { 1749 int l=mpz_sizeinbase((mpz_ptr) a, 10) + 2; 1750 s=(char*)omAlloc(l); 1751 z=mpz_get_str(s,10,(mpz_ptr) a); 1752 StringAppendS(z); 1753 omFreeSize((ADDRESS)s,l); 1754 } 1755 } 1756 } 1757 1758 /*2 1759 * extracts a long integer from s, returns the rest (COPY FROM longrat0.cc) 1760 */ 1761 static const char * nlEatLongC(char *s, mpz_ptr i) 1762 { 1763 const char * start=s; 1764 1765 if (*s<'0' || *s>'9') 1766 { 1767 mpz_set_ui(i,1); 1768 return s; 1769 } 1770 while (*s >= '0' && *s <= '9') s++; 1771 if (*s=='\0') 1772 { 1773 mpz_set_str(i,start,10); 1774 } 1775 else 1776 { 1777 char c=*s; 1778 *s='\0'; 1779 mpz_set_str(i,start,10); 1780 *s=c; 1781 } 1782 return s; 1783 } 1784 1785 static const char * nrzRead (const char *s, number *a, const coeffs) 1786 { 1787 mpz_ptr z = (mpz_ptr) omAllocBin(gmp_nrz_bin); 1788 { 1789 mpz_init(z); 1790 s = nlEatLongC((char *) s, z); 1791 } 1792 *a = nrz_short((number) z); 1793 return s; 1794 } 1795 1796 static void nrzCoeffWrite (const coeffs, BOOLEAN /*details*/) 1797 { 1798 //PrintS("// ZZ\n"); 1799 PrintS("// coeff. ring is : Integers\n"); 1800 } 1801 1802 static char* nrzCoeffString(const coeffs) 1803 { 1804 return omStrDup("integer"); 1805 } 1806 1807 static CanonicalForm nrzConvSingNFactoryN( number n, BOOLEAN setChar, const coeffs /*r*/ ) 1808 { 1809 if (setChar) setCharacteristic( 0 ); 1810 1811 CanonicalForm term; 1812 if ( n_Z_IS_SMALL(n)) 1813 { 1814 term = SR_TO_INT(n); 1815 } 1816 else 1817 { 1818 mpz_t dummy; 1819 mpz_init_set( dummy,n->z ); 1820 term = make_cf( dummy ); 1821 } 1822 return term; 1823 } 1824 1825 static number nrzConvFactoryNSingN( const CanonicalForm n, const coeffs r) 1826 { 1827 if (n.isImm()) 1828 { 1829 return nrzInit(n.intval(),r); 1830 } 1831 else 1832 { 1833 if ( !n.den().isOne() ) 1834 { 1835 WerrorS("rational in conversion to integer"); 1836 return NULL; 1837 } 1838 mpz_ptr z = (mpz_ptr) omAlloc0Bin(gmp_nrz_bin); 1839 gmp_numerator( n,z); 1840 return nrz_short((number)z); 1841 } 1842 } 1843 1844 static void nrzMPZ(mpz_t res, number &a, const coeffs) 1845 { 1846 if (n_Z_IS_SMALL(a)) 1847 mpz_init_set_si(res, SR_TO_INT(a)); 1848 else 1849 mpz_init_set(res, (mpz_ptr) a); 1850 } 1851 1852 static coeffs nrzQuot1(number c, const coeffs r) 1853 { 1854 mpz_t dummy; 1855 if(n_Z_IS_SMALL(c)) 1856 { 1857 long ch = r->cfInt(c, r); 1858 mpz_init_set_ui(dummy, ch); 1859 } 1860 else 1861 { 1862 mpz_init_set(dummy, (mpz_ptr)c); 1863 } 1864 ZnmInfo info; 1865 info.base = dummy; 1866 info.exp = (unsigned long) 1; 1867 coeffs rr = nInitChar(n_Zn, (void*)&info); 1868 mpz_clear(dummy); 1869 return(rr); 1870 } 1871 1872 number nlReadFd(const ssiInfo *d, const coeffs); 1873 void nlWriteFd(number n, const ssiInfo* d, const coeffs); 1874 1875 BOOLEAN nrzInitChar(coeffs r, void *) 1876 { 1877 assume( getCoeffType(r) == n_Z ); 1878 1879 r->is_field=FALSE; 1880 r->is_domain=TRUE; 1881 r->rep=n_rep_gap_gmp; 1882 1883 //r->nCoeffIsEqual = ndCoeffIsEqual; 1884 r->cfCoeffString = nrzCoeffString; 1885 //r->cfKillChar = ndKillChar; 1886 r->cfMult = nrzMult; 1887 r->cfSub = nrzSub; 1888 r->cfAdd = nrzAdd; 1889 r->cfDiv = nrzDiv; 1890 r->cfIntMod= nrzIntMod; 1891 r->cfExactDiv= nrzExactDiv; 1892 r->cfInit = nrzInit; 1893 r->cfInitMPZ = nrzInitMPZ; 1894 r->cfSize = nrzSize; 1895 r->cfInt = nrzInt; 1896 //#ifdef HAVE_RINGS 1897 r->cfDivComp = nrzDivComp; // only for ring stuff 1898 r->cfIsUnit = nrzIsUnit; // only for ring stuff 1899 r->cfGetUnit = nrzGetUnit; // only for ring stuff 1900 r->cfAnn = nrzAnn; 1901 r->cfExtGcd = nrzExtGcd; // only for ring stuff 1902 r->cfXExtGcd = nrzXExtGcd; // only for ring stuff 1903 r->cfEucNorm = nrzEucNorm; 1904 r->cfQuotRem = nrzSmallestQuotRem; 1905 r->cfDivBy = nrzDivBy; // only for ring stuff 1906 //#endif 1907 r->cfInpNeg = nrzNeg; 1908 r->cfInvers= nrzInvers; 1909 r->cfCopy = nrzCopy; 1910 r->cfWriteLong = nrzWrite; 1911 r->cfRead = nrzRead; 1912 r->cfGreater = nrzGreater; 1913 r->cfEqual = nrzEqual; 1914 r->cfIsZero = nrzIsZero; 1915 r->cfIsOne = nrzIsOne; 1916 r->cfIsMOne = nrzIsMOne; 1917 r->cfGreaterZero = nrzGreaterZero; 1918 r->cfPower = nrzPower; 1919 r->cfGcd = nrzGcd; 1920 r->cfLcm = nrzLcm; 1921 r->cfDelete= nrzDelete; 1922 r->cfSetMap = nrzSetMap; 1923 r->cfCoeffWrite = nrzCoeffWrite; 1924 r->convSingNFactoryN = nrzConvSingNFactoryN; 1925 r->convFactoryNSingN = nrzConvFactoryNSingN; 1926 r->cfMPZ = nrzMPZ; 1927 r->cfFarey = nrzFarey; 1928 r->cfWriteFd=nlWriteFd; 1929 r->cfReadFd=nlReadFd; 1930 1931 r->cfQuot1 = nrzQuot1; 1932 // requires conversion to factory: 1933 r->cfChineseRemainder=nlChineseRemainderSym; 1934 // debug stuff 1935 1936 #ifdef LDEBUG 1937 r->cfDBTest=nrzDBTest; 1938 #endif 1939 1940 r->ch = 0; 1941 r->has_simple_Alloc=FALSE; 1942 r->has_simple_Inverse=FALSE; 1943 return FALSE; 1944 } 1945 44 #include "rintegers3.cc" 1946 45 #elif SI_INTEGER_VARIANT == 1 1947 46 BOOLEAN nrzInitChar(coeffs r, void *) -
libpolys/coeffs/rintegers.h
r03ff5f r9e2cf88 10 10 11 11 #ifdef HAVE_RINGS 12 #include "coeffs/coeffs.h" 12 13 13 14 #if SI_INTEGER_VARIANT == 3
Note: See TracChangeset
for help on using the changeset viewer.
