source: git/factory/facFqBivarUtil.cc @ 9fbf02

spielwiese
Last change on this file since 9fbf02 was 9fbf02, checked in by Martin Lee <martinlee84@…>, 10 years ago
fix: bug in computing bounds for factor recombination
  • Property mode set to 100644
File size: 29.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#ifdef HAVE_CONFIG_H
14#include "config.h"
15#endif /* HAVE_CONFIG_H */
16
17#include "timing.h"
18
19#include "cf_map.h"
20#include "algext.h"
21#include "cf_map_ext.h"
22#include "templates/ftmpl_functions.h"
23#include "ExtensionInfo.h"
24#include "cf_algorithm.h"
25#include "cf_factory.h"
26#include "cf_util.h"
27#include "imm.h"
28#include "cf_iter.h"
29#include "facFqBivarUtil.h"
30#include "cfNewtonPolygon.h"
31#include "facHensel.h"
32#include "facMul.h"
33
34TIMING_DEFINE_PRINT(fac_log_deriv_div)
35TIMING_DEFINE_PRINT(fac_log_deriv_mul)
36TIMING_DEFINE_PRINT(fac_log_deriv_pre)
37
38void append (CFList& factors1, const CFList& factors2)
39{
40  for (CFListIterator i= factors2; i.hasItem(); i++)
41  {
42    if (!i.getItem().inCoeffDomain())
43      factors1.append (i.getItem());
44  }
45  return;
46}
47
48void decompress (CFList& factors, const CFMap& N)
49{
50  for (CFListIterator i= factors; i.hasItem(); i++)
51    i.getItem()= N (i.getItem());
52}
53
54void decompress (CFFList& factors, const CFMap& N)
55{
56  for (CFFListIterator i= factors; i.hasItem(); i++)
57    i.getItem()= CFFactor (N (i.getItem().factor()), i.getItem().exp());
58}
59
60void decompress (CFAFList& factors, const CFMap& N)
61{
62  for (CFAFListIterator i= factors; i.hasItem(); i++)
63    i.getItem()= CFAFactor (N (i.getItem().factor()), i.getItem().minpoly(),
64                            i.getItem().exp());
65}
66
67void appendSwapDecompress (CFList& factors1, const CFList& factors2,
68                           const CFList& factors3, const bool swap1,
69                           const bool swap2, const CFMap& N)
70{
71  Variable x= Variable (1);
72  Variable y= Variable (2);
73  for (CFListIterator i= factors1; i.hasItem(); i++)
74  {
75    if (swap1)
76    {
77      if (!swap2)
78        i.getItem()= swapvar (i.getItem(), x, y);
79    }
80    else
81    {
82      if (swap2)
83        i.getItem()= swapvar (i.getItem(), y, x);
84    }
85    i.getItem()= N (i.getItem());
86  }
87  for (CFListIterator i= factors2; i.hasItem(); i++)
88    factors1.append (N (i.getItem()));
89  for (CFListIterator i= factors3; i.hasItem(); i++)
90    factors1.append (N (i.getItem()));
91  return;
92}
93
94void swapDecompress (CFList& factors, const bool swap, const CFMap& N)
95{
96  Variable x= Variable (1);
97  Variable y= Variable (2);
98  for (CFListIterator i= factors; i.hasItem(); i++)
99  {
100    if (swap)
101      i.getItem()= swapvar (i.getItem(), x, y);
102    i.getItem()= N (i.getItem());
103  }
104  return;
105}
106
107static inline
108bool GFInExtensionHelper (const CanonicalForm& F, const int number)
109{
110  if (F.isOne()) return false;
111  InternalCF* buf;
112  int exp;
113  bool result= false;
114  if (F.inBaseDomain())
115  {
116    buf= F.getval();
117    exp= imm2int (buf);
118    if (exp%number != 0)
119      return true;
120    else
121      return result;
122  }
123  else
124  {
125    for (CFIterator i= F; i.hasTerms(); i++)
126    {
127      result= GFInExtensionHelper (i.coeff(), number);
128      if (result == true)
129        return result;
130    }
131  }
132  return result;
133}
134
135static inline
136bool FqInExtensionHelper (const CanonicalForm& F, const CanonicalForm& gamma,
137                          const CanonicalForm& delta, CFList& source,
138                          CFList& dest)
139{
140  bool result= false;
141  if (F.inBaseDomain())
142    return result;
143  else if (F.inCoeffDomain())
144  {
145    if (!fdivides (gamma, F))
146      return true;
147    else
148    {
149      int pos= findItem (source, F);
150      if (pos > 0)
151        return false;
152      Variable a;
153      hasFirstAlgVar (F, a);
154      int bound= ipower (getCharacteristic(), degree (getMipo (a)));
155      CanonicalForm buf= 1;
156      for (int i= 1; i < bound; i++)
157      {
158        buf *= gamma;
159        if (buf == F)
160        {
161          source.append (buf);
162          dest.append (power (delta, i));
163          return false;
164        }
165      }
166      return true;
167    }
168  }
169  else
170  {
171    for (CFIterator i= F; i.hasTerms(); i++)
172    {
173      result= FqInExtensionHelper (i.coeff(), gamma, delta, source, dest);
174      if (result == true)
175        return result;
176    }
177  }
178  return result;
179}
180
181bool isInExtension (const CanonicalForm& F, const CanonicalForm& gamma,
182                    const int k, const CanonicalForm& delta,
183                    CFList& source, CFList& dest)
184{
185  bool result;
186  if (CFFactory::gettype() == GaloisFieldDomain)
187  {
188    int p= getCharacteristic();
189    int orderFieldExtension= ipower (p, getGFDegree()) - 1;
190    int order= ipower (p, k) - 1;
191    int number= orderFieldExtension/order;
192    result= GFInExtensionHelper (F, number);
193    return result;
194  }
195  else
196  {
197    result= FqInExtensionHelper (F, gamma, delta, source, dest);
198    return result;
199  }
200}
201
202CanonicalForm
203mapDown (const CanonicalForm& F, const ExtensionInfo& info, CFList& source,
204         CFList& dest)
205{
206  int k= info.getGFDegree();
207  Variable beta= info.getAlpha();
208  CanonicalForm primElem= info.getGamma();
209  CanonicalForm imPrimElem= info.getDelta();
210  if (k > 1)
211    return GFMapDown (F, k);
212  else if (k == 1)
213    return F;
214  if (/*k==0 &&*/ beta == Variable (1))
215    return F;
216  else /*if (k==0 && beta != Variable (1))*/
217    return mapDown (F, imPrimElem, primElem, beta, source, dest);
218}
219
220void appendTestMapDown (CFList& factors, const CanonicalForm& f,
221                        const ExtensionInfo& info, CFList& source, CFList& dest)
222{
223  int k= info.getGFDegree();
224  Variable beta= info.getBeta();
225  Variable alpha= info.getAlpha();
226  CanonicalForm delta= info.getDelta();
227  CanonicalForm gamma= info.getGamma();
228  CanonicalForm g= f;
229  int degMipoBeta;
230  if (!k && beta.level() == 1)
231    degMipoBeta= 1;
232  else if (!k && beta.level() != 1)
233    degMipoBeta= degree (getMipo (beta));
234  if (k > 1)
235  {
236    if (!isInExtension (g, gamma, k, delta, source, dest))
237    {
238      g= GFMapDown (g, k);
239      factors.append (g);
240    }
241  }
242  else if (k == 1)
243  {
244    if (!isInExtension (g, gamma, k, delta, source, dest))
245      factors.append (g);
246  }
247  else if (!k && beta == Variable (1))
248  {
249    if (degree (g, alpha) < degMipoBeta)
250      factors.append (g);
251  }
252  else if (!k && beta != Variable (1))
253  {
254    if (!isInExtension (g, gamma, k, delta, source, dest))
255    {
256      g= mapDown (g, delta, gamma, alpha, source, dest);
257      factors.append (g);
258    }
259  }
260  return;
261}
262
263void
264appendMapDown (CFList& factors, const CanonicalForm& g,
265               const ExtensionInfo& info, CFList& source, CFList& dest)
266{
267  int k= info.getGFDegree();
268  Variable beta= info.getBeta();
269  Variable alpha= info.getAlpha();
270  CanonicalForm gamma= info.getGamma();
271  CanonicalForm delta= info.getDelta();
272  if (k > 1)
273    factors.append (GFMapDown (g, k));
274  else if (k == 1)
275    factors.append (g);
276  else if (!k && beta == Variable (1))
277    factors.append (g);
278  else if (!k && beta != Variable (1))
279    factors.append (mapDown (g, delta, gamma, alpha, source, dest));
280  return;
281}
282
283void normalize (CFList& factors)
284{
285  CanonicalForm lcinv;
286  for (CFListIterator i= factors; i.hasItem(); i++)
287  {
288    lcinv= 1/Lc (i.getItem());
289    i.getItem() *= lcinv;
290  }
291  return;
292}
293
294void normalize (CFFList& factors)
295{
296  CanonicalForm lcinv;
297  for (CFFListIterator i= factors; i.hasItem(); i++)
298  {
299    lcinv= 1/ Lc (i.getItem().factor());
300    i.getItem()= CFFactor (i.getItem().factor()*lcinv,
301                           i.getItem().exp());
302  }
303  return;
304}
305
306CFList subset (int index [], const int& s, const CFArray& elements,
307               bool& noSubset)
308{
309  int r= elements.size();
310  int i= 0;
311  CFList result;
312  noSubset= false;
313  if (index[s - 1] == 0)
314  {
315    while (i < s)
316    {
317      index[i]= i + 1;
318      result.append (elements[i]);
319      i++;
320    }
321    return result;
322  }
323  int buf;
324  int k;
325  bool found= false;
326  if (index[s - 1] == r)
327  {
328    if (index[0] == r - s + 1)
329    {
330      noSubset= true;
331      return result;
332    }
333    else {
334      while (found == false)
335      {
336        if (index[s - 2 - i] < r - i - 1)
337          found= true;
338        i++;
339      }
340      buf= index[s - i - 1];
341      k= 0;
342      while (s - i - 1 + k < s)
343      {
344        index[s - i - 1 + k]= buf + k + 1;
345        k++;
346      }
347    }
348    for (int j= 0; j < s; j++)
349      result.append (elements[index[j] - 1]);
350    return result;
351  }
352  else
353  {
354    index[s - 1] += 1;
355    for (int j= 0; j < s; j++)
356      result.append (elements[index[j] - 1]);
357    return result;
358  }
359}
360
361CFArray copy (const CFList& list)
362{
363  CFArray array= CFArray (list.length());
364  int j= 0;
365  for (CFListIterator i= list; i.hasItem(); i++, j++)
366    array[j]= i.getItem();
367  return array;
368}
369
370void indexUpdate (int index [], const int& subsetSize, const int& setSize,
371                   bool& noSubset)
372{
373  noSubset= false;
374  if (subsetSize > setSize)
375  {
376    noSubset= true;
377    return;
378  }
379  int * v= new int [setSize];
380  for (int i= 0; i < setSize; i++)
381    v[i]= index[i];
382  if (subsetSize == 1)
383  {
384    v[0]= v[0] - 1;
385    if (v[0] >= setSize)
386    {
387      noSubset= true;
388      delete [] v;
389      return;
390    }
391  }
392  else
393  {
394    if (v[subsetSize - 1] - v[0] + 1 == subsetSize && v[0] > 1)
395    {
396      if (v[0] + subsetSize - 1 > setSize)
397      {
398        noSubset= true;
399        delete [] v;
400        return;
401      }
402      v[0]= v[0] - 1;
403      for (int i= 1; i < subsetSize - 1; i++)
404        v[i]= v[i - 1] + 1;
405      v[subsetSize - 1]= v[subsetSize - 2];
406    }
407    else
408    {
409      if (v[0] + subsetSize - 1 > setSize)
410      {
411        noSubset= true;
412        delete [] v;
413        return;
414      }
415      for (int i= 1; i < subsetSize - 1; i++)
416        v[i]= v[i - 1] + 1;
417      v[subsetSize - 1]= v[subsetSize - 2];
418    }
419  }
420  for (int i= 0; i < setSize; i++)
421    index[i]= v[i];
422  delete [] v;
423}
424
425int subsetDegree (const CFList& S)
426{
427  int result= 0;
428  for (CFListIterator i= S; i.hasItem(); i++)
429    result += degree (i.getItem(), Variable (1));
430  return result;
431}
432
433CFFList multiplicity (CanonicalForm& F, const CFList& factors)
434{
435  if (F.inCoeffDomain())
436    return CFFList (CFFactor (F, 1));
437  CFFList result;
438  int multi= 0;
439  CanonicalForm quot;
440  for (CFListIterator i= factors; i.hasItem(); i++)
441  {
442    while (fdivides (i.getItem(), F, quot))
443    {
444      multi++;
445      F= quot;
446    }
447    if (multi > 0)
448      result.append (CFFactor (i.getItem(), multi));
449    multi= 0;
450  }
451  return result;
452}
453
454#ifdef HAVE_NTL
455CFArray
456logarithmicDerivative (const CanonicalForm& F, const CanonicalForm& G, int l,
457                       CanonicalForm& Q
458                      )
459{
460  Variable x= Variable (2);
461  Variable y= Variable (1);
462  CanonicalForm xToL= power (x, l);
463  CanonicalForm q,r;
464  CanonicalForm logDeriv;
465
466  TIMING_START (fac_log_deriv_div);
467  q= newtonDiv (F, G, xToL);
468  TIMING_END_AND_PRINT (fac_log_deriv_div, "time for division in logderiv1: ");
469
470  TIMING_START (fac_log_deriv_mul);
471  logDeriv= mulMod2 (q, deriv (G, y), xToL);
472  TIMING_END_AND_PRINT (fac_log_deriv_mul, "time to multiply in logderiv1: ");
473
474  if (degree (logDeriv, x) == 0)
475  {
476    Q= q;
477    return CFArray();
478  }
479
480  int j= degree (logDeriv, y) + 1;
481  CFArray result= CFArray (j);
482  CFIterator ii;
483  for (CFIterator i= logDeriv; i.hasTerms() && !logDeriv.isZero(); i++)
484  {
485    if (i.coeff().inCoeffDomain())
486      result[0] += i.coeff()*power (x,i.exp());
487    else
488    {
489      for (ii= i.coeff(); ii.hasTerms(); ii++)
490        result[ii.exp()] += ii.coeff()*power (x,i.exp());
491    }
492  }
493  Q= q;
494  return result;
495}
496
497CFArray
498logarithmicDerivative (const CanonicalForm& F, const CanonicalForm& G, int l,
499                       int oldL, const CanonicalForm& oldQ, CanonicalForm& Q
500                      )
501{
502  Variable x= Variable (2);
503  Variable y= Variable (1);
504  CanonicalForm xToL= power (x, l);
505  CanonicalForm xToOldL= power (x, oldL);
506  CanonicalForm xToLOldL= power (x, l-oldL);
507  CanonicalForm q,r;
508  CanonicalForm logDeriv;
509
510  CanonicalForm bufF;
511  TIMING_START (fac_log_deriv_pre);
512  if ((oldL > 100 && l - oldL < 50) || (oldL < 100 && l - oldL < 30))
513  {
514    bufF= F;
515    CanonicalForm oldF= mulMod2 (G, oldQ, xToL);
516    bufF -= oldF;
517    bufF= div (bufF, xToOldL);
518  }
519  else
520  {
521    //middle product style computation of [G*oldQ]^{l}_{oldL}
522    CanonicalForm G3= div (G, xToOldL);
523    CanonicalForm Up= mulMod2 (G3, oldQ, xToLOldL);
524    CanonicalForm xToOldL2= power (x, (oldL+1)/2);
525    CanonicalForm G2= mod (G, xToOldL);
526    CanonicalForm G1= div (G2, xToOldL2);
527    CanonicalForm G0= mod (G2, xToOldL2);
528    CanonicalForm oldQ1= div (oldQ, xToOldL2);
529    CanonicalForm oldQ0= mod (oldQ, xToOldL2);
530    CanonicalForm Mid;
531    if (oldL % 2 == 1)
532      Mid= mulMod2 (G1, oldQ1*x, xToLOldL);
533    else
534      Mid= mulMod2 (G1, oldQ1, xToLOldL);
535    //computation of Low might be faster using a real middle product?
536    CanonicalForm Low= mulMod2 (G0, oldQ1, xToOldL)+mulMod2 (G1, oldQ0, xToOldL);
537    Low= div (Low, power (x, oldL/2));
538    Low= mod (Low, xToLOldL);
539    Up += Mid + Low;
540    bufF= div (F, xToOldL);
541    bufF -= Up;
542  }
543  TIMING_END_AND_PRINT (fac_log_deriv_pre, "time to preprocess: ");
544
545  TIMING_START (fac_log_deriv_div);
546  if (l-oldL > 0)
547    q= newtonDiv (bufF, G, xToLOldL);
548  else
549    q= 0;
550  q *= xToOldL;
551  q += oldQ;
552  TIMING_END_AND_PRINT (fac_log_deriv_div, "time for div in logderiv2: ");
553
554  TIMING_START (fac_log_deriv_mul);
555  logDeriv= mulMod2 (q, deriv (G, y), xToL);
556  TIMING_END_AND_PRINT (fac_log_deriv_mul, "time for mul in logderiv2: ");
557
558  if (degree (logDeriv, x) == 0)
559  {
560    Q= q;
561    return CFArray();
562  }
563
564  int j= degree (logDeriv,y) + 1;
565  CFArray result= CFArray (j);
566  CFIterator ii;
567  for (CFIterator i= logDeriv; i.hasTerms() && !logDeriv.isZero(); i++)
568  {
569    if (i.coeff().inCoeffDomain())
570      result[0] += i.coeff()*power (x,i.exp());
571    else
572    {
573      for (ii= i.coeff(); ii.hasTerms(); ii++)
574        result[ii.exp()] += ii.coeff()*power (x,i.exp());
575    }
576  }
577  Q= q;
578  return result;
579}
580#endif
581
582void
583writeInMatrix (CFMatrix& M, const CFArray& A, const int column,
584               const int startIndex
585              )
586{
587  ASSERT (A.size () - startIndex >= 0, "wrong starting index");
588  ASSERT (A.size () - startIndex <= M.rows(), "wrong starting index");
589  ASSERT (column > 0 && column <= M.columns(), "wrong column");
590  if (A.size() - startIndex <= 0) return;
591  int j= 1;
592  for (int i= startIndex; i < A.size(); i++, j++)
593    M (j, column)= A [i];
594}
595
596CFArray getCoeffs (const CanonicalForm& F, const int k)
597{
598  ASSERT (F.isUnivariate() || F.inCoeffDomain(), "univariate input expected");
599  if (degree (F, 2) < k)
600    return CFArray();
601
602  CFArray result= CFArray (degree (F) - k + 1);
603  CFIterator j= F;
604  for (int i= degree (F); i >= k; i--)
605  {
606    if (j.exp() == i)
607    {
608      result [i - k]= j.coeff();
609      j++;
610      if (!j.hasTerms())
611        return result;
612    }
613    else
614      result[i - k]= 0;
615  }
616  return result;
617}
618
619CFArray getCoeffs (const CanonicalForm& F, const int k, const Variable& alpha)
620{
621  ASSERT (F.isUnivariate() || F.inCoeffDomain(), "univariate input expected");
622  if (degree (F, 2) < k)
623    return CFArray ();
624
625  int d= degree (getMipo (alpha));
626  CFArray result= CFArray ((degree (F) - k + 1)*d);
627  CFIterator j= F;
628  CanonicalForm buf;
629  CFIterator iter;
630  for (int i= degree (F); i >= k; i--)
631  {
632    if (j.exp() == i)
633    {
634      iter= j.coeff();
635      for (int l= degree (j.coeff(), alpha); l >= 0; l--)
636      {
637        if (iter.exp() == l)
638        {
639          result [(i - k)*d + l]= iter.coeff();
640          iter++;
641          if (!iter.hasTerms())
642            break;
643        }
644      }
645      j++;
646      if (!j.hasTerms())
647        return result;
648    }
649    else
650    {
651      for (int l= 0; l < d; l++)
652        result[(i - k)*d + l]= 0;
653    }
654  }
655  return result;
656}
657
658#ifdef HAVE_NTL
659CFArray
660getCoeffs (const CanonicalForm& G, const int k, const int l, const int degMipo,
661           const Variable& alpha, const CanonicalForm& evaluation,
662           const mat_zz_p& M)
663{
664  ASSERT (G.isUnivariate() || G.inCoeffDomain(), "univariate input expected");
665  CanonicalForm F= G (G.mvar() - evaluation, G.mvar());
666  if (F.isZero())
667    return CFArray ();
668
669  Variable y= Variable (2);
670  F= F (power (y, degMipo), y);
671  F= F (y, alpha);
672  zz_pX NTLF= convertFacCF2NTLzzpX (F);
673  NTLF.rep.SetLength (l*degMipo);
674  NTLF.rep= M*NTLF.rep;
675  NTLF.normalize();
676  F= convertNTLzzpX2CF (NTLF, y);
677
678  if (degree (F, 2) < k)
679    return CFArray();
680
681  CFArray result= CFArray (degree (F) - k + 1);
682
683  CFIterator j= F;
684  for (int i= degree (F); i >= k; i--)
685  {
686    if (j.exp() == i)
687    {
688      result [i - k]= j.coeff();
689      j++;
690      if (!j.hasTerms())
691        return result;
692    }
693    else
694      result[i - k]= 0;
695  }
696  return result;
697}
698#endif
699
700int * computeBounds (const CanonicalForm& F, int& n, bool& isIrreducible)
701{
702  n= degree (F, 1);
703  int* result= new int [n];
704  int sizeOfNewtonPolygon;
705  int** newtonPolyg= newtonPolygon (F, sizeOfNewtonPolygon);
706
707  isIrreducible= false;
708  if (sizeOfNewtonPolygon == 3)
709  {
710    bool check1=
711        (newtonPolyg[0][0]==0 || newtonPolyg[1][0]==0 || newtonPolyg[2][0]==0);
712    if (check1)
713    {
714      bool check2=
715        (newtonPolyg[0][1]==0 || newtonPolyg[1][1]==0 || newtonPolyg[2][0]==0);
716      if (check2)
717      {
718        int p=getCharacteristic();
719        int d=1;
720        char bufGFName='Z';
721        bool GF= (CFFactory::gettype()==GaloisFieldDomain);
722        if (GF)
723        {
724          d= getGFDegree();
725          bufGFName=gf_name;
726        }
727        setCharacteristic(0);
728        CanonicalForm tmp= gcd (newtonPolyg[0][0],newtonPolyg[0][1]);  // maybe it's better to use plain intgcd
729        tmp= gcd (tmp, newtonPolyg[1][0]);
730        tmp= gcd (tmp, newtonPolyg[1][1]);
731        tmp= gcd (tmp, newtonPolyg[2][0]);
732        tmp= gcd (tmp, newtonPolyg[2][1]);
733        isIrreducible= (tmp==1);
734        if (GF)
735          setCharacteristic (p, d, bufGFName);
736        else
737          setCharacteristic(p);
738      }
739    }
740  }
741
742  int minX, minY, maxX, maxY;
743  minX= newtonPolyg [0] [0];
744  minY= newtonPolyg [0] [1];
745  maxX= minX;
746  maxY= minY;
747  int indZero= 0;
748  for (int i= 1; i < sizeOfNewtonPolygon; i++)
749  {
750    if (newtonPolyg[i][1] == 0)
751    {
752      if (newtonPolyg[indZero][1] == 0)
753      {
754        if (newtonPolyg[indZero][0] < newtonPolyg[i][0])
755          indZero= i;
756      }
757      else
758        indZero= i;
759    }
760    if (minX > newtonPolyg [i] [0])
761      minX= newtonPolyg [i] [0];
762    if (maxX < newtonPolyg [i] [0])
763      maxX= newtonPolyg [i] [0];
764    if (minY > newtonPolyg [i] [1])
765      minY= newtonPolyg [i] [1];
766    if (maxY < newtonPolyg [i] [1])
767      maxY= newtonPolyg [i] [1];
768  }
769
770  int slopeNum, slopeDen, constTerm;
771  bool negativeSlope=false;
772  if (indZero != sizeOfNewtonPolygon - 1)
773  {
774    slopeNum= newtonPolyg[indZero+1][0]-newtonPolyg[indZero][0];
775    slopeDen= newtonPolyg[indZero+1][1];
776    constTerm= newtonPolyg[indZero][0];
777  }
778  else
779  {
780    slopeNum= newtonPolyg[0][0]-newtonPolyg[indZero][0];
781    slopeDen= newtonPolyg[0][1];
782    constTerm= newtonPolyg[indZero][0];
783  }
784  if (slopeNum < 0)
785  {
786    slopeNum= -slopeNum;
787    negativeSlope= true;
788  }
789  int k= 0;
790  int* point= new int [2];
791  for (int i= 0; i < n; i++)
792  {
793    if (((indZero+1) < sizeOfNewtonPolygon && (i+1) > newtonPolyg[indZero+1][1])
794        || ((indZero+1) >= sizeOfNewtonPolygon && (i+1) > newtonPolyg[0][1]))
795    {
796      if (indZero + 1 != sizeOfNewtonPolygon)
797        indZero++;
798      else
799        indZero= 0;
800      if (indZero != sizeOfNewtonPolygon - 1)
801      {
802        slopeNum= newtonPolyg[indZero+1][0]-newtonPolyg[indZero][0];
803        slopeDen= newtonPolyg[indZero+1][1]-newtonPolyg[indZero][1];
804        constTerm= newtonPolyg[indZero][0];
805      }
806      else
807      {
808        slopeNum= newtonPolyg[0][0]-newtonPolyg[indZero][0];
809        slopeDen= newtonPolyg[0][1]-newtonPolyg[indZero][1];
810        constTerm= newtonPolyg[indZero][0];
811      }
812      if (slopeNum < 0)
813      {
814        negativeSlope= true;
815        slopeNum= - slopeNum;
816        k= (int) -(((long) slopeNum*((i+1)-newtonPolyg[indZero][1])+slopeDen-1)/
817                   slopeDen) + constTerm;
818      }
819      else
820        k= (int) (((long) slopeNum*((i+1)-newtonPolyg[indZero][1])) / slopeDen)
821                  + constTerm;
822    }
823    else
824    {
825      if (negativeSlope)
826        k= (int) -(((long) slopeNum*((i+1)-newtonPolyg[indZero][1])+slopeDen-1)/
827                   slopeDen) + constTerm;
828      else
829        k= (int) ((long) slopeNum*((i+1)-newtonPolyg[indZero][1])) / slopeDen
830                  + constTerm;
831    }
832    if (i + 1 > maxY || i + 1 < minY)
833    {
834      result [i]= 0;
835      continue;
836    }
837    point [0]= k;
838    point [1]= i + 1;
839    if (!isInPolygon (newtonPolyg, sizeOfNewtonPolygon, point) && k > 0)
840      k= 0;
841    result [i]= k;
842  }
843
844  delete [] point;
845
846  for (int i= 0; i < sizeOfNewtonPolygon; i++)
847    delete [] newtonPolyg[i];
848  delete [] newtonPolyg;
849
850  return result;
851}
852
853int *
854computeBoundsWrtDiffMainvar (const CanonicalForm& F, int& n,
855                             bool& isIrreducible)
856{
857  n= degree (F, 2);
858  int* result= new int [n];
859  int sizeOfNewtonPolygon;
860  int** newtonPolyg= newtonPolygon (F, sizeOfNewtonPolygon);
861
862  isIrreducible= false;
863  if (sizeOfNewtonPolygon == 3)
864  {
865    bool check1=
866        (newtonPolyg[0][0]==0 || newtonPolyg[1][0]==0 || newtonPolyg[2][0]==0);
867    if (check1)
868    {
869      bool check2=
870        (newtonPolyg[0][1]==0 || newtonPolyg[1][1]==0 || newtonPolyg[2][0]==0);
871      if (check2)
872      {
873        int p=getCharacteristic();
874        int d=1;
875        char bufGFName='Z';
876        bool GF= (CFFactory::gettype()==GaloisFieldDomain);
877        if (GF)
878        {
879          d= getGFDegree();
880          bufGFName=gf_name;
881        }
882        setCharacteristic(0);
883        CanonicalForm tmp= gcd (newtonPolyg[0][0],newtonPolyg[0][1]);
884        tmp= gcd (tmp, newtonPolyg[1][0]);
885        tmp= gcd (tmp, newtonPolyg[1][1]);
886        tmp= gcd (tmp, newtonPolyg[2][0]);
887        tmp= gcd (tmp, newtonPolyg[2][1]);
888        isIrreducible= (tmp==1);
889        if (GF)
890          setCharacteristic (p, d, bufGFName);
891        else
892          setCharacteristic(p);
893      }
894    }
895  }
896
897  int swap;
898  for (int i= 0; i < sizeOfNewtonPolygon; i++)
899  {
900    swap= newtonPolyg[i][1];
901    newtonPolyg[i][1]=newtonPolyg[i][0];
902    newtonPolyg[i][0]= swap;
903  }
904
905  sizeOfNewtonPolygon= polygon(newtonPolyg, sizeOfNewtonPolygon);
906
907  int minX, minY, maxX, maxY;
908  minX= newtonPolyg [0] [0];
909  minY= newtonPolyg [0] [1];
910  maxX= minX;
911  maxY= minY;
912  int indZero= 0;
913  for (int i= 1; i < sizeOfNewtonPolygon; i++)
914  {
915    if (newtonPolyg[i][1] == 0)
916    {
917      if (newtonPolyg[indZero][1] == 0)
918      {
919        if (newtonPolyg[indZero][0] < newtonPolyg[i][0])
920          indZero= i;
921      }
922      else
923        indZero= i;
924    }
925    if (minX > newtonPolyg [i] [0])
926      minX= newtonPolyg [i] [0];
927    if (maxX < newtonPolyg [i] [0])
928      maxX= newtonPolyg [i] [0];
929    if (minY > newtonPolyg [i] [1])
930      minY= newtonPolyg [i] [1];
931    if (maxY < newtonPolyg [i] [1])
932      maxY= newtonPolyg [i] [1];
933  }
934
935  int slopeNum, slopeDen, constTerm;
936  bool negativeSlope=false;
937  if (indZero != sizeOfNewtonPolygon - 1)
938  {
939    slopeNum= newtonPolyg[indZero+1][0]-newtonPolyg[indZero][0];
940    slopeDen= newtonPolyg[indZero+1][1];
941    constTerm= newtonPolyg[indZero][0];
942  }
943  else
944  {
945    slopeNum= newtonPolyg[0][0]-newtonPolyg[indZero][0];
946    slopeDen= newtonPolyg[0][1];
947    constTerm= newtonPolyg[indZero][0];
948  }
949  if (slopeNum < 0)
950  {
951    slopeNum= -slopeNum;
952    negativeSlope= true;
953  }
954  int k= 0;
955
956  int* point= new int [2];
957  for (int i= 0; i < n; i++)
958  {
959    if (((indZero+1) < sizeOfNewtonPolygon && (i+1) > newtonPolyg[indZero+1][1])
960        || ((indZero+1) >= sizeOfNewtonPolygon && (i+1) > newtonPolyg[0][1]))
961    {
962      if (indZero + 1 != sizeOfNewtonPolygon)
963        indZero++;
964      else
965        indZero= 0;
966      if (indZero != sizeOfNewtonPolygon - 1)
967      {
968        slopeNum= newtonPolyg[indZero+1][0]-newtonPolyg[indZero][0];
969        slopeDen= newtonPolyg[indZero+1][1]-newtonPolyg[indZero][1];
970        constTerm= newtonPolyg[indZero][0];
971      }
972      else
973      {
974        slopeNum= newtonPolyg[0][0]-newtonPolyg[indZero][0];
975        slopeDen= newtonPolyg[0][1]-newtonPolyg[indZero][1];
976        constTerm= newtonPolyg[indZero][0];
977      }
978      if (slopeNum < 0)
979      {
980        negativeSlope= true;
981        slopeNum= - slopeNum;
982        k= (int) -(((long) slopeNum*((i+1)-newtonPolyg[indZero][1])+slopeDen-1)/
983                   slopeDen) + constTerm;
984      }
985      else
986        k= (int) (((long) slopeNum*((i+1)-newtonPolyg[indZero][1])) / slopeDen)
987                  + constTerm;
988    }
989    else
990    {
991      if (negativeSlope)
992        k= (int) -(((long) slopeNum*((i+1)-newtonPolyg[indZero][1])+slopeDen-1)/
993                   slopeDen) + constTerm;
994      else
995        k= (int) ((long) slopeNum*((i+1)-newtonPolyg[indZero][1])) / slopeDen
996                  + constTerm;
997    }
998    if (i + 1 > maxY || i + 1 < minY)
999    {
1000      result [i]= 0;
1001      continue;
1002    }
1003
1004    point [0]= k;
1005    point [1]= i + 1;
1006    if (!isInPolygon (newtonPolyg, sizeOfNewtonPolygon, point) && k > 0)
1007      k= 0;
1008    result [i]= k;
1009  }
1010
1011  delete [] point;
1012
1013  for (int i= 0; i < sizeOfNewtonPolygon; i++)
1014    delete [] newtonPolyg[i];
1015  delete [] newtonPolyg;
1016
1017  return result;
1018}
1019
1020int
1021substituteCheck (const CanonicalForm& F, const Variable& x)
1022{
1023  if (F.inCoeffDomain())
1024    return 0;
1025  if (degree (F, x) < 0)
1026    return 0;
1027  CanonicalForm f= swapvar (F, F.mvar(), x);
1028  int sizef= 0;
1029  for (CFIterator i= f; i.hasTerms(); i++, sizef++)
1030  {
1031    if (i.exp() == 1)
1032      return 0;
1033  }
1034  int * expf= new int [sizef];
1035  int j= 0;
1036  for (CFIterator i= f; i.hasTerms(); i++, j++)
1037    expf [j]= i.exp();
1038
1039  int indf= sizef - 1;
1040  if (expf[indf] == 0)
1041    indf--;
1042
1043  int result= expf[indf];
1044  for (int i= indf - 1; i >= 0; i--)
1045  {
1046    if (expf [i]%result != 0)
1047    {
1048      delete [] expf;
1049      return 0;
1050    }
1051  }
1052
1053  delete [] expf;
1054  return result;
1055}
1056
1057static int
1058substituteCheck (const CanonicalForm& F, const CanonicalForm& G)
1059{
1060  if (F.inCoeffDomain() || G.inCoeffDomain())
1061    return 0;
1062  Variable x= Variable (1);
1063  if (degree (F, x) <= 1 || degree (G, x) <= 1)
1064    return 0;
1065  CanonicalForm f= swapvar (F, F.mvar(), x);
1066  CanonicalForm g= swapvar (G, G.mvar(), x);
1067  int sizef= 0;
1068  int sizeg= 0;
1069  for (CFIterator i= f; i.hasTerms(); i++, sizef++)
1070  {
1071    if (i.exp() == 1)
1072      return 0;
1073  }
1074  for (CFIterator i= g; i.hasTerms(); i++, sizeg++)
1075  {
1076    if (i.exp() == 1)
1077      return 0;
1078  }
1079  int * expf= new int [sizef];
1080  int * expg= new int [sizeg];
1081  int j= 0;
1082  for (CFIterator i= f; i.hasTerms(); i++, j++)
1083  {
1084    expf [j]= i.exp();
1085  }
1086  j= 0;
1087  for (CFIterator i= g; i.hasTerms(); i++, j++)
1088  {
1089    expg [j]= i.exp();
1090  }
1091
1092  int indf= sizef - 1;
1093  int indg= sizeg - 1;
1094  if (expf[indf] == 0)
1095    indf--;
1096  if (expg[indg] == 0)
1097    indg--;
1098
1099  if ((expg[indg]%expf [indf] != 0 && expf[indf]%expg[indg] != 0) ||
1100      (expg[indg] == 1 && expf[indf] == 1))
1101  {
1102    delete [] expg;
1103    delete [] expf;
1104    return 0;
1105  }
1106
1107  int result;
1108  if (expg [indg]%expf [indf] == 0)
1109    result= expf[indf];
1110  else
1111    result= expg[indg];
1112  for (int i= indf - 1; i >= 0; i--)
1113  {
1114    if (expf [i]%result != 0)
1115    {
1116      delete [] expf;
1117      delete [] expg;
1118      return 0;
1119    }
1120  }
1121
1122  for (int i= indg - 1; i >= 0; i--)
1123  {
1124    if (expg [i]%result != 0)
1125    {
1126      delete [] expf;
1127      delete [] expg;
1128      return 0;
1129    }
1130  }
1131
1132  delete [] expg;
1133  delete [] expf;
1134  return result;
1135}
1136
1137int recSubstituteCheck (const CanonicalForm& F, const int d)
1138{
1139  if (F.inCoeffDomain())
1140    return 0;
1141  Variable x= Variable (1);
1142  if (degree (F, x) <= 1)
1143    return 0;
1144  CanonicalForm f= swapvar (F, F.mvar(), x);
1145  int sizef= 0;
1146  for (CFIterator i= f; i.hasTerms(); i++, sizef++)
1147  {
1148    if (i.exp() == 1)
1149      return 0;
1150  }
1151  int * expf= new int [sizef];
1152  int j= 0;
1153  for (CFIterator i= f; i.hasTerms(); i++, j++)
1154  {
1155    expf [j]= i.exp();
1156  }
1157
1158  int indf= sizef - 1;
1159  if (expf[indf] == 0)
1160    indf--;
1161
1162  if ((d%expf [indf] != 0 && expf[indf]%d != 0) || (expf[indf] == 1))
1163  {
1164    delete [] expf;
1165    return 0;
1166  }
1167
1168  int result;
1169  if (d%expf [indf] == 0)
1170    result= expf[indf];
1171  else
1172    result= d;
1173  for (int i= indf - 1; i >= 0; i--)
1174  {
1175    if (expf [i]%result != 0)
1176    {
1177      delete [] expf;
1178      return 0;
1179    }
1180  }
1181
1182  delete [] expf;
1183  return result;
1184}
1185
1186int substituteCheck (const CFList& L)
1187{
1188  ASSERT (L.length() > 1, "expected a list of at least two elements");
1189  if (L.length() < 2)
1190    return 0;
1191  CFListIterator i= L;
1192  i++;
1193  int result= substituteCheck (L.getFirst(), i.getItem());
1194  if (result <= 1)
1195    return result;
1196  i++;
1197  for (;i.hasItem(); i++)
1198  {
1199    result= recSubstituteCheck (i.getItem(), result);
1200    if (result <= 1)
1201      return result;
1202  }
1203  return result;
1204}
1205
1206void
1207subst (const CanonicalForm& F, CanonicalForm& A, const int d, const Variable& x)
1208{
1209  if (d <= 1)
1210  {
1211    A= F;
1212    return;
1213  }
1214  if (degree (F, x) <= 0)
1215  {
1216    A= F;
1217    return;
1218  }
1219  CanonicalForm C= 0;
1220  CanonicalForm f= swapvar (F, x, F.mvar());
1221  for (CFIterator i= f; i.hasTerms(); i++)
1222    C += i.coeff()*power (f.mvar(), i.exp()/ d);
1223  A= swapvar (C, x, F.mvar());
1224}
1225
1226CanonicalForm
1227reverseSubst (const CanonicalForm& F, const int d, const Variable& x)
1228{
1229  if (d <= 1)
1230    return F;
1231  if (degree (F, x) <= 0)
1232    return F;
1233  CanonicalForm f= swapvar (F, x, F.mvar());
1234  CanonicalForm result= 0;
1235  for (CFIterator i= f; i.hasTerms(); i++)
1236    result += i.coeff()*power (f.mvar(), d*i.exp());
1237  return swapvar (result, x, F.mvar());
1238}
1239
1240void
1241reverseSubst (CFList& L, const int d, const Variable& x)
1242{
1243  for (CFListIterator i= L; i.hasItem(); i++)
1244    i.getItem()= reverseSubst (i.getItem(), d, x);
1245}
1246
Note: See TracBrowser for help on using the repository browser.