source: git/factory/fac_univar.cc @ f0596e

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