source: git/factory/facFqBivarUtil.cc @ 4e2385

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