source: git/factory/cfGcdUtil.cc @ 96ce32

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