PageRenderTime 230ms CodeModel.GetById 20ms app.highlight 194ms RepoModel.GetById 1ms app.codeStats 1ms

/factory/facHensel.cc

https://github.com/hannes14/Singular
C++ | 3225 lines | 2827 code | 299 blank | 99 comment | 745 complexity | 57700a23d586cde3ee6a35c606de1422 MD5 | raw file

Large files files are truncated, but you can click here to view the full file

   1/*****************************************************************************\
   2 * Computer Algebra System SINGULAR
   3\*****************************************************************************/
   4/** @file facHensel.cc
   5 *
   6 * This file implements functions to lift factors via Hensel lifting and
   7 * functions for modular multiplication and division with remainder.
   8 *
   9 * ABSTRACT: Hensel lifting is described in "Efficient Multivariate
  10 * Factorization over Finite Fields" by L. Bernardin & M. Monagon. Division with
  11 * remainder is described in "Fast Recursive Division" by C. Burnikel and
  12 * J. Ziegler. Karatsuba multiplication is described in "Modern Computer
  13 * Algebra" by J. von zur Gathen and J. Gerhard.
  14 *
  15 * @author Martin Lee
  16 *
  17 * @internal @version \$Id$
  18 *
  19 **/
  20/*****************************************************************************/
  21
  22#include "assert.h"
  23#include "debug.h"
  24#include "timing.h"
  25
  26#include "facHensel.h"
  27#include "cf_util.h"
  28#include "fac_util.h"
  29#include "cf_algorithm.h"
  30
  31#ifdef HAVE_NTL
  32#include <NTL/lzz_pEX.h>
  33#include "NTLconvert.h"
  34
  35CanonicalForm
  36mulNTL (const CanonicalForm& F, const CanonicalForm& G)
  37{
  38  if (F.inCoeffDomain() || G.inCoeffDomain() || getCharacteristic() == 0)
  39    return F*G;
  40  ASSERT (F.isUnivariate() && G.isUnivariate(), "expected univariate polys");
  41  ASSERT (F.level() == G.level(), "expected polys of same level");
  42  if (CFFactory::gettype() == GaloisFieldDomain)
  43    return F*G;
  44  zz_p::init (getCharacteristic());
  45  Variable alpha;
  46  CanonicalForm result;
  47  if (hasFirstAlgVar (F, alpha) || hasFirstAlgVar (G, alpha))
  48  {
  49    zz_pX NTLMipo= convertFacCF2NTLzzpX (getMipo (alpha));
  50    zz_pE::init (NTLMipo);
  51    zz_pEX NTLF= convertFacCF2NTLzz_pEX (F, NTLMipo);
  52    zz_pEX NTLG= convertFacCF2NTLzz_pEX (G, NTLMipo);
  53    mul (NTLF, NTLF, NTLG);
  54    result= convertNTLzz_pEX2CF(NTLF, F.mvar(), alpha);
  55  }
  56  else
  57  {
  58    zz_pX NTLF= convertFacCF2NTLzzpX (F);
  59    zz_pX NTLG= convertFacCF2NTLzzpX (G);
  60    mul (NTLF, NTLF, NTLG);
  61    result= convertNTLzzpX2CF(NTLF, F.mvar());
  62  }
  63  return result;
  64}
  65
  66CanonicalForm
  67modNTL (const CanonicalForm& F, const CanonicalForm& G)
  68{
  69  if (F.inCoeffDomain() && G.isUnivariate())
  70    return F;
  71  else if (F.inCoeffDomain() && G.inCoeffDomain())
  72    return mod (F, G);
  73  else if (F.isUnivariate() && G.inCoeffDomain())
  74    return mod (F,G);
  75
  76  if (getCharacteristic() == 0)
  77    return mod (F, G);
  78
  79  ASSERT (F.isUnivariate() && G.isUnivariate(), "expected univariate polys");
  80  ASSERT (F.level() == G.level(), "expected polys of same level");
  81  if (CFFactory::gettype() == GaloisFieldDomain)
  82    return mod (F, G);
  83  zz_p::init (getCharacteristic());
  84  Variable alpha;
  85  CanonicalForm result;
  86  if (hasFirstAlgVar (F, alpha) || hasFirstAlgVar (G, alpha))
  87  {
  88    zz_pX NTLMipo= convertFacCF2NTLzzpX(getMipo (alpha));
  89    zz_pE::init (NTLMipo);
  90    zz_pEX NTLF= convertFacCF2NTLzz_pEX (F, NTLMipo);
  91    zz_pEX NTLG= convertFacCF2NTLzz_pEX (G, NTLMipo);
  92    rem (NTLF, NTLF, NTLG);
  93    result= convertNTLzz_pEX2CF(NTLF, F.mvar(), alpha);
  94  }
  95  else
  96  {
  97    zz_pX NTLF= convertFacCF2NTLzzpX (F);
  98    zz_pX NTLG= convertFacCF2NTLzzpX (G);
  99    rem (NTLF, NTLF, NTLG);
 100    result= convertNTLzzpX2CF(NTLF, F.mvar());
 101  }
 102  return result;
 103}
 104
 105CanonicalForm
 106divNTL (const CanonicalForm& F, const CanonicalForm& G)
 107{
 108  if (F.inCoeffDomain() && G.isUnivariate())
 109    return F;
 110  else if (F.inCoeffDomain() && G.inCoeffDomain())
 111    return div (F, G);
 112  else if (F.isUnivariate() && G.inCoeffDomain())
 113    return div (F,G);
 114
 115  if (getCharacteristic() == 0)
 116    return div (F, G);
 117
 118  ASSERT (F.isUnivariate() && G.isUnivariate(), "expected univariate polys");
 119  ASSERT (F.level() == G.level(), "expected polys of same level");
 120  if (CFFactory::gettype() == GaloisFieldDomain)
 121    return div (F, G);
 122  zz_p::init (getCharacteristic());
 123  Variable alpha;
 124  CanonicalForm result;
 125  if (hasFirstAlgVar (F, alpha) || hasFirstAlgVar (G, alpha))
 126  {
 127    zz_pX NTLMipo= convertFacCF2NTLzzpX(getMipo (alpha));
 128    zz_pE::init (NTLMipo);
 129    zz_pEX NTLF= convertFacCF2NTLzz_pEX (F, NTLMipo);
 130    zz_pEX NTLG= convertFacCF2NTLzz_pEX (G, NTLMipo);
 131    div (NTLF, NTLF, NTLG);
 132    result= convertNTLzz_pEX2CF(NTLF, F.mvar(), alpha);
 133  }
 134  else
 135  {
 136    zz_pX NTLF= convertFacCF2NTLzzpX (F);
 137    zz_pX NTLG= convertFacCF2NTLzzpX (G);
 138    div (NTLF, NTLF, NTLG);
 139    result= convertNTLzzpX2CF(NTLF, F.mvar());
 140  }
 141  return result;
 142}
 143
 144/*
 145void
 146divremNTL (const CanonicalForm& F, const CanonicalForm& G, CanonicalForm& Q,
 147           CanonicalForm& R)
 148{
 149  if (F.inCoeffDomain() && G.isUnivariate())
 150  {
 151    R= F;
 152    Q= 0;
 153  }
 154  else if (F.inCoeffDomain() && G.inCoeffDomain())
 155  {
 156    divrem (F, G, Q, R);
 157    return;
 158  }
 159  else if (F.isUnivariate() && G.inCoeffDomain())
 160  {
 161    divrem (F, G, Q, R);
 162    return;
 163  }
 164
 165  ASSERT (F.isUnivariate() && G.isUnivariate(), "expected univariate polys");
 166  ASSERT (F.level() == G.level(), "expected polys of same level");
 167  if (CFFactory::gettype() == GaloisFieldDomain)
 168  {
 169    divrem (F, G, Q, R);
 170    return;
 171  }
 172  zz_p::init (getCharacteristic());
 173  Variable alpha;
 174  CanonicalForm result;
 175  if (hasFirstAlgVar (F, alpha) || hasFirstAlgVar (G, alpha))
 176  {
 177    zz_pX NTLMipo= convertFacCF2NTLzzpX(getMipo (alpha));
 178    zz_pE::init (NTLMipo);
 179    zz_pEX NTLF= convertFacCF2NTLzz_pEX (F, NTLMipo);
 180    zz_pEX NTLG= convertFacCF2NTLzz_pEX (G, NTLMipo);
 181    zz_pEX NTLQ;
 182    zz_pEX NTLR;
 183    DivRem (NTLQ, NTLR, NTLF, NTLG);
 184    Q= convertNTLzz_pEX2CF(NTLQ, F.mvar(), alpha);
 185    R= convertNTLzz_pEX2CF(NTLR, F.mvar(), alpha);
 186    return;
 187  }
 188  else
 189  {
 190    zz_pX NTLF= convertFacCF2NTLzzpX (F);
 191    zz_pX NTLG= convertFacCF2NTLzzpX (G);
 192    zz_pX NTLQ;
 193    zz_pX NTLR;
 194    DivRem (NTLQ, NTLR, NTLF, NTLG);
 195    Q= convertNTLzzpX2CF(NTLQ, F.mvar());
 196    R= convertNTLzzpX2CF(NTLR, F.mvar());
 197    return;
 198  }
 199}*/
 200
 201CanonicalForm mod (const CanonicalForm& F, const CFList& M)
 202{
 203  CanonicalForm A= F;
 204  for (CFListIterator i= M; i.hasItem(); i++)
 205    A= mod (A, i.getItem());
 206  return A;
 207}
 208
 209zz_pX kronSubFp (const CanonicalForm& A, int d)
 210{
 211  int degAy= degree (A);
 212  zz_pX result;
 213  result.rep.SetLength (d*(degAy + 1));
 214
 215  zz_p *resultp;
 216  resultp= result.rep.elts();
 217  zz_pX buf;
 218  zz_p *bufp;
 219
 220  for (CFIterator i= A; i.hasTerms(); i++)
 221  {
 222    if (i.coeff().inCoeffDomain())
 223      buf= convertFacCF2NTLzzpX (i.coeff());
 224    else
 225      buf= convertFacCF2NTLzzpX (i.coeff());
 226
 227    int k= i.exp()*d;
 228    bufp= buf.rep.elts();
 229    int bufRepLength= (int) buf.rep.length();
 230    for (int j= 0; j < bufRepLength; j++)
 231      resultp [j + k]= bufp [j];
 232  }
 233  result.normalize();
 234
 235  return result;
 236}
 237
 238zz_pEX kronSub (const CanonicalForm& A, int d, const Variable& alpha)
 239{
 240  int degAy= degree (A);
 241  zz_pEX result;
 242  result.rep.SetLength (d*(degAy + 1));
 243
 244  Variable v;
 245  zz_pE *resultp;
 246  resultp= result.rep.elts();
 247  zz_pEX buf1;
 248  zz_pE *buf1p;
 249  zz_pX buf2;
 250  zz_pX NTLMipo= convertFacCF2NTLzzpX (getMipo (alpha));
 251
 252  for (CFIterator i= A; i.hasTerms(); i++)
 253  {
 254    if (i.coeff().inCoeffDomain())
 255    {
 256      buf2= convertFacCF2NTLzzpX (i.coeff());
 257      buf1= to_zz_pEX (to_zz_pE (buf2));
 258    }
 259    else
 260      buf1= convertFacCF2NTLzz_pEX (i.coeff(), NTLMipo);
 261
 262    int k= i.exp()*d;
 263    buf1p= buf1.rep.elts();
 264    int buf1RepLength= (int) buf1.rep.length();
 265    for (int j= 0; j < buf1RepLength; j++)
 266      resultp [j + k]= buf1p [j];
 267  }
 268  result.normalize();
 269
 270  return result;
 271}
 272
 273void
 274kronSubRecipro (zz_pEX& subA1, zz_pEX& subA2,const CanonicalForm& A, int d,
 275                const Variable& alpha)
 276{
 277  int degAy= degree (A);
 278  subA1.rep.SetLength ((long) d*(degAy + 1));
 279  subA2.rep.SetLength ((long) d*(degAy + 1));
 280
 281  Variable v;
 282  zz_pE *subA1p;
 283  zz_pE *subA2p;
 284  subA1p= subA1.rep.elts();
 285  subA2p= subA2.rep.elts();
 286  zz_pEX buf;
 287  zz_pE *bufp;
 288  zz_pX buf2;
 289  zz_pX NTLMipo= convertFacCF2NTLzzpX (getMipo (alpha));
 290
 291  for (CFIterator i= A; i.hasTerms(); i++)
 292  {
 293    if (i.coeff().inCoeffDomain())
 294    {
 295      buf2= convertFacCF2NTLzzpX (i.coeff());
 296      buf= to_zz_pEX (to_zz_pE (buf2));
 297    }
 298    else
 299      buf= convertFacCF2NTLzz_pEX (i.coeff(), NTLMipo);
 300
 301    int k= i.exp()*d;
 302    int kk= (degAy - i.exp())*d;
 303    bufp= buf.rep.elts();
 304    int bufRepLength= (int) buf.rep.length();
 305    for (int j= 0; j < bufRepLength; j++)
 306    {
 307      subA1p [j + k]= bufp [j];
 308      subA2p [j + kk]= bufp [j];
 309    }
 310  }
 311  subA1.normalize();
 312  subA2.normalize();
 313}
 314
 315void
 316kronSubRecipro (zz_pX& subA1, zz_pX& subA2,const CanonicalForm& A, int d)
 317{
 318  int degAy= degree (A);
 319  subA1.rep.SetLength ((long) d*(degAy + 1));
 320  subA2.rep.SetLength ((long) d*(degAy + 1));
 321
 322  Variable v;
 323  zz_p *subA1p;
 324  zz_p *subA2p;
 325  subA1p= subA1.rep.elts();
 326  subA2p= subA2.rep.elts();
 327  zz_pX buf;
 328  zz_p *bufp;
 329
 330  for (CFIterator i= A; i.hasTerms(); i++)
 331  {
 332    buf= convertFacCF2NTLzzpX (i.coeff());
 333
 334    int k= i.exp()*d;
 335    int kk= (degAy - i.exp())*d;
 336    bufp= buf.rep.elts();
 337    int bufRepLength= (int) buf.rep.length();
 338    for (int j= 0; j < bufRepLength; j++)
 339    {
 340      subA1p [j + k]= bufp [j];
 341      subA2p [j + kk]= bufp [j];
 342    }
 343  }
 344  subA1.normalize();
 345  subA2.normalize();
 346}
 347
 348CanonicalForm
 349reverseSubst (const zz_pEX& F, const zz_pEX& G, int d, int k,
 350              const Variable& alpha)
 351{
 352  Variable y= Variable (2);
 353  Variable x= Variable (1);
 354
 355  zz_pEX f= F;
 356  zz_pEX g= G;
 357  int degf= deg(f);
 358  int degg= deg(g);
 359
 360  zz_pEX buf1;
 361  zz_pEX buf2;
 362  zz_pEX buf3;
 363
 364  zz_pE *buf1p;
 365  zz_pE *buf2p;
 366  zz_pE *buf3p;
 367  if (f.rep.length() < (long) d*(k+1)) //zero padding
 368    f.rep.SetLength ((long)d*(k+1));
 369
 370  zz_pE *gp= g.rep.elts();
 371  zz_pE *fp= f.rep.elts();
 372  CanonicalForm result= 0;
 373  int i= 0;
 374  int lf= 0;
 375  int lg= d*k;
 376  int degfSubLf= degf;
 377  int deggSubLg= degg-lg;
 378  int repLengthBuf2;
 379  int repLengthBuf1;
 380  zz_pE zzpEZero= zz_pE();
 381
 382  while (degf >= lf || lg >= 0)
 383  {
 384    if (degfSubLf >= d)
 385      repLengthBuf1= d;
 386    else if (degfSubLf < 0)
 387      repLengthBuf1= 0;
 388    else
 389      repLengthBuf1= degfSubLf + 1;
 390    buf1.rep.SetLength((long) repLengthBuf1);
 391
 392    buf1p= buf1.rep.elts();
 393    for (int ind= 0; ind < repLengthBuf1; ind++)
 394      buf1p [ind]= fp [ind + lf];
 395    buf1.normalize();
 396
 397    repLengthBuf1= buf1.rep.length();
 398
 399    if (deggSubLg >= d - 1)
 400      repLengthBuf2= d - 1;
 401    else if (deggSubLg < 0)
 402      repLengthBuf2= 0;
 403    else
 404      repLengthBuf2= deggSubLg + 1;
 405
 406    buf2.rep.SetLength ((long) repLengthBuf2);
 407    buf2p= buf2.rep.elts();
 408    for (int ind= 0; ind < repLengthBuf2; ind++)
 409    {
 410      buf2p [ind]= gp [ind + lg];
 411    }
 412    buf2.normalize();
 413
 414    repLengthBuf2= buf2.rep.length();
 415
 416    buf3.rep.SetLength((long) repLengthBuf2 + d);
 417    buf3p= buf3.rep.elts();
 418    buf2p= buf2.rep.elts();
 419    buf1p= buf1.rep.elts();
 420    for (int ind= 0; ind < repLengthBuf1; ind++)
 421      buf3p [ind]= buf1p [ind];
 422    for (int ind= repLengthBuf1; ind < d; ind++)
 423      buf3p [ind]= zzpEZero;
 424    for (int ind= 0; ind < repLengthBuf2; ind++)
 425      buf3p [ind + d]= buf2p [ind];
 426    buf3.normalize();
 427
 428    result += convertNTLzz_pEX2CF (buf3, x, alpha)*power (y, i);
 429    i++;
 430
 431
 432    lf= i*d;
 433    degfSubLf= degf - lf;
 434
 435    lg= d*(k-i);
 436    deggSubLg= degg - lg;
 437
 438    buf1p= buf1.rep.elts();
 439
 440    if (lg >= 0 && deggSubLg > 0)
 441    {
 442      if (repLengthBuf2 > degfSubLf + 1)
 443        degfSubLf= repLengthBuf2 - 1;
 444      int tmp= tmin (repLengthBuf1, deggSubLg);
 445      for (int ind= 0; ind < tmp; ind++)
 446        gp [ind + lg] -= buf1p [ind];
 447    }
 448
 449    if (lg < 0)
 450      break;
 451
 452    buf2p= buf2.rep.elts();
 453    if (degfSubLf >= 0)
 454    {
 455      for (int ind= 0; ind < repLengthBuf2; ind++)
 456        fp [ind + lf] -= buf2p [ind];
 457    }
 458  }
 459
 460  return result;
 461}
 462
 463CanonicalForm
 464reverseSubst (const zz_pX& F, const zz_pX& G, int d, int k)
 465{
 466  Variable y= Variable (2);
 467  Variable x= Variable (1);
 468
 469  zz_pX f= F;
 470  zz_pX g= G;
 471  int degf= deg(f);
 472  int degg= deg(g);
 473
 474  zz_pX buf1;
 475  zz_pX buf2;
 476  zz_pX buf3;
 477
 478  zz_p *buf1p;
 479  zz_p *buf2p;
 480  zz_p *buf3p;
 481
 482  if (f.rep.length() < (long) d*(k+1)) //zero padding
 483    f.rep.SetLength ((long)d*(k+1));
 484
 485  zz_p *gp= g.rep.elts();
 486  zz_p *fp= f.rep.elts();
 487  CanonicalForm result= 0;
 488  int i= 0;
 489  int lf= 0;
 490  int lg= d*k;
 491  int degfSubLf= degf;
 492  int deggSubLg= degg-lg;
 493  int repLengthBuf2;
 494  int repLengthBuf1;
 495  zz_p zzpZero= zz_p();
 496  while (degf >= lf || lg >= 0)
 497  {
 498    if (degfSubLf >= d)
 499      repLengthBuf1= d;
 500    else if (degfSubLf < 0)
 501      repLengthBuf1= 0;
 502    else
 503      repLengthBuf1= degfSubLf + 1;
 504    buf1.rep.SetLength((long) repLengthBuf1);
 505
 506    buf1p= buf1.rep.elts();
 507    for (int ind= 0; ind < repLengthBuf1; ind++)
 508      buf1p [ind]= fp [ind + lf];
 509    buf1.normalize();
 510
 511    repLengthBuf1= buf1.rep.length();
 512
 513    if (deggSubLg >= d - 1)
 514      repLengthBuf2= d - 1;
 515    else if (deggSubLg < 0)
 516      repLengthBuf2= 0;
 517    else
 518      repLengthBuf2= deggSubLg + 1;
 519
 520    buf2.rep.SetLength ((long) repLengthBuf2);
 521    buf2p= buf2.rep.elts();
 522    for (int ind= 0; ind < repLengthBuf2; ind++)
 523      buf2p [ind]= gp [ind + lg];
 524
 525    buf2.normalize();
 526
 527    repLengthBuf2= buf2.rep.length();
 528
 529
 530    buf3.rep.SetLength((long) repLengthBuf2 + d);
 531    buf3p= buf3.rep.elts();
 532    buf2p= buf2.rep.elts();
 533    buf1p= buf1.rep.elts();
 534    for (int ind= 0; ind < repLengthBuf1; ind++)
 535      buf3p [ind]= buf1p [ind];
 536    for (int ind= repLengthBuf1; ind < d; ind++)
 537      buf3p [ind]= zzpZero;
 538    for (int ind= 0; ind < repLengthBuf2; ind++)
 539      buf3p [ind + d]= buf2p [ind];
 540    buf3.normalize();
 541
 542    result += convertNTLzzpX2CF (buf3, x)*power (y, i);
 543    i++;
 544
 545
 546    lf= i*d;
 547    degfSubLf= degf - lf;
 548
 549    lg= d*(k-i);
 550    deggSubLg= degg - lg;
 551
 552    buf1p= buf1.rep.elts();
 553
 554    if (lg >= 0 && deggSubLg > 0)
 555    {
 556      if (repLengthBuf2 > degfSubLf + 1)
 557        degfSubLf= repLengthBuf2 - 1;
 558      int tmp= tmin (repLengthBuf1, deggSubLg);
 559      for (int ind= 0; ind < tmp; ind++)
 560        gp [ind + lg] -= buf1p [ind];
 561    }
 562    if (lg < 0)
 563      break;
 564
 565    buf2p= buf2.rep.elts();
 566    if (degfSubLf >= 0)
 567    {
 568      for (int ind= 0; ind < repLengthBuf2; ind++)
 569        fp [ind + lf] -= buf2p [ind];
 570    }
 571  }
 572
 573  return result;
 574}
 575
 576CanonicalForm reverseSubst (const zz_pEX& F, int d, const Variable& alpha)
 577{
 578  Variable y= Variable (2);
 579  Variable x= Variable (1);
 580
 581  zz_pEX f= F;
 582  zz_pE *fp= f.rep.elts();
 583
 584  zz_pEX buf;
 585  zz_pE *bufp;
 586  CanonicalForm result= 0;
 587  int i= 0;
 588  int degf= deg(f);
 589  int k= 0;
 590  int degfSubK;
 591  int repLength;
 592  while (degf >= k)
 593  {
 594    degfSubK= degf - k;
 595    if (degfSubK >= d)
 596      repLength= d;
 597    else
 598      repLength= degfSubK + 1;
 599
 600    buf.rep.SetLength ((long) repLength);
 601    bufp= buf.rep.elts();
 602    for (int j= 0; j < repLength; j++)
 603      bufp [j]= fp [j + k];
 604    buf.normalize();
 605
 606    result += convertNTLzz_pEX2CF (buf, x, alpha)*power (y, i);
 607    i++;
 608    k= d*i;
 609  }
 610
 611  return result;
 612}
 613
 614CanonicalForm reverseSubstFp (const zz_pX& F, int d)
 615{
 616  Variable y= Variable (2);
 617  Variable x= Variable (1);
 618
 619  zz_pX f= F;
 620  zz_p *fp= f.rep.elts();
 621
 622  zz_pX buf;
 623  zz_p *bufp;
 624  CanonicalForm result= 0;
 625  int i= 0;
 626  int degf= deg(f);
 627  int k= 0;
 628  int degfSubK;
 629  int repLength;
 630  while (degf >= k)
 631  {
 632    degfSubK= degf - k;
 633    if (degfSubK >= d)
 634      repLength= d;
 635    else
 636      repLength= degfSubK + 1;
 637
 638    buf.rep.SetLength ((long) repLength);
 639    bufp= buf.rep.elts();
 640    for (int j= 0; j < repLength; j++)
 641      bufp [j]= fp [j + k];
 642    buf.normalize();
 643
 644    result += convertNTLzzpX2CF (buf, x)*power (y, i);
 645    i++;
 646    k= d*i;
 647  }
 648
 649  return result;
 650}
 651
 652// assumes input to be reduced mod M and to be an element of Fq not Fp
 653CanonicalForm
 654mulMod2NTLFpReci (const CanonicalForm& F, const CanonicalForm& G, const
 655                  CanonicalForm& M)
 656{
 657  int d1= tmax (degree (F, 1), degree (G, 1)) + 1;
 658  int d2= tmax (degree (F, 2), degree (G, 2));
 659
 660  zz_pX F1, F2;
 661  kronSubRecipro (F1, F2, F, d1);
 662  zz_pX G1, G2;
 663  kronSubRecipro (G1, G2, G, d1);
 664
 665  int k= d1*degree (M);
 666  MulTrunc (F1, F1, G1, (long) k);
 667
 668  mul (F2, F2, G2);
 669  if (deg (F2) > k)
 670    F2 >>= (k - d1);
 671
 672  return reverseSubst (F1, F2, d1, d2);
 673}
 674
 675//Kronecker substitution
 676CanonicalForm
 677mulMod2NTLFp (const CanonicalForm& F, const CanonicalForm& G, const
 678              CanonicalForm& M)
 679{
 680  CanonicalForm A= F;
 681  CanonicalForm B= G;
 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) && (2*degAy > degree (M)))
 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 k= d1*degree (M);
 718  MulTrunc (F1, F1, G1, (long) k);
 719
 720  mul (F2, F2, G2);
 721  if (deg (F2) > k)
 722    F2 >>= (k-d1);
 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        (2*degAy > degree (M)))
 752      return mulMod2NTLFqReci (A, B, M, alpha);
 753
 754    zz_pEX NTLA= kronSub (A, d1, alpha);
 755    zz_pEX NTLB= kronSub (B, d1, alpha);
 756
 757    int k= d1*degree (M);
 758
 759    MulTrunc (NTLA, NTLA, NTLB, (long) k);
 760
 761    A= reverseSubst (NTLA, d1, alpha);
 762
 763    return A;
 764  }
 765  else
 766    return mulMod2NTLFp (A, B, M);
 767}
 768
 769CanonicalForm mulMod2 (const CanonicalForm& A, const CanonicalForm& B,
 770                       const CanonicalForm& M)
 771{
 772  if (A.isZero() || B.isZero())
 773    return 0;
 774
 775  ASSERT (M.isUnivariate(), "M must be univariate");
 776
 777  CanonicalForm F= mod (A, M);
 778  CanonicalForm G= mod (B, M);
 779  if (F.inCoeffDomain() || G.inCoeffDomain())
 780    return F*G;
 781  Variable y= M.mvar();
 782  int degF= degree (F, y);
 783  int degG= degree (G, y);
 784
 785  if ((degF < 1 && degG < 1) && (F.isUnivariate() && G.isUnivariate()) &&
 786      (F.level() == G.level()))
 787  {
 788    CanonicalForm result= mulNTL (F, G);
 789    return mod (result, M);
 790  }
 791  else if (degF <= 1 && degG <= 1)
 792  {
 793    CanonicalForm result= F*G;
 794    return mod (result, M);
 795  }
 796
 797  int sizeF= size (F);
 798  int sizeG= size (G);
 799
 800  int fallBackToNaive= 50;
 801  if (sizeF < fallBackToNaive || sizeG < fallBackToNaive)
 802    return mod (F*G, M);
 803
 804  if (getCharacteristic() > 0 && CFFactory::gettype() != GaloisFieldDomain &&
 805      (((degF-degG) < 50 && degF > degG) || ((degG-degF) < 50 && degF <= degG)))
 806    return mulMod2NTLFq (F, G, M);
 807
 808  int m= (int) ceil (degree (M)/2.0);
 809  if (degF >= m || degG >= m)
 810  {
 811    CanonicalForm MLo= power (y, m);
 812    CanonicalForm MHi= power (y, degree (M) - m);
 813    CanonicalForm F0= mod (F, MLo);
 814    CanonicalForm F1= div (F, MLo);
 815    CanonicalForm G0= mod (G, MLo);
 816    CanonicalForm G1= div (G, MLo);
 817    CanonicalForm F0G1= mulMod2 (F0, G1, MHi);
 818    CanonicalForm F1G0= mulMod2 (F1, G0, MHi);
 819    CanonicalForm F0G0= mulMod2 (F0, G0, M);
 820    return F0G0 + MLo*(F0G1 + F1G0);
 821  }
 822  else
 823  {
 824    m= (int) ceil (tmax (degF, degG)/2.0);
 825    CanonicalForm yToM= power (y, m);
 826    CanonicalForm F0= mod (F, yToM);
 827    CanonicalForm F1= div (F, yToM);
 828    CanonicalForm G0= mod (G, yToM);
 829    CanonicalForm G1= div (G, yToM);
 830    CanonicalForm H00= mulMod2 (F0, G0, M);
 831    CanonicalForm H11= mulMod2 (F1, G1, M);
 832    CanonicalForm H01= mulMod2 (F0 + F1, G0 + G1, M);
 833    return H11*yToM*yToM + (H01 - H11 - H00)*yToM + H00;
 834  }
 835  DEBOUTLN (cerr, "fatal end in mulMod2");
 836}
 837
 838CanonicalForm mulMod (const CanonicalForm& A, const CanonicalForm& B,
 839                      const CFList& MOD)
 840{
 841  if (A.isZero() || B.isZero())
 842    return 0;
 843
 844  if (MOD.length() == 1)
 845    return mulMod2 (A, B, MOD.getLast());
 846
 847  CanonicalForm M= MOD.getLast();
 848  CanonicalForm F= mod (A, M);
 849  CanonicalForm G= mod (B, M);
 850  if (F.inCoeffDomain() || G.inCoeffDomain())
 851    return F*G;
 852  Variable y= M.mvar();
 853  int degF= degree (F, y);
 854  int degG= degree (G, y);
 855
 856  if ((degF <= 1 && F.level() <= M.level()) &&
 857      (degG <= 1 && G.level() <= M.level()))
 858  {
 859    CFList buf= MOD;
 860    buf.removeLast();
 861    if (degF == 1 && degG == 1)
 862    {
 863      CanonicalForm F0= mod (F, y);
 864      CanonicalForm F1= div (F, y);
 865      CanonicalForm G0= mod (G, y);
 866      CanonicalForm G1= div (G, y);
 867      if (degree (M) > 2)
 868      {
 869        CanonicalForm H00= mulMod (F0, G0, buf);
 870        CanonicalForm H11= mulMod (F1, G1, buf);
 871        CanonicalForm H01= mulMod (F0 + F1, G0 + G1, buf);
 872        return H11*y*y + (H01 - H00 - H11)*y + H00;
 873      }
 874      else //here degree (M) == 2
 875      {
 876        buf.append (y);
 877        CanonicalForm F0G1= mulMod (F0, G1, buf);
 878        CanonicalForm F1G0= mulMod (F1, G0, buf);
 879        CanonicalForm F0G0= mulMod (F0, G0, MOD);
 880        CanonicalForm result= F0G0 + y*(F0G1 + F1G0);
 881        return result;
 882      }
 883    }
 884    else if (degF == 1 && degG == 0)
 885      return mulMod (div (F, y), G, buf)*y + mulMod (mod (F, y), G, buf);
 886    else if (degF == 0 && degG == 1)
 887      return mulMod (div (G, y), F, buf)*y + mulMod (mod (G, y), F, buf);
 888    else
 889      return mulMod (F, G, buf);
 890  }
 891  int m= (int) ceil (degree (M)/2.0);
 892  if (degF >= m || degG >= m)
 893  {
 894    CanonicalForm MLo= power (y, m);
 895    CanonicalForm MHi= power (y, degree (M) - m);
 896    CanonicalForm F0= mod (F, MLo);
 897    CanonicalForm F1= div (F, MLo);
 898    CanonicalForm G0= mod (G, MLo);
 899    CanonicalForm G1= div (G, MLo);
 900    CFList buf= MOD;
 901    buf.removeLast();
 902    buf.append (MHi);
 903    CanonicalForm F0G1= mulMod (F0, G1, buf);
 904    CanonicalForm F1G0= mulMod (F1, G0, buf);
 905    CanonicalForm F0G0= mulMod (F0, G0, MOD);
 906    return F0G0 + MLo*(F0G1 + F1G0);
 907  }
 908  else
 909  {
 910    m= (int) ceil (tmax (degF, degG)/2.0);
 911    CanonicalForm yToM= power (y, m);
 912    CanonicalForm F0= mod (F, yToM);
 913    CanonicalForm F1= div (F, yToM);
 914    CanonicalForm G0= mod (G, yToM);
 915    CanonicalForm G1= div (G, yToM);
 916    CanonicalForm H00= mulMod (F0, G0, MOD);
 917    CanonicalForm H11= mulMod (F1, G1, MOD);
 918    CanonicalForm H01= mulMod (F0 + F1, G0 + G1, MOD);
 919    return H11*yToM*yToM + (H01 - H11 - H00)*yToM + H00;
 920  }
 921  DEBOUTLN (cerr, "fatal end in mulMod");
 922}
 923
 924CanonicalForm prodMod (const CFList& L, const CanonicalForm& M)
 925{
 926  if (L.isEmpty())
 927    return 1;
 928  int l= L.length();
 929  if (l == 1)
 930    return mod (L.getFirst(), M);
 931  else if (l == 2) {
 932    CanonicalForm result= mulMod2 (L.getFirst(), L.getLast(), M);
 933    return result;
 934  }
 935  else
 936  {
 937    l /= 2;
 938    CFList tmp1, tmp2;
 939    CFListIterator i= L;
 940    CanonicalForm buf1, buf2;
 941    for (int j= 1; j <= l; j++, i++)
 942      tmp1.append (i.getItem());
 943    tmp2= Difference (L, tmp1);
 944    buf1= prodMod (tmp1, M);
 945    buf2= prodMod (tmp2, M);
 946    CanonicalForm result= mulMod2 (buf1, buf2, M);
 947    return result;
 948  }
 949}
 950
 951CanonicalForm prodMod (const CFList& L, const CFList& M)
 952{
 953  if (L.isEmpty())
 954    return 1;
 955  else if (L.length() == 1)
 956    return L.getFirst();
 957  else if (L.length() == 2)
 958    return mulMod (L.getFirst(), L.getLast(), M);
 959  else
 960  {
 961    int l= L.length()/2;
 962    CFListIterator i= L;
 963    CFList tmp1, tmp2;
 964    CanonicalForm buf1, buf2;
 965    for (int j= 1; j <= l; j++, i++)
 966      tmp1.append (i.getItem());
 967    tmp2= Difference (L, tmp1);
 968    buf1= prodMod (tmp1, M);
 969    buf2= prodMod (tmp2, M);
 970    return mulMod (buf1, buf2, M);
 971  }
 972}
 973
 974
 975CanonicalForm reverse (const CanonicalForm& F, int d)
 976{
 977  if (d == 0)
 978    return F;
 979  CanonicalForm A= F;
 980  Variable y= Variable (2);
 981  Variable x= Variable (1);
 982  if (degree (A, x) > 0)
 983  {
 984    A= swapvar (A, x, y);
 985    CanonicalForm result= 0;
 986    CFIterator i= A;
 987    while (d - i.exp() < 0)
 988      i++;
 989
 990    for (; i.hasTerms() && (d - i.exp() >= 0); i++)
 991      result += swapvar (i.coeff(),x,y)*power (x, d - i.exp());
 992    return result;
 993  }
 994  else
 995    return A*power (x, d);
 996}
 997
 998CanonicalForm
 999newtonInverse (const CanonicalForm& F, const int n, const CanonicalForm& M)
