source: git/factory/facFqBivarUtil.cc @ 0b447e

jengelh-datetimespielwiese
Last change on this file since 0b447e was 0b447e, checked in by Martin Lee <martinlee84@…>, 10 years ago
fix: bug in log derivative computation over field extensions
  • Property mode set to 100644
File size: 24.3 KB
RevLine 
[7bf145]1/*****************************************************************************\
[806c18]2 * Computer Algebra System SINGULAR
[7bf145]3\*****************************************************************************/
4/** @file facFqBivarUtil.cc
[806c18]5 *
6 * This file provides utility functions for bivariate factorization
[7bf145]7 *
8 * @author Martin Lee
9 *
10 **/
11/*****************************************************************************/
12
[e4fe2b]13#include "config.h"
[7bf145]14
[5e4636]15#include "timing.h"
16
[7bf145]17#include "cf_map.h"
[81d96c]18#include "algext.h"
[7bf145]19#include "cf_map_ext.h"
[6db552]20#include "templates/ftmpl_functions.h"
[7bf145]21#include "ExtensionInfo.h"
[fecc08]22#include "cf_algorithm.h"
23#include "cf_factory.h"
[07718c3]24#include "cf_util.h"
[fecc08]25#include "imm.h"
26#include "cf_iter.h"
[7bf145]27#include "facFqBivarUtil.h"
[f876a66]28#include "cfNewtonPolygon.h"
29#include "facHensel.h"
[0e2e23]30#include "facMul.h"
[7bf145]31
[5e4636]32TIMING_DEFINE_PRINT(fac_log_deriv_div)
33TIMING_DEFINE_PRINT(fac_log_deriv_mul)
34TIMING_DEFINE_PRINT(fac_log_deriv_pre)
[7bf145]35
[806c18]36void append (CFList& factors1, const CFList& factors2)
[7bf145]37{
[806c18]38  for (CFListIterator i= factors2; i.hasItem(); i++)
[7bf145]39  {
40    if (!i.getItem().inCoeffDomain())
41      factors1.append (i.getItem());
42  }
43  return;
44}
45
[806c18]46void decompress (CFList& factors, const CFMap& N)
[7bf145]47{
[806c18]48  for (CFListIterator i= factors; i.hasItem(); i++)
[7bf145]49    i.getItem()= N (i.getItem());
[f876a66]50}
51
52void decompress (CFFList& factors, const CFMap& N)
53{
54  for (CFFListIterator i= factors; i.hasItem(); i++)
55    i.getItem()= CFFactor (N (i.getItem().factor()), i.getItem().exp());
[7bf145]56}
57
[806c18]58void appendSwapDecompress (CFList& factors1, const CFList& factors2,
59                           const CFList& factors3, const bool swap1,
60                           const bool swap2, const CFMap& N)
[7bf145]61{
62  Variable x= Variable (1);
63  Variable y= Variable (2);
[806c18]64  for (CFListIterator i= factors1; i.hasItem(); i++)
[7bf145]65  {
[806c18]66    if (swap1)
[7bf145]67    {
[806c18]68      if (!swap2)
[7bf145]69        i.getItem()= swapvar (i.getItem(), x, y);
70    }
[806c18]71    else
[7bf145]72    {
[806c18]73      if (swap2)
[7bf145]74        i.getItem()= swapvar (i.getItem(), y, x);
75    }
76    i.getItem()= N (i.getItem());
77  }
[806c18]78  for (CFListIterator i= factors2; i.hasItem(); i++)
[7bf145]79    factors1.append (N (i.getItem()));
[806c18]80  for (CFListIterator i= factors3; i.hasItem(); i++)
[7bf145]81    factors1.append (N (i.getItem()));
82  return;
83}
84
[806c18]85void swapDecompress (CFList& factors, const bool swap, const CFMap& N)
[7bf145]86{
87  Variable x= Variable (1);
88  Variable y= Variable (2);
[806c18]89  for (CFListIterator i= factors; i.hasItem(); i++)
[7bf145]90  {
[806c18]91    if (swap)
[7bf145]92      i.getItem()= swapvar (i.getItem(), x, y);
93    i.getItem()= N (i.getItem());
94  }
95  return;
96}
97
[806c18]98static inline
99bool GFInExtensionHelper (const CanonicalForm& F, const int number)
[7bf145]100{
101  if (F.isOne()) return false;
102  InternalCF* buf;
103  int exp;
104  bool result= false;
[806c18]105  if (F.inBaseDomain())
[7bf145]106  {
107    buf= F.getval();
108    exp= imm2int (buf);
[806c18]109    if (exp%number != 0)
[7bf145]110      return true;
[806c18]111    else
[7bf145]112      return result;
113  }
[806c18]114  else
[7bf145]115  {
[806c18]116    for (CFIterator i= F; i.hasTerms(); i++)
[7bf145]117    {
118      result= GFInExtensionHelper (i.coeff(), number);
119      if (result == true)
120        return result;
121    }
122  }
123  return result;
124}
125
[806c18]126static inline
[f876a66]127bool FqInExtensionHelper (const CanonicalForm& F, const CanonicalForm& gamma,
128                          const CanonicalForm& delta, CFList& source,
129                          CFList& dest)
[7bf145]130{
131  bool result= false;
[806c18]132  if (F.inBaseDomain())
[7bf145]133    return result;
[806c18]134  else if (F.inCoeffDomain())
[7bf145]135  {
[806c18]136    if (!fdivides (gamma, F))
[7bf145]137      return true;
[806c18]138    else
[f876a66]139    {
140      int pos= findItem (source, F);
141      if (pos > 0)
142        return false;
143      Variable a;
144      hasFirstAlgVar (F, a);
145      int bound= ipower (getCharacteristic(), degree (getMipo (a)));
146      CanonicalForm buf= 1;
147      for (int i= 1; i < bound; i++)
148      {
149        buf *= gamma;
150        if (buf == F)
151        {
152          source.append (buf);
153          dest.append (power (delta, i));
154          return false;
155        }
156      }
157      return true;
158    }
[7bf145]159  }
[806c18]160  else
[7bf145]161  {
[806c18]162    for (CFIterator i= F; i.hasTerms(); i++)
[7bf145]163    {
[f876a66]164      result= FqInExtensionHelper (i.coeff(), gamma, delta, source, dest);
[7bf145]165      if (result == true)
[806c18]166        return result;
[7bf145]167    }
168  }
169  return result;
170}
171
172bool isInExtension (const CanonicalForm& F, const CanonicalForm& gamma,
[f876a66]173                    const int k, const CanonicalForm& delta,
174                    CFList& source, CFList& dest)
[7bf145]175{
176  bool result;
[806c18]177  if (CFFactory::gettype() == GaloisFieldDomain)
[7bf145]178  {
179    int p= getCharacteristic();
180    int orderFieldExtension= ipower (p, getGFDegree()) - 1;
181    int order= ipower (p, k) - 1;
182    int number= orderFieldExtension/order;
183    result= GFInExtensionHelper (F, number);
184    return result;
185  }
[806c18]186  else
[7bf145]187  {
[f876a66]188    result= FqInExtensionHelper (F, gamma, delta, source, dest);
[7bf145]189    return result;
190  }
191}
192
[806c18]193CanonicalForm
194mapDown (const CanonicalForm& F, const ExtensionInfo& info, CFList& source,
195         CFList& dest)
[7bf145]196{
197  int k= info.getGFDegree();
198  Variable beta= info.getAlpha();
199  CanonicalForm primElem= info.getGamma();
200  CanonicalForm imPrimElem= info.getDelta();
201  if (k > 1)
202    return GFMapDown (F, k);
203  else if (k == 1)
204    return F;
[c1b9927]205  if (/*k==0 &&*/ beta == Variable (1))
[7bf145]206    return F;
[c1b9927]207  else /*if (k==0 && beta != Variable (1))*/
[618da5]208    return mapDown (F, imPrimElem, primElem, beta, source, dest);
[7bf145]209}
210
[806c18]211void appendTestMapDown (CFList& factors, const CanonicalForm& f,
212                        const ExtensionInfo& info, CFList& source, CFList& dest)
[7bf145]213{
214  int k= info.getGFDegree();
215  Variable beta= info.getBeta();
216  Variable alpha= info.getAlpha();
217  CanonicalForm delta= info.getDelta();
218  CanonicalForm gamma= info.getGamma();
219  CanonicalForm g= f;
220  int degMipoBeta;
[1372ae]221  if (!k && beta.level() == 1)
[7bf145]222    degMipoBeta= 1;
[1372ae]223  else if (!k && beta.level() != 1)
[7bf145]224    degMipoBeta= degree (getMipo (beta));
[806c18]225  if (k > 1)
[7bf145]226  {
[f876a66]227    if (!isInExtension (g, gamma, k, delta, source, dest))
[7bf145]228    {
229      g= GFMapDown (g, k);
230      factors.append (g);
231    }
232  }
[806c18]233  else if (k == 1)
234  {
[1709e1]235    if (!isInExtension (g, gamma, k, delta, source, dest))
[7bf145]236      factors.append (g);
237  }
[806c18]238  else if (!k && beta == Variable (1))
[7bf145]239  {
240    if (degree (g, alpha) < degMipoBeta)
241      factors.append (g);
242  }
[806c18]243  else if (!k && beta != Variable (1))
[7bf145]244  {
[f876a66]245    if (!isInExtension (g, gamma, k, delta, source, dest))
[7bf145]246    {
[618da5]247      g= mapDown (g, delta, gamma, alpha, source, dest);
[7bf145]248      factors.append (g);
[806c18]249    }
[7bf145]250  }
251  return;
252}
253
[806c18]254void
255appendMapDown (CFList& factors, const CanonicalForm& g,
256               const ExtensionInfo& info, CFList& source, CFList& dest)
[7bf145]257{
258  int k= info.getGFDegree();
259  Variable beta= info.getBeta();
260  Variable alpha= info.getAlpha();
261  CanonicalForm gamma= info.getGamma();
262  CanonicalForm delta= info.getDelta();
[806c18]263  if (k > 1)
[7bf145]264    factors.append (GFMapDown (g, k));
[806c18]265  else if (k == 1)
[7bf145]266    factors.append (g);
[806c18]267  else if (!k && beta == Variable (1))
[7bf145]268    factors.append (g);
[806c18]269  else if (!k && beta != Variable (1))
[618da5]270    factors.append (mapDown (g, delta, gamma, alpha, source, dest));
[7bf145]271  return;
272}
273
[806c18]274void normalize (CFList& factors)
[7bf145]275{
[168269]276  CanonicalForm lcinv;
[806c18]277  for (CFListIterator i= factors; i.hasItem(); i++)
[168269]278  {
279    lcinv= 1/Lc (i.getItem());
280    i.getItem() *= lcinv;
281  }
[7bf145]282  return;
283}
284
[f876a66]285void normalize (CFFList& factors)
286{
[168269]287  CanonicalForm lcinv;
[c1b9927]288  for (CFFListIterator i= factors; i.hasItem(); i++)
[168269]289  {
290    lcinv= 1/ Lc (i.getItem().factor());
291    i.getItem()= CFFactor (i.getItem().factor()*lcinv,
[f876a66]292                           i.getItem().exp());
[168269]293  }
[f876a66]294  return;
295}
296
[806c18]297CFList subset (int index [], const int& s, const CFArray& elements,
298               bool& noSubset)
[7bf145]299{
300  int r= elements.size();
301  int i= 0;
302  CFList result;
303  noSubset= false;
[806c18]304  if (index[s - 1] == 0)
[7bf145]305  {
[806c18]306    while (i < s)
[7bf145]307    {
308      index[i]= i + 1;
309      result.append (elements[i]);
310      i++;
[806c18]311    }
[7bf145]312    return result;
313  }
314  int buf;
315  int k;
316  bool found= false;
[806c18]317  if (index[s - 1] == r)
[7bf145]318  {
[806c18]319    if (index[0] == r - s + 1)
[7bf145]320    {
321      noSubset= true;
322      return result;
323    }
324    else {
[806c18]325      while (found == false)
[7bf145]326      {
327        if (index[s - 2 - i] < r - i - 1)
328          found= true;
[806c18]329        i++;
[7bf145]330      }
331      buf= index[s - i - 1];
332      k= 0;
[806c18]333      while (s - i - 1 + k < s)
[7bf145]334      {
335        index[s - i - 1 + k]= buf + k + 1;
336        k++;
[806c18]337      }
[7bf145]338    }
[806c18]339    for (int j= 0; j < s; j++)
[7bf145]340      result.append (elements[index[j] - 1]);
341    return result;
342  }
[806c18]343  else
[7bf145]344  {
345    index[s - 1] += 1;
346    for (int j= 0; j < s; j++)
347      result.append (elements[index[j] - 1]);
348    return result;
349  }
[806c18]350}
[7bf145]351
[806c18]352CFArray copy (const CFList& list)
[7bf145]353{
354  CFArray array= CFArray (list.length());
355  int j= 0;
[806c18]356  for (CFListIterator i= list; i.hasItem(); i++, j++)
[7bf145]357    array[j]= i.getItem();
358  return array;
359}
360
[806c18]361void indexUpdate (int index [], const int& subsetSize, const int& setSize,
362                   bool& noSubset)
[7bf145]363{
364  noSubset= false;
[806c18]365  if (subsetSize > setSize)
[7bf145]366  {
367    noSubset= true;
368    return;
369  }
[1673386]370  int * v= new int [setSize];
[806c18]371  for (int i= 0; i < setSize; i++)
[7bf145]372    v[i]= index[i];
[806c18]373  if (subsetSize == 1)
[7bf145]374  {
375    v[0]= v[0] - 1;
[806c18]376    if (v[0] >= setSize)
[7bf145]377    {
378      noSubset= true;
[1673386]379      delete [] v;
[7bf145]380      return;
381    }
382  }
[806c18]383  else
[7bf145]384  {
[806c18]385    if (v[subsetSize - 1] - v[0] + 1 == subsetSize && v[0] > 1)
[7bf145]386    {
[806c18]387      if (v[0] + subsetSize - 1 > setSize)
388      {
[7bf145]389        noSubset= true;
[1673386]390        delete [] v;
[7bf145]391        return;
392      }
393      v[0]= v[0] - 1;
[806c18]394      for (int i= 1; i < subsetSize - 1; i++)
[7bf145]395        v[i]= v[i - 1] + 1;
396      v[subsetSize - 1]= v[subsetSize - 2];
397    }
[806c18]398    else
[7bf145]399    {
[806c18]400      if (v[0] + subsetSize - 1 > setSize)
401      {
[7bf145]402        noSubset= true;
[1673386]403        delete [] v;
[7bf145]404        return;
405      }
406      for (int i= 1; i < subsetSize - 1; i++)
407        v[i]= v[i - 1] + 1;
408      v[subsetSize - 1]= v[subsetSize - 2];
409    }
410  }
411  for (int i= 0; i < setSize; i++)
412    index[i]= v[i];
[1673386]413  delete [] v;
[7bf145]414}
415
[806c18]416int subsetDegree (const CFList& S)
[7bf145]417{
418  int result= 0;
419  for (CFListIterator i= S; i.hasItem(); i++)
420    result += degree (i.getItem(), Variable (1));
421  return result;
422}
423
424CFFList multiplicity (CanonicalForm& F, const CFList& factors)
425{
426  if (F.inCoeffDomain())
427    return CFFList (CFFactor (F, 1));
428  CFFList result;
429  int multi= 0;
[21b8f4]430  CanonicalForm quot;
[7bf145]431  for (CFListIterator i= factors; i.hasItem(); i++)
432  {
[21b8f4]433    while (fdivides (i.getItem(), F, quot))
[7bf145]434    {
435      multi++;
[21b8f4]436      F= quot;
[7bf145]437    }
[575e1d]438    if (multi > 0)
439      result.append (CFFactor (i.getItem(), multi));
[7bf145]440    multi= 0;
441  }
442  return result;
443}
444
[615ca8]445#ifdef HAVE_NTL
[f876a66]446CFArray
447logarithmicDerivative (const CanonicalForm& F, const CanonicalForm& G, int l,
448                       CanonicalForm& Q
449                      )
450{
451  Variable x= Variable (2);
452  Variable y= Variable (1);
453  CanonicalForm xToL= power (x, l);
454  CanonicalForm q,r;
455  CanonicalForm logDeriv;
456
[5e4636]457  TIMING_START (fac_log_deriv_div);
[f876a66]458  q= newtonDiv (F, G, xToL);
[5e4636]459  TIMING_END_AND_PRINT (fac_log_deriv_div, "time for division in logderiv1: ");
[f876a66]460
[5e4636]461  TIMING_START (fac_log_deriv_mul);
[f876a66]462  logDeriv= mulMod2 (q, deriv (G, y), xToL);
[5e4636]463  TIMING_END_AND_PRINT (fac_log_deriv_mul, "time to multiply in logderiv1: ");
[f876a66]464
[9614bb]465  int j= degree (logDeriv, y) + 1;
[f876a66]466  CFArray result= CFArray (j);
[9614bb]467  CFIterator ii;
[f876a66]468  for (CFIterator i= logDeriv; i.hasTerms() && !logDeriv.isZero(); i++)
469  {
[0b447e]470    if (i.coeff().inCoeffDomain())
471      result[0] += i.coeff()*power (x,i.exp());
472    else
473    {
474      for (ii= i.coeff(); ii.hasTerms(); ii++)
475        result[ii.exp()] += ii.coeff()*power (x,i.exp());
476    }
[f876a66]477  }
478  Q= q;
479  return result;
480}
481
482CFArray
483logarithmicDerivative (const CanonicalForm& F, const CanonicalForm& G, int l,
484                       int oldL, const CanonicalForm& oldQ, CanonicalForm& Q
485                      )
486{
487  Variable x= Variable (2);
488  Variable y= Variable (1);
489  CanonicalForm xToL= power (x, l);
490  CanonicalForm xToOldL= power (x, oldL);
491  CanonicalForm xToLOldL= power (x, l-oldL);
492  CanonicalForm q,r;
493  CanonicalForm logDeriv;
494
[ac3fcca]495  CanonicalForm bufF;
[5e4636]496  TIMING_START (fac_log_deriv_pre);
[ac3fcca]497  if ((oldL > 100 && l - oldL < 50) || (oldL < 100 && l - oldL < 30))
498  {
499    bufF= F;
500    CanonicalForm oldF= mulMod2 (G, oldQ, xToL);
501    bufF -= oldF;
502    bufF= div (bufF, xToOldL);
503  }
504  else
505  {
506    //middle product style computation of [G*oldQ]^{l}_{oldL}
507    CanonicalForm G3= div (G, xToOldL);
508    CanonicalForm Up= mulMod2 (G3, oldQ, xToLOldL);
[fd80670]509    CanonicalForm xToOldL2= power (x, (oldL+1)/2);
[ac3fcca]510    CanonicalForm G2= mod (G, xToOldL);
511    CanonicalForm G1= div (G2, xToOldL2);
512    CanonicalForm G0= mod (G2, xToOldL2);
513    CanonicalForm oldQ1= div (oldQ, xToOldL2);
514    CanonicalForm oldQ0= mod (oldQ, xToOldL2);
[fd80670]515    CanonicalForm Mid;
516    if (oldL % 2 == 1)
517      Mid= mulMod2 (G1, oldQ1*x, xToLOldL);
518    else
519      Mid= mulMod2 (G1, oldQ1, xToLOldL);
[ac3fcca]520    //computation of Low might be faster using a real middle product?
521    CanonicalForm Low= mulMod2 (G0, oldQ1, xToOldL)+mulMod2 (G1, oldQ0, xToOldL);
[fd80670]522    Low= div (Low, power (x, oldL/2));
523    Low= mod (Low, xToLOldL);
[ac3fcca]524    Up += Mid + Low;
525    bufF= div (F, xToOldL);
526    bufF -= Up;
527  }
[5e4636]528  TIMING_END_AND_PRINT (fac_log_deriv_pre, "time to preprocess: ");
[f876a66]529
[5e4636]530  TIMING_START (fac_log_deriv_div);
[72f1e4b]531  if (l-oldL > 0)
532    q= newtonDiv (bufF, G, xToLOldL);
533  else
534    q= 0;
[f876a66]535  q *= xToOldL;
536  q += oldQ;
[5e4636]537  TIMING_END_AND_PRINT (fac_log_deriv_div, "time for div in logderiv2: ");
[f876a66]538
[5e4636]539  TIMING_START (fac_log_deriv_mul);
[f876a66]540  logDeriv= mulMod2 (q, deriv (G, y), xToL);
[5e4636]541  TIMING_END_AND_PRINT (fac_log_deriv_mul, "time for mul in logderiv2: ");
[f876a66]542
[9614bb]543  int j= degree (logDeriv,y) + 1;
[f876a66]544  CFArray result= CFArray (j);
[9614bb]545  CFIterator ii;
[f876a66]546  for (CFIterator i= logDeriv; i.hasTerms() && !logDeriv.isZero(); i++)
547  {
[0b447e]548    if (i.coeff().inCoeffDomain())
549      result[0] += i.coeff()*power (x,i.exp());
550    else
551    {
552      for (ii= i.coeff(); ii.hasTerms(); ii++)
553        result[ii.exp()] += ii.coeff()*power (x,i.exp());
554    }
[f876a66]555  }
556  Q= q;
557  return result;
558}
[615ca8]559#endif
[f876a66]560
561void
562writeInMatrix (CFMatrix& M, const CFArray& A, const int column,
563               const int startIndex
564              )
565{
[050d1b]566  ASSERT (A.size () - startIndex >= 0, "wrong starting index");
567  ASSERT (A.size () - startIndex <= M.rows(), "wrong starting index");
[f876a66]568  ASSERT (column > 0 && column <= M.columns(), "wrong column");
569  if (A.size() - startIndex <= 0) return;
570  int j= 1;
571  for (int i= startIndex; i < A.size(); i++, j++)
572    M (j, column)= A [i];
573}
574
575CFArray getCoeffs (const CanonicalForm& F, const int k)
576{
[050d1b]577  ASSERT (F.isUnivariate() || F.inCoeffDomain(), "univariate input expected");
[f876a66]578  if (degree (F, 2) < k)
579    return CFArray();
580
581  CFArray result= CFArray (degree (F) - k + 1);
582  CFIterator j= F;
583  for (int i= degree (F); i >= k; i--)
584  {
585    if (j.exp() == i)
586    {
587      result [i - k]= j.coeff();
[c1b9927]588      j++;
[f876a66]589      if (!j.hasTerms())
590        return result;
591    }
592    else
593      result[i - k]= 0;
594  }
595  return result;
596}
597
598CFArray getCoeffs (const CanonicalForm& F, const int k, const Variable& alpha)
599{
[050d1b]600  ASSERT (F.isUnivariate() || F.inCoeffDomain(), "univariate input expected");
[f876a66]601  if (degree (F, 2) < k)
602    return CFArray ();
603
604  int d= degree (getMipo (alpha));
605  CFArray result= CFArray ((degree (F) - k + 1)*d);
606  CFIterator j= F;
607  CanonicalForm buf;
608  CFIterator iter;
609  for (int i= degree (F); i >= k; i--)
610  {
611    if (j.exp() == i)
612    {
613      iter= j.coeff();
614      for (int l= degree (j.coeff(), alpha); l >= 0; l--)
615      {
616        if (iter.exp() == l)
617        {
618          result [(i - k)*d + l]= iter.coeff();
619          iter++;
620          if (!iter.hasTerms())
621            break;
622        }
623      }
624      j++;
625      if (!j.hasTerms())
626        return result;
627    }
628    else
629    {
630      for (int l= 0; l < d; l++)
631        result[(i - k)*d + l]= 0;
632    }
633  }
634  return result;
635}
636
637#ifdef HAVE_NTL
638CFArray
639getCoeffs (const CanonicalForm& G, const int k, const int l, const int degMipo,
640           const Variable& alpha, const CanonicalForm& evaluation,
641           const mat_zz_p& M)
642{
[050d1b]643  ASSERT (G.isUnivariate() || G.inCoeffDomain(), "univariate input expected");
[f876a66]644  CanonicalForm F= G (G.mvar() - evaluation, G.mvar());
645  if (F.isZero())
646    return CFArray ();
647
648  Variable y= Variable (2);
649  F= F (power (y, degMipo), y);
650  F= F (y, alpha);
651  zz_pX NTLF= convertFacCF2NTLzzpX (F);
652  NTLF.rep.SetLength (l*degMipo);
653  NTLF.rep= M*NTLF.rep;
654  NTLF.normalize();
655  F= convertNTLzzpX2CF (NTLF, y);
656
657  if (degree (F, 2) < k)
658    return CFArray();
659
660  CFArray result= CFArray (degree (F) - k + 1);
661
662  CFIterator j= F;
663  for (int i= degree (F); i >= k; i--)
664  {
665    if (j.exp() == i)
666    {
667      result [i - k]= j.coeff();
[c1b9927]668      j++;
[f876a66]669      if (!j.hasTerms())
670        return result;
671    }
672    else
673      result[i - k]= 0;
674  }
675  return result;
676}
677#endif
678
[542864]679int * computeBounds (const CanonicalForm& F, int& n, bool& isIrreducible)
[f876a66]680{
681  n= degree (F, 1);
682  int* result= new int [n];
683  int sizeOfNewtonPolygon;
684  int** newtonPolyg= newtonPolygon (F, sizeOfNewtonPolygon);
685
[542864]686  isIrreducible= false;
687  if (sizeOfNewtonPolygon == 3)
688  {
689    bool check1=
690        (newtonPolyg[0][0]==0 || newtonPolyg[1][0]==0 || newtonPolyg[2][0]==0);
691    if (check1)
692    {
693      bool check2=
694        (newtonPolyg[0][1]==0 || newtonPolyg[1][1]==0 || newtonPolyg[2][0]==0);
695      if (check2)
696      {
697        int p=getCharacteristic();
698        int d=1;
699        char bufGFName='Z';
700        bool GF= (CFFactory::gettype()==GaloisFieldDomain);
701        if (GF)
702        {
703          d= getGFDegree();
704          bufGFName=gf_name;
705        }
706        setCharacteristic(0);
707        CanonicalForm tmp= gcd (newtonPolyg[0][0],newtonPolyg[0][1]);
708        tmp= gcd (tmp, newtonPolyg[1][0]);
709        tmp= gcd (tmp, newtonPolyg[1][1]);
710        tmp= gcd (tmp, newtonPolyg[2][0]);
711        tmp= gcd (tmp, newtonPolyg[2][1]);
712        isIrreducible= (tmp==1);
713        if (GF)
714          setCharacteristic (p, d, bufGFName);
715        else
716          setCharacteristic(p);
717      }
718    }
719  }
720
[f876a66]721  int minX, minY, maxX, maxY;
722  minX= newtonPolyg [0] [0];
723  minY= newtonPolyg [0] [1];
724  maxX= minX;
725  maxY= minY;
726  for (int i= 1; i < sizeOfNewtonPolygon; i++)
727  {
728    if (minX > newtonPolyg [i] [0])
729      minX= newtonPolyg [i] [0];
730    if (maxX < newtonPolyg [i] [0])
731      maxX= newtonPolyg [i] [0];
732    if (minY > newtonPolyg [i] [1])
733      minY= newtonPolyg [i] [1];
734    if (maxY < newtonPolyg [i] [1])
735      maxY= newtonPolyg [i] [1];
736  }
737
738  int k= maxX;
739  for (int i= 0; i < n; i++)
740  {
741    if (i + 1 > maxY || i + 1 < minY)
742    {
743      result [i]= 0;
744      continue;
745    }
746    int* point= new int [2];
747    point [0]= k;
748    point [1]= i + 1;
749    while (!isInPolygon (newtonPolyg, sizeOfNewtonPolygon, point) && k > 0)
750    {
751      k--;
752      point [0]= k;
753    }
754    result [i]= k;
755    k= maxX;
756    delete [] point;
757  }
758
759  return result;
760}
761
[d8a7da]762int *
763computeBoundsWrtDiffMainvar (const CanonicalForm& F, int& n,
764                             bool& isIrreducible)
765{
766  n= degree (F, 2);
767  int* result= new int [n];
768  int sizeOfNewtonPolygon;
769  int** newtonPolyg= newtonPolygon (F, sizeOfNewtonPolygon);
770
771  isIrreducible= false;
772  if (sizeOfNewtonPolygon == 3)
773  {
774    bool check1=
775        (newtonPolyg[0][0]==0 || newtonPolyg[1][0]==0 || newtonPolyg[2][0]==0);
776    if (check1)
777    {
778      bool check2=
779        (newtonPolyg[0][1]==0 || newtonPolyg[1][1]==0 || newtonPolyg[2][0]==0);
780      if (check2)
781      {
782        int p=getCharacteristic();
783        int d=1;
784        char bufGFName='Z';
785        bool GF= (CFFactory::gettype()==GaloisFieldDomain);
786        if (GF)
787        {
788          d= getGFDegree();
789          bufGFName=gf_name;
790        }
791        setCharacteristic(0);
792        CanonicalForm tmp= gcd (newtonPolyg[0][0],newtonPolyg[0][1]);
793        tmp= gcd (tmp, newtonPolyg[1][0]);
794        tmp= gcd (tmp, newtonPolyg[1][1]);
795        tmp= gcd (tmp, newtonPolyg[2][0]);
796        tmp= gcd (tmp, newtonPolyg[2][1]);
797        isIrreducible= (tmp==1);
798        if (GF)
799          setCharacteristic (p, d, bufGFName);
800        else
801          setCharacteristic(p);
802      }
803    }
804  }
805
806  int minX, minY, maxX, maxY;
807  minX= newtonPolyg [0] [0];
808  minY= newtonPolyg [0] [1];
809  maxX= minX;
810  maxY= minY;
811  for (int i= 1; i < sizeOfNewtonPolygon; i++)
812  {
813    if (minX > newtonPolyg [i] [0])
814      minX= newtonPolyg [i] [0];
815    if (maxX < newtonPolyg [i] [0])
816      maxX= newtonPolyg [i] [0];
817    if (minY > newtonPolyg [i] [1])
818      minY= newtonPolyg [i] [1];
819    if (maxY < newtonPolyg [i] [1])
820      maxY= newtonPolyg [i] [1];
821  }
822
823  int k= maxY;
824  for (int i= 0; i < n; i++)
825  {
826    if (i + 1 > maxX || i + 1 < minX)
827    {
828      result [i]= 0;
829      continue;
830    }
831    int* point= new int [2];
832    point [0]= i + 1;
833    point [1]= k;
834    while (!isInPolygon (newtonPolyg, sizeOfNewtonPolygon, point) && k > 0)
835    {
836      k--;
837      point [1]= k;
838    }
839    result [i]= k;
840    k= maxY;
841    delete [] point;
842  }
843
844  return result;
845}
846
[dd8047]847int
848substituteCheck (const CanonicalForm& F, const Variable& x)
849{
850  if (F.inCoeffDomain())
851    return 0;
852  if (degree (F, x) < 0)
853    return 0;
854  CanonicalForm f= swapvar (F, F.mvar(), x);
855  int sizef= 0;
856  for (CFIterator i= f; i.hasTerms(); i++, sizef++)
857  {
858    if (i.exp() == 1)
859      return 0;
860  }
861  int * expf= new int [sizef];
862  int j= 0;
863  for (CFIterator i= f; i.hasTerms(); i++, j++)
864    expf [j]= i.exp();
865
866  int indf= sizef - 1;
867  if (expf[indf] == 0)
868    indf--;
869
870  int result= expf[indf];
871  for (int i= indf - 1; i >= 0; i--)
872  {
873    if (expf [i]%result != 0)
874    {
875      delete [] expf;
876      return 0;
877    }
878  }
879
880  delete [] expf;
881  return result;
882}
883
[686ce3]884static int
885substituteCheck (const CanonicalForm& F, const CanonicalForm& G)
886{
887  if (F.inCoeffDomain() || G.inCoeffDomain())
888    return 0;
889  Variable x= Variable (1);
890  if (degree (F, x) <= 1 || degree (G, x) <= 1)
891    return 0;
892  CanonicalForm f= swapvar (F, F.mvar(), x);
893  CanonicalForm g= swapvar (G, G.mvar(), x);
894  int sizef= 0;
895  int sizeg= 0;
896  for (CFIterator i= f; i.hasTerms(); i++, sizef++)
897  {
898    if (i.exp() == 1)
899      return 0;
900  }
901  for (CFIterator i= g; i.hasTerms(); i++, sizeg++)
902  {
903    if (i.exp() == 1)
904      return 0;
905  }
906  int * expf= new int [sizef];
907  int * expg= new int [sizeg];
908  int j= 0;
909  for (CFIterator i= f; i.hasTerms(); i++, j++)
910  {
911    expf [j]= i.exp();
912  }
913  j= 0;
914  for (CFIterator i= g; i.hasTerms(); i++, j++)
915  {
916    expg [j]= i.exp();
917  }
918
919  int indf= sizef - 1;
920  int indg= sizeg - 1;
921  if (expf[indf] == 0)
922    indf--;
923  if (expg[indg] == 0)
924    indg--;
925
926  if ((expg[indg]%expf [indf] != 0 && expf[indf]%expg[indg] != 0) ||
927      (expg[indg] == 1 && expf[indf] == 1))
928  {
929    delete [] expg;
930    delete [] expf;
931    return 0;
932  }
933
934  int result;
935  if (expg [indg]%expf [indf] == 0)
936    result= expf[indf];
937  else
938    result= expg[indg];
939  for (int i= indf - 1; i >= 0; i--)
940  {
941    if (expf [i]%result != 0)
942    {
943      delete [] expf;
944      delete [] expg;
945      return 0;
946    }
947  }
948
949  for (int i= indg - 1; i >= 0; i--)
950  {
951    if (expg [i]%result != 0)
952    {
953      delete [] expf;
954      delete [] expg;
955      return 0;
956    }
957  }
958
959  delete [] expg;
960  delete [] expf;
961  return result;
962}
963
964int recSubstituteCheck (const CanonicalForm& F, const int d)
965{
966  if (F.inCoeffDomain())
967    return 0;
968  Variable x= Variable (1);
969  if (degree (F, x) <= 1)
970    return 0;
971  CanonicalForm f= swapvar (F, F.mvar(), x);
972  int sizef= 0;
973  for (CFIterator i= f; i.hasTerms(); i++, sizef++)
974  {
975    if (i.exp() == 1)
976      return 0;
977  }
978  int * expf= new int [sizef];
979  int j= 0;
980  for (CFIterator i= f; i.hasTerms(); i++, j++)
981  {
982    expf [j]= i.exp();
983  }
984
985  int indf= sizef - 1;
986  if (expf[indf] == 0)
987    indf--;
988
989  if ((d%expf [indf] != 0 && expf[indf]%d != 0) || (expf[indf] == 1))
990  {
991    delete [] expf;
992    return 0;
993  }
994
995  int result;
996  if (d%expf [indf] == 0)
997    result= expf[indf];
998  else
999    result= d;
1000  for (int i= indf - 1; i >= 0; i--)
1001  {
1002    if (expf [i]%result != 0)
1003    {
1004      delete [] expf;
1005      return 0;
1006    }
1007  }
1008
1009  delete [] expf;
1010  return result;
1011}
1012
1013int substituteCheck (const CFList& L)
1014{
1015  ASSERT (L.length() > 1, "expected a list of at least two elements");
1016  if (L.length() < 2)
1017    return 0;
1018  CFListIterator i= L;
1019  i++;
1020  int result= substituteCheck (L.getFirst(), i.getItem());
1021  if (result <= 1)
1022    return result;
1023  i++;
1024  for (;i.hasItem(); i++)
1025  {
1026    result= recSubstituteCheck (i.getItem(), result);
1027    if (result <= 1)
1028      return result;
1029  }
1030  return result;
1031}
1032
1033void
1034subst (const CanonicalForm& F, CanonicalForm& A, const int d, const Variable& x)
1035{
1036  if (d <= 1)
1037  {
1038    A= F;
1039    return;
1040  }
1041  if (degree (F, x) <= 0)
1042  {
1043    A= F;
1044    return;
1045  }
1046  CanonicalForm C= 0;
1047  CanonicalForm f= swapvar (F, x, F.mvar());
1048  for (CFIterator i= f; i.hasTerms(); i++)
1049    C += i.coeff()*power (f.mvar(), i.exp()/ d);
1050  A= swapvar (C, x, F.mvar());
1051}
1052
1053CanonicalForm
1054reverseSubst (const CanonicalForm& F, const int d, const Variable& x)
1055{
1056  if (d <= 1)
1057    return F;
1058  if (degree (F, x) <= 0)
1059    return F;
1060  CanonicalForm f= swapvar (F, x, F.mvar());
1061  CanonicalForm result= 0;
1062  for (CFIterator i= f; i.hasTerms(); i++)
1063    result += i.coeff()*power (f.mvar(), d*i.exp());
1064  return swapvar (result, x, F.mvar());
1065}
1066
1067void
1068reverseSubst (CFList& L, const int d, const Variable& x)
1069{
1070  for (CFListIterator i= L; i.hasItem(); i++)
1071    i.getItem()= reverseSubst (i.getItem(), d, x);
1072}
1073
Note: See TracBrowser for help on using the repository browser.