source:git/factory/cfGcdUtil.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.1 KB
Line
1#include "config.h"
2
3#include "canonicalform.h"
4#include "cf_factory.h"
5#include "cf_reval.h"
6#include "cf_util.h"
7#include "cf_iter.h"
8#include "gfops.h"
9#include "cf_map_ext.h"
10#include "templates/ftmpl_functions.h"
11
12#ifdef HAVE_NTL
13#include "NTLconvert.h"
14#endif
15
16#ifdef HAVE_FLINT
17#include "FLINTconvert.h"
18#endif
19
20/// Coprimality Check. f and g are assumed to have the same level. If swap is
21/// true, the main variables of f and g are swapped with Variable(1). If the
22/// result is false, d is set to the degree of the gcd of f and g evaluated at a
23/// random point in K^n-1. This gcd is a gcd of univariate polynomials.
24#ifdef HAVE_NTL // mapPrimElem
25bool
26gcd_test_one ( const CanonicalForm & f, const CanonicalForm & g, bool swap, int & d )
27{
28    d= 0;
29    int count = 0;
30    // assume polys have same level;
31
32    Variable v= Variable (1);
33    bool algExtension= (hasFirstAlgVar (f, v) || hasFirstAlgVar (g, v));
34    CanonicalForm lcf, lcg;
35    if ( swap )
36    {
37        lcf = swapvar( LC( f ), Variable(1), f.mvar() );
38        lcg = swapvar( LC( g ), Variable(1), f.mvar() );
39    }
40    else
41    {
42        lcf = LC( f, Variable(1) );
43        lcg = LC( g, Variable(1) );
44    }
45
46    CanonicalForm F, G;
47    if ( swap )
48    {
49        F=swapvar( f, Variable(1), f.mvar() );
50        G=swapvar( g, Variable(1), g.mvar() );
51    }
52    else
53    {
54        F = f;
55        G = g;
56    }
57
58    #define TEST_ONE_MAX 50
59    int p= getCharacteristic();
60    bool passToGF= false;
61    int k= 1;
62    bool extOfExt= false;
63    Variable v3;
64    if (p > 0 && p < TEST_ONE_MAX && CFFactory::gettype() != GaloisFieldDomain && !algExtension)
65    {
66      if (p == 2)
67        setCharacteristic (2, 6, 'Z');
68      else if (p == 3)
69        setCharacteristic (3, 4, 'Z');
70      else if (p == 5 || p == 7)
71        setCharacteristic (p, 3, 'Z');
72      else
73        setCharacteristic (p, 2, 'Z');
74      passToGF= true;
75    }
76    else if (p > 0 && CFFactory::gettype() == GaloisFieldDomain && ipower (p , getGFDegree()) < TEST_ONE_MAX)
77    {
78      k= getGFDegree();
79      if (ipower (p, 2*k) > TEST_ONE_MAX)
80        setCharacteristic (p, 2*k, gf_name);
81      else
82        setCharacteristic (p, 3*k, gf_name);
83      F= GFMapUp (F, k);
84      G= GFMapUp (G, k);
85      lcf= GFMapUp (lcf, k);
86      lcg= GFMapUp (lcg, k);
87    }
88    else if (p > 0 && p < TEST_ONE_MAX && algExtension)
89    {
90#if defined(HAVE_NTL) || defined(HAVE_FLINT)
91      int d= degree (getMipo (v));
92      CFList source, dest;
93      Variable v2;
94      CanonicalForm primElem, imPrimElem;
95      if (p == 2 && d < 6)
96      {
97        bool primFail= false;
98        Variable vBuf;
99        primElem= primitiveElement (v, vBuf, primFail);
100        ASSERT (!primFail, "failure in integer factorizer");
101        if (d < 3)
102        {
103          #ifdef HAVE_FLINT
104          nmod_poly_t Irredpoly;
105          nmod_poly_init(Irredpoly,p);
106          nmod_poly_randtest_monic_irreducible(Irredpoly, FLINTrandom, 3*d+1);
107          CanonicalForm newMipo=convertnmod_poly_t2FacCF(Irredpoly,Variable(1));
108          nmod_poly_clear(Irredpoly);
109          #elif defined(HAVE_NTL)
110          if (fac_NTL_char != 2)
111          {
112            fac_NTL_char= 2;
113            zz_p::init (p);
114          }
115          zz_pX NTLIrredpoly;
116          BuildIrred (NTLIrredpoly, d*3);
117          CanonicalForm newMipo= convertNTLzzpX2CF (NTLIrredpoly, Variable (1));
118          #endif
119          v2= rootOf (newMipo);
120        }
121        else
122        {
123          #ifdef HAVE_FLINT
124          nmod_poly_t Irredpoly;
125          nmod_poly_init(Irredpoly,p);
126          nmod_poly_randtest_monic_irreducible(Irredpoly, FLINTrandom, 3*d+1);
127          CanonicalForm newMipo=convertnmod_poly_t2FacCF(Irredpoly,Variable(1));
128          nmod_poly_clear(Irredpoly);
129          #elif defined(HAVE_NTL)
130          if (fac_NTL_char != 2)
131          {
132            fac_NTL_char= 2;
133            zz_p::init (p);
134          }
135          zz_pX NTLIrredpoly;
136          BuildIrred (NTLIrredpoly, d*2);
137          CanonicalForm newMipo= convertNTLzzpX2CF (NTLIrredpoly, Variable (1));
138          #endif
139          v2= rootOf (newMipo);
140        }
141        imPrimElem= mapPrimElem (primElem, v, v2);
142        extOfExt= true;
143      }
144      else if ((p == 3 && d < 4) || ((p == 5 || p == 7) && d < 3))
145      {
146        bool primFail= false;
147        Variable vBuf;
148        primElem= primitiveElement (v, vBuf, primFail);
149        ASSERT (!primFail, "failure in integer factorizer");
150        #ifdef HAVE_FLINT
151        nmod_poly_t Irredpoly;
152        nmod_poly_init(Irredpoly,p);
153        nmod_poly_randtest_monic_irreducible(Irredpoly, FLINTrandom, 2*d+1);
154        CanonicalForm newMipo=convertnmod_poly_t2FacCF(Irredpoly,Variable(1));
155        nmod_poly_clear(Irredpoly);
156        #elif defined(HAVE_NTL)
157        if (fac_NTL_char != p)
158        {
159          fac_NTL_char= p;
160          zz_p::init (p);
161        }
162        zz_pX NTLIrredpoly;
163        BuildIrred (NTLIrredpoly, d*2);
164        CanonicalForm newMipo= convertNTLzzpX2CF (NTLIrredpoly, Variable (1));
165        #endif
166        v2= rootOf (newMipo);
167        imPrimElem= mapPrimElem (primElem, v, v2);
168        extOfExt= true;
169      }
170      if (extOfExt)
171      {
172        v3= v;
173        F= mapUp (F, v, v2, primElem, imPrimElem, source, dest);
174        G= mapUp (G, v, v2, primElem, imPrimElem, source, dest);
175        lcf= mapUp (lcf, v, v2, primElem, imPrimElem, source, dest);
176        lcg= mapUp (lcg, v, v2, primElem, imPrimElem, source, dest);
177        v= v2;
178      }
179#endif
180    }
181
182    CFRandom * sample;
183    if ((!algExtension && p > 0) || p == 0)
184      sample = CFRandomFactory::generate();
185    else
186      sample = AlgExtRandomF (v).clone();
187
188    REvaluation e( 2, tmax( f.level(), g.level() ), *sample );
189    delete sample;
190
191    if (passToGF)
192    {
193      lcf= lcf.mapinto();
194      lcg= lcg.mapinto();
195    }
196
197    CanonicalForm eval1, eval2;
198    if (passToGF)
199    {
200      eval1= e (lcf);
201      eval2= e (lcg);
202    }
203    else
204    {
205      eval1= e (lcf);
206      eval2= e (lcg);
207    }
208
209    while ( ( eval1.isZero() || eval2.isZero() ) && count < TEST_ONE_MAX )
210    {
211        e.nextpoint();
212        count++;
213        eval1= e (lcf);
214        eval2= e (lcg);
215    }
216    if ( count >= TEST_ONE_MAX )
217    {
218        if (passToGF)
219          setCharacteristic (p);
220        if (k > 1)
221          setCharacteristic (p, k, gf_name);
222        if (extOfExt)
223          prune1 (v3);
224        return false;
225    }
226
227
228    if (passToGF)
229    {
230      F= F.mapinto();
231      G= G.mapinto();
232      eval1= e (F);
233      eval2= e (G);
234    }
235    else
236    {
237      eval1= e (F);
238      eval2= e (G);
239    }
240
241    CanonicalForm c= gcd (eval1, eval2);
242    d= c.degree();
243    bool result= d < 1;
244    if (d < 0)
245      d= 0;
246
247    if (passToGF)
248      setCharacteristic (p);
249    if (k > 1)
250      setCharacteristic (p, k, gf_name);
251    if (extOfExt)
252      prune1 (v3);
253    return result;
254}
255#endif
256
257/**
258 * same as balance_p ( const CanonicalForm & f, const CanonicalForm & q )
259 * but qh= q/2 is provided, too.
260**/
261CanonicalForm
262balance_p ( const CanonicalForm & f, const CanonicalForm & q, const CanonicalForm & qh )
263{
264    Variable x = f.mvar();
265    CanonicalForm result = 0;
266    CanonicalForm c;
267    CFIterator i;
268    for ( i = f; i.hasTerms(); i++ )
269    {
270        c = i.coeff();
271        if ( c.inCoeffDomain())
272        {
273          if ( c > qh )
274            result += power( x, i.exp() ) * (c - q);
275          else
276            result += power( x, i.exp() ) * c;
277        }
278        else
279          result += power( x, i.exp() ) * balance_p(c,q,qh);
280    }
281    return result;
282}
283
284/** static CanonicalForm balance_p ( const CanonicalForm & f, const CanonicalForm & q )
285 *
286 * balance_p() - map f from positive to symmetric representation
287 *   mod q.
288 *
289 * This makes sense for polynomials over Z only.
290 * q should be an integer.
291 *
292**/
293CanonicalForm
294balance_p ( const CanonicalForm & f, const CanonicalForm & q )
295{
296    CanonicalForm qh = q / 2;
297    return balance_p (f, q, qh);
298}
299
300
301/*static CanonicalForm
302balance ( const CanonicalForm & f, const CanonicalForm & q )
303{
304    Variable x = f.mvar();
305    CanonicalForm result = 0, qh = q / 2;
306    CanonicalForm c;
307    CFIterator i;
308    for ( i = f; i.hasTerms(); i++ ) {
309        c = mod( i.coeff(), q );
310        if ( c > qh )
311            result += power( x, i.exp() ) * (c - q);
312        else
313            result += power( x, i.exp() ) * c;
314    }
315    return result;
316}*/
Note: See TracBrowser for help on using the repository browser.