Changeset 3ef2d6 in git


Ignore:
Timestamp:
Feb 2, 2012, 2:34:22 PM (12 years ago)
Author:
Martin Lee <martinlee84@…>
Branches:
(u'fieker-DuVal', '117eb8c30fc9e991c4decca4832b1d19036c4c65')(u'spielwiese', 'b21a664aa22dc6e196223af8a74ad4885e83547c')
Children:
deb0f602cd605c4cc65110a3abe2c5eeb3b47945
Parents:
935e8385a885664de6c3ca79f6d8a8373d118f85
git-author:
Martin Lee <martinlee84@web.de>2012-02-02 14:34:22+01:00
git-committer:
Martin Lee <martinlee84@web.de>2012-04-04 14:42:25+02:00
Message:
chg: use Newton inversion for division over Q(a)
File:
1 edited

Legend:

Unmodified
Added
Removed
  • factory/facHensel.cc

    r935e83 r3ef2d6  
    239239  fmpz_poly_clear (FLINTB);
    240240  return A;
     241}
     242
     243CanonicalForm
     244mulFLINTQTrunc (const CanonicalForm& F, const CanonicalForm& G, int m)
     245{
     246  if (F.inCoeffDomain() || G.inCoeffDomain())
     247    return mod (F*G, power (Variable (1), m));
     248  Variable alpha;
     249  if (hasFirstAlgVar (F, alpha) || hasFirstAlgVar (G, alpha))
     250    return mulFLINTQaTrunc (F, G, alpha, m);
     251
     252  CanonicalForm A= F;
     253  CanonicalForm B= G;
     254
     255  CanonicalForm denA= bCommonDen (A);
     256  CanonicalForm denB= bCommonDen (B);
     257
     258  A *= denA;
     259  B *= denB;
     260  fmpz_poly_t FLINTA,FLINTB;
     261  convertFacCF2Fmpz_poly_t (FLINTA, A);
     262  convertFacCF2Fmpz_poly_t (FLINTB, B);
     263  fmpz_poly_mullow (FLINTA, FLINTA, FLINTB, m);
     264  denA *= denB;
     265  A= convertFmpz_poly_t2FacCF (FLINTA, F.mvar());
     266  A /= denA;
     267  fmpz_poly_clear (FLINTA);
     268  fmpz_poly_clear (FLINTB);
     269
     270  return A;
     271}
     272
     273CanonicalForm uniReverse (const CanonicalForm& F, int d)
     274{
     275  if (d == 0)
     276    return F;
     277  if (F.inCoeffDomain())
     278    return F*power (Variable (1),d);
     279  Variable x= Variable (1);
     280  CanonicalForm result= 0;
     281  CFIterator i= F;
     282  while (d - i.exp() < 0)
     283    i++;
     284
     285  for (; i.hasTerms() && (d - i.exp() >= 0); i++)
     286    result += i.coeff()*power (x, d - i.exp());
     287  return result;
     288}
     289
     290CanonicalForm
     291newtonInverse (const CanonicalForm& F, const int n)
     292{
     293  int l= ilog2(n);
     294
     295  CanonicalForm g= F [0];
     296
     297  ASSERT (!g.isZero(), "expected a unit");
     298
     299  if (!g.isOne())
     300    g = 1/g;
     301  Variable x= Variable (1);
     302  CanonicalForm result;
     303  int exp= 0;
     304  if (n & 1)
     305  {
     306    result= g;
     307    exp= 1;
     308  }
     309  CanonicalForm h;
     310
     311  for (int i= 1; i <= l; i++)
     312  {
     313    h= mulNTL (g, mod (F, power (x, (1 << i))));
     314    h= mod (h, power (x, (1 << i)) - 1);
     315    h= div (h, power (x, (1 << (i - 1))));
     316    g -= power (x, (1 << (i - 1)))*
     317         mulFLINTQTrunc (g, h, 1 << (i-1));
     318
     319    if (n & (1 << i))
     320    {
     321      if (exp)
     322      {
     323        h= mulNTL (result, mod (F, power (x, exp + (1 << i))));
     324        h= mod (h, power (x, exp + (1 << i)) - 1);
     325        h= div (h, power (x, exp));
     326        result -= power(x, exp)*mulFLINTQTrunc (g, h, 1 << i);
     327        exp += (1 << i);
     328      }
     329      else
     330      {
     331        exp= (1 << i);
     332        result= g;
     333      }
     334    }
     335  }
     336
     337  return result;
     338}
     339
     340void
     341newtonDivrem (const CanonicalForm& F, const CanonicalForm& G, CanonicalForm& Q,
     342              CanonicalForm& R)
     343{
     344  CanonicalForm A= F;
     345  CanonicalForm B= G;
     346  Variable x= Variable (1);
     347  int degA= degree (A, x);
     348  int degB= degree (B, x);
     349  int m= degA - degB;
     350
     351  if (m < 0)
     352  {
     353    R= A;
     354    Q= 0;
     355    return;
     356  }
     357
     358  if (degB <= 1)
     359    divrem (A, B, Q, R);
     360  else
     361  {
     362    R= uniReverse (A, degA);
     363
     364    CanonicalForm revB= uniReverse (B, degB);
     365    CanonicalForm buf= revB;
     366    revB= newtonInverse (revB, m + 1);
     367    Q= mulFLINTQTrunc (R, revB, m + 1);
     368    Q= uniReverse (Q, m);
     369
     370    R= A - mulNTL (Q, B);
     371  }
     372}
     373
     374void
     375newtonDiv (const CanonicalForm& F, const CanonicalForm& G, CanonicalForm& Q)
     376{
     377  CanonicalForm A= F;
     378  CanonicalForm B= G;
     379  Variable x= Variable (1);
     380  int degA= degree (A, x);
     381  int degB= degree (B, x);
     382  int m= degA - degB;
     383
     384  if (m < 0)
     385  {
     386    Q= 0;
     387    return;
     388  }
     389
     390  if (degB <= 1)
     391    Q= div (A, B);
     392  else
     393  {
     394    CanonicalForm R= uniReverse (A, degA);
     395
     396    CanonicalForm revB= uniReverse (B, degB);
     397    revB= newtonInverse (revB, m + 1);
     398    Q= mulFLINTQTrunc (R, revB, m + 1);
     399    Q= uniReverse (Q, m);
     400  }
    241401}
    242402
     
    610770      return modFLINTQ (F, G);
    611771    else
    612       //TODO newtonDivrem
    613      return mod (F, G);
     772    {
     773      CanonicalForm Q, R;
     774      newtonDivrem (F, G, Q, R);
     775      return R;
     776    }
    614777#else
    615778    return mod (F, G);
     
    670833      return divFLINTQ (F,G);
    671834    else
    672       //TODO newtonDivrem
    673       return div (F, G);
     835    {
     836      CanonicalForm Q;
     837      newtonDiv (F, G, Q);
     838      return Q;
     839    }
    674840#else
    675841    return div (F, G);
Note: See TracChangeset for help on using the changeset viewer.