source: git/factory/facFqBivarUtil.cc @ fd80670

jengelh-datetimespielwiese
Last change on this file since fd80670 was fd80670, checked in by Martin Lee <martinlee84@…>, 10 years ago
fix: bug in divrem2 and logDeriv
  • Property mode set to 100644
File size: 21.7 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  CanonicalForm lcinv;
272  for (CFListIterator i= factors; i.hasItem(); i++)
273  {
274    lcinv= 1/Lc (i.getItem());
275    i.getItem() *= lcinv;
276  }
277  return;
278}
279
280void normalize (CFFList& factors)
281{
282  CanonicalForm lcinv;
283  for (CFFListIterator i= factors; i.hasItem(); i++)
284  {
285    lcinv= 1/ Lc (i.getItem().factor());
286    i.getItem()= CFFactor (i.getItem().factor()*lcinv,
287                           i.getItem().exp());
288  }
289  return;
290}
291
292CFList subset (int index [], const int& s, const CFArray& elements,
293               bool& noSubset)
294{
295  int r= elements.size();
296  int i= 0;
297  CFList result;
298  noSubset= false;
299  if (index[s - 1] == 0)
300  {
301    while (i < s)
302    {
303      index[i]= i + 1;
304      result.append (elements[i]);
305      i++;
306    }
307    return result;
308  }
309  int buf;
310  int k;
311  bool found= false;
312  if (index[s - 1] == r)
313  {
314    if (index[0] == r - s + 1)
315    {
316      noSubset= true;
317      return result;
318    }
319    else {
320      while (found == false)
321      {
322        if (index[s - 2 - i] < r - i - 1)
323          found= true;
324        i++;
325      }
326      buf= index[s - i - 1];
327      k= 0;
328      while (s - i - 1 + k < s)
329      {
330        index[s - i - 1 + k]= buf + k + 1;
331        k++;
332      }
333    }
334    for (int j= 0; j < s; j++)
335      result.append (elements[index[j] - 1]);
336    return result;
337  }
338  else
339  {
340    index[s - 1] += 1;
341    for (int j= 0; j < s; j++)
342      result.append (elements[index[j] - 1]);
343    return result;
344  }
345}
346
347CFArray copy (const CFList& list)
348{
349  CFArray array= CFArray (list.length());
350  int j= 0;
351  for (CFListIterator i= list; i.hasItem(); i++, j++)
352    array[j]= i.getItem();
353  return array;
354}
355
356void indexUpdate (int index [], const int& subsetSize, const int& setSize,
357                   bool& noSubset)
358{
359  noSubset= false;
360  if (subsetSize > setSize)
361  {
362    noSubset= true;
363    return;
364  }
365  int * v= new int [setSize];
366  for (int i= 0; i < setSize; i++)
367    v[i]= index[i];
368  if (subsetSize == 1)
369  {
370    v[0]= v[0] - 1;
371    if (v[0] >= setSize)
372    {
373      noSubset= true;
374      delete [] v;
375      return;
376    }
377  }
378  else
379  {
380    if (v[subsetSize - 1] - v[0] + 1 == subsetSize && v[0] > 1)
381    {
382      if (v[0] + subsetSize - 1 > setSize)
383      {
384        noSubset= true;
385        delete [] v;
386        return;
387      }
388      v[0]= v[0] - 1;
389      for (int i= 1; i < subsetSize - 1; i++)
390        v[i]= v[i - 1] + 1;
391      v[subsetSize - 1]= v[subsetSize - 2];
392    }
393    else
394    {
395      if (v[0] + subsetSize - 1 > setSize)
396      {
397        noSubset= true;
398        delete [] v;
399        return;
400      }
401      for (int i= 1; i < subsetSize - 1; i++)
402        v[i]= v[i - 1] + 1;
403      v[subsetSize - 1]= v[subsetSize - 2];
404    }
405  }
406  for (int i= 0; i < setSize; i++)
407    index[i]= v[i];
408  delete [] v;
409}
410
411int subsetDegree (const CFList& S)
412{
413  int result= 0;
414  for (CFListIterator i= S; i.hasItem(); i++)
415    result += degree (i.getItem(), Variable (1));
416  return result;
417}
418
419CFFList multiplicity (CanonicalForm& F, const CFList& factors)
420{
421  if (F.inCoeffDomain())
422    return CFFList (CFFactor (F, 1));
423  CFFList result;
424  int multi= 0;
425  CanonicalForm quot;
426  for (CFListIterator i= factors; i.hasItem(); i++)
427  {
428    while (fdivides (i.getItem(), F, quot))
429    {
430      multi++;
431      F= quot;
432    }
433    if (multi > 0)
434      result.append (CFFactor (i.getItem(), multi));
435    multi= 0;
436  }
437  return result;
438}
439
440#ifdef HAVE_NTL
441CFArray
442logarithmicDerivative (const CanonicalForm& F, const CanonicalForm& G, int l,
443                       CanonicalForm& Q
444                      )
445{
446  Variable x= Variable (2);
447  Variable y= Variable (1);
448  CanonicalForm xToL= power (x, l);
449  CanonicalForm q,r;
450  CanonicalForm logDeriv;
451
452  q= newtonDiv (F, G, xToL);
453
454  logDeriv= mulMod2 (q, deriv (G, y), xToL);
455
456  logDeriv= swapvar (logDeriv, x, y);
457  int j= degree (logDeriv) + 1;
458  CFArray result= CFArray (j);
459  for (CFIterator i= logDeriv; i.hasTerms() && !logDeriv.isZero(); i++)
460  {
461    if (i.exp() == j - 1)
462    {
463      result [j - 1]= swapvar (i.coeff(), x, y);
464      j--;
465    }
466    else
467    {
468      for (; i.exp() < j - 1; j--)
469        result [j - 1]= 0;
470      result [j - 1]= swapvar (i.coeff(), x, y);
471      j--;
472    }
473  }
474  for (; j > 0; j--)
475    result [j - 1]= 0;
476  Q= q;
477  return result;
478}
479
480CFArray
481logarithmicDerivative (const CanonicalForm& F, const CanonicalForm& G, int l,
482                       int oldL, const CanonicalForm& oldQ, CanonicalForm& Q
483                      )
484{
485  Variable x= Variable (2);
486  Variable y= Variable (1);
487  CanonicalForm xToL= power (x, l);
488  CanonicalForm xToOldL= power (x, oldL);
489  CanonicalForm xToLOldL= power (x, l-oldL);
490  CanonicalForm q,r;
491  CanonicalForm logDeriv;
492
493  CanonicalForm bufF;
494  if ((oldL > 100 && l - oldL < 50) || (oldL < 100 && l - oldL < 30))
495  {
496    bufF= F;
497    CanonicalForm oldF= mulMod2 (G, oldQ, xToL);
498    bufF -= oldF;
499    bufF= div (bufF, xToOldL);
500  }
501  else
502  {
503    //middle product style computation of [G*oldQ]^{l}_{oldL}
504    CanonicalForm G3= div (G, xToOldL);
505    CanonicalForm Up= mulMod2 (G3, oldQ, xToLOldL);
506    CanonicalForm xToOldL2= power (x, (oldL+1)/2);
507    CanonicalForm G2= mod (G, xToOldL);
508    CanonicalForm G1= div (G2, xToOldL2);
509    CanonicalForm G0= mod (G2, xToOldL2);
510    CanonicalForm oldQ1= div (oldQ, xToOldL2);
511    CanonicalForm oldQ0= mod (oldQ, xToOldL2);
512    CanonicalForm Mid;
513    if (oldL % 2 == 1)
514      Mid= mulMod2 (G1, oldQ1*x, xToLOldL);
515    else
516      Mid= mulMod2 (G1, oldQ1, xToLOldL);
517    //computation of Low might be faster using a real middle product?
518    CanonicalForm Low= mulMod2 (G0, oldQ1, xToOldL)+mulMod2 (G1, oldQ0, xToOldL);
519    Low= div (Low, power (x, oldL/2));
520    Low= mod (Low, xToLOldL);
521    Up += Mid + Low;
522    bufF= div (F, xToOldL);
523    bufF -= Up;
524  }
525
526  if (l-oldL > 0)
527    q= newtonDiv (bufF, G, xToLOldL);
528  else
529    q= 0;
530  q *= xToOldL;
531  q += oldQ;
532
533  logDeriv= mulMod2 (q, deriv (G, y), xToL);
534
535  logDeriv= swapvar (logDeriv, x, y);
536  int j= degree (logDeriv) + 1;
537  CFArray result= CFArray (j);
538  for (CFIterator i= logDeriv; i.hasTerms() && !logDeriv.isZero(); i++)
539  {
540    if (i.exp() == j - 1)
541    {
542      result [j - 1]= swapvar (i.coeff(), x, y);
543      j--;
544    }
545    else
546    {
547      for (; i.exp() < j - 1; j--)
548        result [j - 1]= 0;
549      result [j - 1]= swapvar (i.coeff(), x, y);
550      j--;
551    }
552  }
553  for (; j > 0; j--)
554    result [j - 1]= 0;
555  Q= q;
556  return result;
557}
558#endif
559
560void
561writeInMatrix (CFMatrix& M, const CFArray& A, const int column,
562               const int startIndex
563              )
564{
565  ASSERT (A.size () - startIndex >= 0, "wrong starting index");
566  ASSERT (A.size () - startIndex <= M.rows(), "wrong starting index");
567  ASSERT (column > 0 && column <= M.columns(), "wrong column");
568  if (A.size() - startIndex <= 0) return;
569  int j= 1;
570  for (int i= startIndex; i < A.size(); i++, j++)
571    M (j, column)= A [i];
572}
573
574CFArray getCoeffs (const CanonicalForm& F, const int k)
575{
576  ASSERT (F.isUnivariate() || F.inCoeffDomain(), "univariate input expected");
577  if (degree (F, 2) < k)
578    return CFArray();
579
580  CFArray result= CFArray (degree (F) - k + 1);
581  CFIterator j= F;
582  for (int i= degree (F); i >= k; i--)
583  {
584    if (j.exp() == i)
585    {
586      result [i - k]= j.coeff();
587      j++;
588      if (!j.hasTerms())
589        return result;
590    }
591    else
592      result[i - k]= 0;
593  }
594  return result;
595}
596
597CFArray getCoeffs (const CanonicalForm& F, const int k, const Variable& alpha)
598{
599  ASSERT (F.isUnivariate() || F.inCoeffDomain(), "univariate input expected");
600  if (degree (F, 2) < k)
601    return CFArray ();
602
603  int d= degree (getMipo (alpha));
604  CFArray result= CFArray ((degree (F) - k + 1)*d);
605  CFIterator j= F;
606  CanonicalForm buf;
607  CFIterator iter;
608  for (int i= degree (F); i >= k; i--)
609  {
610    if (j.exp() == i)
611    {
612      iter= j.coeff();
613      for (int l= degree (j.coeff(), alpha); l >= 0; l--)
614      {
615        if (iter.exp() == l)
616        {
617          result [(i - k)*d + l]= iter.coeff();
618          iter++;
619          if (!iter.hasTerms())
620            break;
621        }
622      }
623      j++;
624      if (!j.hasTerms())
625        return result;
626    }
627    else
628    {
629      for (int l= 0; l < d; l++)
630        result[(i - k)*d + l]= 0;
631    }
632  }
633  return result;
634}
635
636#ifdef HAVE_NTL
637CFArray
638getCoeffs (const CanonicalForm& G, const int k, const int l, const int degMipo,
639           const Variable& alpha, const CanonicalForm& evaluation,
640           const mat_zz_p& M)
641{
642  ASSERT (G.isUnivariate() || G.inCoeffDomain(), "univariate input expected");
643  CanonicalForm F= G (G.mvar() - evaluation, G.mvar());
644  if (F.isZero())
645    return CFArray ();
646
647  Variable y= Variable (2);
648  F= F (power (y, degMipo), y);
649  F= F (y, alpha);
650  zz_pX NTLF= convertFacCF2NTLzzpX (F);
651  NTLF.rep.SetLength (l*degMipo);
652  NTLF.rep= M*NTLF.rep;
653  NTLF.normalize();
654  F= convertNTLzzpX2CF (NTLF, y);
655
656  if (degree (F, 2) < k)
657    return CFArray();
658
659  CFArray result= CFArray (degree (F) - k + 1);
660
661  CFIterator j= F;
662  for (int i= degree (F); i >= k; i--)
663  {
664    if (j.exp() == i)
665    {
666      result [i - k]= j.coeff();
667      j++;
668      if (!j.hasTerms())
669        return result;
670    }
671    else
672      result[i - k]= 0;
673  }
674  return result;
675}
676#endif
677
678int * computeBounds (const CanonicalForm& F, int& n, bool& isIrreducible)
679{
680  n= degree (F, 1);
681  int* result= new int [n];
682  int sizeOfNewtonPolygon;
683  int** newtonPolyg= newtonPolygon (F, sizeOfNewtonPolygon);
684
685  isIrreducible= false;
686  if (sizeOfNewtonPolygon == 3)
687  {
688    bool check1=
689        (newtonPolyg[0][0]==0 || newtonPolyg[1][0]==0 || newtonPolyg[2][0]==0);
690    if (check1)
691    {
692      bool check2=
693        (newtonPolyg[0][1]==0 || newtonPolyg[1][1]==0 || newtonPolyg[2][0]==0);
694      if (check2)
695      {
696        int p=getCharacteristic();
697        int d=1;
698        char bufGFName='Z';
699        bool GF= (CFFactory::gettype()==GaloisFieldDomain);
700        if (GF)
701        {
702          d= getGFDegree();
703          bufGFName=gf_name;
704        }
705        setCharacteristic(0);
706        CanonicalForm tmp= gcd (newtonPolyg[0][0],newtonPolyg[0][1]);
707        tmp= gcd (tmp, newtonPolyg[1][0]);
708        tmp= gcd (tmp, newtonPolyg[1][1]);
709        tmp= gcd (tmp, newtonPolyg[2][0]);
710        tmp= gcd (tmp, newtonPolyg[2][1]);
711        isIrreducible= (tmp==1);
712        if (GF)
713          setCharacteristic (p, d, bufGFName);
714        else
715          setCharacteristic(p);
716      }
717    }
718  }
719
720  int minX, minY, maxX, maxY;
721  minX= newtonPolyg [0] [0];
722  minY= newtonPolyg [0] [1];
723  maxX= minX;
724  maxY= minY;
725  for (int i= 1; i < sizeOfNewtonPolygon; i++)
726  {
727    if (minX > newtonPolyg [i] [0])
728      minX= newtonPolyg [i] [0];
729    if (maxX < newtonPolyg [i] [0])
730      maxX= newtonPolyg [i] [0];
731    if (minY > newtonPolyg [i] [1])
732      minY= newtonPolyg [i] [1];
733    if (maxY < newtonPolyg [i] [1])
734      maxY= newtonPolyg [i] [1];
735  }
736
737  int k= maxX;
738  for (int i= 0; i < n; i++)
739  {
740    if (i + 1 > maxY || i + 1 < minY)
741    {
742      result [i]= 0;
743      continue;
744    }
745    int* point= new int [2];
746    point [0]= k;
747    point [1]= i + 1;
748    while (!isInPolygon (newtonPolyg, sizeOfNewtonPolygon, point) && k > 0)
749    {
750      k--;
751      point [0]= k;
752    }
753    result [i]= k;
754    k= maxX;
755    delete [] point;
756  }
757
758  return result;
759}
760
761int
762substituteCheck (const CanonicalForm& F, const Variable& x)
763{
764  if (F.inCoeffDomain())
765    return 0;
766  if (degree (F, x) < 0)
767    return 0;
768  CanonicalForm f= swapvar (F, F.mvar(), x);
769  int sizef= 0;
770  for (CFIterator i= f; i.hasTerms(); i++, sizef++)
771  {
772    if (i.exp() == 1)
773      return 0;
774  }
775  int * expf= new int [sizef];
776  int j= 0;
777  for (CFIterator i= f; i.hasTerms(); i++, j++)
778    expf [j]= i.exp();
779
780  int indf= sizef - 1;
781  if (expf[indf] == 0)
782    indf--;
783
784  int result= expf[indf];
785  for (int i= indf - 1; i >= 0; i--)
786  {
787    if (expf [i]%result != 0)
788    {
789      delete [] expf;
790      return 0;
791    }
792  }
793
794  delete [] expf;
795  return result;
796}
797
798static int
799substituteCheck (const CanonicalForm& F, const CanonicalForm& G)
800{
801  if (F.inCoeffDomain() || G.inCoeffDomain())
802    return 0;
803  Variable x= Variable (1);
804  if (degree (F, x) <= 1 || degree (G, x) <= 1)
805    return 0;
806  CanonicalForm f= swapvar (F, F.mvar(), x);
807  CanonicalForm g= swapvar (G, G.mvar(), x);
808  int sizef= 0;
809  int sizeg= 0;
810  for (CFIterator i= f; i.hasTerms(); i++, sizef++)
811  {
812    if (i.exp() == 1)
813      return 0;
814  }
815  for (CFIterator i= g; i.hasTerms(); i++, sizeg++)
816  {
817    if (i.exp() == 1)
818      return 0;
819  }
820  int * expf= new int [sizef];
821  int * expg= new int [sizeg];
822  int j= 0;
823  for (CFIterator i= f; i.hasTerms(); i++, j++)
824  {
825    expf [j]= i.exp();
826  }
827  j= 0;
828  for (CFIterator i= g; i.hasTerms(); i++, j++)
829  {
830    expg [j]= i.exp();
831  }
832
833  int indf= sizef - 1;
834  int indg= sizeg - 1;
835  if (expf[indf] == 0)
836    indf--;
837  if (expg[indg] == 0)
838    indg--;
839
840  if ((expg[indg]%expf [indf] != 0 && expf[indf]%expg[indg] != 0) ||
841      (expg[indg] == 1 && expf[indf] == 1))
842  {
843    delete [] expg;
844    delete [] expf;
845    return 0;
846  }
847
848  int result;
849  if (expg [indg]%expf [indf] == 0)
850    result= expf[indf];
851  else
852    result= expg[indg];
853  for (int i= indf - 1; i >= 0; i--)
854  {
855    if (expf [i]%result != 0)
856    {
857      delete [] expf;
858      delete [] expg;
859      return 0;
860    }
861  }
862
863  for (int i= indg - 1; i >= 0; i--)
864  {
865    if (expg [i]%result != 0)
866    {
867      delete [] expf;
868      delete [] expg;
869      return 0;
870    }
871  }
872
873  delete [] expg;
874  delete [] expf;
875  return result;
876}
877
878int recSubstituteCheck (const CanonicalForm& F, const int d)
879{
880  if (F.inCoeffDomain())
881    return 0;
882  Variable x= Variable (1);
883  if (degree (F, x) <= 1)
884    return 0;
885  CanonicalForm f= swapvar (F, F.mvar(), x);
886  int sizef= 0;
887  for (CFIterator i= f; i.hasTerms(); i++, sizef++)
888  {
889    if (i.exp() == 1)
890      return 0;
891  }
892  int * expf= new int [sizef];
893  int j= 0;
894  for (CFIterator i= f; i.hasTerms(); i++, j++)
895  {
896    expf [j]= i.exp();
897  }
898
899  int indf= sizef - 1;
900  if (expf[indf] == 0)
901    indf--;
902
903  if ((d%expf [indf] != 0 && expf[indf]%d != 0) || (expf[indf] == 1))
904  {
905    delete [] expf;
906    return 0;
907  }
908
909  int result;
910  if (d%expf [indf] == 0)
911    result= expf[indf];
912  else
913    result= d;
914  for (int i= indf - 1; i >= 0; i--)
915  {
916    if (expf [i]%result != 0)
917    {
918      delete [] expf;
919      return 0;
920    }
921  }
922
923  delete [] expf;
924  return result;
925}
926
927int substituteCheck (const CFList& L)
928{
929  ASSERT (L.length() > 1, "expected a list of at least two elements");
930  if (L.length() < 2)
931    return 0;
932  CFListIterator i= L;
933  i++;
934  int result= substituteCheck (L.getFirst(), i.getItem());
935  if (result <= 1)
936    return result;
937  i++;
938  for (;i.hasItem(); i++)
939  {
940    result= recSubstituteCheck (i.getItem(), result);
941    if (result <= 1)
942      return result;
943  }
944  return result;
945}
946
947void
948subst (const CanonicalForm& F, CanonicalForm& A, const int d, const Variable& x)
949{
950  if (d <= 1)
951  {
952    A= F;
953    return;
954  }
955  if (degree (F, x) <= 0)
956  {
957    A= F;
958    return;
959  }
960  CanonicalForm C= 0;
961  CanonicalForm f= swapvar (F, x, F.mvar());
962  for (CFIterator i= f; i.hasTerms(); i++)
963    C += i.coeff()*power (f.mvar(), i.exp()/ d);
964  A= swapvar (C, x, F.mvar());
965}
966
967CanonicalForm
968reverseSubst (const CanonicalForm& F, const int d, const Variable& x)
969{
970  if (d <= 1)
971    return F;
972  if (degree (F, x) <= 0)
973    return F;
974  CanonicalForm f= swapvar (F, x, F.mvar());
975  CanonicalForm result= 0;
976  for (CFIterator i= f; i.hasTerms(); i++)
977    result += i.coeff()*power (f.mvar(), d*i.exp());
978  return swapvar (result, x, F.mvar());
979}
980
981void
982reverseSubst (CFList& L, const int d, const Variable& x)
983{
984  for (CFListIterator i= L; i.hasItem(); i++)
985    i.getItem()= reverseSubst (i.getItem(), d, x);
986}
987
Note: See TracBrowser for help on using the repository browser.