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
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 "timing.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
32TIMING_DEFINE_PRINT(fac_log_deriv_div)
33TIMING_DEFINE_PRINT(fac_log_deriv_mul)
34TIMING_DEFINE_PRINT(fac_log_deriv_pre)
35
36void append (CFList& factors1, const CFList& factors2)
37{
38  for (CFListIterator i= factors2; i.hasItem(); i++)
39  {
40    if (!i.getItem().inCoeffDomain())
41      factors1.append (i.getItem());
42  }
43  return;
44}
45
46void decompress (CFList& factors, const CFMap& N)
47{
48  for (CFListIterator i= factors; i.hasItem(); i++)
49    i.getItem()= N (i.getItem());
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());
56}
57
58void appendSwapDecompress (CFList& factors1, const CFList& factors2,
59                           const CFList& factors3, const bool swap1,
60                           const bool swap2, const CFMap& N)
61{
62  Variable x= Variable (1);
63  Variable y= Variable (2);
64  for (CFListIterator i= factors1; i.hasItem(); i++)
65  {
66    if (swap1)
67    {
68      if (!swap2)
69        i.getItem()= swapvar (i.getItem(), x, y);
70    }
71    else
72    {
73      if (swap2)
74        i.getItem()= swapvar (i.getItem(), y, x);
75    }
76    i.getItem()= N (i.getItem());
77  }
78  for (CFListIterator i= factors2; i.hasItem(); i++)
79    factors1.append (N (i.getItem()));
80  for (CFListIterator i= factors3; i.hasItem(); i++)
81    factors1.append (N (i.getItem()));
82  return;
83}
84
85void swapDecompress (CFList& factors, const bool swap, const CFMap& N)
86{
87  Variable x= Variable (1);
88  Variable y= Variable (2);
89  for (CFListIterator i= factors; i.hasItem(); i++)
90  {
91    if (swap)
92      i.getItem()= swapvar (i.getItem(), x, y);
93    i.getItem()= N (i.getItem());
94  }
95  return;
96}
97
98static inline
99bool GFInExtensionHelper (const CanonicalForm& F, const int number)
100{
101  if (F.isOne()) return false;
102  InternalCF* buf;
103  int exp;
104  bool result= false;
105  if (F.inBaseDomain())
106  {
107    buf= F.getval();
108    exp= imm2int (buf);
109    if (exp%number != 0)
110      return true;
111    else
112      return result;
113  }
114  else
115  {
116    for (CFIterator i= F; i.hasTerms(); i++)
117    {
118      result= GFInExtensionHelper (i.coeff(), number);
119      if (result == true)
120        return result;
121    }
122  }
123  return result;
124}
125
126static inline
127bool FqInExtensionHelper (const CanonicalForm& F, const CanonicalForm& gamma,
128                          const CanonicalForm& delta, CFList& source,
129                          CFList& dest)
130{
131  bool result= false;
132  if (F.inBaseDomain())
133    return result;
134  else if (F.inCoeffDomain())
135  {
136    if (!fdivides (gamma, F))
137      return true;
138    else
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    }
159  }
160  else
161  {
162    for (CFIterator i= F; i.hasTerms(); i++)
163    {
164      result= FqInExtensionHelper (i.coeff(), gamma, delta, source, dest);
165      if (result == true)
166        return result;
167    }
168  }
169  return result;
170}
171
172bool isInExtension (const CanonicalForm& F, const CanonicalForm& gamma,
173                    const int k, const CanonicalForm& delta,
174                    CFList& source, CFList& dest)
175{
176  bool result;
177  if (CFFactory::gettype() == GaloisFieldDomain)
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  }
186  else
187  {
188    result= FqInExtensionHelper (F, gamma, delta, source, dest);
189    return result;
190  }
191}
192
193CanonicalForm
194mapDown (const CanonicalForm& F, const ExtensionInfo& info, CFList& source,
195         CFList& dest)
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;
205  if (/*k==0 &&*/ beta == Variable (1))
206    return F;
207  else /*if (k==0 && beta != Variable (1))*/
208    return mapDown (F, imPrimElem, primElem, beta, source, dest);
209}
210
211void appendTestMapDown (CFList& factors, const CanonicalForm& f,
212                        const ExtensionInfo& info, CFList& source, CFList& dest)
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;
221  if (!k && beta.level() == 1)
222    degMipoBeta= 1;
223  else if (!k && beta.level() != 1)
224    degMipoBeta= degree (getMipo (beta));
225  if (k > 1)
226  {
227    if (!isInExtension (g, gamma, k, delta, source, dest))
228    {
229      g= GFMapDown (g, k);
230      factors.append (g);
231    }
232  }
233  else if (k == 1)
234  {
235    if (!isInExtension (g, gamma, k, delta, source, dest))
236      factors.append (g);
237  }
238  else if (!k && beta == Variable (1))
239  {
240    if (degree (g, alpha) < degMipoBeta)
241      factors.append (g);
242  }
243  else if (!k && beta != Variable (1))
244  {
245    if (!isInExtension (g, gamma, k, delta, source, dest))
246    {
247      g= mapDown (g, delta, gamma, alpha, source, dest);
248      factors.append (g);
249    }
250  }
251  return;
252}
253
254void
255appendMapDown (CFList& factors, const CanonicalForm& g,
256               const ExtensionInfo& info, CFList& source, CFList& dest)
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();
263  if (k > 1)
264    factors.append (GFMapDown (g, k));
265  else if (k == 1)
266    factors.append (g);
267  else if (!k && beta == Variable (1))
268    factors.append (g);
269  else if (!k && beta != Variable (1))
270    factors.append (mapDown (g, delta, gamma, alpha, source, dest));
271  return;
272}
273
274void normalize (CFList& factors)
275{
276  CanonicalForm lcinv;
277  for (CFListIterator i= factors; i.hasItem(); i++)
278  {
279    lcinv= 1/Lc (i.getItem());
280    i.getItem() *= lcinv;
281  }
282  return;
283}
284
285void normalize (CFFList& factors)
286{
287  CanonicalForm lcinv;
288  for (CFFListIterator i= factors; i.hasItem(); i++)
289  {
290    lcinv= 1/ Lc (i.getItem().factor());
291    i.getItem()= CFFactor (i.getItem().factor()*lcinv,
292                           i.getItem().exp());
293  }
294  return;
295}
296
297CFList subset (int index [], const int& s, const CFArray& elements,
298               bool& noSubset)
299{
300  int r= elements.size();
301  int i= 0;
302  CFList result;
303  noSubset= false;
304  if (index[s - 1] == 0)
305  {
306    while (i < s)
307    {
308      index[i]= i + 1;
309      result.append (elements[i]);
310      i++;
311    }
312    return result;
313  }
314  int buf;
315  int k;
316  bool found= false;
317  if (index[s - 1] == r)
318  {
319    if (index[0] == r - s + 1)
320    {
321      noSubset= true;
322      return result;
323    }
324    else {
325      while (found == false)
326      {
327        if (index[s - 2 - i] < r - i - 1)
328          found= true;
329        i++;
330      }
331      buf= index[s - i - 1];
332      k= 0;
333      while (s - i - 1 + k < s)
334      {
335        index[s - i - 1 + k]= buf + k + 1;
336        k++;
337      }
338    }
339    for (int j= 0; j < s; j++)
340      result.append (elements[index[j] - 1]);
341    return result;
342  }
343  else
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  }
350}
351
352CFArray copy (const CFList& list)
353{
354  CFArray array= CFArray (list.length());
355  int j= 0;
356  for (CFListIterator i= list; i.hasItem(); i++, j++)
357    array[j]= i.getItem();
358  return array;
359}
360
361void indexUpdate (int index [], const int& subsetSize, const int& setSize,
362                   bool& noSubset)
363{
364  noSubset= false;
365  if (subsetSize > setSize)
366  {
367    noSubset= true;
368    return;
369  }
370  int * v= new int [setSize];
371  for (int i= 0; i < setSize; i++)
372    v[i]= index[i];
373  if (subsetSize == 1)
374  {
375    v[0]= v[0] - 1;
376    if (v[0] >= setSize)
377    {
378      noSubset= true;
379      delete [] v;
380      return;
381    }
382  }
383  else
384  {
385    if (v[subsetSize - 1] - v[0] + 1 == subsetSize && v[0] > 1)
386    {
387      if (v[0] + subsetSize - 1 > setSize)
388      {
389        noSubset= true;
390        delete [] v;
391        return;
392      }
393      v[0]= v[0] - 1;
394      for (int i= 1; i < subsetSize - 1; i++)
395        v[i]= v[i - 1] + 1;
396      v[subsetSize - 1]= v[subsetSize - 2];
397    }
398    else
399    {
400      if (v[0] + subsetSize - 1 > setSize)
401      {
402        noSubset= true;
403        delete [] v;
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];
413  delete [] v;
414}
415
416int subsetDegree (const CFList& S)
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;
430  CanonicalForm quot;
431  for (CFListIterator i= factors; i.hasItem(); i++)
432  {
433    while (fdivides (i.getItem(), F, quot))
434    {
435      multi++;
436      F= quot;
437    }
438    if (multi > 0)
439      result.append (CFFactor (i.getItem(), multi));
440    multi= 0;
441  }
442  return result;
443}
444
445#ifdef HAVE_NTL
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
457  TIMING_START (fac_log_deriv_div);
458  q= newtonDiv (F, G, xToL);
459  TIMING_END_AND_PRINT (fac_log_deriv_div, "time for division in logderiv1: ");
460
461  TIMING_START (fac_log_deriv_mul);
462  logDeriv= mulMod2 (q, deriv (G, y), xToL);
463  TIMING_END_AND_PRINT (fac_log_deriv_mul, "time to multiply in logderiv1: ");
464
465  int j= degree (logDeriv, y) + 1;
466  CFArray result= CFArray (j);
467  CFIterator ii;
468  for (CFIterator i= logDeriv; i.hasTerms() && !logDeriv.isZero(); i++)
469  {
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    }
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
495  CanonicalForm bufF;
496  TIMING_START (fac_log_deriv_pre);
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);
509    CanonicalForm xToOldL2= power (x, (oldL+1)/2);
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);
515    CanonicalForm Mid;
516    if (oldL % 2 == 1)
517      Mid= mulMod2 (G1, oldQ1*x, xToLOldL);
518    else
519      Mid= mulMod2 (G1, oldQ1, xToLOldL);
520    //computation of Low might be faster using a real middle product?
521    CanonicalForm Low= mulMod2 (G0, oldQ1, xToOldL)+mulMod2 (G1, oldQ0, xToOldL);
522    Low= div (Low, power (x, oldL/2));
523    Low= mod (Low, xToLOldL);
524    Up += Mid + Low;
525    bufF= div (F, xToOldL);
526    bufF -= Up;
527  }
528  TIMING_END_AND_PRINT (fac_log_deriv_pre, "time to preprocess: ");
529
530  TIMING_START (fac_log_deriv_div);
531  if (l-oldL > 0)
532    q= newtonDiv (bufF, G, xToLOldL);
533  else
534    q= 0;
535  q *= xToOldL;
536  q += oldQ;
537  TIMING_END_AND_PRINT (fac_log_deriv_div, "time for div in logderiv2: ");
538
539  TIMING_START (fac_log_deriv_mul);
540  logDeriv= mulMod2 (q, deriv (G, y), xToL);
541  TIMING_END_AND_PRINT (fac_log_deriv_mul, "time for mul in logderiv2: ");
542
543  int j= degree (logDeriv,y) + 1;
544  CFArray result= CFArray (j);
545  CFIterator ii;
546  for (CFIterator i= logDeriv; i.hasTerms() && !logDeriv.isZero(); i++)
547  {
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    }
555  }
556  Q= q;
557  return result;
558}
559#endif
560
561void
562writeInMatrix (CFMatrix& M, const CFArray& A, const int column,
563               const int startIndex
564              )
565{
566  ASSERT (A.size () - startIndex >= 0, "wrong starting index");
567  ASSERT (A.size () - startIndex <= M.rows(), "wrong starting index");
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{
577  ASSERT (F.isUnivariate() || F.inCoeffDomain(), "univariate input expected");
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();
588      j++;
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{
600  ASSERT (F.isUnivariate() || F.inCoeffDomain(), "univariate input expected");
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{
643  ASSERT (G.isUnivariate() || G.inCoeffDomain(), "univariate input expected");
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();
668      j++;
669      if (!j.hasTerms())
670        return result;
671    }
672    else
673      result[i - k]= 0;
674  }
675  return result;
676}
677#endif
678
679int * computeBounds (const CanonicalForm& F, int& n, bool& isIrreducible)
680{
681  n= degree (F, 1);
682  int* result= new int [n];
683  int sizeOfNewtonPolygon;
684  int** newtonPolyg= newtonPolygon (F, sizeOfNewtonPolygon);
685
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
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
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
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
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.