source: git/factory/fac_univar.cc @ 346edc8

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