source: git/factory/fac_univar.cc @ 3ace5b6

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