source: git/factory/fac_univar.cc @ 398b15

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