Changeset eb3abbb in git


Ignore:
Timestamp:
May 23, 2013, 7:20:14 PM (11 years ago)
Author:
Oleksandr Motsak <motsak@…>
Branches:
(u'spielwiese', 'fe61d9c35bf7c61f2b6cbf1b56e25e2f08d536cc')
Children:
02d3d641de579bc515b51a2769fe96f9149cd06755150e5fbf40e9e0fa5d2c2fe5efd309f4ab2e12
Parents:
e2c9b2fd1bc99d892a2a72d1c9e5046a794ac330aba90e858d998d23e9f343db0210fff2a167df2f
Message:
Merge pull request #329 from mmklee/factory_newtonpolygon_sw

Factory newtonpolygon sw
Location:
factory
Files:
5 edited

Legend:

Unmodified
Added
Removed

  • factory/cfNewtonPolygon.cc

    re2c9b2 reb3abbb  
    499499}
    500500
    501 #ifdef HAVE_NTL
    502 void convexDense(int** points, int sizePoints, mat_ZZ& M, vec_ZZ& A)
     501void mpz_mat_mul (const mpz_t* N, mpz_t*& M)
     502{
     503  mpz_t * tmp= new mpz_t[4];
     504
     505  mpz_init_set (tmp[0], N[0]);
     506  mpz_mul (tmp[0], tmp[0], M[0]);
     507  mpz_addmul (tmp[0], N[1], M[2]);
     508
     509  mpz_init_set (tmp[1], N[0]);
     510  mpz_mul (tmp[1], tmp[1], M[1]);
     511  mpz_addmul (tmp[1], N[1], M[3]);
     512
     513  mpz_init_set (tmp[2], N[2]);
     514  mpz_mul (tmp[2], tmp[2], M[0]);
     515  mpz_addmul (tmp[2], N[3], M[2]);
     516
     517  mpz_init_set (tmp[3], N[2]);
     518  mpz_mul (tmp[3], tmp[3], M[1]);
     519  mpz_addmul (tmp[3], N[3], M[3]);
     520
     521  mpz_set (M[0], tmp[0]);
     522  mpz_set (M[1], tmp[1]);
     523  mpz_set (M[2], tmp[2]);
     524  mpz_set (M[3], tmp[3]);
     525
     526  mpz_clear (tmp[0]);
     527  mpz_clear (tmp[1]);
     528  mpz_clear (tmp[2]);
     529  mpz_clear (tmp[3]);
     530
     531  delete [] tmp;
     532}
     533
     534void mpz_mat_inv (mpz_t*& M)
     535{
     536  mpz_t det;
     537  mpz_init_set (det, M[0]);
     538  mpz_mul (det, det, M[3]);
     539  mpz_submul (det, M[1], M[2]);
     540
     541  mpz_t tmp;
     542  mpz_init_set (tmp, M[0]);
     543  mpz_divexact (tmp, tmp, det);
     544  mpz_set (M[0], M[3]);
     545  mpz_divexact (M[0], M[0], det);
     546  mpz_set (M[3], tmp);
     547
     548  mpz_neg (M[1], M[1]);
     549  mpz_divexact (M[1], M[1], det);
     550  mpz_neg (M[2], M[2]);
     551  mpz_divexact (M[2], M[2], det);
     552
     553  mpz_clear (det);
     554  mpz_clear (tmp);
     555}
     556
     557void convexDense(int** points, int sizePoints, mpz_t*& M, mpz_t*& A)
    503558{
    504559  if (sizePoints < 3)
    … …  
    506561    if (sizePoints == 2)
    507562    {
    508       int maxX= (points [1] [1] < points [0] [1])?points [0] [1]:points [1] [1];
    509       int maxY= (points [1] [0] < points [0] [0])?points [0] [0]:points [1] [0];
    510       long u,v,g;
    511       XGCD (g, u, v, maxX, maxY);
    512       M.SetDims (2,2);
    513       A.SetLength (2);
     563      mpz_t u,v,g,maxX,maxY;
     564      mpz_init (u);
     565      mpz_init (v);
     566      mpz_init (g);
     567      mpz_init_set_si (maxX,
     568                       (points[1][1] < points[0][1])?points[0][1]:points[1][1]);
     569      mpz_init_set_si (maxY,
     570                       (points[1][0] < points[0][0])?points[0][0]:points[1][0]);
     571      mpz_gcdext (g, u, v, maxX, maxY);
    514572      if (points [0] [1] != points [0] [0] && points [1] [0] != points [1] [1])
    515573      {
    516         M (1,1)= -u;
    517         M (1,2)= v;
    518         M (2,1)= maxY/g;
    519         M (2,2)= maxX/g;
    520         A (1)= u*maxX;
    521         A (2)= -(maxY/g)*maxX;
     574        mpz_set (A[0], u);
     575        mpz_mul (A[0], A[0], maxX);
     576        mpz_set (M[2], maxY);
     577        mpz_divexact (M[2], M[2], g);
     578        mpz_set (A[1], M[2]);
     579        mpz_neg (A[1], A[1]);
     580        mpz_mul (A[1], A[1], maxX);
     581        mpz_neg (u, u);
     582        mpz_set (M[0], u);
     583        mpz_set (M[1], v);
     584        mpz_set (M[3], maxX);
     585        mpz_divexact (M[3], M[3], g);
    522586      }
    523587      else
    524588      {
    525         M (1,1)= u;
    526         M (1,2)= v;
    527         M (2,1)= -maxY/g;
    528         M (2,2)= maxX/g;
    529         A (1)= to_ZZ (0);
    530         A (2)= to_ZZ (0);
     589        mpz_set (M[0], u);
     590        mpz_set (M[1], v);
     591        mpz_set (M[2], maxY);
     592        mpz_divexact (M[2], M[2], g);
     593        mpz_neg (M[2], M[2]);
     594        mpz_set (M[3], maxX);
     595        mpz_divexact (M[3], M[3], g);
    531596      }
     597      mpz_clear (u);
     598      mpz_clear (v);
     599      mpz_clear (g);
     600      mpz_clear (maxX);
     601      mpz_clear (maxY);
    532602    }
    533603    else if (sizePoints == 1)
    534604    {
    535       ident (M, 2);
    536       A.SetLength (2);
    537       A (1)= to_ZZ (0);
    538       A (2)= to_ZZ (0);
     605      mpz_set_si (M[0], 1);
     606      mpz_set_si (M[3], 1);
    539607    }
    540608    return;
    541609  }
    542   A.SetLength (2);
    543   A (1)= to_ZZ (0);
    544   A (2)= to_ZZ (0);
    545   ident (M, 2);
    546   mat_ZZ Mu;
    547   Mu.SetDims (2, 2);
    548   Mu (2,1)= to_ZZ (1);
    549   Mu (1,2)= to_ZZ (1);
    550   Mu (1,1)= to_ZZ (0);
    551   Mu (2,2)= to_ZZ (0);
    552   mat_ZZ Lambda;
    553   Lambda.SetDims (2, 2);
    554   Lambda (1,1)= to_ZZ (1);
    555   Lambda (1,2)= to_ZZ (-1);
    556   Lambda (2,2)= to_ZZ (1);
    557   Lambda (2,1)= to_ZZ (0);
    558   mat_ZZ InverseLambda;
    559   InverseLambda.SetDims (2,2);
    560   InverseLambda (1,1)= to_ZZ (1);
    561   InverseLambda (1,2)= to_ZZ (1);
    562   InverseLambda (2,2)= to_ZZ (1);
    563   InverseLambda (2,1)= to_ZZ (0);
    564   ZZ tmp;
     610  mpz_set_si (M[0], 1);
     611  mpz_set_si (M[3], 1);
     612
     613  mpz_t * Mu= new mpz_t[4];
     614  mpz_init_set_si (Mu[1], 1);
     615  mpz_init_set_si (Mu[2], 1);
     616  mpz_init (Mu[0]);
     617  mpz_init (Mu[3]);
     618
     619  mpz_t * Lambda= new mpz_t[4];
     620  mpz_init_set_si (Lambda[0], 1);
     621  mpz_init_set_si (Lambda[1], -1);
     622  mpz_init_set_si (Lambda[3], 1);
     623  mpz_init (Lambda[2]);
     624
     625  mpz_t * InverseLambda= new mpz_t[4];
     626  mpz_init_set_si (InverseLambda[0], 1);
     627  mpz_init_set_si (InverseLambda[1], 1);
     628  mpz_init_set_si (InverseLambda[3], 1);
     629  mpz_init (InverseLambda[2]);
     630
     631  mpz_t tmp;
     632  mpz_init (tmp);
    565633  int minDiff, minSum, maxDiff, maxSum, maxX, maxY, b, d, f, h;
    566   getMaxMin (points, sizePoints, minDiff, minSum, maxDiff, maxSum, maxX, maxY);
     634  getMaxMin(points, sizePoints, minDiff, minSum, maxDiff, maxSum, maxX, maxY);
    567635  do
    568636  {
    … …  
    570638    {
    571639      mu (points, sizePoints);
    572       M= Mu*M;
    573       tmp= A (1);
    574       A (1)= A (2);
    575       A (2)= tmp;
    576     }
    577     getMaxMin (points, sizePoints, minDiff, minSum, maxDiff, maxSum, maxX, maxY);
     640
     641      mpz_mat_mul (Mu, M);
     642
     643      mpz_set (tmp, A[0]);
     644      mpz_set (A[0], A[1]);
     645      mpz_set (A[1], tmp);
     646    }
     647    getMaxMin(points, sizePoints, minDiff, minSum, maxDiff, maxSum, maxX, maxY);
    578648    b= maxX - maxDiff;
    579649    d= maxX + maxY - maxSum;
    … …  
    584654      lambda (points, sizePoints);
    585655      tau (points, sizePoints, maxY - f);
    586       M= Lambda*M;
    587       A [0] += (maxY-f);
     656
     657      mpz_mat_mul (Lambda, M);
     658
     659      if (maxY-f > 0)
     660        mpz_add_ui (A[0], A[0], maxY-f);
     661      else
     662        mpz_add_ui (A[0], A[0], f-maxY);
    588663      maxX= maxX + maxY - b - f;
    589664    }
    … …  
    592667      lambdaInverse (points, sizePoints);
    593668      tau (points, sizePoints, -h);
    594       M= InverseLambda*M;
    595       A [0] += (-h);
     669
     670      mpz_mat_mul (InverseLambda, M);
     671
     672      if (h < 0)
     673        mpz_add_ui (A[0], A[0], -h);
     674      else
     675        mpz_sub_ui (A[0], A[0], h);
    596676      maxX= maxX + maxY - d - h;
    597677    }
    598678    else
     679    {
     680      mpz_clear (tmp);
     681      mpz_clear (Mu[0]);
     682      mpz_clear (Mu[1]);
     683      mpz_clear (Mu[2]);
     684      mpz_clear (Mu[3]);
     685      delete [] Mu;
     686
     687      mpz_clear (Lambda[0]);
     688      mpz_clear (Lambda[1]);
     689      mpz_clear (Lambda[2]);
     690      mpz_clear (Lambda[3]);
     691      delete [] Lambda;
     692
     693      mpz_clear (InverseLambda[0]);
     694      mpz_clear (InverseLambda[1]);
     695      mpz_clear (InverseLambda[2]);
     696      mpz_clear (InverseLambda[3]);
     697      delete [] InverseLambda;
     698
    599699      return;
     700    }
    600701  } while (1);
    601702}
    602703
    603704CanonicalForm
    604 compress (const CanonicalForm& F, mat_ZZ& M, vec_ZZ& A, bool computeMA)
     705compress (const CanonicalForm& F, mpz_t*& M, mpz_t*& A, bool computeMA)
    605706{
    606707  int n;
    … …  
    611712    convexDense (newtonPolyg, n, M, A);
    612713  }
     714
    613715  CanonicalForm result= 0;
    614   ZZ expX, expY;
    615716  Variable x= Variable (1);
    616717  Variable y= Variable (2);
    617718
    618   ZZ minExpX, minExpY;
     719  mpz_t expX, expY, minExpX, minExpY;
     720  mpz_init (expX);
     721  mpz_init (expY);
     722  mpz_init (minExpX);
     723  mpz_init (minExpY);
    619724
    620725  int k= 0;
    … …  
    624729    if (i.coeff().inCoeffDomain() && hasFirstAlgVar (i.coeff(), alpha))
    625730    {
    626       expX= i.exp()*M (1,2) + A (1);
    627       expY= i.exp()*M (2,2) + A (2);
     731      mpz_set (expX, A[0]);
     732      mpz_set (expY, A[1]);
     733      mpz_addmul_ui (expX, M[1], i.exp());
     734      mpz_addmul_ui (expY, M[3], i.exp());
     735
    628736      if (k == 0)
    629737      {
    630         minExpY= expY;
    631         minExpX= expX;
     738        mpz_set (minExpX, expX);
     739        mpz_set (minExpY, expY);
    632740        k= 1;
    633741      }
    634742      else
    635743      {
    636         minExpY= (minExpY > expY) ? expY : minExpY;
    637         minExpX= (minExpX > expX) ? expX : minExpX;
     744        if (mpz_cmp (minExpY, expY) > 0)
     745          mpz_set (minExpY, expY);
     746        if (mpz_cmp (minExpX, expX) > 0)
     747          mpz_set (minExpX, expX);
    638748      }
    639       result += i.coeff()*power (x, to_long (expX))*power (y, to_long (expY));
     749      result += i.coeff()*power (x, mpz_get_si (expX))*
     750                power (y, mpz_get_si (expY));
    640751      continue;
    641752    }
    … …  
    643754    if (k == 0)
    644755    {
    645       expX= j.exp()*M (1,1) + i.exp()*M (1,2) + A (1);
    646       expY= j.exp()*M (2,1) + i.exp()*M (2,2) + A (2);
    647       minExpX= expX;
    648       minExpY= expY;
    649       result += j.coeff()*power (x, to_long (expX))*power (y, to_long (expY));
     756      mpz_set (expX, A[0]);
     757      mpz_addmul_ui (expX, M[1], i.exp());
     758      mpz_addmul_ui (expX, M[0], j.exp());
     759
     760      mpz_set (expY, A[1]);
     761      mpz_addmul_ui (expY, M[3], i.exp());
     762      mpz_addmul_ui (expY, M[2], j.exp());
     763
     764      mpz_set (minExpX, expX);
     765      mpz_set (minExpY, expY);
     766      result += j.coeff()*power (x, mpz_get_si (expX))*
     767                power (y, mpz_get_si (expY));
    650768      j++;
    651769      k= 1;
    … …  
    654772    for (; j.hasTerms(); j++)
    655773    {
    656       expX= j.exp()*M (1,1) + i.exp()*M (1,2) + A (1);
    657       expY= j.exp()*M (2,1) + i.exp()*M (2,2) + A (2);
    658       result += j.coeff()*power (x, to_long (expX))*power (y, to_long (expY));
    659       minExpY= (minExpY > expY) ? expY : minExpY;
    660       minExpX= (minExpX > expX) ? expX : minExpX;
    661     }
    662   }
    663 
    664   if (to_long (minExpX) < 0)
    665   {
    666     result *= power (x,-to_long(minExpX));
     774      mpz_set (expX, A[0]);
     775      mpz_addmul_ui (expX, M[1], i.exp());
     776      mpz_addmul_ui (expX, M[0], j.exp());
     777
     778      mpz_set (expY, A[1]);
     779      mpz_addmul_ui (expY, M[3], i.exp());
     780      mpz_addmul_ui (expY, M[2], j.exp());
     781
     782      result += j.coeff()*power (x, mpz_get_si (expX))*
     783                power (y, mpz_get_si (expY));
     784      if (mpz_cmp (minExpY, expY) > 0)
     785        mpz_set (minExpY, expY);
     786      if (mpz_cmp (minExpX, expX) > 0)
     787        mpz_set (minExpX, expX);
     788    }
     789  }
     790
     791  if (mpz_sgn (minExpX) < 0)
     792  {
     793    result *= power (x,-mpz_get_si(minExpX));
    667794    result /= CanonicalForm (x, 0);
    668795  }
    669796  else
    670     result /= power (x,to_long(minExpX));
    671 
    672   if (to_long (minExpY) < 0)
    673   {
    674     result *= power (y,-to_long(minExpY));
     797    result /= power (x,mpz_get_si(minExpX));
     798
     799  if (mpz_sgn (minExpY) < 0)
     800  {
     801    result *= power (y,-mpz_get_si(minExpY));
    675802    result /= CanonicalForm (y, 0);
    676803  }
    677804  else
    678     result /= power (y,to_long(minExpY));
     805    result /= power (y,mpz_get_si(minExpY));
    679806
    680807  CanonicalForm tmp= LC (result);
    … …  
    692819      delete [] newtonPolyg [i];
    693820    delete [] newtonPolyg;
    694     M= inv (M);
    695   }
     821    mpz_mat_inv (M);
     822  }
     823
     824  mpz_clear (expX);
     825  mpz_clear (expY);
     826  mpz_clear (minExpX);
     827  mpz_clear (minExpY);
    696828
    697829  return result;
    … …  
    699831
    700832CanonicalForm
    701 decompress (const CanonicalForm& F, const mat_ZZ& inverseM, const vec_ZZ& A)
     833decompress (const CanonicalForm& F, const mpz_t* inverseM, const mpz_t * A)
    702834{
    703835  CanonicalForm result= 0;
    704   ZZ expX, expY;
    705836  Variable x= Variable (1);
    706837  Variable y= Variable (2);
    707   ZZ minExpX, minExpY;
     838
     839  mpz_t expX, expY, minExpX, minExpY;
     840  mpz_init (expX);
     841  mpz_init (expY);
     842  mpz_init (minExpX);
     843  mpz_init (minExpY);
     844
    708845  if (F.isUnivariate() && F.level() == 1)
    709846  {
    710847    CFIterator i= F;
    711     expX= (i.exp() - A (1))*inverseM (1,1) + (-A (2))*inverseM (1,2);
    712     expY= (i.exp() - A (1))*inverseM (2,1) + (-A (2))*inverseM (2,2);
    713     minExpX= expX;
    714     minExpY= expY;
    715     result += i.coeff()*power (x, to_long (expX))*power (y, to_long (expY));
     848
     849    mpz_set_si (expX, i.exp());
     850    mpz_sub (expX, expX, A[0]);
     851    mpz_mul (expX, expX, inverseM[0]);
     852    mpz_submul (expX, inverseM[1], A[1]);
     853
     854    mpz_set_si (expY, i.exp());
     855    mpz_sub (expY, expY, A[0]);
     856    mpz_mul (expY, expY, inverseM[2]);
     857    mpz_submul (expY, inverseM[3], A[1]);
     858
     859    mpz_set (minExpX, expX);
     860    mpz_set (minExpY, expY);
     861    result += i.coeff()*power (x, mpz_get_si (expX))*
     862              power (y, mpz_get_si (expY));
    716863    i++;
    717864    for (; i.hasTerms(); i++)
    718865    {
    719       expX= (i.exp() - A (1))*inverseM (1,1) + (-A (2))*inverseM (1,2);
    720       expY= (i.exp() - A (1))*inverseM (2,1) + (-A (2))*inverseM (2,2);
    721       result += i.coeff()*power (x, to_long (expX))*power (y, to_long (expY));
    722       minExpY= (minExpY > expY) ? expY : minExpY;
    723       minExpX= (minExpX > expX) ? expX : minExpX;
    724     }
    725 
    726     if (to_long (minExpX) < 0)
    727     {
    728       result *= power (x,-to_long(minExpX));
     866      mpz_set_si (expX, i.exp());
     867      mpz_sub (expX, expX, A[0]);
     868      mpz_mul (expX, expX, inverseM[0]);
     869      mpz_submul (expX, inverseM[1], A[1]);
     870
     871      mpz_set_si (expY, i.exp());
     872      mpz_sub (expY, expY, A[0]);
     873      mpz_mul (expY, expY, inverseM[2]);
     874      mpz_submul (expY, inverseM[3], A[1]);
     875
     876      result += i.coeff()*power (x, mpz_get_si (expX))*
     877                power (y, mpz_get_si (expY));
     878      if (mpz_cmp (minExpY, expY) > 0)
     879        mpz_set (minExpY, expY);
     880      if (mpz_cmp (minExpX, expX) > 0)
     881        mpz_set (minExpX, expX);
     882    }
     883
     884    if (mpz_sgn (minExpX) < 0)
     885    {
     886      result *= power (x,-mpz_get_si(minExpX));
    729887      result /= CanonicalForm (x, 0);
    730888    }
    731889    else
    732       result /= power (x,to_long(minExpX));
    733 
    734     if (to_long (minExpY) < 0)
    735     {
    736       result *= power (y,-to_long(minExpY));
     890      result /= power (x,mpz_get_si(minExpX));
     891
     892    if (mpz_sgn (minExpY) < 0)
     893    {
     894      result *= power (y,-mpz_get_si(minExpY));
    737895      result /= CanonicalForm (y, 0);
    738896    }
    739897    else
    740       result /= power (y,to_long(minExpY));
     898      result /= power (y,mpz_get_si(minExpY));
     899
     900    mpz_clear (expX);
     901    mpz_clear (expY);
     902    mpz_clear (minExpX);
     903    mpz_clear (minExpY);
    741904
    742905    return result/ Lc (result); //normalize
    743906  }
    744907
     908  mpz_t tmp;
     909  mpz_init (tmp);
    745910  int k= 0;
    746911  Variable alpha;
    … …  
    749914    if (i.coeff().inCoeffDomain() && hasFirstAlgVar (i.coeff(), alpha))
    750915    {
    751       expX= -A(1)*inverseM (1,1) + (i.exp() - A (2))*inverseM (1,2);
    752       expY= -A(1)*inverseM (2,1) + (i.exp() - A (2))*inverseM (2,2);
     916      mpz_set_si (expX, i.exp());
     917      mpz_sub (expX, expX, A[1]);
     918      mpz_mul (expX, expX, inverseM[1]);
     919      mpz_submul (expX, A[0], inverseM[0]);
     920
     921      mpz_set_si (expY, i.exp());
     922      mpz_sub (expY, expY, A[1]);
     923      mpz_mul (expY, expY, inverseM[3]);
     924      mpz_submul (expY, A[0], inverseM[2]);
     925
    753926      if (k == 0)
    754927      {
    755         minExpY= expY;
    756         minExpX= expX;
     928        mpz_set (minExpX, expX);
     929        mpz_set (minExpY, expY);
    757930        k= 1;
    758931      }
    759932      else
    760933      {
    761         minExpY= (minExpY > expY) ? expY : minExpY;
    762         minExpX= (minExpX > expX) ? expX : minExpX;
     934        if (mpz_cmp (minExpY, expY) > 0)
     935          mpz_set (minExpY, expY);
     936        if (mpz_cmp (minExpX, expX) > 0)
     937          mpz_set (minExpX, expX);
    763938      }
    764       result += i.coeff()*power (x, to_long (expX))*power (y, to_long (expY));
     939      result += i.coeff()*power (x, mpz_get_si (expX))*
     940                power (y, mpz_get_si (expY));
    765941      continue;
    766942    }
    … …  
    768944    if (k == 0)
    769945    {
    770       expX= (j.exp() - A (1))*inverseM (1,1) + (i.exp() - A (2))*inverseM (1,2);
    771       expY= (j.exp() - A (1))*inverseM (2,1) + (i.exp() - A (2))*inverseM (2,2);
    772       minExpX= expX;
    773       minExpY= expY;
    774       result += j.coeff()*power (x, to_long (expX))*power (y, to_long (expY));
     946      mpz_set_si (expX, j.exp());
     947      mpz_sub (expX, expX, A[0]);
     948      mpz_mul (expX, expX, inverseM[0]);
     949      mpz_set_si (tmp, i.exp());
     950      mpz_sub (tmp, tmp, A[1]);
     951      mpz_addmul (expX, tmp, inverseM[1]);
     952
     953      mpz_set_si (expY, j.exp());
     954      mpz_sub (expY, expY, A[0]);
     955      mpz_mul (expY, expY, inverseM[2]);
     956      mpz_set_si (tmp, i.exp());
     957      mpz_sub (tmp, tmp, A[1]);
     958      mpz_addmul (expY, tmp, inverseM[3]);
     959
     960      mpz_set (minExpX, expX);
     961      mpz_set (minExpY, expY);
     962      result += j.coeff()*power (x, mpz_get_si (expX))*
     963                power (y, mpz_get_si (expY));
    775964      j++;
    776965      k= 1;
    … …  
    779968    for (; j.hasTerms(); j++)
    780969    {
    781       expX= (j.exp() - A (1))*inverseM (1,1) + (i.exp() - A (2))*inverseM (1,2);
    782       expY= (j.exp() - A (1))*inverseM (2,1) + (i.exp() - A (2))*inverseM (2,2);
    783       result += j.coeff()*power (x, to_long (expX))*power (y, to_long (expY));
    784       minExpY= (minExpY > expY) ? expY : minExpY;
    785       minExpX= (minExpX > expX) ? expX : minExpX;
    786     }
    787   }
    788 
    789   if (to_long (minExpX) < 0)
    790   {
    791     result *= power (x,-to_long(minExpX));
     970      mpz_set_si (expX, j.exp());
     971      mpz_sub (expX, expX, A[0]);
     972      mpz_mul (expX, expX, inverseM[0]);
     973      mpz_set_si (tmp, i.exp());
     974      mpz_sub (tmp, tmp, A[1]);
     975      mpz_addmul (expX, tmp, inverseM[1]);
     976
     977      mpz_set_si (expY, j.exp());
     978      mpz_sub (expY, expY, A[0]);
     979      mpz_mul (expY, expY, inverseM[2]);
     980      mpz_set_si (tmp, i.exp());
     981      mpz_sub (tmp, tmp, A[1]);
     982      mpz_addmul (expY, tmp, inverseM[3]);
     983
     984      result += j.coeff()*power (x, mpz_get_si (expX))*
     985                power (y, mpz_get_si (expY));
     986      if (mpz_cmp (minExpY, expY) > 0)
     987        mpz_set (minExpY, expY);
     988      if (mpz_cmp (minExpX, expX) > 0)
     989        mpz_set (minExpX, expX);
     990    }
     991  }
     992
     993  if (mpz_sgn (minExpX) < 0)
     994  {
     995    result *= power (x,-mpz_get_si(minExpX));
    792996    result /= CanonicalForm (x, 0);
    793997  }
    794998  else
    795     result /= power (x,to_long(minExpX));
    796 
    797   if (to_long (minExpY) < 0)
    798   {
    799     result *= power (y,-to_long(minExpY));
     999    result /= power (x,mpz_get_si(minExpX));
     1000
     1001  if (mpz_sgn (minExpY) < 0)
     1002  {
     1003    result *= power (y,-mpz_get_si(minExpY));
    8001004    result /= CanonicalForm (y, 0);
    8011005  }
    8021006  else
    803     result /= power (y,to_long(minExpY));
     1007    result /= power (y,mpz_get_si(minExpY));
     1008
     1009  mpz_clear (expX);
     1010  mpz_clear (expY);
     1011  mpz_clear (minExpX);
     1012  mpz_clear (minExpY);
     1013  mpz_clear (tmp);
    8041014
    8051015  return result/Lc (result); //normalize
    8061016}
    807 #endif
    8081017
    8091018//assumes the input is a Newton polygon of a bivariate polynomial which is

  • factory/cfNewtonPolygon.h

    re2c9b2 reb3abbb  
    1616
    1717// #include "config.h"
    18 
    19 #ifdef HAVE_NTL
    20 #include "NTLconvert.h"
    21 #endif
    2218
    2319/// compute a polygon
    … …  
    104100                          );
    105101
    106 
    107 #ifdef HAVE_NTL
    108102/// Algorithm 5 as described in Convex-Dense Bivariate Polynomial Factorization
    109103/// by Berthomieu, Lecerf
    … …  
    111105                                 ///< M (points)+A
    112106                  int sizePoints,///< [in] size of points
    113                   mat_ZZ& M,     ///< [in,out] returns an invertible matrix
    114                   vec_ZZ& A      ///< [in,out] returns translation
     107                  mpz_t*& M,     ///< [in,out] returns an invertible 2x2 matrix
     108                  mpz_t*& A      ///< [in,out] returns translation
    115109                 );
    116110
    … …  
    122116compress (const CanonicalForm& F, ///< [in] compressed, i.e. F.level()==2,
    123117                                  ///< bivariate poly
    124           mat_ZZ& inverseM,       ///< [in,out] returns the inverse of M,
     118          mpz_t*& inverseM,       ///< [in,out] returns the inverse of M,
    125119                                  ///< if computeMA==true, M otherwise
    126           vec_ZZ& A,              ///< [in,out] returns translation
     120          mpz_t*& A,              ///< [in,out] returns translation
    127121          bool computeMA= true    ///< [in] whether to compute M and A
    128122         );
    … …  
    135129decompress (const CanonicalForm& F,///< [in] compressed, i.e. F.level()<= 2,
    136130                                   ///< uni- or bivariate poly
    137             const mat_ZZ& M,       ///< [in] matrix M obtained from compress
    138             const vec_ZZ& A        ///< [in] vector A obtained from compress
     131            const mpz_t* M,       ///< [in] matrix M obtained from compress
     132            const mpz_t* A        ///< [in] vector A obtained from compress
    139133           );
    140 #endif
    141134
    142135#endif

  • factory/cf_gcd_smallp.cc

    re2c9b2 reb3abbb  
    537537  ppB= B/cB;
    538538
    539   int sizeNewtonPolyg;
    540   int ** newtonPolyg= NULL;
    541   mat_ZZ MM;
    542   vec_ZZ V;
    543   bool compressConvexDense= false; //(ppA.level() == 2 && ppB.level() == 2);
    544   if (compressConvexDense)
    545   {
    546     CanonicalForm bufcA= cA;
    547     CanonicalForm bufcB= cB;
    548     cA= content (ppA, 1);
    549     cB= content (ppB, 1);
    550     ppA /= cA;
    551     ppB /= cB;
    552     gcdcAcB *= gcd (cA, cB);
    553     cA *= bufcA;
    554     cB *= bufcB;
    555     if (ppA.isUnivariate() || ppB.isUnivariate())
    556     {
    557       if (ppA.level() == ppB.level())
    558       {
    559         CanonicalForm result= gcd (ppA, ppB);
    560         coF= N ((ppA/result)*(cA/gcdcAcB));
    561         coG= N ((ppB/result)*(cB/gcdcAcB));
    562         return N (result*gcdcAcB);
    563       }
    564       else
    565       {
    566         coF= N (ppA*(cA/gcdcAcB));
    567         coG= N (ppB*(cB/gcdcAcB));
    568         return N (gcdcAcB);
    569       }
    570     }
    571 
    572     newtonPolyg= newtonPolygon (ppA,ppB, sizeNewtonPolyg);
    573     convexDense (newtonPolyg, sizeNewtonPolyg, MM, V);
    574 
    575     for (int i= 0; i < sizeNewtonPolyg; i++)
    576       delete [] newtonPolyg[i];
    577     delete [] newtonPolyg;
    578 
    579     ppA= compress (ppA, MM, V, false);
    580     ppB= compress (ppB, MM, V, false);
    581     MM= inv (MM);
    582 
    583     if (ppA.isUnivariate() && ppB.isUnivariate())
    584     {
    585       if (ppA.level() == ppB.level())
    586       {
    587         CanonicalForm result= gcd (ppA, ppB);
    588         coF= N (decompress ((ppA/result), MM, V)*(cA/gcdcAcB));
    589         coG= N (decompress ((ppB/result), MM, V)*(cB/gcdcAcB));
    590         return N (decompress (result, MM, V)*gcdcAcB);
    591       }
    592       else
    593       {
    594         coF= N (decompress (ppA, MM, V));
    595         coG= N (decompress (ppB, MM, V));
    596         return N (gcdcAcB);
    597       }
    598     }
    599   }
    600 
    601539  CanonicalForm lcA, lcB;  // leading coefficients of A and B
    602540  CanonicalForm gcdlcAlcB;
    … …  
    604542  lcA= uni_lcoeff (ppA);
    605543  lcB= uni_lcoeff (ppB);
    606 
    607   /*if (fdivides (lcA, lcB))
    608   {
    609     if (fdivides (A, B))
    610       return F/Lc(F);
    611   }
    612   if (fdivides (lcB, lcA))
    613   {
    614     if (fdivides (B, A))
    615       return G/Lc(G);
    616   }*/
    617544
    618545  gcdlcAlcB= gcd (lcA, lcB);
    … …  
    817744          ppCoG= mapDown (ppCoG, prim_elem, im_prim_elem, alpha, u, v);
    818745          DEBOUTLN (cerr, "ppH after mapDown= " << ppH);
    819           if (compressConvexDense)
    820           {
    821             ppH= decompress (ppH, MM, V);
    822             ppCoF= decompress (ppCoF, MM, V);
    823             ppCoG= decompress (ppCoG, MM, V);
    824           }
    825746          coF= N ((cA/gcdcAcB)*ppCoF);
    826747          coG= N ((cB/gcdcAcB)*ppCoG);
    … …  
    835756                (fdivides (ppH, ppA, ppCoF) && fdivides (ppH, ppB, ppCoG)))
    836757      {
    837         if (compressConvexDense)
    838         {
    839           ppH= decompress (ppH, MM, V);
    840           ppCoF= decompress (ppCoF, MM, V);
    841           ppCoG= decompress (ppCoG, MM, V);
    842         }
    843758        coF= N ((cA/gcdcAcB)*ppCoF);
    844759        coG= N ((cB/gcdcAcB)*ppCoG);
    … …  
    998913  ppB= B/cB;
    999914
    1000   int sizeNewtonPolyg;
    1001   int ** newtonPolyg= NULL;
    1002   mat_ZZ MM;
    1003   vec_ZZ V;
    1004   bool compressConvexDense= false; //(ppA.level() == 2 && ppB.level() == 2);
    1005   if (compressConvexDense)
    1006   {
    1007     CanonicalForm bufcA= cA;
    1008     CanonicalForm bufcB= cB;
    1009     cA= content (ppA, 1);
    1010     cB= content (ppB, 1);
    1011     ppA /= cA;
    1012     ppB /= cB;
    1013     gcdcAcB *= gcd (cA, cB);
    1014     cA *= bufcA;
    1015     cB *= bufcB;
    1016     if (ppA.isUnivariate() || ppB.isUnivariate())
    1017     {
    1018       if (ppA.level() == ppB.level())
    1019       {
    1020         CanonicalForm result= gcd (ppA, ppB);
    1021         coF= N ((ppA/result)*(cA/gcdcAcB));
    1022         coG= N ((ppB/result)*(cB/gcdcAcB));
    1023         return N (result*gcdcAcB);
    1024       }
    1025       else
    1026       {
    1027         coF= N (ppA*(cA/gcdcAcB));
    1028         coG= N (ppB*(cB/gcdcAcB));
    1029         return N (gcdcAcB);
    1030       }
    1031     }
    1032 
    1033     newtonPolyg= newtonPolygon (ppA,ppB, sizeNewtonPolyg);
    1034     convexDense (newtonPolyg, sizeNewtonPolyg, MM, V);
    1035 
    1036     for (int i= 0; i < sizeNewtonPolyg; i++)
    1037       delete [] newtonPolyg[i];
    1038     delete [] newtonPolyg;
    1039 
    1040     ppA= compress (ppA, MM, V, false);
    1041     ppB= compress (ppB, MM, V, false);
    1042     MM= inv (MM);
    1043 
    1044     if (ppA.isUnivariate() && ppB.isUnivariate())
    1045     {
    1046       if (ppA.level() == ppB.level())
    1047       {
    1048         CanonicalForm result= gcd (ppA, ppB);
    1049         coF= N (decompress ((ppA/result), MM, V)*(cA/gcdcAcB));
    1050         coG= N (decompress ((ppB/result), MM, V)*(cB/gcdcAcB));
    1051         return N (decompress (result, MM, V)*gcdcAcB);
    1052       }
    1053       else
    1054       {
    1055         coF= N (decompress (ppA, MM, V));
    1056         coG= N (decompress (ppB, MM, V));
    1057         return N (gcdcAcB);
    1058       }
    1059     }
    1060   }
    1061 
    1062915  CanonicalForm lcA, lcB;  // leading coefficients of A and B
    1063916  CanonicalForm gcdlcAlcB;
    … …  
    1065918  lcA= uni_lcoeff (ppA);
    1066919  lcB= uni_lcoeff (ppB);
    1067 
    1068   /*if (fdivides (lcA, lcB))
    1069   {
    1070     if (fdivides (ppA, ppB, coG))
    1071     {
    1072       coF= 1;
    1073       if (compressConvexDense)
    1074         coG= decompress (coG, MM, V);
    1075       coG= N (coG*(cB/gcdcAcB));
    1076       return F;
    1077     }
    1078   }
    1079   if (fdivides (lcB, lcA))
    1080   {
    1081     if (fdivides (ppB, ppA, coF))
    1082     {
    1083       coG= 1;
    1084       if (compressConvexDense)
    1085         coF= decompress (coF, MM, V);
    1086       coF= N (coF*(cA/gcdcAcB));
    1087       return G;
    1088     }
    1089   }*/
    1090920
    1091921  gcdlcAlcB= gcd (lcA, lcB);
    … …  
    1119949  H= 0;
    1120950  bool fail= false;
    1121   //topLevel= false;
     951  topLevel= false;
    1122952  bool inextension= false;
    1123953  int p=-1;
    … …  
    12651095          ppCoG= GFMapDown (ppCoG, k);
    12661096          DEBOUTLN (cerr, "ppH after mapDown= " << ppH);
    1267           if (compressConvexDense)
    1268           {
    1269             ppH= decompress (ppH, MM, V);
    1270             ppCoF= decompress (ppCoF, MM, V);
    1271             ppCoG= decompress (ppCoG, MM, V);
    1272           }
    12731097          coF= N ((cA/gcdcAcB)*ppCoF);
    12741098          coG= N ((cB/gcdcAcB)*ppCoG);
    … …  
    12861110           (fdivides (ppH, ppA, ppCoF) && fdivides (ppH, ppB, ppCoG)))
    12871111        {
    1288           if (compressConvexDense)
    1289           {
    1290             ppH= decompress (ppH, MM, V);
    1291             ppCoF= decompress (ppCoF, MM, V);
    1292             ppCoG= decompress (ppCoG, MM, V);
    1293           }
    12941112          coF= N ((cA/gcdcAcB)*ppCoF);
    12951113          coG= N ((cB/gcdcAcB)*ppCoG);
    … …  
    14591277  ppB= B/cB;
    14601278
    1461   int sizeNewtonPolyg;
    1462   int ** newtonPolyg= NULL;
    1463   mat_ZZ MM;
    1464   vec_ZZ V;
    1465   bool compressConvexDense= false; //(ppA.level() == 2 && ppB.level() == 2);
    1466   if (compressConvexDense)
    1467   {
    1468     CanonicalForm bufcA= cA;
    1469     CanonicalForm bufcB= cB;
    1470     cA= content (ppA, 1);
    1471     cB= content (ppB, 1);
    1472     ppA /= cA;
    1473     ppB /= cB;
    1474     gcdcAcB *= gcd (cA, cB);
    1475     cA *= bufcA;
    1476     cB *= bufcB;
    1477     if (ppA.isUnivariate() || ppB.isUnivariate())
    1478     {
    1479       if (ppA.level() == ppB.level())
    1480       {
    1481         CanonicalForm result= gcd (ppA, ppB);
    1482         coF= N ((ppA/result)*(cA/gcdcAcB));
    1483         coG= N ((ppB/result)*(cB/gcdcAcB));
    1484         return N (result*gcdcAcB);
    1485       }
    1486       else
    1487       {
    1488         coF= N (ppA*(cA/gcdcAcB));
    1489         coG= N (ppB*(cB/gcdcAcB));
    1490         return N (gcdcAcB);
    1491       }
    1492     }
    1493 
    1494     newtonPolyg= newtonPolygon (ppA,ppB, sizeNewtonPolyg);
    1495     convexDense (newtonPolyg, sizeNewtonPolyg, MM, V);
    1496 
    1497     for (int i= 0; i < sizeNewtonPolyg; i++)
    1498       delete [] newtonPolyg[i];
    1499     delete [] newtonPolyg;
    1500 
    1501     ppA= compress (ppA, MM, V, false);
    1502     ppB= compress (ppB, MM, V, false);
    1503     MM= inv (MM);
    1504 
    1505     if (ppA.isUnivariate() && ppB.isUnivariate())
    1506     {
    1507       if (ppA.level() == ppB.level())
    1508       {
    1509         CanonicalForm result= gcd (ppA, ppB);
    1510         coF= N (decompress ((ppA/result), MM, V)*(cA/gcdcAcB));
    1511         coG= N (decompress ((ppB/result), MM, V)*(cB/gcdcAcB));
    1512         return N (decompress (result, MM, V)*gcdcAcB);
    1513       }
    1514       else
    1515       {
    1516         coF= N (decompress (ppA, MM, V));
    1517         coG= N (decompress (ppB, MM, V));
    1518         return N (gcdcAcB);
    1519       }
    1520     }
    1521   }
    1522 
    15231279  CanonicalForm lcA, lcB;  // leading coefficients of A and B
    15241280  CanonicalForm gcdlcAlcB;
    15251281  lcA= uni_lcoeff (ppA);
    15261282  lcB= uni_lcoeff (ppB);
    1527 
    1528   /*if (fdivides (lcA, lcB))
    1529   {
    1530     if (fdivides (A, B))
    1531       return F/Lc(F);
    1532   }
    1533   if (fdivides (lcB, lcA))
    1534   {
    1535     if (fdivides (B, A))
    1536       return G/Lc(G);
    1537   }*/
    15381283
    15391284  gcdlcAlcB= gcd (lcA, lcB);
    … …  
    17761521           (fdivides (ppH, ppA, ppCoF) && fdivides (ppH, ppB, ppCoG)))
    17771522      {
    1778         if (compressConvexDense)
    1779         {
    1780           ppH= decompress (ppH, MM, V);
    1781           ppCoF= decompress (ppCoF, MM, V);
    1782           ppCoG= decompress (ppCoG, MM, V);
    1783         }
    17841523        coF= N ((cA/gcdcAcB)*ppCoF);
    17851524        coG= N ((cB/gcdcAcB)*ppCoG);

  • factory/facBivar.h

    re2c9b2 reb3abbb  
    8383    return result;
    8484  }
    85   mat_ZZ M;
    86   vec_ZZ S;
     85
     86  mpz_t * M=new mpz_t [4];
     87  mpz_init (M[0]);
     88  mpz_init (M[1]);
     89  mpz_init (M[2]);
     90  mpz_init (M[3]);
     91
     92  mpz_t * S=new mpz_t [2];
     93  mpz_init (S[0]);
     94  mpz_init (S[1]);
     95
    8796  F= compress (F, M, S);
    8897  CFList result= biFactorize (F, v);
    … …  
    98107    result.insert (Lc (G));
    99108  }
     109
     110  mpz_clear (M[0]);
     111  mpz_clear (M[1]);
     112  mpz_clear (M[2]);
     113  mpz_clear (M[3]);
     114  delete [] M;
     115
     116  mpz_clear (S[0]);
     117  mpz_clear (S[1]);
     118  delete [] S;
     119
    100120  return result;
    101121}
    … …  
    197217    return result;
    198218  }
    199   mat_ZZ M;
    200   vec_ZZ S;
     219
     220  mpz_t * M=new mpz_t [4];
     221  mpz_init (M[0]);
     222  mpz_init (M[1]);
     223  mpz_init (M[2]);
     224  mpz_init (M[3]);
     225
     226  mpz_t * S=new mpz_t [2];
     227  mpz_init (S[0]);
     228  mpz_init (S[1]);
     229
    201230  F= compress (F, M, S);
    202231  TIMING_START (fac_bi_sqrf);
    … …  
    233262    result.insert (CFFactor (LcF, 1));
    234263  }
     264
     265  mpz_clear (M[0]);
     266  mpz_clear (M[1]);
     267  mpz_clear (M[2]);
     268  mpz_clear (M[3]);
     269  delete [] M;
     270
     271  mpz_clear (S[0]);
     272  mpz_clear (S[1]);
     273  delete [] S;
     274
    235275  return result;
    236276}

  • factory/facFqBivar.h

    re2c9b2 reb3abbb  
    9494    return result;
    9595  }
    96   mat_ZZ M;
    97   vec_ZZ S;
     96  mpz_t * M=new mpz_t [4];
     97  mpz_init (M[0]);
     98  mpz_init (M[1]);
     99  mpz_init (M[2]);
     100  mpz_init (M[3]);
     101
     102  mpz_t * S=new mpz_t [2];
     103  mpz_init (S[0]);
     104  mpz_init (S[1]);
     105
    98106  F= compress (F, M, S);
    99107  CFList result= biFactorize (F, info);
    … …  
    106114  normalize (result);
    107115  result.insert (Lc(G));
     116
     117  mpz_clear (M[0]);
     118  mpz_clear (M[1]);
     119  mpz_clear (M[2]);
     120  mpz_clear (M[3]);
     121  delete [] M;
     122
     123  mpz_clear (S[0]);
     124  mpz_clear (S[1]);
     125  delete [] S;
     126
    108127  return result;
    109128}
    … …  
    228247    return result;
    229248  }
    230   mat_ZZ M;
    231   vec_ZZ S;
     249  mpz_t * M=new mpz_t [4];
     250  mpz_init (M[0]);
     251  mpz_init (M[1]);
     252  mpz_init (M[2]);
     253  mpz_init (M[3]);
     254
     255  mpz_t * S=new mpz_t [2];
     256  mpz_init (S[0]);
     257  mpz_init (S[1]);
     258
    232259  F= compress (F, M, S);
    233260
    … …  
    254281  normalize (result);
    255282  result.insert (CFFactor (LcF, 1));
     283
     284  mpz_clear (M[0]);
     285  mpz_clear (M[1]);
     286  mpz_clear (M[2]);
     287  mpz_clear (M[3]);
     288  delete [] M;
     289
     290  mpz_clear (S[0]);
     291  mpz_clear (S[1]);
     292  delete [] S;
     293
    256294  return result;
    257295}
    … …  
    335373    return result;
    336374  }
    337   mat_ZZ M;
    338   vec_ZZ S;
     375
     376  mpz_t * M=new mpz_t [4];
     377  mpz_init (M[0]);
     378  mpz_init (M[1]);
     379  mpz_init (M[2]);
     380  mpz_init (M[3]);
     381
     382  mpz_t * S=new mpz_t [2];
     383  mpz_init (S[0]);
     384  mpz_init (S[1]);
     385
    339386  F= compress (F, M, S);
    340387
    … …  
    361408  normalize (result);
    362409  result.insert (CFFactor (LcF, 1));
     410
     411  mpz_clear (M[0]);
     412  mpz_clear (M[1]);
     413  mpz_clear (M[2]);
     414  mpz_clear (M[3]);
     415  delete [] M;
     416
     417  mpz_clear (S[0]);
     418  mpz_clear (S[1]);
     419  delete [] S;
     420
    363421  return result;
    364422}
    … …  
    443501    return result;
    444502  }
    445   mat_ZZ M;
    446   vec_ZZ S;
     503
     504  mpz_t * M=new mpz_t [4];
     505  mpz_init (M[0]);
     506  mpz_init (M[1]);
     507  mpz_init (M[2]);
     508  mpz_init (M[3]);
     509
     510  mpz_t * S=new mpz_t [2];
     511  mpz_init (S[0]);
     512  mpz_init (S[1]);
     513
    447514  F= compress (F, M, S);
    448515
    … …  
    469536  normalize (result);
    470537  result.insert (CFFactor (LcF, 1));
     538
     539  mpz_clear (M[0]);
     540  mpz_clear (M[1]);
     541  mpz_clear (M[2]);
     542  mpz_clear (M[3]);
     543  delete [] M;
     544
     545  mpz_clear (S[0]);
     546  mpz_clear (S[1]);
     547  delete [] S;
     548
    471549  return result;
    472550}
Note: See TracChangeset for help on using the changeset viewer.