source: git/factory/facFqBivarUtil.cc @ 542864

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