source: git/factory/facFqBivarUtil.cc @ 6f1268

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