Changeset dd2a9b3 in git


Ignore:
Timestamp:
May 13, 2011, 9:44:24 AM (13 years ago)
Author:
Martin Lee <martinlee84@…>
Branches:
(u'fieker-DuVal', '117eb8c30fc9e991c4decca4832b1d19036c4c65')(u'spielwiese', '79dfb9a6d258bfeb991428bdb25b8f55e9e809b2')
Children:
a2dd9b278a186a85e5955f8e7adad56483fddffd
Parents:
4775bdf78134ad29d5eb334ec54f2e1fa7f460ff
Message:
Kronecker substitution multiplication for bivariate multiplication


git-svn-id: file:///usr/local/Singular/svn/trunk@14214 2c84dea3-7e68-4137-9b89-c4e89433aadc
File:
1 edited

Legend:

Unmodified
Added
Removed
  • factory/facHensel.cc

    r4775bd rdd2a9b3  
    204204}
    205205
     206zz_pX kronSubFp (const CanonicalForm& A, int d)
     207{
     208  int degAy= degree (A);
     209  zz_pX result;
     210  result.rep.SetLength (d*(degAy + 1));
     211
     212  zz_p *resultp;
     213  resultp= result.rep.elts();
     214  zz_pX buf;
     215  zz_p *bufp;
     216
     217  for (CFIterator i= A; i.hasTerms(); i++)
     218  {
     219    if (i.coeff().inCoeffDomain())
     220      buf= convertFacCF2NTLzzpX (i.coeff());
     221    else
     222      buf= convertFacCF2NTLzzpX (i.coeff());
     223
     224    int k= i.exp()*d;
     225    bufp= buf.rep.elts();
     226    int bufRepLength= (int) buf.rep.length();
     227    for (int j= 0; j < bufRepLength; j++)
     228      resultp [j + k]= bufp [j];
     229  }
     230  result.normalize();
     231
     232  return result;
     233}
     234
     235zz_pEX kronSub (const CanonicalForm& A, int d, const Variable& alpha)
     236{
     237  int degAy= degree (A);
     238  zz_pEX result;
     239  result.rep.SetLength (d*(degAy + 1));
     240
     241  Variable v;
     242  zz_pE *resultp;
     243  resultp= result.rep.elts();
     244  zz_pEX buf1;
     245  zz_pE *buf1p;
     246  zz_pX buf2;
     247  zz_pX NTLMipo= convertFacCF2NTLzzpX (getMipo (alpha));
     248
     249  for (CFIterator i= A; i.hasTerms(); i++)
     250  {
     251    if (i.coeff().inCoeffDomain())
     252    {
     253      buf2= convertFacCF2NTLzzpX (i.coeff());
     254      buf1= to_zz_pEX (to_zz_pE (buf2));
     255    }
     256    else
     257      buf1= convertFacCF2NTLzz_pEX (i.coeff(), NTLMipo);
     258
     259    int k= i.exp()*d;
     260    buf1p= buf1.rep.elts();
     261    int buf1RepLength= (int) buf1.rep.length();
     262    for (int j= 0; j < buf1RepLength; j++)
     263      resultp [j + k]= buf1p [j];
     264  }
     265  result.normalize();
     266
     267  return result;
     268}
     269
     270void
     271kronSubRecipro (zz_pEX& subA1, zz_pEX& subA2,const CanonicalForm& A, int d,
     272                const Variable& alpha)
     273{
     274  int degAy= degree (A);
     275  subA1.rep.SetLength ((long) d*(degAy + 1));
     276  subA2.rep.SetLength ((long) d*(degAy + 1));
     277
     278  Variable v;
     279  zz_pE *subA1p;
     280  zz_pE *subA2p;
     281  subA1p= subA1.rep.elts();
     282  subA2p= subA2.rep.elts();
     283  zz_pEX buf;
     284  zz_pE *bufp;
     285  zz_pX buf2;
     286  zz_pX NTLMipo= convertFacCF2NTLzzpX (getMipo (alpha));
     287
     288  for (CFIterator i= A; i.hasTerms(); i++)
     289  {
     290    if (i.coeff().inCoeffDomain())
     291    {
     292      buf2= convertFacCF2NTLzzpX (i.coeff());
     293      buf= to_zz_pEX (to_zz_pE (buf2));
     294    }
     295    else
     296      buf= convertFacCF2NTLzz_pEX (i.coeff(), NTLMipo);
     297
     298    int k= i.exp()*d;
     299    int kk= (degAy - i.exp())*d;
     300    bufp= buf.rep.elts();
     301    int bufRepLength= (int) buf.rep.length();
     302    for (int j= 0; j < bufRepLength; j++)
     303    {
     304      subA1p [j + k]= bufp [j];
     305      subA2p [j + kk]= bufp [j];
     306    }
     307  }
     308  subA1.normalize();
     309  subA2.normalize();
     310}
     311
     312void
     313kronSubRecipro (zz_pX& subA1, zz_pX& subA2,const CanonicalForm& A, int d)
     314{
     315  int degAy= degree (A);
     316  subA1.rep.SetLength ((long) d*(degAy + 1));
     317  subA2.rep.SetLength ((long) d*(degAy + 1));
     318
     319  Variable v;
     320  zz_p *subA1p;
     321  zz_p *subA2p;
     322  subA1p= subA1.rep.elts();
     323  subA2p= subA2.rep.elts();
     324  zz_pX buf;
     325  zz_p *bufp;
     326
     327  for (CFIterator i= A; i.hasTerms(); i++)
     328  {
     329    buf= convertFacCF2NTLzzpX (i.coeff());
     330
     331    int k= i.exp()*d;
     332    int kk= (degAy - i.exp())*d;
     333    bufp= buf.rep.elts();
     334    int bufRepLength= (int) buf.rep.length();
     335    for (int j= 0; j < bufRepLength; j++)
     336    {
     337      subA1p [j + k]= bufp [j];
     338      subA2p [j + kk]= bufp [j];
     339    }
     340  }
     341  subA1.normalize();
     342  subA2.normalize();
     343}
     344
     345CanonicalForm
     346reverseSubst (const zz_pEX& F, const zz_pEX& G, int d, int k,
     347              const Variable& alpha)
     348{
     349  Variable y= Variable (2);
     350  Variable x= Variable (1);
     351
     352  zz_pEX f= F;
     353  zz_pEX g= G;
     354  int degf= deg(f);
     355  int degg= deg(g);
     356
     357  zz_pEX buf1;
     358  zz_pEX buf2;
     359  zz_pEX buf3;
     360
     361  zz_pE *buf1p;
     362  zz_pE *buf2p;
     363  zz_pE *buf3p;
     364  if (f.rep.length() < (long) d*(k+1)) //zero padding
     365    f.rep.SetLength ((long)d*(k+1));
     366
     367  zz_pE *gp= g.rep.elts();
     368  zz_pE *fp= f.rep.elts();
     369  CanonicalForm result= 0;
     370  int i= 0;
     371  int lf= 0;
     372  int lg= d*k;
     373  int degfSubLf= degf;
     374  int deggSubLg= degg-lg;
     375  int repLengthBuf2;
     376  int repLengthBuf1;
     377  int ii;
     378  zz_pE zzpEZero= zz_pE();
     379
     380  while (degf >= lf || lg >= 0)
     381  {
     382    if (degfSubLf >= d)
     383      repLengthBuf1= d;
     384    else if (degfSubLf < 0)
     385      repLengthBuf1= 0;
     386    else
     387      repLengthBuf1= degfSubLf + 1;
     388    buf1.rep.SetLength((long) repLengthBuf1);
     389
     390    buf1p= buf1.rep.elts();
     391    for (int ind= 0; ind < repLengthBuf1; ind++)
     392      buf1p [ind]= fp [ind + lf];
     393    buf1.normalize();
     394
     395    repLengthBuf1= buf1.rep.length();
     396
     397    if (deggSubLg >= d - 1)
     398      repLengthBuf2= d - 1;
     399    else if (deggSubLg < 0)
     400      repLengthBuf2= 0;
     401    else
     402      repLengthBuf2= deggSubLg + 1;
     403
     404    buf2.rep.SetLength ((long) repLengthBuf2);
     405    buf2p= buf2.rep.elts();
     406    for (int ind= 0; ind < repLengthBuf2; ind++)
     407    {
     408      buf2p [ind]= gp [ind + lg];
     409    }
     410    buf2.normalize();
     411
     412    repLengthBuf2= buf2.rep.length();
     413
     414    buf3.rep.SetLength((long) repLengthBuf2 + d);
     415    buf3p= buf3.rep.elts();
     416    buf2p= buf2.rep.elts();
     417    buf1p= buf1.rep.elts();
     418    for (int ind= 0; ind < repLengthBuf1; ind++)
     419      buf3p [ind]= buf1p [ind];
     420    for (int ind= repLengthBuf1; ind < d; ind++)
     421      buf3p [ind]= zzpEZero;
     422    for (int ind= 0; ind < repLengthBuf2; ind++)
     423      buf3p [ind + d]= buf2p [ind];
     424    buf3.normalize();
     425
     426    result += convertNTLzz_pEX2CF (buf3, x, alpha)*power (y, i);
     427    i++;
     428
     429
     430    lf= i*d;
     431    degfSubLf= degf - lf;
     432
     433    lg= d*(k-i);
     434    deggSubLg= degg - lg;
     435
     436    buf1p= buf1.rep.elts();
     437
     438    if (lg >= 0 && deggSubLg > 0)
     439    {
     440      if (repLengthBuf2 > degfSubLf + 1)
     441        degfSubLf= repLengthBuf2 - 1;
     442      int tmp= tmin (repLengthBuf1, deggSubLg);
     443      for (int ind= 0; ind < tmp; ind++)
     444        gp [ind + lg] -= buf1p [ind];
     445    }
     446
     447    if (lg < 0)
     448      break;
     449
     450    buf2p= buf2.rep.elts();
     451    if (degfSubLf >= 0)
     452    {
     453      for (int ind= 0; ind < repLengthBuf2; ind++)
     454        fp [ind + lf] -= buf2p [ind];
     455    }
     456  }
     457
     458  return result;
     459}
     460
     461CanonicalForm
     462reverseSubst (const zz_pX& F, const zz_pX& G, int d, int k)
     463{
     464  Variable y= Variable (2);
     465  Variable x= Variable (1);
     466
     467  zz_pX f= F;
     468  zz_pX g= G;
     469  int degf= deg(f);
     470  int degg= deg(g);
     471
     472  zz_pX buf1;
     473  zz_pX buf2;
     474  zz_pX buf3;
     475
     476  zz_p *buf1p;
     477  zz_p *buf2p;
     478  zz_p *buf3p;
     479
     480  if (f.rep.length() < (long) d*(k+1)) //zero padding
     481    f.rep.SetLength ((long)d*(k+1));
     482
     483  zz_p *gp= g.rep.elts();
     484  zz_p *fp= f.rep.elts();
     485  CanonicalForm result= 0;
     486  int i= 0;
     487  int lf= 0;
     488  int lg= d*k;
     489  int degfSubLf= degf;
     490  int deggSubLg= degg-lg;
     491  int repLengthBuf2;
     492  int repLengthBuf1;
     493  int ii;
     494  zz_p zzpZero= zz_p();
     495  while (degf >= lf || lg >= 0)
     496  {
     497    if (degfSubLf >= d)
     498      repLengthBuf1= d;
     499    else if (degfSubLf < 0)
     500      repLengthBuf1= 0;
     501    else
     502      repLengthBuf1= degfSubLf + 1;
     503    buf1.rep.SetLength((long) repLengthBuf1);
     504
     505    buf1p= buf1.rep.elts();
     506    for (int ind= 0; ind < repLengthBuf1; ind++)
     507      buf1p [ind]= fp [ind + lf];
     508    buf1.normalize();
     509
     510    repLengthBuf1= buf1.rep.length();
     511   
     512    if (deggSubLg >= d - 1)
     513      repLengthBuf2= d - 1;
     514    else if (deggSubLg < 0)
     515      repLengthBuf2= 0;
     516    else
     517      repLengthBuf2= deggSubLg + 1;
     518
     519    buf2.rep.SetLength ((long) repLengthBuf2);
     520    buf2p= buf2.rep.elts();
     521    for (int ind= 0; ind < repLengthBuf2; ind++)
     522      buf2p [ind]= gp [ind + lg];
     523
     524    buf2.normalize();
     525
     526    repLengthBuf2= buf2.rep.length();
     527
     528
     529    buf3.rep.SetLength((long) repLengthBuf2 + d);
     530    buf3p= buf3.rep.elts();
     531    buf2p= buf2.rep.elts();
     532    buf1p= buf1.rep.elts();
     533    for (int ind= 0; ind < repLengthBuf1; ind++)
     534      buf3p [ind]= buf1p [ind];
     535    for (int ind= repLengthBuf1; ind < d; ind++)
     536      buf3p [ind]= zzpZero;
     537    for (int ind= 0; ind < repLengthBuf2; ind++)
     538      buf3p [ind + d]= buf2p [ind];
     539    buf3.normalize();
     540
     541    result += convertNTLzzpX2CF (buf3, x)*power (y, i);
     542    i++;
     543
     544
     545    lf= i*d;
     546    degfSubLf= degf - lf;
     547
     548    lg= d*(k-i);
     549    deggSubLg= degg - lg;
     550
     551    buf1p= buf1.rep.elts();
     552
     553    if (lg >= 0 && deggSubLg > 0)
     554    {
     555      if (repLengthBuf2 > degfSubLf + 1)
     556        degfSubLf= repLengthBuf2 - 1;
     557      int tmp= tmin (repLengthBuf1, deggSubLg);
     558      for (int ind= 0; ind < tmp; ind++)
     559        gp [ind + lg] -= buf1p [ind];
     560    }
     561    if (lg < 0)
     562      break;
     563
     564    buf2p= buf2.rep.elts();
     565    if (degfSubLf >= 0)
     566    {
     567      for (int ind= 0; ind < repLengthBuf2; ind++)
     568        fp [ind + lf] -= buf2p [ind];
     569    }
     570  }
     571
     572  return result;
     573}
     574
     575CanonicalForm reverseSubst (const zz_pEX& F, int d, const Variable& alpha)
     576{
     577  Variable y= Variable (2);
     578  Variable x= Variable (1);
     579
     580  zz_pEX f= F;
     581  zz_pE *fp= f.rep.elts();
     582
     583  zz_pEX buf;
     584  zz_pE *bufp;
     585  CanonicalForm result= 0;
     586  int i= 0;
     587  int degf= deg(f);
     588  int k= 0;
     589  int degfSubK;
     590  int repLength;
     591  while (degf >= k)
     592  {
     593    degfSubK= degf - k;
     594    if (degfSubK >= d)
     595      repLength= d;
     596    else
     597      repLength= degfSubK + 1;
     598
     599    buf.rep.SetLength ((long) repLength);
     600    bufp= buf.rep.elts();
     601    for (int j= 0; j < repLength; j++)
     602      bufp [j]= fp [j + k];
     603    buf.normalize();
     604
     605    result += convertNTLzz_pEX2CF (buf, x, alpha)*power (y, i);
     606    i++;
     607    k= d*i;
     608  }
     609
     610  return result;
     611}
     612
     613CanonicalForm reverseSubstFp (const zz_pX& F, int d)
     614{
     615  Variable y= Variable (2);
     616  Variable x= Variable (1);
     617
     618  zz_pX f= F;
     619  zz_p *fp= f.rep.elts();
     620
     621  zz_pX buf;
     622  zz_p *bufp;
     623  CanonicalForm result= 0;
     624  int i= 0;
     625  int degf= deg(f);
     626  int k= 0;
     627  int degfSubK;
     628  int repLength;
     629  while (degf >= k)
     630  {
     631    degfSubK= degf - k;
     632    if (degfSubK >= d)
     633      repLength= d;
     634    else
     635      repLength= degfSubK + 1;
     636
     637    buf.rep.SetLength ((long) repLength);
     638    bufp= buf.rep.elts();
     639    for (int j= 0; j < repLength; j++)
     640      bufp [j]= fp [j + k];
     641    buf.normalize();
     642
     643    result += convertNTLzzpX2CF (buf, x)*power (y, i);
     644    i++;
     645    k= d*i;
     646  }
     647
     648  return result;
     649}
     650
     651// assumes input to be reduced mod M and to be an element of Fq not Fp
     652CanonicalForm
     653mulMod2NTLFpReci (const CanonicalForm& F, const CanonicalForm& G, const
     654                  CanonicalForm& M)
     655{
     656  int d1= tmax (degree (F, 1), degree (G, 1)) + 1;
     657  int d2= tmax (degree (F, 2), degree (G, 2));
     658
     659  zz_pX F1, F2;
     660  kronSubRecipro (F1, F2, F, d1);
     661  zz_pX G1, G2;
     662  kronSubRecipro (G1, G2, G, d1);
     663
     664  int k= d1*degree (M);
     665  MulTrunc (F1, F1, G1, (long) k);
     666
     667  mul (F2, F2, G2);
     668  F2 >>= k;
     669
     670  return reverseSubst (F1, F2, d1, d2);
     671}
     672
     673//Kronecker substitution
     674CanonicalForm
     675mulMod2NTLFp (const CanonicalForm& F, const CanonicalForm& G, const
     676              CanonicalForm& M)
     677{
     678  CanonicalForm A= F;
     679  CanonicalForm B= G;
     680
     681  int p= getCharacteristic();
     682
     683  int degAx= degree (A, 1);
     684  int degAy= degree (A, 2);
     685  int degBx= degree (B, 1);
     686  int degBy= degree (B, 2);
     687  int d1= degAx + 1 + degBx;
     688  int d2= tmax (degree (A, 2), degree (B, 2));
     689
     690  if (d1 > 128 && d2 > 160 && (degAy == degBy))
     691    return mulMod2NTLFpReci (A, B, M);
     692
     693  zz_pX NTLA= kronSubFp (A, d1);
     694  zz_pX NTLB= kronSubFp (B, d1);
     695
     696  int k= d1*degree (M);
     697  MulTrunc (NTLA, NTLA, NTLB, (long) k);
     698
     699  A= reverseSubstFp (NTLA, d1);
     700
     701  return A;
     702}
     703
     704// assumes input to be reduced mod M and to be an element of Fq not Fp
     705CanonicalForm
     706mulMod2NTLFqReci (const CanonicalForm& F, const CanonicalForm& G, const
     707                  CanonicalForm& M, const Variable& alpha)
     708{
     709  int d1= tmax (degree (F, 1), degree (G, 1)) + 1;
     710  int d2= tmax (degree (F, 2), degree (G, 2));
     711
     712  zz_pEX F1, F2;
     713  kronSubRecipro (F1, F2, F, d1, alpha);
     714  zz_pEX G1, G2;
     715  kronSubRecipro (G1, G2, G, d1, alpha);
     716
     717  int k1= d1*degree (M);
     718  MulTrunc (F1, F1, G1, (long) k1);
     719
     720  mul (F2, F2, G2);
     721
     722  F2 >>= k1;
     723
     724  CanonicalForm result= reverseSubst (F1, F2, d1, d2, alpha);
     725
     726  return result;
     727}
     728
     729CanonicalForm
     730mulMod2NTLFq (const CanonicalForm& F, const CanonicalForm& G, const
     731              CanonicalForm& M)
     732{
     733  Variable alpha;
     734  CanonicalForm A= F;
     735  CanonicalForm B= G;
     736
     737  if (hasFirstAlgVar (A, alpha) || hasFirstAlgVar (B, alpha))
     738  {
     739    int degAx= degree (A, 1);
     740    int degAy= degree (A, 2);
     741    int degBx= degree (B, 1);
     742    int degBy= degree (B, 2);
     743    int d1= degAx + degBx + 1;
     744    int d2= tmax (degree (A, 2), degree (B, 2));
     745    zz_p::init (getCharacteristic());
     746    zz_pX NTLMipo= convertFacCF2NTLzzpX (getMipo (alpha));
     747    zz_pE::init (NTLMipo);
     748
     749    int degMipo= degree (getMipo (alpha));
     750    if ((d1 > 128/degMipo) && (d2 > 160/degMipo) && (degAy == degBy))
     751      return mulMod2NTLFqReci (A, B, M, alpha);
     752
     753    zz_pEX NTLA= kronSub (A, d1, alpha);
     754    zz_pEX NTLB= kronSub (B, d1, alpha);
     755
     756    int k= d1*degree (M);
     757
     758    MulTrunc (NTLA, NTLA, NTLB, (long) k);
     759
     760    A= reverseSubst (NTLA, d1, alpha);
     761
     762    return A;
     763  }
     764  else
     765    return mulMod2NTLFp (A, B, M);
     766}
     767
    206768CanonicalForm mulMod2 (const CanonicalForm& A, const CanonicalForm& B,
    207769                       const CanonicalForm& M)
     
    220782  int degG= degree (G, y);
    221783
    222   if ((degF < 1 && degG < 1) && (F.isUnivariate() && G.isUnivariate()) &&
     784  if ((degF < 1 && degG < 1) && (F.isUnivariate() && G.isUnivariate()) && 
    223785      (F.level() == G.level()))
    224   {
     786  { 
    225787    CanonicalForm result= mulNTL (F, G);
    226788    return mod (result, M);
    227789  }
    228790  else if (degF <= 1 && degG <= 1)
    229   {
     791  { 
    230792    CanonicalForm result= F*G;
    231793    return mod (result, M);
    232794  }
     795
     796  int sizeF= size (F);
     797  int sizeG= size (G);
     798 
     799  int fallBackToNaive= 50;
     800  if (sizeF < fallBackToNaive || sizeG < fallBackToNaive)
     801    return mod (F*G, M);
     802
     803  if (CFFactory::gettype() != GaloisFieldDomain &&
     804      (((degF-degG) < 50 && degF > degG) || ((degG-degF) < 50 && degF <= degG)))
     805    return mulMod2NTLFq (F, G, M);
    233806
    234807  int m= (int) ceil (degree (M)/2.0);
     
    404977    return F;
    405978  CanonicalForm A= F;
    406   Variable y= F.mvar();
     979  Variable y= Variable (2);
    407980  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;
     981  if (degree (A, x) > 0)
     982  {
     983    A= swapvar (A, x, y);
     984    CanonicalForm result= 0;
     985    CFIterator i= A;
     986    while (d - i.exp() < 0)
     987      i++;
     988
     989    for (; i.hasTerms() && (d - i.exp() >= 0); i++)
     990      result += swapvar (i.coeff(),x,y)*power (x, d - i.exp());
     991    return result;
     992  }
     993  else
     994    return A*power (x, d);
    417995}
    418996
     
    4901068
    4911069  CanonicalForm Q;
    492   if (degB <= 1)
     1070  if (degB <= 1 || CFFactory::gettype() == GaloisFieldDomain)
    4931071  {
    4941072    CanonicalForm R;
     
    5261104  }
    5271105
    528   if (degB <= 1)
     1106  if (degB <= 1 || CFFactory::gettype() == GaloisFieldDomain)
    5291107  {
    5301108     divrem2 (A, B, Q, R, M);
Note: See TracChangeset for help on using the changeset viewer.