source: git/factory/fac_univar.cc @ 17b1f3

spielwiese
Last change on this file since 17b1f3 was c1b9927, checked in by Hans Schoenemann <hannes@…>, 13 years ago
- removed some unsed variables - never put static inline routine without a body in a .h file git-svn-id: file:///usr/local/Singular/svn/trunk@14265 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 16.4 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 "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        {
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        i++;
308    }
309    if ( e.isZero() ) {
310        tmp = content( g );
311        gk = g / tmp; hk = h / ( lf / tmp );
312    }
313    else {
314        gk = pk(g); hk = pk(h);
315    }
316    setCharacteristic( save );
317    return e.isZero();
318}
319
320static int
321UnivariateLinearLift ( const CanonicalForm &F, const  CanonicalForm & G, const CanonicalForm &H, const modpk & pk, const CanonicalForm & Gamma, CanonicalForm & gk, CanonicalForm & hk )
322{
323    CanonicalForm lf, f, gamma;
324    CanonicalForm a, b, c, g, h, g1, h1, e, modulus, tmp, q, r;
325    int i, save;
326    int p = pk.getp(), k = pk.getk();
327    save = getCharacteristic(); setCharacteristic( 0 );
328
329    lf = lc( F );
330    f = lf * F;
331    {
332        setCharacteristic( p );
333        g1 = mapinto( lf ) / lc( G ) * G;
334        h1 = mapinto( lf ) / lc( H ) * H;
335        (void)extgcd( g1, h1, a, b );
336        setCharacteristic( 0 );
337    }
338    g = mapinto( g1 ); h = mapinto( h1 );
339    g = replaceLc( g, lf ); h = replaceLc( h, lf );
340    e = f - g * h;
341    modulus = p;
342    i = 1;
343
344    while ( ! e.isZero() && i <= k ) {
345        c = e / modulus;
346        {
347            setCharacteristic( p );
348            c = mapinto( c );
349            divrem( a * c, h1, q, r );
350            tmp = b * c + q * g1;
351            setCharacteristic( 0 );
352        }
353        g += mapinto( tmp ) * modulus;
354        h += mapinto( r ) * modulus;
355//        g = replaceLc( g, lf ); h = replaceLc( h, lf );
356        e = f - g * h;
357        modulus *= p;
358        ASSERT( mod( f - g * h, modulus ) == 0, "error at lift stage" );
359        i++;
360    }
361    if ( e.isZero() ) {
362        tmp = content( g );
363        gk = g / tmp; hk = h / ( lf / tmp );
364    }
365    else {
366        gk = pk(g); hk = pk(h);
367    }
368    setCharacteristic( save );
369//    return e.isZero();
370    return (F-gk*hk).isZero();
371}
372
373
374CFFList
375ZFactorizeUnivariate( const CanonicalForm& ff, bool issqrfree )
376{
377    bool symmsave = isOn( SW_SYMMETRIC_FF );
378    CanonicalForm cont = content( ff );
379    CanonicalForm lf, recip_lf, fp, f, g = ff / cont, dummy1, dummy2;
380    int i, k, exp, n;
381    bool ok;
382    CFFList H, F[max_fp_fac];
383    CFFList ZF;
384    int * p = new int [max_fp_fac];
385    int * D = 0;
386    int * Dh = 0;
387    ListIterator<CFFactor> J, I;
388
389    DEBINCLEVEL( cerr, "ZFactorizeUnivariate" );
390    On( SW_SYMMETRIC_FF );
391
392    // get squarefree decomposition of f
393    if ( issqrfree )
394        H.append( CFFactor( g, 1 ) );
395    else
396        H = sqrFree( g );
397
398    DEBOUTLN( cerr, "H = " << H );
399
400    // cycle through squarefree factors of f
401    for ( J = H; J.hasItem(); ++J ) {
402        f = J.getItem().factor();
403        if ( f.inCoeffDomain() ) continue;
404        n = f.degree() / 2 + 1;
405        delete [] D;
406        delete [] Dh;
407        D = new int [n]; D[0] = n;
408        Dh = new int [n]; Dh[0] = n;
409        exp = J.getItem().exp();
410
411        // choose primes to factor f
412        TIMING_START(fac_choosePrimes);
413        ok = choosePrimes( p, f );
414        TIMING_END_AND_PRINT(fac_choosePrimes, "time to choose the primes: ");
415        if ( ! ok ) {
416            DEBOUTLN( cerr, "warning: no good prime found to factorize " << f );
417            STICKYWARN( ok, "there occured an error.  We went out of primes p\n"
418                        "to factorize mod p.  Please send the example which caused\n"
419                        "this error to the authors.  Nonetheless we will go on with the\n"
420                        "calculations hoping the result will be correct.  Thank you.");
421            ZF.append( CFFactor( f, exp ) );
422            continue;
423        }
424
425        // factorize f modulo certain primes
426        TIMING_START(fac_facModPrimes);
427        for ( i = 0; i < max_fp_fac; i++ ) {
428            setCharacteristic( p[i] );
429            fp = mapinto( f );
430            F[i] = FpFactorizeUnivariateCZ( fp, true, 0, Variable(), Variable() );
431//                if ( p[i] < 23 && fp.degree() < 10 )
432//                    F[i] = FpFactorizeUnivariateB( fp, true );
433//                else
434//                    F[i] = FpFactorizeUnivariateCZ( fp, true, 0, Variable, Variable() );
435            DEBOUTLN( cerr, "F[i] = " << F[i] << ", p = " << p[i] );
436        }
437        TIMING_END_AND_PRINT(fac_facModPrimes, "time to factorize mod primes: ");
438        setCharacteristic( 0 );
439
440        // do some strange things with the D's
441        initHG( D, F[0] );
442        hgroup( D );
443        DEBOUTHPRINT( cerr, "D = ", D );
444        for ( i = 1; i < max_fp_fac; i++ ) {
445            initHG( Dh, F[i] );
446            hgroup( Dh );
447            DEBOUTHPRINT( cerr, "Dh = ", Dh );
448            hintersect( D, Dh );
449            DEBOUTHPRINT( cerr, "D = ", D );
450        }
451
452        // look which p gives the shortest factorization of f mod p
453        // j: index of that p in p[]
454        int min, j;
455        min = F[0].length(), j = 0;
456        for ( i = 1; i < max_fp_fac; i++ ) {
457            if ( min >= F[i].length() ) {
458                j = i; min = F[i].length();
459            }
460        }
461        k = kBound( f, p[j] );
462        CFArray theFactors( F[j].length() );
463//            pk = power( CanonicalForm( p[j] ), k );
464//            pkhalf = pk / 2;
465        modpk pk( p[j], k );
466        DEBOUTLN( cerr, "coeff bound = " << pk.getpk() );
467        theModulus = pk;
468        setCharacteristic( p[j] );
469        fp = mapinto( f );
470        F[j].sort( cmpFactor );
471        I = F[j]; i = 0;
472        TIMING_START(fac_liftFactors);
473        while ( I.hasItem() ) {
474            DEBOUTLN( cerr, "factor to lift = " << I.getItem().factor() );
475            if ( isOn( SW_FAC_QUADRATICLIFT ) )
476                ok = UnivariateQuadraticLift( f, I.getItem().factor(), fp / I.getItem().factor(), pk, lc( f ), dummy1, dummy2 );
477            else
478                ok = UnivariateLinearLift( f, I.getItem().factor(), fp / I.getItem().factor(), pk, lc( f ), dummy1, dummy2 );
479            if ( ok ) {
480                // should be done in a more efficient way
481                DEBOUTLN( cerr, "dummy1 = " << dummy1 );
482                DEBOUTLN( cerr, "dummy2 = " << dummy2 );
483                f = dummy2;
484                fp /= I.getItem().factor();
485                ZF.append( CFFactor( dummy1, exp ) );
486                I.remove( 0 );
487                I = F[j];
488                i = 0;
489                DEBOUTLN( cerr, "F[j] = " << F[j] );
490            }
491            else {
492                DEBOUTLN( cerr, "i = " << i );
493                DEBOUTLN( cerr, "dummy1 = " << dummy1 );
494                setCharacteristic( 0 );
495//                    theFactors[i] = pk( dummy1 * pk.inverse( lc( dummy1 ) ) );
496                theFactors[i] = pk( dummy1 );
497                setCharacteristic( p[j] );
498                i++;
499                I++;
500            }
501        }
502        TIMING_END_AND_PRINT(fac_liftFactors, "time to lift the factors: ");
503        DEBOUTLN( cerr, "ZF = " << ZF );
504        initHG( Dh, theFactors );
505        hgroup( Dh );
506        DEBOUTHPRINT( cerr, "Dh = ", Dh );
507        hintersect( D, Dh );
508        setCharacteristic( 0 );
509        for ( int l = i; l < F[j].length(); l++ )
510            theFactors[l] = 1;
511        DEBOUTLN( cerr, "theFactors = " << theFactors );
512        DEBOUTLN( cerr, "f = " << f );
513        DEBOUTLN( cerr, "p = " << pk.getp() << ", k = " << pk.getk() );
514        DEBOUTHPRINT( cerr, "D = ", D );
515        lf = lc( f );
516        (void)bextgcd( pk.getpk(), lf, dummy1, recip_lf );
517        DEBOUTLN( cerr, "recip_lf = " << recip_lf );
518        TIMING_START(fac_combineFactors);
519        for ( i = 1; i < D[0]; i++ )
520            if ( D[i] != 0 )
521                while ( liftDegreeFactRec( theFactors, f, recip_lf, lf, pk, 0, i, ZF, exp ) );
522        if ( degree( f ) > 0 )
523            ZF.append( CFFactor( f, exp ) );
524        TIMING_END_AND_PRINT(fac_combineFactors, "time to combine the factors: ");
525    }
526
527    // brush up our result
528    if ( ZF.getFirst().factor().inCoeffDomain() )
529        ZF.removeFirst();
530    if ( lc( ff ).sign() < 0 )
531        ZF.insert( CFFactor( -cont, 1 ) );
532    else
533        ZF.insert( CFFactor( cont, 1 ) );
534    delete [] D;
535    delete [] Dh;
536    if ( ! symmsave )
537        Off( SW_SYMMETRIC_FF );
538
539    DEBDECLEVEL( cerr, "ZFactorizeUnivariate" );
540    return ZF;
541}
Note: See TracBrowser for help on using the repository browser.