source: git/factory/facFqBivarUtil.cc

spielwiese
Last change on this file was 88faa1, checked in by Hans Schoenemann <hannes@…>, 3 months ago
compiler warning
  • 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=0;
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  nmod_poly_clear(FLINTF);
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.