source: git/factory/facFqBivarUtil.cc @ 5e4636

jengelh-datetimespielwiese
Last change on this file since 5e4636 was 5e4636, checked in by Martin Lee <martinlee84@…>, 10 years ago
chg: more timing info
  • Property mode set to 100644
File size: 24.1 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    for (ii= i.coeff(); ii.hasTerms(); ii++)
471      result[ii.exp()] += ii.coeff()*power (x,i.exp());
472  }
473  Q= q;
474  return result;
475}
476
477CFArray
478logarithmicDerivative (const CanonicalForm& F, const CanonicalForm& G, int l,
479                       int oldL, const CanonicalForm& oldQ, CanonicalForm& Q
480                      )
481{
482  Variable x= Variable (2);
483  Variable y= Variable (1);
484  CanonicalForm xToL= power (x, l);
485  CanonicalForm xToOldL= power (x, oldL);
486  CanonicalForm xToLOldL= power (x, l-oldL);
487  CanonicalForm q,r;
488  CanonicalForm logDeriv;
489
490  CanonicalForm bufF;
491  TIMING_START (fac_log_deriv_pre);
492  if ((oldL > 100 && l - oldL < 50) || (oldL < 100 && l - oldL < 30))
493  {
494    bufF= F;
495    CanonicalForm oldF= mulMod2 (G, oldQ, xToL);
496    bufF -= oldF;
497    bufF= div (bufF, xToOldL);
498  }
499  else
500  {
501    //middle product style computation of [G*oldQ]^{l}_{oldL}
502    CanonicalForm G3= div (G, xToOldL);
503    CanonicalForm Up= mulMod2 (G3, oldQ, xToLOldL);
504    CanonicalForm xToOldL2= power (x, (oldL+1)/2);
505    CanonicalForm G2= mod (G, xToOldL);
506    CanonicalForm G1= div (G2, xToOldL2);
507    CanonicalForm G0= mod (G2, xToOldL2);
508    CanonicalForm oldQ1= div (oldQ, xToOldL2);
509    CanonicalForm oldQ0= mod (oldQ, xToOldL2);
510    CanonicalForm Mid;
511    if (oldL % 2 == 1)
512      Mid= mulMod2 (G1, oldQ1*x, xToLOldL);
513    else
514      Mid= mulMod2 (G1, oldQ1, xToLOldL);
515    //computation of Low might be faster using a real middle product?
516    CanonicalForm Low= mulMod2 (G0, oldQ1, xToOldL)+mulMod2 (G1, oldQ0, xToOldL);
517    Low= div (Low, power (x, oldL/2));
518    Low= mod (Low, xToLOldL);
519    Up += Mid + Low;
520    bufF= div (F, xToOldL);
521    bufF -= Up;
522  }
523  TIMING_END_AND_PRINT (fac_log_deriv_pre, "time to preprocess: ");
524
525  TIMING_START (fac_log_deriv_div);
526  if (l-oldL > 0)
527    q= newtonDiv (bufF, G, xToLOldL);
528  else
529    q= 0;
530  q *= xToOldL;
531  q += oldQ;
532  TIMING_END_AND_PRINT (fac_log_deriv_div, "time for div in logderiv2: ");
533
534  TIMING_START (fac_log_deriv_mul);
535  logDeriv= mulMod2 (q, deriv (G, y), xToL);
536  TIMING_END_AND_PRINT (fac_log_deriv_mul, "time for mul in logderiv2: ");
537
538  int j= degree (logDeriv,y) + 1;
539  CFArray result= CFArray (j);
540  CFIterator ii;
541  for (CFIterator i= logDeriv; i.hasTerms() && !logDeriv.isZero(); i++)
542  {
543    for (ii= i.coeff(); ii.hasTerms(); ii++)
544      result[ii.exp()] += ii.coeff()*power (x,i.exp());
545  }
546  Q= q;
547  return result;
548}
549#endif
550
551void
552writeInMatrix (CFMatrix& M, const CFArray& A, const int column,
553               const int startIndex
554              )
555{
556  ASSERT (A.size () - startIndex >= 0, "wrong starting index");
557  ASSERT (A.size () - startIndex <= M.rows(), "wrong starting index");
558  ASSERT (column > 0 && column <= M.columns(), "wrong column");
559  if (A.size() - startIndex <= 0) return;
560  int j= 1;
561  for (int i= startIndex; i < A.size(); i++, j++)
562    M (j, column)= A [i];
563}
564
565CFArray getCoeffs (const CanonicalForm& F, const int k)
566{
567  ASSERT (F.isUnivariate() || F.inCoeffDomain(), "univariate input expected");
568  if (degree (F, 2) < k)
569    return CFArray();
570
571  CFArray result= CFArray (degree (F) - k + 1);
572  CFIterator j= F;
573  for (int i= degree (F); i >= k; i--)
574  {
575    if (j.exp() == i)
576    {
577      result [i - k]= j.coeff();
578      j++;
579      if (!j.hasTerms())
580        return result;
581    }
582    else
583      result[i - k]= 0;
584  }
585  return result;
586}
587
588CFArray getCoeffs (const CanonicalForm& F, const int k, const Variable& alpha)
589{
590  ASSERT (F.isUnivariate() || F.inCoeffDomain(), "univariate input expected");
591  if (degree (F, 2) < k)
592    return CFArray ();
593
594  int d= degree (getMipo (alpha));
595  CFArray result= CFArray ((degree (F) - k + 1)*d);
596  CFIterator j= F;
597  CanonicalForm buf;
598  CFIterator iter;
599  for (int i= degree (F); i >= k; i--)
600  {
601    if (j.exp() == i)
602    {
603      iter= j.coeff();
604      for (int l= degree (j.coeff(), alpha); l >= 0; l--)
605      {
606        if (iter.exp() == l)
607        {
608          result [(i - k)*d + l]= iter.coeff();
609          iter++;
610          if (!iter.hasTerms())
611            break;
612        }
613      }
614      j++;
615      if (!j.hasTerms())
616        return result;
617    }
618    else
619    {
620      for (int l= 0; l < d; l++)
621        result[(i - k)*d + l]= 0;
622    }
623  }
624  return result;
625}
626
627#ifdef HAVE_NTL
628CFArray
629getCoeffs (const CanonicalForm& G, const int k, const int l, const int degMipo,
630           const Variable& alpha, const CanonicalForm& evaluation,
631           const mat_zz_p& M)
632{
633  ASSERT (G.isUnivariate() || G.inCoeffDomain(), "univariate input expected");
634  CanonicalForm F= G (G.mvar() - evaluation, G.mvar());
635  if (F.isZero())
636    return CFArray ();
637
638  Variable y= Variable (2);
639  F= F (power (y, degMipo), y);
640  F= F (y, alpha);
641  zz_pX NTLF= convertFacCF2NTLzzpX (F);
642  NTLF.rep.SetLength (l*degMipo);
643  NTLF.rep= M*NTLF.rep;
644  NTLF.normalize();
645  F= convertNTLzzpX2CF (NTLF, y);
646
647  if (degree (F, 2) < k)
648    return CFArray();
649
650  CFArray result= CFArray (degree (F) - k + 1);
651
652  CFIterator j= F;
653  for (int i= degree (F); i >= k; i--)
654  {
655    if (j.exp() == i)
656    {
657      result [i - k]= j.coeff();
658      j++;
659      if (!j.hasTerms())
660        return result;
661    }
662    else
663      result[i - k]= 0;
664  }
665  return result;
666}
667#endif
668
669int * computeBounds (const CanonicalForm& F, int& n, bool& isIrreducible)
670{
671  n= degree (F, 1);
672  int* result= new int [n];
673  int sizeOfNewtonPolygon;
674  int** newtonPolyg= newtonPolygon (F, sizeOfNewtonPolygon);
675
676  isIrreducible= false;
677  if (sizeOfNewtonPolygon == 3)
678  {
679    bool check1=
680        (newtonPolyg[0][0]==0 || newtonPolyg[1][0]==0 || newtonPolyg[2][0]==0);
681    if (check1)
682    {
683      bool check2=
684        (newtonPolyg[0][1]==0 || newtonPolyg[1][1]==0 || newtonPolyg[2][0]==0);
685      if (check2)
686      {
687        int p=getCharacteristic();
688        int d=1;
689        char bufGFName='Z';
690        bool GF= (CFFactory::gettype()==GaloisFieldDomain);
691        if (GF)
692        {
693          d= getGFDegree();
694          bufGFName=gf_name;
695        }
696        setCharacteristic(0);
697        CanonicalForm tmp= gcd (newtonPolyg[0][0],newtonPolyg[0][1]);
698        tmp= gcd (tmp, newtonPolyg[1][0]);
699        tmp= gcd (tmp, newtonPolyg[1][1]);
700        tmp= gcd (tmp, newtonPolyg[2][0]);
701        tmp= gcd (tmp, newtonPolyg[2][1]);
702        isIrreducible= (tmp==1);
703        if (GF)
704          setCharacteristic (p, d, bufGFName);
705        else
706          setCharacteristic(p);
707      }
708    }
709  }
710
711  int minX, minY, maxX, maxY;
712  minX= newtonPolyg [0] [0];
713  minY= newtonPolyg [0] [1];
714  maxX= minX;
715  maxY= minY;
716  for (int i= 1; i < sizeOfNewtonPolygon; i++)
717  {
718    if (minX > newtonPolyg [i] [0])
719      minX= newtonPolyg [i] [0];
720    if (maxX < newtonPolyg [i] [0])
721      maxX= newtonPolyg [i] [0];
722    if (minY > newtonPolyg [i] [1])
723      minY= newtonPolyg [i] [1];
724    if (maxY < newtonPolyg [i] [1])
725      maxY= newtonPolyg [i] [1];
726  }
727
728  int k= maxX;
729  for (int i= 0; i < n; i++)
730  {
731    if (i + 1 > maxY || i + 1 < minY)
732    {
733      result [i]= 0;
734      continue;
735    }
736    int* point= new int [2];
737    point [0]= k;
738    point [1]= i + 1;
739    while (!isInPolygon (newtonPolyg, sizeOfNewtonPolygon, point) && k > 0)
740    {
741      k--;
742      point [0]= k;
743    }
744    result [i]= k;
745    k= maxX;
746    delete [] point;
747  }
748
749  return result;
750}
751
752int *
753computeBoundsWrtDiffMainvar (const CanonicalForm& F, int& n,
754                             bool& isIrreducible)
755{
756  n= degree (F, 2);
757  int* result= new int [n];
758  int sizeOfNewtonPolygon;
759  int** newtonPolyg= newtonPolygon (F, sizeOfNewtonPolygon);
760
761  isIrreducible= false;
762  if (sizeOfNewtonPolygon == 3)
763  {
764    bool check1=
765        (newtonPolyg[0][0]==0 || newtonPolyg[1][0]==0 || newtonPolyg[2][0]==0);
766    if (check1)
767    {
768      bool check2=
769        (newtonPolyg[0][1]==0 || newtonPolyg[1][1]==0 || newtonPolyg[2][0]==0);
770      if (check2)
771      {
772        int p=getCharacteristic();
773        int d=1;
774        char bufGFName='Z';
775        bool GF= (CFFactory::gettype()==GaloisFieldDomain);
776        if (GF)
777        {
778          d= getGFDegree();
779          bufGFName=gf_name;
780        }
781        setCharacteristic(0);
782        CanonicalForm tmp= gcd (newtonPolyg[0][0],newtonPolyg[0][1]);
783        tmp= gcd (tmp, newtonPolyg[1][0]);
784        tmp= gcd (tmp, newtonPolyg[1][1]);
785        tmp= gcd (tmp, newtonPolyg[2][0]);
786        tmp= gcd (tmp, newtonPolyg[2][1]);
787        isIrreducible= (tmp==1);
788        if (GF)
789          setCharacteristic (p, d, bufGFName);
790        else
791          setCharacteristic(p);
792      }
793    }
794  }
795
796  int minX, minY, maxX, maxY;
797  minX= newtonPolyg [0] [0];
798  minY= newtonPolyg [0] [1];
799  maxX= minX;
800  maxY= minY;
801  for (int i= 1; i < sizeOfNewtonPolygon; i++)
802  {
803    if (minX > newtonPolyg [i] [0])
804      minX= newtonPolyg [i] [0];
805    if (maxX < newtonPolyg [i] [0])
806      maxX= newtonPolyg [i] [0];
807    if (minY > newtonPolyg [i] [1])
808      minY= newtonPolyg [i] [1];
809    if (maxY < newtonPolyg [i] [1])
810      maxY= newtonPolyg [i] [1];
811  }
812
813  int k= maxY;
814  for (int i= 0; i < n; i++)
815  {
816    if (i + 1 > maxX || i + 1 < minX)
817    {
818      result [i]= 0;
819      continue;
820    }
821    int* point= new int [2];
822    point [0]= i + 1;
823    point [1]= k;
824    while (!isInPolygon (newtonPolyg, sizeOfNewtonPolygon, point) && k > 0)
825    {
826      k--;
827      point [1]= k;
828    }
829    result [i]= k;
830    k= maxY;
831    delete [] point;
832  }
833
834  return result;
835}
836
837int
838substituteCheck (const CanonicalForm& F, const Variable& x)
839{
840  if (F.inCoeffDomain())
841    return 0;
842  if (degree (F, x) < 0)
843    return 0;
844  CanonicalForm f= swapvar (F, F.mvar(), x);
845  int sizef= 0;
846  for (CFIterator i= f; i.hasTerms(); i++, sizef++)
847  {
848    if (i.exp() == 1)
849      return 0;
850  }
851  int * expf= new int [sizef];
852  int j= 0;
853  for (CFIterator i= f; i.hasTerms(); i++, j++)
854    expf [j]= i.exp();
855
856  int indf= sizef - 1;
857  if (expf[indf] == 0)
858    indf--;
859
860  int result= expf[indf];
861  for (int i= indf - 1; i >= 0; i--)
862  {
863    if (expf [i]%result != 0)
864    {
865      delete [] expf;
866      return 0;
867    }
868  }
869
870  delete [] expf;
871  return result;
872}
873
874static int
875substituteCheck (const CanonicalForm& F, const CanonicalForm& G)
876{
877  if (F.inCoeffDomain() || G.inCoeffDomain())
878    return 0;
879  Variable x= Variable (1);
880  if (degree (F, x) <= 1 || degree (G, x) <= 1)
881    return 0;
882  CanonicalForm f= swapvar (F, F.mvar(), x);
883  CanonicalForm g= swapvar (G, G.mvar(), x);
884  int sizef= 0;
885  int sizeg= 0;
886  for (CFIterator i= f; i.hasTerms(); i++, sizef++)
887  {
888    if (i.exp() == 1)
889      return 0;
890  }
891  for (CFIterator i= g; i.hasTerms(); i++, sizeg++)
892  {
893    if (i.exp() == 1)
894      return 0;
895  }
896  int * expf= new int [sizef];
897  int * expg= new int [sizeg];
898  int j= 0;
899  for (CFIterator i= f; i.hasTerms(); i++, j++)
900  {
901    expf [j]= i.exp();
902  }
903  j= 0;
904  for (CFIterator i= g; i.hasTerms(); i++, j++)
905  {
906    expg [j]= i.exp();
907  }
908
909  int indf= sizef - 1;
910  int indg= sizeg - 1;
911  if (expf[indf] == 0)
912    indf--;
913  if (expg[indg] == 0)
914    indg--;
915
916  if ((expg[indg]%expf [indf] != 0 && expf[indf]%expg[indg] != 0) ||
917      (expg[indg] == 1 && expf[indf] == 1))
918  {
919    delete [] expg;
920    delete [] expf;
921    return 0;
922  }
923
924  int result;
925  if (expg [indg]%expf [indf] == 0)
926    result= expf[indf];
927  else
928    result= expg[indg];
929  for (int i= indf - 1; i >= 0; i--)
930  {
931    if (expf [i]%result != 0)
932    {
933      delete [] expf;
934      delete [] expg;
935      return 0;
936    }
937  }
938
939  for (int i= indg - 1; i >= 0; i--)
940  {
941    if (expg [i]%result != 0)
942    {
943      delete [] expf;
944      delete [] expg;
945      return 0;
946    }
947  }
948
949  delete [] expg;
950  delete [] expf;
951  return result;
952}
953
954int recSubstituteCheck (const CanonicalForm& F, const int d)
955{
956  if (F.inCoeffDomain())
957    return 0;
958  Variable x= Variable (1);
959  if (degree (F, x) <= 1)
960    return 0;
961  CanonicalForm f= swapvar (F, F.mvar(), x);
962  int sizef= 0;
963  for (CFIterator i= f; i.hasTerms(); i++, sizef++)
964  {
965    if (i.exp() == 1)
966      return 0;
967  }
968  int * expf= new int [sizef];
969  int j= 0;
970  for (CFIterator i= f; i.hasTerms(); i++, j++)
971  {
972    expf [j]= i.exp();
973  }
974
975  int indf= sizef - 1;
976  if (expf[indf] == 0)
977    indf--;
978
979  if ((d%expf [indf] != 0 && expf[indf]%d != 0) || (expf[indf] == 1))
980  {
981    delete [] expf;
982    return 0;
983  }
984
985  int result;
986  if (d%expf [indf] == 0)
987    result= expf[indf];
988  else
989    result= d;
990  for (int i= indf - 1; i >= 0; i--)
991  {
992    if (expf [i]%result != 0)
993    {
994      delete [] expf;
995      return 0;
996    }
997  }
998
999  delete [] expf;
1000  return result;
1001}
1002
1003int substituteCheck (const CFList& L)
1004{
1005  ASSERT (L.length() > 1, "expected a list of at least two elements");
1006  if (L.length() < 2)
1007    return 0;
1008  CFListIterator i= L;
1009  i++;
1010  int result= substituteCheck (L.getFirst(), i.getItem());
1011  if (result <= 1)
1012    return result;
1013  i++;
1014  for (;i.hasItem(); i++)
1015  {
1016    result= recSubstituteCheck (i.getItem(), result);
1017    if (result <= 1)
1018      return result;
1019  }
1020  return result;
1021}
1022
1023void
1024subst (const CanonicalForm& F, CanonicalForm& A, const int d, const Variable& x)
1025{
1026  if (d <= 1)
1027  {
1028    A= F;
1029    return;
1030  }
1031  if (degree (F, x) <= 0)
1032  {
1033    A= F;
1034    return;
1035  }
1036  CanonicalForm C= 0;
1037  CanonicalForm f= swapvar (F, x, F.mvar());
1038  for (CFIterator i= f; i.hasTerms(); i++)
1039    C += i.coeff()*power (f.mvar(), i.exp()/ d);
1040  A= swapvar (C, x, F.mvar());
1041}
1042
1043CanonicalForm
1044reverseSubst (const CanonicalForm& F, const int d, const Variable& x)
1045{
1046  if (d <= 1)
1047    return F;
1048  if (degree (F, x) <= 0)
1049    return F;
1050  CanonicalForm f= swapvar (F, x, F.mvar());
1051  CanonicalForm result= 0;
1052  for (CFIterator i= f; i.hasTerms(); i++)
1053    result += i.coeff()*power (f.mvar(), d*i.exp());
1054  return swapvar (result, x, F.mvar());
1055}
1056
1057void
1058reverseSubst (CFList& L, const int d, const Variable& x)
1059{
1060  for (CFListIterator i= L; i.hasItem(); i++)
1061    i.getItem()= reverseSubst (i.getItem(), d, x);
1062}
1063
Note: See TracBrowser for help on using the repository browser.