source: git/factory/fac_univar.cc @ a8fbae

fieker-DuValspielwiese
Last change on this file since a8fbae was 564dd8, checked in by Hans Schoenemann <hannes@…>, 4 years ago
factorize in Z[x,..] w/o NTL
  • Property mode set to 100644
File size: 17.1 KB
Line 
1/* emacs edit mode for this file is -*- C++ -*- */
2
3
4#include "config.h"
5
6
7#include "cf_assert.h"
8#include "debug.h"
9#include "timing.h"
10
11#include "cf_defs.h"
12#include "cf_algorithm.h"
13#include "fac_util.h"
14#include "fac_univar.h"
15#include "fac_cantzass.h"
16#include "fac_berlekamp.h"
17#include "cf_iter.h"
18#include "cf_primes.h"
19#include "fac_sqrfree.h"
20#include "cfUnivarGcd.h"
21
22#ifdef HAVE_FLINT
23#include "FLINTconvert.h" // for __FLINT_RELEASE
24#endif
25
26#if !defined(HAVE_NTL)
27TIMING_DEFINE_PRINT(fac_choosePrimes)
28TIMING_DEFINE_PRINT(fac_facModPrimes)
29TIMING_DEFINE_PRINT(fac_liftFactors)
30TIMING_DEFINE_PRINT(fac_combineFactors)
31
32
33const int max_fp_fac = 3;
34
35STATIC_INST_VAR modpk theModulus;
36
37#ifdef DEBUGOUTPUT
38#define DEBOUTHPRINT(stream, msg, hg) \
39{std::stream << deb_level_msg << msg, std::stream.flush(); hprint( hg ); std::stream << std::endl;}
40
41static void
42hprint ( int * a )
43{
44    int n = a[0];
45    std::cerr << "( " << n << ": ";
46    int i = 1;
47    while ( i < n ) {
48        if ( a[i] != 0 )
49            std::cerr << i << " ";
50        i++;
51    }
52    std::cerr << ")";
53}
54#else /* DEBUGOUTPUT */
55#define DEBOUTHPRINT(stream, msg, hg)
56#endif /* DEBUGOUTPUT */
57
58static void
59hgroup ( int * a )
60{
61    int n = a[0];
62    int i, j, k;
63    for ( i = 1; i < n; i++ )
64        if ( a[i] != 0 )
65            for ( j = 1; j <= i; j++ )
66                if ( a[j] != 0 )
67                    for ( k = i; k < n; k += j )
68                        a[k] = 1;
69}
70
71static void
72hintersect( int * a, const int * const b )
73{
74    int i, n, na = a[0], nb = b[0];
75    if ( nb < na )
76        n = nb;
77    else
78        n = na;
79    for ( i = 1; i < n; i++ )
80        if ( b[i] == 0 )
81            a[i] = 0;
82    a[0] = n;
83}
84
85/*
86static int
87hcount ( const int * const a )
88{
89    int n = a[0], sum = 0, i;
90    for ( i = 1; i < n; i++ )
91        if ( a[i] != 0 ) sum++;
92    return sum;
93}
94*/
95
96static void
97initHG ( int * a, const CFFList & F )
98{
99    ListIterator<CFFactor> i;
100
101    int n = a[0], k;
102    for ( int j = 1; j < n; j++ ) a[j] = 0;
103    for ( i = F; i.hasItem(); i++ )
104        if ( (k = i.getItem().factor().degree()) < n )
105        {
106            if ( k == -1 ) {
107                STICKYWARN( k == -1, "there occured an error.  factory was not able to factorize\n"
108                            "correctly mod p.  Please send the example which caused\n"
109                            "this error to the authors.  Nonetheless we will go on with the\n"
110                            "calculations hoping the result will be correct.  Thank you." );
111            }
112            else if ( k != 0 )
113                a[k] = 1;
114        }
115}
116
117static void
118initHG ( int * a, const Array<CanonicalForm> & F )
119{
120    int i, n = a[0], m = F.size(), k;
121    for ( i = 1; i < n; i++ ) a[i] = 0;
122    for ( i = 1; i < m; i++ )
123        if ( (k = F[i].degree()) < n )
124        {
125            if ( k == -1 )
126            {
127                STICKYWARN( k == -1, "there occured an error.  factory was not able to factorize\n"
128                            "correctly mod p.  Please send the example which caused\n"
129                            "this error to the authors.  Nonetheless we will go on with the\n"
130                            "calculations hoping the result will be correct.  Thank you." );
131            }
132            else if ( k != 0 )
133                a[k] = 1;
134        }
135}
136
137static int cmpFactor( const CFFactor & a, const CFFactor & b )
138{
139    CFFactor A( a ), B( b );
140    return degree( A.factor() ) > degree( B.factor() );
141}
142
143//static double cf2double ( const CanonicalForm & f )
144//{
145//    CanonicalForm a = f, q, r;
146//    double m = 1, res = 0;
147//    if ( a.sign() < 0 ) a = -a;
148//    while ( ! a.isZero() ) {
149//        divrem( a, 10, q, r );
150//        res += m * (double)(r.intval());
151//        m *= 10;
152//        a = q;
153//    }
154//    if ( f.sign() < 0 ) res = -res;
155//    return res;
156//}
157
158//{{{ static int kBound ( const CanonicalForm & f, int p )
159//{{{ docu
160//
161// kBound() - return bound of coefficients of factors of f.
162//
163// The bound is returned as an integer k such that p^k is larger
164// than all coefficients of all possible factors of f.  f should
165// be an univariate polynomial over Z.
166//
167// For a discussion of the formula, see the article Mignotte -
168// 'Some Usefull Bounds' in Buchberger, Collins, Loos (eds.) -
169// 'Computer Algebra: Symbolic and Algebraic Computation', 2nd
170// ed.
171//
172// Use by ZFactorizeUnivariate().
173//
174//}}}
175static int
176kBound ( const CanonicalForm & f, int p )
177{
178    return (int)(f.degree() + (double)(ilog2( euclideanNorm(f)+1 ) + 1) / (double)ilog2(p)) + 1;
179}
180//}}}
181
182modpk
183getZFacModulus()
184{
185    return theModulus;
186}
187
188static bool
189liftDegreeFactRec( CFArray & theFactors, CanonicalForm & F, const CanonicalForm & recip_lf, const CanonicalForm & f, const modpk & pk, int i, int d, CFFList & ZF, int exp )
190{
191    if ( i >= theFactors.size() )
192        return false;
193    else  if ( degree( f ) + degree( theFactors[i] ) == d ) {
194        DEBOUTLN( cerr, "ldfr (f) = " << f );
195        DEBOUTLN( cerr, "ldfr (g) = " << theFactors[i] );
196        CanonicalForm g = pp( pk( recip_lf * f * theFactors[i] ) );
197        DEBOUTLN( cerr, "ldfr (pk(f*g)) = " << g );
198        CanonicalForm gg, hh;
199        DEBOUTLN( cerr, "F = " << F );
200        DEBOUTLN( cerr, "g = " << g );
201        if ( divremt( F, g, gg, hh ) && hh.isZero() ) {
202            ZF.append( CFFactor( g, exp ) );
203            F = gg;
204            theFactors[i] = 1;
205            return true;
206        }
207        else {
208            return liftDegreeFactRec( theFactors, F, recip_lf, f, pk, i+1, d, ZF, exp );
209        }
210    }
211    else  if ( degree( f ) + degree( theFactors[i] ) > d )
212        return false;
213    else {
214        bool ok = liftDegreeFactRec( theFactors, F, recip_lf, pk( recip_lf * f * theFactors[i] ), pk, i+1, d, ZF, exp );
215        if ( ok )
216            theFactors[i] = 1;
217        else
218            ok = liftDegreeFactRec( theFactors, F, recip_lf, f, pk, i+1, d, ZF, exp );
219        return ok;
220    }
221}
222
223bool isSqrFreeZ ( const CanonicalForm & f )
224{
225    return gcd( f, f.deriv() ).degree() == 0;
226}
227
228bool isSqrFreeFp( const CanonicalForm & f )
229{
230  CFFList F = sqrFreeFp( f );
231  return ( F.length() == 1 && F.getFirst().exp() == 1 );
232}
233
234bool isSqrFree ( const CanonicalForm & f )
235{
236//    ASSERT( f.isUnivariate(), "multivariate factorization not implemented" );
237    if ( getCharacteristic() == 0 )
238        return isSqrFreeZ( f );
239    else
240        return isSqrFreeFp( f );
241}
242
243
244static int choosePrimes ( int * p, const CanonicalForm & f )
245{
246    int ptr = 0;
247    int i = 0;
248    int maxp = cf_getNumPrimes();
249    int prime;
250
251    while ( ptr < maxp && i < max_fp_fac )
252    {
253        prime = cf_getPrime( ptr );
254        if ( mod( lc( f ), prime ) != 0 )
255        {
256            setCharacteristic( prime );
257            if ( isSqrFree( mapinto( f ) ) )
258            {
259                p[i] = prime;
260                i++;
261            }
262            setCharacteristic( 0 );
263        }
264        ptr++;
265    }
266    return ( i == max_fp_fac );
267}
268
269
270static int
271UnivariateQuadraticLift ( const CanonicalForm &F, const  CanonicalForm & G, const CanonicalForm &H, const modpk & pk, const CanonicalForm & /*Gamma*/, CanonicalForm & gk, CanonicalForm & hk )
272{
273    CanonicalForm lf, f, gamma;
274    CanonicalForm a, b, aa, bb, c, g, h, g1, h1, e, modulus, tmp, q, r;
275    int i, j, save;
276    int p = pk.getp(), k = pk.getk();
277    int no_iter = SI_LOG2(k)+2;
278    int * kvals = new int[no_iter];
279
280    DEBOUTLN( cerr, "quadratic lift called with p = " << p << "  and k = " << k );
281    for ( j = 0, i = k; i > 1; i = ( i+1 ) / 2, j++ ) kvals[j] = i;
282    kvals[j] = 1;
283
284    save = getCharacteristic(); setCharacteristic( 0 );
285
286    lf = lc( F );
287    f = lf * F;
288    {
289        setCharacteristic( p );
290        g1 = mapinto( lf ) / lc( G ) * G;
291        h1 = mapinto( lf ) / lc( H ) * H;
292        (void)extgcd( g1, h1, a, b );
293        setCharacteristic( 0 );
294    }
295    a = mapinto( a ); b = mapinto( b );
296    g = mapinto( g1 ); h = mapinto( h1 );
297    g = replaceLc( g, lf ); h = replaceLc( h, lf );
298    e = f - g * h;
299    modulus = p;
300    i = 1;
301
302    while ( ! e.isZero() && j > 0 ) {
303        c = e / modulus;
304        {
305            j--;
306            setCharacteristic( p, kvals[j+1]);
307            DEBOUTLN( cerr, "lifting from p^" << kvals[j+1] << " to p^" << kvals[j] );
308            c = mapinto( c );
309            DEBOUTLN( cerr, " !!! g = " << mapinto( g ) );
310            g1 = mapinto( lf ) / mapinto( lc( g ) ) * mapinto( g );
311            h1 = mapinto( lf ) / mapinto( lc( h ) ) * mapinto( h );
312//            (void)extgcd( g1, h1, a, b );
313//            DEBOUTLN( cerr, " a = " << aa );
314//            DEBOUTLN( cerr, " b = " << bb );
315            a = mapinto( a ); b = mapinto( b );
316            a += ( ( 1 - a * g1 ) *  a ) % h1;
317            b += ( ( 1 - b * h1 ) *  b ) % g1;
318            DEBOUTLN( cerr, " a = " << a );
319            DEBOUTLN( cerr, " b = " << b );
320            divrem( a * c, h1, q, r );
321            tmp = b * c + q * g1;
322            setCharacteristic( 0 );
323        }
324        a = mapinto( a ); b = mapinto( b );
325        g += mapinto( tmp ) * modulus;
326        h += mapinto( r ) * modulus;
327//        g = replaceLc( g, lf ); h = replaceLc( h, lf );
328        e = f - g * h;
329        modulus = power( CanonicalForm(p), kvals[j] );
330        if ( mod( f - g * h, modulus ) != 0 ) {
331            DEBOUTLN( cerr, "error at lift stage " << i );
332        }
333        i++;
334    }
335    if ( e.isZero() ) {
336        tmp = content( g );
337        gk = g / tmp; hk = h / ( lf / tmp );
338    }
339    else {
340        gk = pk(g); hk = pk(h);
341    }
342    setCharacteristic( save );
343    delete [] kvals;
344    return e.isZero();
345}
346
347static int
348UnivariateLinearLift ( const CanonicalForm &F, const  CanonicalForm & G, const CanonicalForm &H, const modpk & pk, const CanonicalForm & /*Gamma*/, CanonicalForm & gk, CanonicalForm & hk )
349{
350    CanonicalForm lf, f, gamma;
351    CanonicalForm a, b, c, g, h, g1, h1, e, modulus, tmp, q, r;
352    int i, save;
353    int p = pk.getp(), k = pk.getk();
354    save = getCharacteristic(); setCharacteristic( 0 );
355
356    lf = lc( F );
357    f = lf * F;
358    {
359        setCharacteristic( p );
360        g1 = mapinto( lf ) / lc( G ) * G;
361        h1 = mapinto( lf ) / lc( H ) * H;
362        (void)extgcd( g1, h1, a, b );
363        setCharacteristic( 0 );
364    }
365    g = mapinto( g1 ); h = mapinto( h1 );
366    g = replaceLc( g, lf ); h = replaceLc( h, lf );
367    e = f - g * h;
368    modulus = p;
369    i = 1;
370
371    while ( ! e.isZero() && i <= k ) {
372        c = e / modulus;
373        {
374            setCharacteristic( p );
375            c = mapinto( c );
376            divrem( a * c, h1, q, r );
377            tmp = b * c + q * g1;
378            setCharacteristic( 0 );
379        }
380        g += mapinto( tmp ) * modulus;
381        h += mapinto( r ) * modulus;
382//        g = replaceLc( g, lf ); h = replaceLc( h, lf );
383        e = f - g * h;
384        modulus *= p;
385        ASSERT( mod( f - g * h, modulus ) == 0, "error at lift stage" );
386        i++;
387    }
388    if ( e.isZero() ) {
389        tmp = content( g );
390        gk = g / tmp; hk = h / ( lf / tmp );
391    }
392    else {
393        gk = pk(g); hk = pk(h);
394    }
395    setCharacteristic( save );
396//    return e.isZero();
397    return (F-gk*hk).isZero();
398}
399
400
401CFFList
402ZFactorizeUnivariate( const CanonicalForm& ff, bool issqrfree )
403{
404    bool symmsave = isOn( SW_SYMMETRIC_FF );
405    CanonicalForm cont = content( ff );
406    CanonicalForm lf, recip_lf, fp, f, g = ff / cont, dummy1, dummy2;
407    int i, k, exp, n;
408    bool ok;
409    CFFList H, F[max_fp_fac];
410    CFFList ZF;
411    int * p = new int [max_fp_fac];
412    int * D = 0;
413    int * Dh = 0;
414    ListIterator<CFFactor> J, I;
415
416    DEBINCLEVEL( cerr, "ZFactorizeUnivariate" );
417    On( SW_SYMMETRIC_FF );
418
419    // get squarefree decomposition of f
420    if ( issqrfree )
421        H.append( CFFactor( g, 1 ) );
422    else
423        H = sqrFree( g );
424
425    DEBOUTLN( cerr, "H = " << H );
426
427    // cycle through squarefree factors of f
428    for ( J = H; J.hasItem(); ++J ) {
429        f = J.getItem().factor();
430        if ( f.inCoeffDomain() ) continue;
431        n = f.degree() / 2 + 1;
432        delete [] D;
433        delete [] Dh;
434        D = new int [n]; D[0] = n;
435        Dh = new int [n]; Dh[0] = n;
436        exp = J.getItem().exp();
437
438        // choose primes to factor f
439        TIMING_START(fac_choosePrimes);
440        ok = choosePrimes( p, f );
441        TIMING_END_AND_PRINT(fac_choosePrimes, "time to choose the primes: ");
442        if ( ! ok ) {
443            DEBOUTLN( cerr, "warning: no good prime found to factorize " << f );
444            STICKYWARN( ok, "there occured an error.  We went out of primes p\n"
445                        "to factorize mod p.  Please send the example which caused\n"
446                        "this error to the authors.  Nonetheless we will go on with the\n"
447                        "calculations hoping the result will be correct.  Thank you.");
448            ZF.append( CFFactor( f, exp ) );
449            continue;
450        }
451
452        // factorize f modulo certain primes
453        TIMING_START(fac_facModPrimes);
454        for ( i = 0; i < max_fp_fac; i++ ) {
455            setCharacteristic( p[i] );
456            fp = mapinto( f );
457            F[i] = FpFactorizeUnivariateCZ( fp, true, 0, Variable(), Variable() );
458//                if ( p[i] < 23 && fp.degree() < 10 )
459//                    F[i] = FpFactorizeUnivariateB( fp, true );
460//                else
461//                    F[i] = FpFactorizeUnivariateCZ( fp, true, 0, Variable, Variable() );
462            DEBOUTLN( cerr, "F[i] = " << F[i] << ", p = " << p[i] );
463        }
464        TIMING_END_AND_PRINT(fac_facModPrimes, "time to factorize mod primes: ");
465        setCharacteristic( 0 );
466
467        // do some strange things with the D's
468        initHG( D, F[0] );
469        hgroup( D );
470        DEBOUTHPRINT( cerr, "D = ", D );
471        for ( i = 1; i < max_fp_fac; i++ ) {
472            initHG( Dh, F[i] );
473            hgroup( Dh );
474            DEBOUTHPRINT( cerr, "Dh = ", Dh );
475            hintersect( D, Dh );
476            DEBOUTHPRINT( cerr, "D = ", D );
477        }
478
479        // look which p gives the shortest factorization of f mod p
480        // j: index of that p in p[]
481        int min, j;
482        min = F[0].length(), j = 0;
483        for ( i = 1; i < max_fp_fac; i++ ) {
484            if ( min >= F[i].length() ) {
485                j = i; min = F[i].length();
486            }
487        }
488        k = kBound( f, p[j] );
489        CFArray theFactors( F[j].length() );
490//            pk = power( CanonicalForm( p[j] ), k );
491//            pkhalf = pk / 2;
492        modpk pk( p[j], k );
493        DEBOUTLN( cerr, "coeff bound = " << pk.getpk() );
494        theModulus = pk;
495        setCharacteristic( p[j] );
496        fp = mapinto( f );
497        F[j].sort( cmpFactor );
498        I = F[j]; i = 0;
499        TIMING_START(fac_liftFactors);
500        while ( I.hasItem() ) {
501            DEBOUTLN( cerr, "factor to lift = " << I.getItem().factor() );
502            if ( isOn( SW_FAC_QUADRATICLIFT ) )
503                ok = UnivariateQuadraticLift( f, I.getItem().factor(), fp / I.getItem().factor(), pk, lc( f ), dummy1, dummy2 );
504            else
505                ok = UnivariateLinearLift( f, I.getItem().factor(), fp / I.getItem().factor(), pk, lc( f ), dummy1, dummy2 );
506            if ( ok ) {
507                // should be done in a more efficient way
508                DEBOUTLN( cerr, "dummy1 = " << dummy1 );
509                DEBOUTLN( cerr, "dummy2 = " << dummy2 );
510                f = dummy2;
511                fp /= I.getItem().factor();
512                ZF.append( CFFactor( dummy1, exp ) );
513                I.remove( 0 );
514                I = F[j];
515                i = 0;
516                DEBOUTLN( cerr, "F[j] = " << F[j] );
517            }
518            else {
519                DEBOUTLN( cerr, "i = " << i );
520                DEBOUTLN( cerr, "dummy1 = " << dummy1 );
521                setCharacteristic( 0 );
522//                    theFactors[i] = pk( dummy1 * pk.inverse( lc( dummy1 ) ) );
523                theFactors[i] = pk( dummy1 );
524                setCharacteristic( p[j] );
525                i++;
526                I++;
527            }
528        }
529        TIMING_END_AND_PRINT(fac_liftFactors, "time to lift the factors: ");
530        DEBOUTLN( cerr, "ZF = " << ZF );
531        initHG( Dh, theFactors );
532        hgroup( Dh );
533        DEBOUTHPRINT( cerr, "Dh = ", Dh );
534        hintersect( D, Dh );
535        setCharacteristic( 0 );
536        for ( int l = i; l < F[j].length(); l++ )
537            theFactors[l] = 1;
538        DEBOUTLN( cerr, "theFactors = " << theFactors );
539        DEBOUTLN( cerr, "f = " << f );
540        DEBOUTLN( cerr, "p = " << pk.getp() << ", k = " << pk.getk() );
541        DEBOUTHPRINT( cerr, "D = ", D );
542        lf = lc( f );
543        (void)bextgcd( pk.getpk(), lf, dummy1, recip_lf );
544        DEBOUTLN( cerr, "recip_lf = " << recip_lf );
545        TIMING_START(fac_combineFactors);
546        for ( i = 1; i < D[0]; i++ )
547            if ( D[i] != 0 )
548                while ( liftDegreeFactRec( theFactors, f, recip_lf, lf, pk, 0, i, ZF, exp ) );
549        if ( degree( f ) > 0 )
550            ZF.append( CFFactor( f, exp ) );
551        TIMING_END_AND_PRINT(fac_combineFactors, "time to combine the factors: ");
552    }
553
554    // brush up our result
555    if ( ZF.getFirst().factor().inCoeffDomain() )
556        ZF.removeFirst();
557    if ( lc( ff ).sign() < 0 )
558        ZF.insert( CFFactor( -cont, 1 ) );
559    else
560        ZF.insert( CFFactor( cont, 1 ) );
561    delete [] D;
562    delete [] Dh;
563    delete [] p;
564    if ( ! symmsave )
565        Off( SW_SYMMETRIC_FF );
566
567    DEBDECLEVEL( cerr, "ZFactorizeUnivariate" );
568    return ZF;
569}
570#endif
Note: See TracBrowser for help on using the repository browser.