Changeset 34636b in git


Ignore:
Timestamp:
Apr 29, 2013, 12:35:53 PM (10 years ago)
Author:
Martin Lee <martinlee84@…>
Branches:
(u'jengelh-datetime', 'ceac47cbc86fe4a15902392bdbb9bd2ae0ea02c6')(u'spielwiese', 'a800fe4b3e9d37a38c5a10cc0ae9dfa0c15a4ee6')
Children:
dd5ab893f9075772f5a92bae7a22a2e986ce4f10
Parents:
e2c9b2fd1bc99d892a2a72d1c9e5046a794ac330
git-author:
Martin Lee <martinlee84@web.de>2013-04-29 12:35:53+02:00
git-committer:
Martin Lee <martinlee84@web.de>2013-05-23 18:18:13+02:00
Message:
chg: rewrote convex dense compression with mpz_t
Location:
factory
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • factory/cfNewtonPolygon.cc

    re2c9b2 r34636b  
    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 r34636b  
    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
Note: See TracChangeset for help on using the changeset viewer.