1000{
1001  int l= ilog2(n);
1002
1003  CanonicalForm g= mod (F, M)[0] [0];
1004
1005  ASSERT (!g.isZero(), "expected a unit");
1006
1007  Variable alpha;
1008
1009  if (!g.isOne())
1010    g = 1/g;
1011  Variable x= Variable (1);
1012  CanonicalForm result;
1013  int exp= 0;
1014  if (n & 1)
1015  {
1016    result= g;
1017    exp= 1;
1018  }
1019  CanonicalForm h;
1020
1021  for (int i= 1; i <= l; i++)
1022  {
1023    h= mulMod2 (g, mod (F, power (x, (1 << i))), M);
1024    h= mod (h, power (x, (1 << i)) - 1);
1025    h= div (h, power (x, (1 << (i - 1))));
1026    h= mod (h, M);
1027    g -= power (x, (1 << (i - 1)))*
1028         mod (mulMod2 (g, h, M), power (x, (1 << (i - 1))));
1029
1030    if (n & (1 << i))
1031    {
1032      if (exp)
1033      {
1034        h= mulMod2 (result, mod (F, power (x, exp + (1 << i))), M);
1035        h= mod (h, power (x, exp + (1 << i)) - 1);
1036        h= div (h, power (x, exp));
1037        h= mod (h, M);
1038        result -= power(x, exp)*mod (mulMod2 (g, h, M),
1039                                       power (x, (1 << i)));
1040        exp += (1 << i);
1041      }
1042      else
1043      {
1044        exp= (1 << i);
1045        result= g;
1046      }
1047    }
1048  }
1049
1050  return result;
1051}
1052
1053CanonicalForm
1054newtonDiv (const CanonicalForm& F, const CanonicalForm& G, const CanonicalForm&
1055           M)
1056{
1057  ASSERT (getCharacteristic() > 0, "positive characteristic expected");
1058  ASSERT (CFFactory::gettype() != GaloisFieldDomain, "no GF expected");
1059
1060  CanonicalForm A= mod (F, M);
1061  CanonicalForm B= mod (G, M);
1062
1063  Variable x= Variable (1);
1064  int degA= degree (A, x);
1065  int degB= degree (B, x);
1066  int m= degA - degB;
1067  if (m < 0)
1068    return 0;
1069
1070  CanonicalForm Q;
1071  if (degB <= 1 || CFFactory::gettype() == GaloisFieldDomain)
1072  {
1073    CanonicalForm R;
1074    divrem2 (A, B, Q, R, M);
1075  }
1076  else
1077  {
1078    CanonicalForm R= reverse (A, degA);
1079    CanonicalForm revB= reverse (B, degB);
1080    revB= newtonInverse (revB, m + 1, M);
1081    Q= mulMod2 (R, revB, M);
1082    Q= mod (Q, power (x, m + 1));
1083    Q= reverse (Q, m);
1084  }
1085
1086  return Q;
1087}
1088
1089void
1090newtonDivrem (const CanonicalForm& F, const CanonicalForm& G, CanonicalForm& Q,
1091              CanonicalForm& R, const CanonicalForm& M)
1092{
1093  CanonicalForm A= mod (F, M);
1094  CanonicalForm B= mod (G, M);
1095  Variable x= Variable (1);
1096  int degA= degree (A, x);
1097  int degB= degree (B, x);
1098  int m= degA - degB;
1099
1100  if (m < 0)
1101  {
1102    R= A;
1103    Q= 0;
1104    return;
1105  }
1106
1107  if (degB <= 1 || CFFactory::gettype() == GaloisFieldDomain)
1108  {
1109     divrem2 (A, B, Q, R, M);
1110  }
1111  else
1112  {
1113    R= reverse (A, degA);
1114
1115    CanonicalForm revB= reverse (B, degB);
1116    revB= newtonInverse (revB, m + 1, M);
1117    Q= mulMod2 (R, revB, M);
1118
1119    Q= mod (Q, power (x, m + 1));
1120    Q= reverse (Q, m);
1121
1122    R= A - mulMod2 (Q, B, M);
1123  }
1124}
1125
1126static inline
1127CFList split (const CanonicalForm& F, const int m, const Variable& x)
1128{
1129  CanonicalForm A= F;
1130  CanonicalForm buf= 0;
1131  bool swap= false;
1132  if (degree (A, x) <= 0)
1133    return CFList(A);
1134  else if (x.level() != A.level())
1135  {
1136    swap= true;
1137    A= swapvar (A, x, A.mvar());
1138  }
1139
1140  int j= (int) floor ((double) degree (A)/ m);
1141  CFList result;
1142  CFIterator i= A;
1143  for (; j >= 0; j--)
1144  {
1145    while (i.hasTerms() && i.exp() - j*m >= 0)
1146    {
1147      if (swap)
1148        buf += i.coeff()*power (A.mvar(), i.exp() - j*m);
1149      else
1150        buf += i.coeff()*power (x, i.exp() - j*m);
1151      i++;
1152    }
1153    if (swap)
1154      result.append (swapvar (buf, x, F.mvar()));
1155    else
1156      result.append (buf);
1157    buf= 0;
1158  }
1159  return result;
1160}
1161
1162static inline
1163void divrem32 (const CanonicalForm& F, const CanonicalForm& G, CanonicalForm& Q,
1164               CanonicalForm& R, const CFList& M);
1165
1166static inline
1167void divrem21 (const CanonicalForm& F, const CanonicalForm& G, CanonicalForm& Q,
1168               CanonicalForm& R, const CFList& M)
1169{
1170  CanonicalForm A= mod (F, M);
1171  CanonicalForm B= mod (G, M);
1172  Variable x= Variable (1);
1173  int degB= degree (B, x);
1174  int degA= degree (A, x);
1175  if (degA < degB)
1176  {
1177    Q= 0;
1178    R= A;
1179    return;
1180  }
1181  ASSERT (2*degB > degA, "expected degree (F, 1) < 2*degree (G, 1)");
1182  if (degB < 1)
1183  {
1184    divrem (A, B, Q, R);
1185    Q= mod (Q, M);
1186    R= mod (R, M);
1187    return;
1188  }
1189
1190  int m= (int) ceil ((double) (degB + 1)/2.0) + 1;
1191  CFList splitA= split (A, m, x);
1192  if (splitA.length() == 3)
1193    splitA.insert (0);
1194  if (splitA.length() == 2)
1195  {
1196    splitA.insert (0);
1197    splitA.insert (0);
1198  }
1199  if (splitA.length() == 1)
1200  {
1201    splitA.insert (0);
1202    splitA.insert (0);
1203    splitA.insert (0);
1204  }
1205
1206  CanonicalForm xToM= power (x, m);
1207
1208  CFListIterator i= splitA;
1209  CanonicalForm H= i.getItem();
1210  i++;
1211  H *= xToM;
1212  H += i.getItem();
1213  i++;
1214  H *= xToM;
1215  H += i.getItem();
1216  i++;
1217
1218  divrem32 (H, B, Q, R, M);
1219
1220  CFList splitR= split (R, m, x);
1221  if (splitR.length() == 1)
1222    splitR.insert (0);
1223
1224  H= splitR.getFirst();
1225  H *= xToM;
1226  H += splitR.getLast();
1227  H *= xToM;
1228  H += i.getItem();
1229
1230  CanonicalForm bufQ;
1231  divrem32 (H, B, bufQ, R, M);
1232
1233  Q *= xToM;
1234  Q += bufQ;
1235  return;
1236}
1237
1238static inline
1239void divrem32 (const CanonicalForm& F, const CanonicalForm& G, CanonicalForm& Q,
1240               CanonicalForm& R, const CFList& M)
1241{
1242  CanonicalForm A= mod (F, M);
1243  CanonicalForm B= mod (G, M);
1244  Variable x= Variable (1);
1245  int degB= degree (B, x);
1246  int degA= degree (A, x);
1247  if (degA < degB)
1248  {
1249    Q= 0;
1250    R= A;
1251    return;
1252  }
1253  ASSERT (3*(degB/2) > degA, "expected degree (F, 1) < 3*(degree (G, 1)/2)");
1254  if (degB < 1)
1255  {
1256    divrem (A, B, Q, R);
1257    Q= mod (Q, M);
1258    R= mod (R, M);
1259    return;
1260  }
1261  int m= (int) ceil ((double) (degB + 1)/ 2.0);
1262
1263  CFList splitA= split (A, m, x);
1264  CFList splitB= split (B, m, x);
1265
1266  if (splitA.length() == 2)
1267  {
1268    splitA.insert (0);
1269  }
1270  if (splitA.length() == 1)
1271  {
1272    splitA.insert (0);
1273    splitA.insert (0);
1274  }
1275  CanonicalForm xToM= power (x, m);
1276
1277  CanonicalForm H;
1278  CFListIterator i= splitA;
1279  i++;
1280
1281  if (degree (splitA.getFirst(), x) < degree (splitB.getFirst(), x))
1282  {
1283    H= splitA.getFirst()*xToM + i.getItem();
1284    divrem21 (H, splitB.getFirst(), Q, R, M);
1285  }
1286  else
1287  {
1288    R= splitA.getFirst()*xToM + i.getItem() + splitB.getFirst() -
1289       splitB.getFirst()*xToM;
1290    Q= xToM - 1;
1291  }
1292
1293  H= mulMod (Q, splitB.getLast(), M);
1294
1295  R= R*xToM + splitA.getLast() - H;
1296
1297  while (degree (R, x) >= degB)
1298  {
1299    xToM= power (x, degree (R, x) - degB);
1300    Q += LC (R, x)*xToM;
1301    R -= mulMod (LC (R, x), B, M)*xToM;
1302    Q= mod (Q, M);
1303    R= mod (R, M);
1304  }
1305
1306  return;
1307}
1308
1309void divrem2 (const CanonicalForm& F, const CanonicalForm& G, CanonicalForm& Q,
1310              CanonicalForm& R, const CanonicalForm& M)
1311{
1312  CanonicalForm A= mod (F, M);
1313  CanonicalForm B= mod (G, M);
1314
1315  if (B.inCoeffDomain())
1316  {
1317    divrem (A, B, Q, R);
1318    return;
1319  }
1320  if (A.inCoeffDomain() && !B.inCoeffDomain())
1321  {
1322    Q= 0;
1323    R= A;
1324    return;
1325  }
1326
1327  if (B.level() < A.level())
1328  {
1329    divrem (A, B, Q, R);
1330    return;
1331  }
1332  if (A.level() > B.level())
1333  {
1334    R= A;
1335    Q= 0;
1336    return;
1337  }
1338  if (B.level() == 1 && B.isUnivariate())
1339  {
1340    divrem (A, B, Q, R);
1341    return;
1342  }
1343  if (!(B.level() == 1 && B.isUnivariate()) && (A.level() == 1 && A.isUnivariate()))
1344  {
1345    Q= 0;
1346    R= A;
1347    return;
1348  }
1349
1350  Variable x= Variable (1);
1351  int degB= degree (B, x);
1352  if (degB > degree (A, x))
1353  {
1354    Q= 0;
1355    R= A;
1356    return;
1357  }
1358
1359  CFList splitA= split (A, degB, x);
1360
1361  CanonicalForm xToDegB= power (x, degB);
1362  CanonicalForm H, bufQ;
1363  Q= 0;
1364  CFListIterator i= splitA;
1365  H= i.getItem()*xToDegB;
1366  i++;
1367  H += i.getItem();
1368  CFList buf;
1369  while (i.hasItem())
1370  {
1371    buf= CFList (M);
1372    divrem21 (H, B, bufQ, R, buf);
1373    i++;
1374    if (i.hasItem())
1375      H= R*xToDegB + i.getItem();
1376    Q *= xToDegB;
1377    Q += bufQ;
1378  }
1379  return;
1380}
1381
1382void divrem (const CanonicalForm& F, const CanonicalForm& G, CanonicalForm& Q,
1383             CanonicalForm& R, const CFList& MOD)
1384{
1385  CanonicalForm A= mod (F, MOD);
1386  CanonicalForm B= mod (G, MOD);
1387  Variable x= Variable (1);
1388  int degB= degree (B, x);
1389  if (degB > degree (A, x))
1390  {
1391    Q= 0;
1392    R= A;
1393    return;
1394  }
1395
1396  if (degB <= 0)
1397  {
1398    divrem (A, B, Q, R);
1399    Q= mod (Q, MOD);
1400    R= mod (R, MOD);
1401    return;
1402  }
1403  CFList splitA= split (A, degB, x);
1404
1405  CanonicalForm xToDegB= power (x, degB);
1406  CanonicalForm H, bufQ;
1407  Q= 0;
1408  CFListIterator i= splitA;
1409  H= i.getItem()*xToDegB;
1410  i++;
1411  H += i.getItem();
1412  while (i.hasItem())
1413  {
1414    divrem21 (H, B, bufQ, R, MOD);
1415    i++;
1416    if (i.hasItem())
1417      H= R*xToDegB + i.getItem();
1418    Q *= xToDegB;
1419    Q += bufQ;
1420  }
1421  return;
1422}
1423
1424void sortList (CFList& list, const Variable& x)
1425{
1426  int l= 1;
1427  int k= 1;
1428  CanonicalForm buf;
1429  CFListIterator m;
1430  for (CFListIterator i= list; l <= list.length(); i++, l++)
1431  {
1432    for (CFListIterator j= list; k <= list.length() - l; k++)
1433    {
1434      m= j;
1435      m++;
1436      if (degree (j.getItem(), x) > degree (m.getItem(), x))
1437      {
1438        buf= m.getItem();
1439        m.getItem()= j.getItem();
1440        j.getItem()= buf;
1441        j++;
1442        j.getItem()= m.getItem();
1443      }
1444      else
1445        j++;
1446    }
1447    k= 1;
1448  }
1449}
1450
1451static inline
1452CFList diophantine (const CanonicalForm& F, const CFList& factors)
1453{
1454  CanonicalForm buf1, buf2, buf3, S, T;
1455  CFListIterator i= factors;
1456  CFList result;
1457  if (i.hasItem())
1458    i++;
1459  buf1= F/factors.getFirst();
1460  buf2= divNTL (F, i.getItem());
1461  buf3= extgcd (buf1, buf2, S, T);
1462  result.append (S);
1463  result.append (T);
1464  if (i.hasItem())
1465    i++;
1466  for (; i.hasItem(); i++)
1467  {
1468    buf1= divNTL (F, i.getItem());
1469    buf3= extgcd (buf3, buf1, S, T);
1470    CFListIterator k= factors;
1471    for (CFListIterator j= result; j.hasItem(); j++, k++)
1472    {
1473      j.getItem()= mulNTL (j.getItem(), S);
1474      j.getItem()= modNTL (j.getItem(), k.getItem());
1475    }
1476    result.append (T);
1477  }
1478  return result;
1479}
1480
1481void
1482henselStep12 (const CanonicalForm& F, const CFList& factors,
1483              CFArray& bufFactors, const CFList& diophant, CFMatrix& M,
1484              CFArray& Pi, int j)
1485{
1486  CanonicalForm E;
1487  CanonicalForm xToJ= power (F.mvar(), j);
1488  Variable x= F.mvar();
1489  // compute the error
1490  if (j == 1)
1491    E= F[j];
1492  else
1493  {
1494    if (degree (Pi [factors.length() - 2], x) > 0)
1495      E= F[j] - Pi [factors.length() - 2] [j];
1496    else
1497      E= F[j];
1498  }
1499
1500  CFArray buf= CFArray (diophant.length());
1501  bufFactors[0]= mod (factors.getFirst(), power (F.mvar(), j + 1));
1502  int k= 0;
1503  CanonicalForm remainder;
1504  // actual lifting
1505  for (CFListIterator i= diophant; i.hasItem(); i++, k++)
1506  {
1507    if (degree (bufFactors[k], x) > 0)
1508    {
1509      if (k > 0)
1510        remainder= modNTL (E, bufFactors[k] [0]);
1511      else
1512        remainder= E;
1513    }
1514    else
1515      remainder= modNTL (E, bufFactors[k]);
1516
1517    buf[k]= mulNTL (i.getItem(), remainder);
1518    if (degree (bufFactors[k], x) > 0)
1519      buf[k]= modNTL (buf[k], bufFactors[k] [0]);
1520    else
1521      buf[k]= modNTL (buf[k], bufFactors[k]);
1522  }
1523  for (k= 1; k < factors.length(); k++)
1524    bufFactors[k] += xToJ*buf[k];
1525
1526  // update Pi [0]
1527  int degBuf0= degree (bufFactors[0], x);
1528  int degBuf1= degree (bufFactors[1], x);
1529  if (degBuf0 > 0 && degBuf1 > 0)
1530    M (j + 1, 1)= mulNTL (bufFactors[0] [j], bufFactors[1] [j]);
1531  CanonicalForm uIZeroJ;
1532  if (j == 1)
1533  {
1534    if (degBuf0 > 0 && degBuf1 > 0)
1535      uIZeroJ= mulNTL ((bufFactors[0] [0] + bufFactors[0] [j]),
1536                  (bufFactors[1] [0] + buf[1])) - M(1, 1) - M(j + 1, 1);
1537    else if (degBuf0 > 0)
1538      uIZeroJ= mulNTL (bufFactors[0] [j], bufFactors[1]);
1539    else if (degBuf1 > 0)
1540      uIZeroJ= mulNTL (bufFactors[0], buf[1]);
1541    else
1542      uIZeroJ= 0;
1543    Pi [0] += xToJ*uIZeroJ;
1544  }
1545  else
1546  {
1547    if (degBuf0 > 0 && degBuf1 > 0)
1548      uIZeroJ= mulNTL ((bufFactors[0] [0] + bufFactors[0] [j]),
1549                  (bufFactors[1] [0] + buf[1])) - M(1, 1) - M(j + 1, 1);
1550    else if (degBuf0 > 0)
1551      uIZeroJ= mulNTL (bufFactors[0] [j], bufFactors[1]);
1552    else if (degBuf1 > 0)
1553      uIZeroJ= mulNTL (bufFactors[0], buf[1]);
1554    else
1555      uIZeroJ= 0;
1556    Pi [0] += xToJ*uIZeroJ;
1557  }
1558  CFArray tmp= CFArray (factors.length() - 1);
1559  for (k= 0; k < factors.length() - 1; k++)
1560    tmp[k]= 0;
1561  CFIterator one, two;
1562  one= bufFactors [0];
1563  two= bufFactors [1];
1564  if (degBuf0 > 0 && degBuf1 > 0)
1565  {
1566    for (k= 1; k <= (int) ceil (j/2.0); k++)
1567    {
1568      if (k != j - k + 1)
1569      {
1570        if ((one.hasTerms() && one.exp() == j - k + 1) && (two.hasTerms() && two.exp() == j - k + 1))
1571        {
1572          tmp[0] += mulNTL ((bufFactors[0] [k] + one.coeff()), (bufFactors[1] [k] +
1573                     two.coeff())) - M (k + 1, 1) - M (j - k + 2, 1);
1574          one++;
1575          two++;
1576        }
1577        else if (one.hasTerms() && one.exp() == j - k + 1)
1578        {
1579          tmp[0] += mulNTL ((bufFactors[0] [k] + one.coeff()), bufFactors[1] [k]) -
1580                     M (k + 1, 1);
1581          one++;
1582        }
1583        else if (two.hasTerms() && two.exp() == j - k + 1)
1584        {
1585          tmp[0] += mulNTL (bufFactors[0] [k], (bufFactors[1] [k] + two.coeff())) -
1586                    M (k + 1, 1);
1587          two++;
1588        }
1589      }
1590      else
1591      {
1592        tmp[0] += M (k + 1, 1);
1593      }
1594    }
1595  }
1596  Pi [0] += tmp[0]*xToJ*F.mvar();
1597
1598  // update Pi [l]
1599  int degPi, degBuf;
1600  for (int l= 1; l < factors.length() - 1; l++)
1601  {
1602    degPi= degree (Pi [l - 1], x);
1603    degBuf= degree (bufFactors[l + 1], x);
1604    if (degPi > 0 && degBuf > 0)
1605      M (j + 1, l + 1)= mulNTL (Pi [l - 1] [j], bufFactors[l + 1] [j]);
1606    if (j == 1)
1607    {
1608      if (degPi > 0 && degBuf > 0)
1609        Pi [l] += xToJ*(mulNTL (Pi [l - 1] [0] + Pi [l - 1] [j],
1610                  bufFactors[l + 1] [0] + buf[l + 1]) - M (j + 1, l +1) -
1611                  M (1, l + 1));
1612      else if (degPi > 0)
1613        Pi [l] += xToJ*(mulNTL (Pi [l - 1] [j], bufFactors[l + 1]));
1614      else if (degBuf > 0)
1615        Pi [l] += xToJ*(mulNTL (Pi [l - 1], buf[l + 1]));
1616    }
1617    else
1618    {
1619      if (degPi > 0 && degBuf > 0)
1620      {
1621        uIZeroJ= mulNTL (uIZeroJ, bufFactors [l + 1] [0]);
1622        uIZeroJ += mulNTL (Pi [l - 1] [0], buf [l + 1]);
1623      }
1624      else if (degPi > 0)
1625        uIZeroJ= mulNTL (uIZeroJ, bufFactors [l + 1]);
1626      else if (degBuf > 0)
1627      {
1628        uIZeroJ= mulNTL (uIZeroJ, bufFactors [l + 1] [0]);
1629        uIZeroJ += mulNTL (Pi [l - 1], buf[l + 1]);
1630      }
1631      Pi[l] += xToJ*uIZeroJ;
1632    }
1633    one= bufFactors [l + 1];
1634    two= Pi [l - 1];
1635    if (two.hasTerms() && two.exp() == j + 1)
1636    {
1637      if (degBuf > 0 && degPi > 0)
1638      {
1639          tmp[l] += mulNTL (two.coeff(), bufFactors[l + 1][0]);
1640          two++;
1641      }
1642      else if (degPi > 0)
1643      {
1644          tmp[l] += mulNTL (two.coeff(), bufFactors[l + 1]);
1645          two++;
1646      }
1647    }
1648    if (degBuf > 0 && degPi > 0)
1649    {
1650      for (k= 1; k <= (int) ceil (j/2.0); k++)
1651      {
1652        if (k != j - k + 1)
1653        {
1654          if ((one.hasTerms() && one.exp() == j - k + 1) &&
1655              (two.hasTerms() && two.exp() == j - k + 1))
1656          {
1657            tmp[l] += mulNTL ((bufFactors[l + 1] [k] + one.coeff()), (Pi[l - 1] [k] +
1658                       two.coeff())) - M (k + 1, l + 1) - M (j - k + 2, l + 1);
1659            one++;
1660            two++;
1661          }
1662          else if (one.hasTerms() && one.exp() == j - k + 1)
1663          {
1664            tmp[l] += mulNTL ((bufFactors[l + 1] [k] + one.coeff()), Pi[l - 1] [k]) -
1665                       M (k + 1, l + 1);
1666            one++;
1667          }
1668          else if (two.hasTerms() && two.exp() == j - k + 1)
1669          {
1670            tmp[l] += mulNTL (bufFactors[l + 1] [k], (Pi[l - 1] [k] + two.coeff())) -
1671                      M (k + 1, l + 1);
1672            two++;
1673          }
1674        }
1675        else
1676          tmp[l] += M (k + 1, l + 1);
1677      }
1678    }
1679    Pi[l] += tmp[l]*xToJ*F.mvar();
1680  }
1681  return;
1682}
1683
1684void
1685henselLift12 (const CanonicalForm& F, CFList& factors, int l, CFArray& Pi,
1686              CFList& diophant, CFMatrix& M, bool sort)
1687{
1688  if (sort)
1689    sortList (factors, Variable (1));
1690  Pi= CFArray (factors.length() - 1);
1691  CFListIterator j= factors;
1692  diophant= diophantine (F[0], factors);
1693  DEBOUTLN (cerr, "diophant= " << diophant);
1694  j++;
1695  Pi [0]= mulNTL (j.getItem(), mod (factors.getFirst(), F.mvar()));
1696  M (1, 1)= Pi [0];
1697  int i= 1;
1698  if (j.hasItem())
1699    j++;
1700  for (; j.hasItem(); j++, i++)
1701  {
1702    Pi [i]= mulNTL (Pi [i - 1], j.getItem());
1703    M (1, i + 1)= Pi [i];
1704  }
1705  CFArray bufFactors= CFArray (factors.length());
1706  i= 0;
1707  for (CFListIterator k= factors; k.hasItem(); i++, k++)
1708  {
1709    if (i == 0)
1710      bufFactors[i]= mod (k.getItem(), F.mvar());
1711    else
1712      bufFactors[i]= k.getItem();
1713  }
1714  for (i= 1; i < l; i++)
1715    henselStep12 (F, factors, bufFactors, diophant, M, Pi, i);
1716
1717  CFListIterator k= factors;
1718  for (i= 0; i < factors.length (); i++, k++)
1719    k.getItem()= bufFactors[i];
1720  factors.removeFirst();
1721  return;
1722}
1723
1724void
1725henselLiftResume12 (const CanonicalForm& F, CFList& factors, int start, int
1726                    end, CFArray& Pi, const CFList& diophant, CFMatrix& M)
1727{
1728  CFArray bufFactors= CFArray (factors.length());
1729  int i= 0;
1730  CanonicalForm xToStart= power (F.mvar(), start);
1731  for (CFListIterator k= factors; k.hasItem(); k++, i++)
1732  {
1733    if (i == 0)
1734      bufFactors[i]= mod (k.getItem(), xToStart);
1735    else
1736      bufFactors[i]= k.getItem();
1737  }
1738  for (i= start; i < end; i++)
1739    henselStep12 (F, factors, bufFactors, diophant, M, Pi, i);
1740
1741  CFListIterator k= factors;
1742  for (i= 0; i < factors.length(); k++, i++)
1743    k.getItem()= bufFactors [i];
1744  factors.removeFirst();
1745  return;
1746}
1747
1748static inline
1749CFList
1750biDiophantine (const CanonicalForm& F, const CFList& factors, const int d)
1751{
1752  Variable y= F.mvar();
1753  CFList result;
1754  if (y.level() == 1)
1755  {
1756    result= diophantine (F, factors);
1757    return result;
1758  }
1759  else
1760  {
1761    CFList buf= factors;
1762    for (CFListIterator i= buf; i.hasItem(); i++)
1763      i.getItem()= mod (i.getItem(), y);
1764    CanonicalForm A= mod (F, y);
1765    int bufD= 1;
1766    CFList recResult= biDiophantine (A, buf, bufD);
1767    CanonicalForm e= 1;
1768    CFList p;
1769    CFArray bufFactors= CFArray (factors.length());
1770    CanonicalForm yToD= power (y, d);
1771    int k= 0;
1772    for (CFListIterator i= factors; i.hasItem(); i++, k++)
1773    {
1774      bufFactors [k]= i.getItem();
1775    }
1776    CanonicalForm b;
1777    for (k= 0; k < factors.length(); k++) //TODO compute b's faster
1778    {
1779      b= 1;
1780      if (fdivides (bufFactors[k], F))
1781        b= F/bufFactors[k];
1782      else
1783      {
1784        for (int l= 0; l < factors.length(); l++)
1785        {
1786          if (l == k)
1787            continue;
1788          else
1789          {
1790            b= mulMod2 (b, bufFactors[l], yToD);
1791          }
1792        }
1793      }
1794      p.append (b);
1795    }
1796
1797    CFListIterator j= p;
1798    for (CFListIterator i= recResult; i.hasItem(); i++, j++)
1799      e -= i.getItem()*j.getItem();
1800
1801    if (e.isZero())
1802      return recResult;
1803    CanonicalForm coeffE;
1804    CFList s;
1805    result= recResult;
1806    CanonicalForm g;
1807    for (int i= 1; i < d; i++)
1808    {
1809      if (degree (e, y) > 0)
1810        coeffE= e[i];
1811      else
1812        coeffE= 0;
1813      if (!coeffE.isZero())
1814      {
1815        CFListIterator k= result;
1816        CFListIterator l= p;
1817        int ii= 0;
1818        j= recResult;
1819        for (; j.hasItem(); j++, k++, l++, ii++)
1820        {
1821          g= coeffE*j.getItem();
1822          if (degree (bufFactors[ii], y) <= 0)
1823            g= mod (g, bufFactors[ii]);
1824          else
1825            g= mod (g, bufFactors[ii][0]);
1826          k.getItem() += g*power (y, i);
1827          e -= mulMod2 (g*power(y, i), l.getItem(), yToD);
1828          DEBOUTLN (cerr, "mod (e, power (y, i + 1))= " <<
1829                    mod (e, power (y, i + 1)));
1830        }
1831      }
1832      if (e.isZero())
1833        break;
1834    }
1835
1836    DEBOUTLN (cerr, "mod (e, y)= " << mod (e, y));
1837
1838#ifdef DEBUGOUTPUT
1839    CanonicalForm test= 0;
1840    j= p;
1841    for (CFListIterator i= result; i.hasItem(); i++, j++)
1842      test += mod (i.getItem()*j.getItem(), power (y, d));
1843    DEBOUTLN (cerr, "test= " << test);
1844#endif
1845    return result;
1846  }
1847}
1848
1849static inline
1850CFList
1851multiRecDiophantine (const CanonicalForm& F, const CFList& factors,
1852                     const CFList& recResult, const CFList& M, const int d)
1853{
1854  Variable y= F.mvar();
1855  CFList result;
1856  CFListIterator i;
1857  CanonicalForm e= 1;
1858  CFListIterator j= factors;
1859  CFList p;
1860  CFArray bufFactors= CFArray (factors.length());
1861  CanonicalForm yToD= power (y, d);
1862  int k= 0;
1863  for (CFListIterator i= factors; i.hasItem(); i++, k++)
1864    bufFactors [k]= i.getItem();
1865  CanonicalForm b;
1866  CFList buf= M;
1867  buf.removeLast();
1868  buf.append (yToD);
1869  for (k= 0; k < factors.length(); k++) //TODO compute b's faster
1870  {
1871    b= 1;
1872    if (fdivides (bufFactors[k], F))
1873      b= F/bufFactors[k];
1874    else
1875    {
1876      for (int l= 0; l < factors.length(); l++)
1877      {
1878        if (l == k)
1879          continue;
1880        else
1881        {
1882          b= mulMod (b, bufFactors[l], buf);
1883        }
1884      }
1885    }
1886    p.append (b);
1887  }
1888  j= p;
1889  for (CFListIterator i= recResult; i.hasItem(); i++, j++)
1890    e -= mulMod (i.getItem(), j.getItem(), M);
1891
1892  if (e.isZero())
1893    return recResult;
1894  CanonicalForm coeffE;
1895  CFList s;
1896  result= recResult;
1897  CanonicalForm g;
1898  for (int i= 1; i < d; i++)
1899  {
1900    if (degree (e, y) > 0)
1901      coeffE= e[i];
1902    else
1903      coeffE= 0;
1904    if (!coeffE.isZero())
1905    {
1906      CFListIterator k= result;
1907      CFListIterator l= p;
1908      j= recResult;
1909      int ii= 0;
1910      CanonicalForm dummy;
1911      for (; j.hasItem(); j++, k++, l++, ii++)
1912      {
1913        g= mulMod (coeffE, j.getItem(), M);
1914        if (degree (bufFactors[ii], y) <= 0)
1915          divrem (g, mod (bufFactors[ii], Variable (y.level() - 1)), dummy,
1916                  g, M);
1917        else
1918          divrem (g, bufFactors[ii][0], dummy, g, M);
1919        k.getItem() += g*power (y, i);
1920        e -= mulMod (g*power (y, i), l.getItem(), M);
1921      }
1922    }
1923
1924    if (e.isZero())
1925      break;
1926  }
1927
1928#ifdef DEBUGOUTPUT
1929    CanonicalForm test= 0;
1930    j= p;
1931    for (CFListIterator i= result; i.hasItem(); i++, j++)
1932      test += mod (i.getItem()*j.getItem(), power (y, d));
1933    DEBOUTLN (cerr, "test= " << test);
1934#endif
1935  return result;
1936}
1937
1938static inline
1939void
1940henselStep (const CanonicalForm& F, const CFList& factors, CFArray& bufFactors,
1941            const CFList& diophant, CFMatrix& M, CFArray& Pi, int j,
1942            const CFList& MOD)
1943{
1944  CanonicalForm E;
1945  CanonicalForm xToJ= power (F.mvar(), j);
1946  Variable x= F.mvar();
1947  // compute the error
1948  if (j == 1)
1949  {
1950    E= F[j];
1951#ifdef DEBUGOUTPUT
1952    CanonicalForm test= 1;
1953    for (int i= 0; i < factors.length(); i++)
1954    {
1955      if (i == 0)
1956        test= mulMod (test, mod (bufFactors [i], xToJ), MOD);
1957      else
1958        test= mulMod (test, bufFactors[i], MOD);
1959    }
1960    CanonicalForm test2= mod (F-test, xToJ);
1961
1962    test2= mod (test2, MOD);
1963    DEBOUTLN (cerr, "test= " << test2);
1964#endif
1965  }
1966  else
1967  {
1968#ifdef DEBUGOUTPUT
1969    CanonicalForm test= 1;
1970    for (int i= 0; i < factors.length(); i++)
1971    {
1972      if (i == 0)
1973        test *= mod (bufFactors [i], power (x, j));
1974      else
1975        test *= bufFactors[i];
1976    }
1977    test= mod (test, power (x, j));
1978    test= mod (test, MOD);
1979    CanonicalForm test2= mod (F, power (x, j - 1)) - mod (test, power (x, j-1));
1980    DEBOUTLN (cerr, "test= " << test2);
1981#endif
1982
1983    if (degree (Pi [factors.length() - 2], x) > 0)
1984      E= F[j] - Pi [factors.length() - 2] [j];
1985    else
1986      E= F[j];
1987  }
1988
1989  CFArray buf= CFArray (diophant.length());
1990  bufFactors[0]= mod (factors.getFirst(), power (F.mvar(), j + 1));
1991  int k= 0;
1992  // actual lifting
1993  CanonicalForm dummy, rest1;
1994  for (CFListIterator i= diophant; i.hasItem(); i++, k++)
1995  {
1996    if (degree (bufFactors[k], x) > 0)
1997    {
1998      if (k > 0)
1999        divrem (E, bufFactors[k] [0], dummy, rest1, MOD);
2000      else
2001        rest1= E;
2002    }
2003    else
2004      divrem (E, bufFactors[k], dummy, rest1, MOD);
2005
2006    buf[k]= mulMod (i.getItem(), rest1, MOD);
2007
2008    if (degree (bufFactors[k], x) > 0)
2009      divrem (buf[k], bufFactors[k] [0], dummy, buf[k], MOD);
2010    else
2011      divrem (buf[k], bufFactors[k], dummy, buf[k], MOD);
2012  }
2013  for (k= 1; k < factors.length(); k++)
2014    bufFactors[k] += xToJ*buf[k];
2015
2016  // update Pi [0]
2017  int degBuf0= degree (bufFactors[0], x);
2018  int degBuf1= degree (bufFactors[1], x);
2019  if (degBuf0 > 0 && degBuf1 > 0)
2020    M (j + 1, 1)= mulMod (bufFactors[0] [j], bufFactors[1] [j], MOD);
2021  CanonicalForm uIZeroJ;
2022  if (j == 1)
2023  {
2024    if (degBuf0 > 0 && degBuf1 > 0)
2025      uIZeroJ= mulMod ((bufFactors[0] [0] + bufFactors[0] [j]),
2026                  (bufFactors[1] [0] + buf[1]), MOD) - M(1, 1) - M(j + 1, 1);
2027    else if (degBuf0 > 0)
2028      uIZeroJ= mulMod (bufFactors[0] [j], bufFactors[1], MOD);
2029    else if (degBuf1 > 0)
2030      uIZeroJ= mulMod (bufFactors[0], buf[1], MOD);
2031    else
2032      uIZeroJ= 0;
2033    Pi [0] += xToJ*uIZeroJ;
2034  }
2035  else
2036  {
2037    if (degBuf0 > 0 && degBuf1 > 0)
2038      uIZeroJ= mulMod ((bufFactors[0] [0] + bufFactors[0] [j]),
2039                  (bufFactors[1] [0] + buf[1]), MOD) - M(1, 1) - M(j + 1, 1);
2040    else if (degBuf0 > 0)
2041      uIZeroJ= mulMod (bufFactors[0] [j], bufFactors[1], MOD);
2042    else if (degBuf1 > 0)
2043      uIZeroJ= mulMod (bufFactors[0], buf[1], MOD);
2044    else
2045      uIZeroJ= 0;
2046    Pi [0] += xToJ*uIZeroJ;
2047  }
2048
2049  CFArray tmp= CFArray (factors.length() - 1);
2050  for (k= 0; k < factors.length() - 1; k++)
2051    tmp[k]= 0;
2052  CFIterator one, two;
2053  one= bufFactors [0];
2054  two= bufFactors [1];
2055  if (degBuf0 > 0 && degBuf1 > 0)
2056  {
2057    for (k= 1; k <= (int) ceil (j/2.0); k++)
2058    {
2059      if (k != j - k + 1)
2060      {
2061        if ((one.hasTerms() && one.exp() == j - k + 1) &&
2062            (two.hasTerms() && two.exp() == j - k + 1))
2063        {
2064          tmp[0] += mulMod 

Large files files are truncated, but you can click here to view the full file