source: git/factory/facFqBivarUtil.cc @ 72bfc8

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