Ignore:
Timestamp:
Aug 10, 2012, 2:49:42 PM (12 years ago)
Author:
Hans Schoenemann <hannes@…>
Branches:
(u'spielwiese', '5b153614cbc72bfa198d75b1e9e33dab2645d9fe')
Children:
73260639099e75033ef78f48ad34a0ee10e1c15e
Parents:
e11ed0be7258f8a5789019dc86eb3edff31cfaf5
Message:
fix: div/mod stuff for paraPlanCurve: aDiv, aMod
File:
1 edited

Legend:

Unmodified
Added
Removed
  • Singular/LIB/paraplanecurves.lib

    re11ed0 re8eba7  
    11///////////////////////////////////////////////////////////////////////////////
    2 version="$Id$";
     2version="$Id: paraplanecurves.lib 14461 2011-12-14 11:00:17Z motsak $";
    33category="Algebraic Geometry";
    44info="
     
    2626library, also van Hoeij's algorithm for computing the integral basis will
    2727be available. @*
    28 From the integral basis, the adjoint ideal is obtained by linear algebra.
     28>From the integral basis, the adjoint ideal is obtained by linear algebra.
    2929Alternatively, the algorithm starts with a local analysis of the singular
    3030locus of C. Then, for each  primary component of the singular locus which
     
    255255see @ref{integralbasis_lib}. In a future edition of the library, also van Hoeij's
    256256algorithm for computing the integral basis will be available. @*
    257 From the integral basis, the adjoint ideal is obtained by linear algebra.
     257>From the integral basis, the adjoint ideal is obtained by linear algebra.
    258258Alternatively, the algorithm starts with a local analysis of the singular
    259259locus of C. Then, for each  primary component of the singular locus which
     
    486486from @ref{normal_lib}: see @ref{integralbasis_lib}. In a future edition of the library,
    487487also van Hoeij's algorithm for computing the integral basis will be available.@*
    488 From the integral basis, the adjoint ideal is obtained by linear algebra.
     488>From the integral basis, the adjoint ideal is obtained by linear algebra.
    489489Alternatively, the algorithm starts with a local analysis of the singular
    490490locus of C. Then, for each  primary component of the singular locus which
     
    19411941
    19421942///////////////////////////////////////////////////////////////////////////////
     1943static proc aMod(bigint a, bigint b)
     1944"USAGE:  aMod(a,b); a, b bigint
     1945RETURN:  r bigint
     1946THEORY:  The asymmetric residue r of the division with remainder a mod b.
     1947"
     1948{
     1949  bigint c = a mod b;
     1950  if (c<0)
     1951  {
     1952    return(c+b);
     1953  }
     1954return(c);
     1955}
     1956
     1957///////////////////////////////////////////////////////////////////////////////
     1958static proc aDiv(bigint a, bigint b)
     1959"USAGE:  aDiv(a,b); a, b bigint
     1960RETURN:  q bigint
     1961THEORY:  Quotient with remainder q = a div b with asymmetric residue.
     1962"
     1963{
     1964  bigint q = a div b;
     1965  if ((a mod b)<0)
     1966  {
     1967    return(q-1);
     1968  }
     1969return(q);
     1970}
     1971
     1972
     1973
     1974///////////////////////////////////////////////////////////////////////////////
    19431975static proc polyModP(poly q, bigint p)
    19441976"USAGE:  polyModP(q, p); q poly, p bigint
     
    20152047      b = inverseModP(bigint(leadcoef(q)), p);
    20162048      q = q - lead(q);
    2017       root = (bigint(q) * b) mod p;
     2049      root = aMod((bigint(q) * b),p);
    20182050      if (((root * root - r) mod p) != 0) { root = 0; }
    20192051    }
     
    20352067  list L = extendedEuclid(r, p);
    20362068  if (L[1] != 1) { "ERROR: GCD of", r, "and", p, "should be 1."; }
    2037   L[2] = L[2] mod p;
     2069  L[2] = aMod(L[2],p);
    20382070  return (L[2]);
    20392071}
     
    20832115    for (i = 1; i <= size(L[1]); i++)
    20842116    {
    2085       roots = insert(roots, rootModP(r mod L[1][i], L[1][i]), size(roots));
     2117      roots = insert(roots, rootModP(aMod(r,L[1][i]), L[1][i]), size(roots));
    20862118    }
    20872119
     
    20952127        c = bigint(roots[i]); l = pPower;
    20962128        temp = r - c * c; l = bigint(2) * c * l; c = temp;
    2097         c = c div pPower; l = l div pPower;
    2098         c = c mod L[1][i]; l = l mod L[1][i];
    2099         c = (c * bigint(inverseModP(l, L[1][i]))) mod L[1][i];
     2129        c = aDiv(c,pPower); l = aDiv(l,pPower);
     2130        c = aMod(c,L[1][i]); l = aMod(l,L[1][i]);
     2131        c = aMod((c * bigint(inverseModP(l, L[1][i]))), L[1][i]);
    21002132        c = bigint(roots[i]) + c * pPower;
    21012133        pPower = pPower * L[1][i]; roots[i] = c;
     
    21142146      mm = insert(mm, z, size(mm));
    21152147    }
    2116     return (chineseRemainder(roots, mm) mod m);
     2148    return (aMod(chineseRemainder(roots, mm) , m));
    21172149  }
    21182150}
     
    21342166  for (i = 1; i <= size(mm); i++)
    21352167  {
    2136     N = M div mm[i];
     2168    N = aDiv(M,mm[i]);
    21372169    l = extendedEuclid(mm[i], N);
    21382170    x = x + rr[i]*l[3]*N;
     
    22092241    cc = c; if (cc <= 0) { cc = -cc; }
    22102242    R1 = squareRoot(-a*b, cc, 0);
    2211     r1 = (R1 * inverseModP(a, cc)) mod cc;
     2243    r1 = aMod((R1 * inverseModP(a, cc)), cc);
    22122244    if (r1*bigint(2) > cc) { r1 = r1 - cc; }
    22132245    Q = (a*r1*r1 + b)/c;
     
    22722304  if (a == 0)         { l = b, 0, 1; if (b < 0) { l = -b, 0, -1; } }
    22732305  if (b == 0)         { l = a, 1, 0; if (a < 0) { l = -a, -1, 0; } }
    2274   if ((a mod b) == 0) { l = b, 0, 1; if (b < 0) { l = -b, 0, -1; } }
    2275   if ((b mod a) == 0) { l = a, 1, 0; if (a < 0) { l = -a, -1, 0; } }
     2306  if (aMod(a , b) == 0) { l = b, 0, 1; if (b < 0) { l = -b, 0, -1; } }
     2307  if (aMod(b , a) == 0) { l = a, 1, 0; if (a < 0) { l = -a, -1, 0; } }
    22762308  if (size(l) > 1) { return (l); }
    22772309
    2278   temp = a mod b;
     2310  temp = aMod(a , b);
    22792311  l = extendedEuclid(b, temp);
    22802312  temp = (a - temp) / b;
     
    22922324"
    22932325{
    2294   bigint rr = r mod p;
     2326  bigint rr = aMod(r , p);
    22952327  if (rr == 0) { return (0) }
    22962328  if (rr == 1) { return (1) }
     
    23012333  while (e > 0)
    23022334  {
    2303     if ((e mod 2) == 1) { result = (result * power) mod p; }
     2335    if ((e mod 2) == 1) { result = aMod((result * power), p); }
    23042336    e = e / bigint(2);
    2305     power = (power * power) mod p;
     2337    power = aMod((power * power), p);
    23062338  }
    23072339  if (result > 1) { result = result - p; /* should be -1 */ }
Note: See TracChangeset for help on using the changeset viewer.