source: git/factory/fac_univar.cc @ 19f74fe

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