Changeset fe99ed in git


Ignore:
Timestamp:
Jul 22, 2014, 12:34:15 PM (10 years ago)
Author:
Claus Fieker <fieker@…>
Branches:
(u'fieker-DuVal', '117eb8c30fc9e991c4decca4832b1d19036c4c65')(u'spielwiese', 'b52fc4b2495505785981d640dcf7eb3e456778ef')
Children:
4e65f09322d32da3c0c1b29ba2d27d1aaf5a50ff
Parents:
9a6f8c15b898f962115ce36f094aaa7f9dc0be1f
git-author:
Claus Fieker <fieker@mathematik.uni-kl.de>2014-07-22 12:34:15+02:00
git-committer:
Hans Schoenemann <hannes@mathematik.uni-kl.de>2014-07-23 13:12:50+02:00
Message:
"fix" the SI==3 case in rintegers and add the missing bits in rmodulon
again
Try a short patch to keep it simple
Location:
libpolys/coeffs
Files:
4 edited

Legend:

Unmodified
Added
Removed
  • libpolys/coeffs/rintegers.cc

    r9a6f8c1 rfe99ed  
    10051005  return c;
    10061006}
    1007 number _nrzAdd (number a, number b, const coeffs R)
     1007number _nrzAdd (number a, number b, const coeffs )
    10081008#else
    1009 number nrzAdd (number a, number b, const coeffs R)
     1009number nrzAdd (number a, number b, const coeffs )
    10101010#endif
    10111011{
     
    10501050}
    10511051
    1052 number nrzSub (number a, number b, const coeffs)
     1052number nrzSub (number a, number b,  const coeffs )
    10531053{
    10541054  if (n_Z_IS_SMALL(a) && n_Z_IS_SMALL(b))
     
    10801080    int_number erg = (int_number) omAllocBin(gmp_nrz_bin);
    10811081    mpz_init(erg);
    1082     if (INT_TO_SR(b)>0)
     1082    if (SR_TO_INT(b)>0)
    10831083      mpz_sub_ui(erg, (int_number) a, (unsigned long)SR_TO_INT(b));
    10841084    else
     
    10971097number  nrzGetUnit (number n, const coeffs r)
    10981098{
    1099   //if (n_Z_IS_SMALL(n))
    1100   //{
    1101   //  if (SR_TO_INT(n)<0)
    1102   //    return INT_TO_SR(-1);
    1103   //  else
    1104   //    return INT_TO_SR(1);
    1105   //}
    11061099  if (nrzGreaterZero(n, r))
    11071100    return INT_TO_SR(1);
     
    11931186}
    11941187
    1195 number nrzDiv (number a,number b, const coeffs R)
     1188number nrzDiv (number a,number b, const coeffs)
    11961189{
    11971190  assume(SR_TO_INT(b));
     
    14081401}
    14091402
    1410 number nrzModNMap(number from, const coeffs src, const coeffs /*dst*/)
     1403number nrzModNMap(number from, const coeffs /* src */, const coeffs /*dst*/)
    14111404{
    14121405  int_number erg = (int_number) omAllocBin(gmp_nrz_bin);
     
    14151408}
    14161409
    1417 number nrzMapQ(number from, const coeffs src, const coeffs dst)
     1410number nrzMapQ(number from, const coeffs /* src */, const coeffs dst)
    14181411{
    14191412  if (SR_HDL(from) & SR_INT)
     
    15621555  return omStrDup("integer");
    15631556}
     1557
     1558#include "factory/canonicalform.h"
     1559#include <factory/cf_gmp.h>
     1560#include "factory/singext.h"
    15641561
    15651562static CanonicalForm nrzConvSingNFactoryN( number n, BOOLEAN setChar, const coeffs /*r*/ )
  • libpolys/coeffs/rintegers.h

    r9a6f8c1 rfe99ed  
    1010#include <coeffs/coeffs.h>
    1111
    12 extern int nrzExp;
    13 extern NATNUMBER nrzModul;
     12#if SI_INTEGER_VARIANT == 3
     13#define SR_HDL(A) ((long)(A))
     14#define SR_INT    1L
     15#define INT_TO_SR(INT)  ((number) (((long)INT << 2) + SR_INT))
     16#define SR_TO_INT(SR)   (((long)SR) >> 2)
     17#define n_Z_IS_SMALL(A)     (SR_HDL(A) & SR_INT)
     18#define INT_IS_SMALL(A) ( ((A << 1) >> 1) == A )
     19#endif
     20
     21//extern int nrzExp;
     22//extern NATNUMBER nrzModul;
    1423
    1524BOOLEAN nrzInitChar    (coeffs r,  void * parameter);
  • libpolys/coeffs/rmodulon.cc

    r9a6f8c1 rfe99ed  
    5757  omFreeSize(b,l);
    5858  return s;
     59}
     60
     61static void nrnKillChar(coeffs r)
     62{
     63  mpz_clear(r->modNumber);
     64  mpz_clear(r->modBase);
     65  omFreeBin((void *) r->modBase, gmp_nrz_bin);
     66  omFreeBin((void *) r->modNumber, gmp_nrz_bin);
    5967}
    6068
     
    104112}
    105113
    106 static number nrnAnn(number b, const coeffs r);
    107114/* for initializing function pointers */
    108115BOOLEAN nrnInitChar (coeffs r, void* p)
     
    110117  assume( (getCoeffType(r) == ID) || (getCoeffType (r) == ID2) );
    111118  ZnmInfo * info= (ZnmInfo *) p;
    112   r->modBase= info->base;
     119  r->modBase= (int_number)nrnCopy((number)info->base, r); //this circumvents the problem
     120  //in bigintmat.cc where we cannot create a "legal" nrn that can be freed.
     121  //If we take a copy, we can do whatever we want.
    113122
    114123  nrnInitExp (info->exp, r);
     
    157166  r->cfGetUnit     = nrnGetUnit;
    158167  r->cfExtGcd      = nrnExtGcd;
     168  r->cfXExtGcd     = nrnXExtGcd;
     169  r->cfQuotRem     = nrnQuotRem;
    159170  r->cfName        = ndName;
    160171  r->cfCoeffWrite  = nrnCoeffWrite;
    161172  r->nCoeffIsEqual = nrnCoeffsEqual;
    162   r->cfKillChar    = ndKillChar;
     173  r->cfKillChar    = nrnKillChar;
    163174  r->cfQuot1       = nrnQuot1;
    164175#ifdef LDEBUG
     
    307318 * Give the largest k, such that a = x * k, b = y * k has
    308319 * a solution and r, s, s.t. k = s*a + t*b
     320 * CF: careful: ExtGcd is wrong as implemented (or at least may not
     321 * give you what you want:
     322 * ExtGcd(5, 10 modulo 12):
     323 * the gcdext will return 5 = 1*5 + 0*10
     324 * however, mod 12, the gcd should be 1
    309325 */
    310326number nrnExtGcd(number a, number b, number *s, number *t, const coeffs r)
     
    323339  return (number)erg;
    324340}
     341/* XExtGcd  returns a unimodular matrix ((s,t)(u,v)) sth.
     342 * (a,b)^t ((st)(uv)) = (g,0)^t
     343 * Beware, the ExtGcd will not necessaairly do this.
     344 * Problem: if g = as+bt then (in Z/nZ) it follows NOT that
     345 *             1 = (a/g)s + (b/g) t
     346 * due to the zero divisors.
     347 */
     348
     349//#define CF_DEB;
     350number nrnXExtGcd(number a, number b, number *s, number *t, number *u, number *v, const coeffs r)
     351{
     352  number xx;
     353#ifdef CF_DEB
     354  StringSetS("XExtGcd of ");
     355  nrnWrite(a, r);
     356  StringAppendS("\t");
     357  nrnWrite(b, r);
     358  StringAppendS(" modulo ");
     359  nrnWrite(xx = (number)r->modNumber, r);
     360  Print("%s\n", StringEndS());
     361#endif
     362
     363  int_number one = (int_number)omAllocBin(gmp_nrz_bin);
     364  int_number erg = (int_number)omAllocBin(gmp_nrz_bin);
     365  int_number bs  = (int_number)omAllocBin(gmp_nrz_bin);
     366  int_number bt  = (int_number)omAllocBin(gmp_nrz_bin);
     367  int_number bu  = (int_number)omAllocBin(gmp_nrz_bin);
     368  int_number bv  = (int_number)omAllocBin(gmp_nrz_bin);
     369  mpz_init(erg);
     370  mpz_init(one);
     371  mpz_init_set(bs, (int_number) a);
     372  mpz_init_set(bt, (int_number) b);
     373  mpz_init(bu);
     374  mpz_init(bv);
     375  mpz_gcd(erg, bs, bt);
     376
     377#ifdef CF_DEB
     378  StringSetS("1st gcd:");
     379  nrnWrite(xx= (number)erg, r);
     380#endif
     381
     382  mpz_gcd(erg, erg, r->modNumber);
     383
     384  mpz_div(bs, bs, erg);
     385  mpz_div(bt, bt, erg);
     386
     387#ifdef CF_DEB
     388  Print("%s\n", StringEndS());
     389  StringSetS("xgcd: ");
     390#endif
     391
     392  mpz_gcdext(one, bu, bv, bs, bt);
     393  number ui = nrnGetUnit(xx = (number) one, r);
     394#ifdef CF_DEB
     395  n_Write(xx, r);
     396  StringAppendS("\t");
     397  n_Write(ui, r);
     398  Print("%s\n", StringEndS());
     399#endif
     400  nrnDelete(&xx, r);
     401  if (!nrnIsOne(ui, r)) {
     402#ifdef CF_DEB
     403    Print("Scaling\n");
     404#endif
     405    number uii = nrnInvers(ui, r);
     406    nrnDelete(&ui, r);
     407    ui = uii;
     408    int_number uu = (int_number)omAllocBin(gmp_nrz_bin);
     409    mpz_init_set(uu, (int_number)ui);
     410    mpz_mul(bu, bu, uu);
     411    mpz_mul(bv, bv, uu);
     412    mpz_clear(uu);
     413    omFreeBin(uu, gmp_nrz_bin);
     414  }
     415  nrnDelete(&ui, r);
     416#ifdef CF_DEB
     417  StringSetS("xgcd");
     418  nrnWrite(xx= (number)bs, r);
     419  StringAppendS("*");
     420  nrnWrite(xx= (number)bu, r);
     421  StringAppendS(" + ");
     422  nrnWrite(xx= (number)bt, r);
     423  StringAppendS("*");
     424  nrnWrite(xx= (number)bv, r);
     425  Print("%s\n", StringEndS());
     426#endif
     427
     428  mpz_mod(bs, bs, r->modNumber);
     429  mpz_mod(bt, bt, r->modNumber);
     430  mpz_mod(bu, bu, r->modNumber);
     431  mpz_mod(bv, bv, r->modNumber);
     432  *s = (number)bu;
     433  *t = (number)bv;
     434  *u = (number)bt;
     435  *u = nrnNeg(*u, r);
     436  *v = (number)bs;
     437  return (number)erg;
     438}
     439
    325440
    326441BOOLEAN nrnIsZero(number a, const coeffs)
     
    410525  nrnDelete((number*) &gcd, NULL);
    411526  return (number)unit;
     527}
     528
     529number nrnAnn(number k, const coeffs r)
     530{
     531  int_number tmp = (int_number) omAllocBin(gmp_nrz_bin);
     532  mpz_init(tmp);
     533  mpz_gcd(tmp, (int_number) k, r->modNumber);
     534  if (mpz_cmp_si(tmp, 1)==0) {
     535    mpz_set_si(tmp, 0);
     536    return (number) tmp;
     537  }
     538  mpz_divexact(tmp, r->modNumber, tmp);
     539  return (number) tmp;
    412540}
    413541
     
    508636}
    509637
    510 static number nrnAnn(number b, const coeffs r)
    511 {
    512   int_number erg = (int_number)omAllocBin(gmp_nrz_bin);
    513   mpz_init(erg);
    514   mpz_tdiv_q(erg, (int_number)r->modNumber, (int_number)b);
    515   return (number)erg;
     638/* CF: note that Z/nZ has (at least) two distinct euclidean structures
     639 * 1st phi(a) := (a mod n) which is just the structure directly
     640 *     inherited from Z
     641 * 2nd phi(a) := gcd(a, n)
     642 * The 1st version is probably faster as everything just comes from Z,
     643 * but the 2nd version behaves nicely wrt. to quotient operations
     644 * and HNF and such. In agreement with nrnMod we imlement the 2nd here
     645 *
     646 * For quotrem note that if b exactly divides a, then
     647 *   min(v_p(a), v_p(n))  >= min(v_p(b), v_p(n))
     648 * so if we divide  a and b by g:= gcd(a,b,n), then   b becomes a
     649 * unit mod n/g.
     650 * Thus we 1st compute the remainder (similar to nrnMod) and then
     651 * the exact quotient.
     652 */
     653number nrnQuotRem(number a, number b, number  * rem, const coeffs r)
     654{
     655  mpz_t g, aa, bb;
     656  int_number qq = (int_number)omAllocBin(gmp_nrz_bin);
     657  int_number rr = (int_number)omAllocBin(gmp_nrz_bin);
     658  mpz_init(qq);
     659  mpz_init(rr);
     660  mpz_init(g);
     661  mpz_init_set(aa, (int_number)a);
     662  mpz_init_set(bb, (int_number)b);
     663
     664  mpz_gcd(g, bb, r->modNumber);
     665  mpz_mod(rr, aa, g);
     666  mpz_sub(aa, aa, rr);
     667  mpz_gcd(g, aa, g);
     668  mpz_div(aa, aa, g);
     669  mpz_div(bb, bb, g);
     670  mpz_div(g, r->modNumber, g);
     671  mpz_invert(g, bb, g);
     672  mpz_mul(qq, aa, g);
     673  if (rem)
     674    *rem = (number)rr;
     675  else {
     676    mpz_clear(rr);
     677    omFreeBin(rr, gmp_nrz_bin);
     678  }
     679  mpz_clear(g);
     680  mpz_clear(aa);
     681  mpz_clear(bb);
     682  return (number) qq;
    516683}
    517684
     
    554721}
    555722
     723#if SI_INTEGER_VARIANT==3
     724number nrnMapZ(number from, const coeffs /*src*/, const coeffs dst)
     725{
     726  int_number erg = (int_number)omAllocBin(gmp_nrz_bin);
     727  if (n_Z_IS_SMALL(from))
     728    mpz_init_set_si(erg, SR_TO_INT(from));
     729  else
     730    mpz_init_set(erg, (int_number) from);
     731  mpz_mod(erg, erg, dst->modNumber);
     732  return (number)erg;
     733}
     734#elif SI_INTEGER_VARIANT==2
     735
    556736number nrnMapZ(number from, const coeffs src, const coeffs dst)
    557737{
     
    563743  return nrnMapGMP(from,src,dst);
    564744}
     745#endif
    565746
    566747number nrnMapQ(number from, const coeffs src, const coeffs dst)
     
    578759  if ((src->rep==n_rep_gmp) && nCoeff_is_Ring_Z(src))
    579760  {
    580     return nrnMapGMP;
     761    return nrnMapZ;
    581762  }
    582763  if ((src->rep==n_rep_gap_gmp) /*&& nCoeff_is_Ring_Z(src)*/)
     
    674855  nrnSetExp(m, r);
    675856  assume (r->modNumber != NULL);
    676   if (mpz_cmp_ui(r->modNumber,2) <= 0)
    677     WarnS("nrnInitExp failed (m in Z/m too small)");
     857//CF: in geenral, the modulus is computed somewhere. I don't want to
     858//  check it's size before I construct the best ring.
     859//  if (mpz_cmp_ui(r->modNumber,2) <= 0)
     860//    WarnS("nrnInitExp failed (m in Z/m too small)");
    678861}
    679862
  • libpolys/coeffs/rmodulon.h

    r9a6f8c1 rfe99ed  
    3333BOOLEAN nrnIsUnit      (number a, const coeffs r);
    3434number  nrnGetUnit     (number a, const coeffs r);
     35number  nrnAnn         (number a, const coeffs r);
    3536number  nrnDiv         (number a, number b, const coeffs r);
    3637number  nrnMod         (number a,number b, const coeffs r);
     
    4546number  nrnGcd         (number a,number b, const coeffs r);
    4647number  nrnExtGcd      (number a, number b, number *s, number *t, const coeffs r);
     48number  nrnXExtGcd      (number a, number b, number *s, number *t, number *u, number *v, const coeffs r);
     49number  nrnQuotRem      (number a, number b, number *s, const coeffs r);
    4750nMapFunc nrnSetMap     (const coeffs src, const coeffs dst);
    4851#define  nrnWrite      nrzWrite
Note: See TracChangeset for help on using the changeset viewer.