source: git/factory/fac_univar.cc @ 1101a8

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