Changeset c9b1d66 in git


Ignore:
Timestamp:
Dec 16, 2010, 2:22:47 PM (13 years ago)
Author:
Martin Lee <martinlee84@…>
Branches:
(u'spielwiese', '5b153614cbc72bfa198d75b1e9e33dab2645d9fe')
Children:
8af49d09040b6ced3fa697fae0b4338a908944a5
Parents:
d08ed84a213b4110fe1dbe0d38276dce7695a481
Message:
added division with remainder via Newton inversion


git-svn-id: file:///usr/local/Singular/svn/trunk@13774 2c84dea3-7e68-4137-9b89-c4e89433aadc
Location:
factory
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • factory/facHensel.cc

    rd08ed8 rc9b1d66  
    395395    buf2= prodMod (tmp2, M);
    396396    return mulMod (buf1, buf2, M);
     397  }
     398}
     399
     400
     401CanonicalForm reverse (const CanonicalForm& F, int d)
     402{
     403  if (d == 0)
     404    return F;
     405  CanonicalForm A= F;
     406  Variable y= F.mvar();
     407  Variable x= Variable (1);
     408  A= swapvar (A, x, y);
     409  CanonicalForm result= 0;
     410  CFIterator i= A;
     411  while (d - i.exp() < 0)
     412    i++;
     413
     414  for (; i.hasTerms() && (d - i.exp() >= 0); i++)
     415    result += swapvar (i.coeff(),x,y)*power (x, d - i.exp());
     416  return result;
     417}
     418
     419CanonicalForm
     420newtonInverse (const CanonicalForm& F, const int n, const CanonicalForm& M)
     421{
     422  int l= ilog2(n);
     423
     424  CanonicalForm g= mod (F, M)[0] [0];
     425
     426  ASSERT (!g.isZero(), "expected a unit");
     427
     428  Variable alpha;
     429
     430  if (!g.isOne())
     431    g = 1/g;
     432  Variable x= Variable (1);
     433  CanonicalForm result;
     434  int exp= 0;
     435  if (n & 1)
     436  {
     437    result= g;
     438    exp= 1;
     439  }
     440  CanonicalForm h;
     441
     442  for (int i= 1; i <= l; i++)
     443  {
     444    h= mulMod2 (g, mod (F, power (x, (1 << i))), M);
     445    h= mod (h, power (x, (1 << i)) - 1);
     446    h= div (h, power (x, (1 << (i - 1))));
     447    h= mod (h, M);
     448    g -= power (x, (1 << (i - 1)))*
     449         mod (mulMod2 (g, h, M), power (x, (1 << (i - 1))));
     450
     451    if (n & (1 << i))
     452    {
     453      if (exp)
     454      {
     455        h= mulMod2 (result, mod (F, power (x, exp + (1 << i))), M);
     456        h= mod (h, power (x, exp + (1 << i)) - 1);
     457        h= div (h, power (x, exp));
     458        h= mod (h, M);
     459        result -= power(x, exp)*mod (mulMod2 (g, h, M),
     460                                       power (x, (1 << i)));
     461        exp += (1 << i);
     462      }
     463      else
     464      {
     465        exp= (1 << i);
     466        result= g;
     467      }
     468    }
     469  }
     470
     471  return result;
     472}
     473
     474CanonicalForm
     475newtonDiv (const CanonicalForm& F, const CanonicalForm& G, const CanonicalForm&
     476           M)
     477{
     478  ASSERT (getCharacteristic() > 0, "positive characteristic expected");
     479  ASSERT (CFFactory::gettype() != GaloisFieldDomain, "no GF expected");
     480
     481  CanonicalForm A= mod (F, M);
     482  CanonicalForm B= mod (G, M);
     483
     484  Variable x= Variable (1);
     485  int degA= degree (A, x);
     486  int degB= degree (B, x);
     487  int m= degA - degB;
     488  if (m < 0)
     489    return 0;
     490
     491  CanonicalForm Q;
     492  if (degB <= 1)
     493  {
     494    CanonicalForm R;
     495    divrem2 (A, B, Q, R, M);
     496  }
     497  else
     498  {
     499    CanonicalForm R= reverse (A, degA);
     500    CanonicalForm revB= reverse (B, degB);
     501    revB= newtonInverse (revB, m + 1, M);
     502    Q= mulMod2 (R, revB, M);
     503    Q= mod (Q, power (x, m + 1));
     504    Q= reverse (Q, m);
     505  }
     506
     507  return Q;
     508}
     509
     510void
     511newtonDivrem (const CanonicalForm& F, const CanonicalForm& G, CanonicalForm& Q,
     512              CanonicalForm& R, const CanonicalForm& M)
     513{
     514  CanonicalForm A= mod (F, M);
     515  CanonicalForm B= mod (G, M);
     516  Variable x= Variable (1);
     517  int degA= degree (A, x);
     518  int degB= degree (B, x);
     519  int m= degA - degB;
     520
     521  if (m < 0)
     522  {
     523    R= A;
     524    Q= 0;
     525    return;
     526  }
     527
     528  if (degB <= 1)
     529  {
     530     divrem2 (A, B, Q, R, M);
     531  }
     532  else
     533  {
     534    R= reverse (A, degA);
     535
     536    CanonicalForm revB= reverse (B, degB);
     537    revB= newtonInverse (revB, m + 1, M);
     538    Q= mulMod2 (R, revB, M);
     539
     540    Q= mod (Q, power (x, m + 1));
     541    Q= reverse (Q, m);
     542
     543    R= A - mulMod2 (Q, B, M);
    397544  }
    398545}
  • factory/facHensel.h

    rd08ed8 rc9b1d66  
    2929#include "algext.h"
    3030
    31 CanonicalForm mul (const CanonicalForm& F, const CanonicalForm& G);
    32 CanonicalForm mulMod3 (const CanonicalForm& F, const CanonicalForm& G, const CFList& MOD);
    3331/// multiplication of univariate polys over a finite field using NTL, if we are
    3432/// in GF factory's default multiplication is used.
     
    139137                                   ///< Variables of level higher than 1
    140138            );
     139
     140
     141/// division with remainder of @a F by
     142/// @a G wrt Variable (1) modulo @a M using Newton inversion
     143///
     144/// @return @a Q returns the dividend, @a R returns the remainder.
     145/// @sa divrem2(), newtonDiv()
     146void
     147newtonDivrem (const CanonicalForm& F, ///< [in] bivariate, compressed polynomial
     148              const CanonicalForm& G, ///< [in] bivariate, compressed polynomial
     149                                      ///< which is monic in Variable (1)
     150              CanonicalForm& Q,       ///< [in,out] dividend
     151              CanonicalForm& R,       ///< [in,out] remainder, degree (R, 1) <
     152                                      ///< degree (G, 1)
     153              const CanonicalForm& M  ///< [in] power of Variable (2)
     154             );
     155
     156/// division of @a F by
     157/// @a G wrt Variable (1) modulo @a M using Newton inversion
     158///
     159/// @return @a newtonDiv returns the dividend
     160/// @sa divrem2(), newtonDivrem()
     161CanonicalForm
     162newtonDiv (const CanonicalForm& F, ///< [in] bivariate, compressed polynomial
     163           const CanonicalForm& G, ///< [in] bivariate, compressed polynomial
     164                                   ///< which is monic in Variable (1)
     165           const CanonicalForm& M  ///< [in] power of Variable (2)
     166          );
    141167
    142168/// reduce @a F modulo elements in @a M.
Note: See TracChangeset for help on using the changeset viewer.