source: git/factory/cf_gcd.cc @ b7cc2b

spielwiese
Last change on this file since b7cc2b was b7cc2b, checked in by Hans Schoenemann <hannes@…>, 3 years ago
factory: FindRoot via FLINT
  • Property mode set to 100644
File size: 8.9 KB
Line 
1/* emacs edit mode for this file is -*- C++ -*- */
2
3/**
4 * @file cf_gcd.cc
5 *
6 * gcd/content/lcm of polynomials
7 *
8 * To compute the GCD different variants are chosen automatically
9**/
10
11#include "config.h"
12
13
14#include "timing.h"
15#include "cf_assert.h"
16#include "debug.h"
17
18#include "cf_defs.h"
19#include "canonicalform.h"
20#include "cf_iter.h"
21#include "cf_reval.h"
22#include "cf_primes.h"
23#include "cf_algorithm.h"
24#include "cfEzgcd.h"
25#include "cfGcdAlgExt.h"
26#include "cfSubResGcd.h"
27#include "cfModGcd.h"
28#include "FLINTconvert.h"
29#include "facAlgFuncUtil.h"
30
31#ifdef HAVE_NTL
32#include <NTL/ZZX.h>
33#include "NTLconvert.h"
34bool isPurePoly(const CanonicalForm & );
35#endif
36
37void out_cf(const char *s1,const CanonicalForm &f,const char *s2);
38
39/** static CanonicalForm icontent ( const CanonicalForm & f, const CanonicalForm & c )
40 *
41 * icontent() - return gcd of c and all coefficients of f which
42 *   are in a coefficient domain.
43 *
44 * @sa icontent().
45 *
46**/
47static CanonicalForm
48icontent ( const CanonicalForm & f, const CanonicalForm & c )
49{
50    if ( f.inBaseDomain() )
51    {
52      if (c.isZero()) return abs(f);
53      return bgcd( f, c );
54    }
55    //else if ( f.inCoeffDomain() )
56    //   return gcd(f,c);
57    else
58    {
59        CanonicalForm g = c;
60        for ( CFIterator i = f; i.hasTerms() && ! g.isOne(); i++ )
61            g = icontent( i.coeff(), g );
62        return g;
63    }
64}
65
66/** CanonicalForm icontent ( const CanonicalForm & f )
67 *
68 * icontent() - return gcd over all coefficients of f which are
69 *   in a coefficient domain.
70 *
71**/
72CanonicalForm
73icontent ( const CanonicalForm & f )
74{
75    return icontent( f, 0 );
76}
77
78/** CanonicalForm gcd_poly ( const CanonicalForm & f, const CanonicalForm & g )
79 *
80 * gcd_poly() - calculate polynomial gcd.
81 *
82 * This is the dispatcher for polynomial gcd calculation.
83 * Different gcd variants get called depending the input, characteristic, and
84 * on switches (cf_defs.h)
85 *
86 * With the current settings from Singular (i.e. SW_USE_EZGCD= on,
87 * SW_USE_EZGCD_P= on, SW_USE_CHINREM_GCD= on, the EZ GCD variants are the
88 * default algorithms for multivariate polynomial GCD computations)
89 *
90 * @sa gcd(), cf_defs.h
91 *
92**/
93CanonicalForm gcd_poly ( const CanonicalForm & f, const CanonicalForm & g )
94{
95  CanonicalForm fc, gc;
96  bool fc_isUnivariate=f.isUnivariate();
97  bool gc_isUnivariate=g.isUnivariate();
98  bool fc_and_gc_Univariate=fc_isUnivariate && gc_isUnivariate;
99  fc = f;
100  gc = g;
101  if ( getCharacteristic() != 0 )
102  {
103    #if defined(HAVE_FLINT) && ( __FLINT_RELEASE >= 20503)
104    if ( isOn( SW_USE_FL_GCD_P)
105    && (CFFactory::gettype() != GaloisFieldDomain)
106    && (getCharacteristic()>10)
107    &&(!hasAlgVar(fc)) && (!hasAlgVar(gc)))
108    {
109      return gcdFlintMP_Zp(fc,gc);
110    }
111    #endif
112    #ifdef HAVE_NTL
113    if ((!fc_and_gc_Univariate) && (isOn( SW_USE_EZGCD_P )))
114    {
115      fc= EZGCD_P (fc, gc);
116    }
117    #endif
118    #ifdef HAVE_NTL
119    else if (isOn(SW_USE_FF_MOD_GCD) && !fc_and_gc_Univariate)
120    {
121      Variable a;
122      if (hasFirstAlgVar (fc, a) || hasFirstAlgVar (gc, a))
123        fc=modGCDFq (fc, gc, a);
124      else if (CFFactory::gettype() == GaloisFieldDomain)
125        fc=modGCDGF (fc, gc);
126      else
127        fc=modGCDFp (fc, gc);
128    }
129    #endif
130    else
131    fc = subResGCD_p( fc, gc );
132  }
133  else if (!fc_and_gc_Univariate) /* && char==0*/
134  {
135    #if defined(HAVE_FLINT) && ( __FLINT_RELEASE >= 20503)
136    if (( isOn( SW_USE_FL_GCD_0) )
137    &&(!hasAlgVar(fc)) && (!hasAlgVar(gc)))
138    {
139      return gcdFlintMP_QQ(fc,gc);
140    }
141    else
142    #endif
143    #ifdef HAVE_NTL
144    if ( isOn( SW_USE_EZGCD ) )
145      fc= ezgcd (fc, gc);
146    #endif 
147    #ifdef HAVE_NTL
148    else if (isOn(SW_USE_CHINREM_GCD))
149      fc = modGCDZ( fc, gc);
150    else
151    #endif
152    {
153       fc = subResGCD_0( fc, gc );
154    }
155  }
156  else
157  {
158    fc = subResGCD_0( fc, gc );
159  }
160  if ((getCharacteristic()>0)&&(!hasAlgVar(fc))) fc/=fc.lc();
161  return fc;
162}
163
164/** static CanonicalForm cf_content ( const CanonicalForm & f, const CanonicalForm & g )
165 *
166 * cf_content() - return gcd(g, content(f)).
167 *
168 * content(f) is calculated with respect to f's main variable.
169 *
170 * @sa gcd(), content(), content( CF, Variable ).
171 *
172**/
173static CanonicalForm
174cf_content ( const CanonicalForm & f, const CanonicalForm & g )
175{
176    if ( f.inPolyDomain() || ( f.inExtension() && ! getReduce( f.mvar() ) ) )
177    {
178        CFIterator i = f;
179        CanonicalForm result = g;
180        while ( i.hasTerms() && ! result.isOne() )
181        {
182            result = gcd( i.coeff(), result );
183            i++;
184        }
185        return result;
186    }
187    else
188        return abs( f );
189}
190
191/** CanonicalForm content ( const CanonicalForm & f )
192 *
193 * content() - return content(f) with respect to main variable.
194 *
195 * Normalizes result.
196 *
197**/
198CanonicalForm
199content ( const CanonicalForm & f )
200{
201    if ( f.inPolyDomain() || ( f.inExtension() && ! getReduce( f.mvar() ) ) )
202    {
203        CFIterator i = f;
204        CanonicalForm result = abs( i.coeff() );
205        i++;
206        while ( i.hasTerms() && ! result.isOne() )
207        {
208            result = gcd( i.coeff(), result );
209            i++;
210        }
211        return result;
212    }
213    else
214        return abs( f );
215}
216
217/** CanonicalForm content ( const CanonicalForm & f, const Variable & x )
218 *
219 * content() - return content(f) with respect to x.
220 *
221 * x should be a polynomial variable.
222 *
223**/
224CanonicalForm
225content ( const CanonicalForm & f, const Variable & x )
226{
227    if (f.inBaseDomain()) return f;
228    ASSERT( x.level() > 0, "cannot calculate content with respect to algebraic variable" );
229    Variable y = f.mvar();
230
231    if ( y == x )
232        return cf_content( f, 0 );
233    else  if ( y < x )
234        return f;
235    else
236        return swapvar( content( swapvar( f, y, x ), y ), y, x );
237}
238
239/** CanonicalForm vcontent ( const CanonicalForm & f, const Variable & x )
240 *
241 * vcontent() - return content of f with repect to variables >= x.
242 *
243 * The content is recursively calculated over all coefficients in
244 * f having level less than x.  x should be a polynomial
245 * variable.
246 *
247**/
248CanonicalForm
249vcontent ( const CanonicalForm & f, const Variable & x )
250{
251    ASSERT( x.level() > 0, "cannot calculate vcontent with respect to algebraic variable" );
252
253    if ( f.mvar() <= x )
254        return content( f, x );
255    else {
256        CFIterator i;
257        CanonicalForm d = 0;
258        for ( i = f; i.hasTerms() && ! d.isOne(); i++ )
259            d = gcd( d, vcontent( i.coeff(), x ) );
260        return d;
261    }
262}
263
264/** CanonicalForm pp ( const CanonicalForm & f )
265 *
266 * pp() - return primitive part of f.
267 *
268 * Returns zero if f equals zero, otherwise f / content(f).
269 *
270**/
271CanonicalForm
272pp ( const CanonicalForm & f )
273{
274    if ( f.isZero() )
275        return f;
276    else
277        return f / content( f );
278}
279
280CanonicalForm
281gcd ( const CanonicalForm & f, const CanonicalForm & g )
282{
283    bool b = f.isZero();
284    if ( b || g.isZero() )
285    {
286        if ( b )
287            return abs( g );
288        else
289            return abs( f );
290    }
291    if ( f.inPolyDomain() || g.inPolyDomain() )
292    {
293        if ( f.mvar() != g.mvar() )
294        {
295            if ( f.mvar() > g.mvar() )
296                return cf_content( f, g );
297            else
298                return cf_content( g, f );
299        }
300        if (isOn(SW_USE_QGCD))
301        {
302          Variable m;
303          if (
304          (getCharacteristic() == 0) &&
305          (hasFirstAlgVar(f,m) || hasFirstAlgVar(g,m))
306          )
307          {
308            bool on_rational = isOn(SW_RATIONAL);
309            CanonicalForm r=QGCD(f,g);
310            On(SW_RATIONAL);
311            CanonicalForm cdF = bCommonDen( r );
312            if (!on_rational) Off(SW_RATIONAL);
313            return cdF*r;
314          }
315        }
316
317        if ( f.inExtension() && getReduce( f.mvar() ) )
318            return CanonicalForm(1);
319        else
320        {
321            if ( fdivides( f, g ) )
322                return abs( f );
323            else  if ( fdivides( g, f ) )
324                return abs( g );
325            if ( !( getCharacteristic() == 0 && isOn( SW_RATIONAL ) ) )
326            {
327                CanonicalForm d;
328                d = gcd_poly( f, g );
329                return abs( d );
330            }
331            else
332            {
333                CanonicalForm cdF = bCommonDen( f );
334                CanonicalForm cdG = bCommonDen( g );
335                CanonicalForm F = f * cdF, G = g * cdG;
336                Off( SW_RATIONAL );
337                CanonicalForm l = gcd_poly( F, G );
338                On( SW_RATIONAL );
339                return abs( l );
340            }
341        }
342    }
343    if ( f.inBaseDomain() && g.inBaseDomain() )
344        return bgcd( f, g );
345    else
346        return 1;
347}
348
349/** CanonicalForm lcm ( const CanonicalForm & f, const CanonicalForm & g )
350 *
351 * lcm() - return least common multiple of f and g.
352 *
353 * The lcm is calculated using the formula lcm(f, g) = f * g / gcd(f, g).
354 *
355 * Returns zero if one of f or g equals zero.
356 *
357**/
358CanonicalForm
359lcm ( const CanonicalForm & f, const CanonicalForm & g )
360{
361    if ( f.isZero() || g.isZero() )
362        return 0;
363    else
364        return ( f / gcd( f, g ) ) * g;
365}
366
Note: See TracBrowser for help on using the repository browser.