source: git/factory/fac_univar.cc @ b1dfaf

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