source: git/factory/fac_univar.cc @ af4c56

spielwiese
Last change on this file since af4c56 was af4c56, checked in by Rüdiger Stobbe <stobbe@…>, 27 years ago
"ZFactorizeUnivariate: now handles the sign of the argument right. " git-svn-id: file:///usr/local/Singular/svn/trunk@30 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 13.7 KB
Line 
1// emacs edit mode for this file is -*- C++ -*-
2// $Id: fac_univar.cc,v 1.3 1996-06-26 13:17:03 stobbe Exp $
3
4/*
5$Log: not supported by cvs2svn $
6Revision 1.2  1996/06/13 10:43:49  stobbe
7"ZFactorizeUnivariate: fix to last bug fix (no assignment)
8"
9
10Revision 1.1  1996/06/13 10:34:04  stobbe
11"ZFactorizeUnivariate: do not use Berlekamp-Algorithm since there is a
12                      bug in the Factory-Implementation of Berlekamp
13"
14
15Revision 1.0  1996/05/17 10:59:45  stobbe
16Initial revision
17
18*/
19
20//#define TIMING
21
22#ifndef NDEBUG
23//#define DEBUGOUTPUT
24#endif
25
26#ifdef TIMING
27#include <sys/times.h>
28#endif
29
30#include <math.h>
31#include "assert.h"
32#include "cf_defs.h"
33#include "fac_util.h"
34#include "fac_univar.h"
35#include "fac_cantzass.h"
36#include "fac_berlekamp.h"
37#include "cf_iter.h"
38#include "cf_primes.h"
39#include "fac_sqrfree.h"
40
41#define MAX_FP_FAC 3
42
43static modpk theModulus;
44
45CanonicalForm
46iextgcd ( const CanonicalForm & f, const CanonicalForm & g, CanonicalForm & a, CanonicalForm & b );
47
48#ifndef NDEBUG
49
50static void
51hprint ( int * a )
52{
53    int n = a[0];
54    cerr << "( " << n << ": ";
55    int i = 1;
56    while ( i < n ) {
57        if ( a[i] != 0 )
58            cerr << i << " ";
59        i++;
60    }
61    cerr << ")" << endl;
62}
63
64#endif
65
66static void
67hgroup ( int * a )
68{
69    int n = a[0];
70    int i, j, k;
71    for ( i = 1; i < n; i++ )
72        if ( a[i] != 0 )
73            for ( j = 1; j <= i; j++ )
74                if ( a[j] != 0 )
75                    for ( k = i; k < n; k += j )
76                        a[k] = 1;
77}
78
79static void
80hintersect( int * a, const int * const b )
81{
82    int i, n, na = a[0], nb = b[0];
83    if ( nb < na )
84        n = nb;
85    else
86        n = na;
87    for ( i = 1; i < n; i++ )
88        if ( b[i] == 0 )
89            a[i] = 0;
90    a[0] = n;
91}
92
93/*
94static int
95hcount ( const int * const a )
96{
97    int n = a[0], sum = 0, i;
98    for ( i = 1; i < n; i++ )
99        if ( a[i] != 0 ) sum++;
100    return sum;
101}
102*/
103
104static void
105initHG ( int * a, const CFFList & F )
106{
107    ListIterator<CFFactor> i;
108
109    int n = a[0], k;
110    for ( int j = 1; j < n; j++ ) a[j] = 0;
111    for ( i = F; i.hasItem(); i++ )
112        if ( (k = i.getItem().factor().degree()) < n )
113            if ( k != 0 )
114                a[k] = 1;
115}
116
117static void
118initHG ( int * a, const Array<CanonicalForm> & F )
119{
120    int i, n = a[0], m = F.size(), k;
121    for ( i = 1; i < n; i++ ) a[i] = 0;
122    for ( i = 1; i < m; i++ )
123        if ( (k = F[i].degree()) < n )
124            if ( k != 0 )
125                a[k] = 1;
126}
127
128static int
129cmpFactor( const CFFactor & a, const CFFactor & b )
130{
131    CFFactor A( a ), B( b );
132    return degree( A.factor() ) > degree( B.factor() );
133}
134
135static double
136cf2double ( const CanonicalForm & f )
137{
138    CanonicalForm a = f, q, r;
139    double m = 1, res = 0;
140    if ( a.sign() < 0 ) a = -a;
141    while ( ! a.isZero() ) {
142        divrem( a, 10, q, r );
143        res += m * (double)(r.intval());
144        m *= 10;
145        a = q;
146    }
147    if ( f.sign() < 0 ) res = -res;
148    return res;
149}
150
151static double
152dnorm ( const CanonicalForm & f )
153{
154    CFIterator i;
155    CanonicalForm sum = 0;
156    for ( i = f; i.hasTerms(); i++ ) sum += i.coeff() * i.coeff();
157    DEBOUTLN( cerr, "sum = ", sum );
158    return sqrt( cf2double( sum ) );
159}
160
161static int
162kBound ( const CanonicalForm & f, int p )
163{
164    DEBOUTLN( cerr, "lc(f) = ", lc(f) );
165    return (int)((f.degree()*log(2)+log( fabs(cf2double(lc(f))) )+log( dnorm( f ) )) / log( p ) + 0.5) + 1;
166}
167
168modpk
169getZFacModulus()
170{
171    return theModulus;
172}
173
174static bool
175liftDegreeFactRec( CFArray & theFactors, CanonicalForm & F, const CanonicalForm & recip_lf, const CanonicalForm & f, const modpk & pk, int i, int d, CFFList & ZF, int exp )
176{
177    if ( i >= theFactors.size() )
178        return false;
179    else  if ( degree( f ) + degree( theFactors[i] ) == d ) {
180        DEBOUTLN( cerr, "ldfr (f) = ", f );
181        DEBOUTLN( cerr, "ldfr (g) = ", theFactors[i] );
182        CanonicalForm g = pp( pk( recip_lf * f * theFactors[i] ) );
183        DEBOUTLN( cerr, "ldfr (pk(f*g)) = ", g );
184        CanonicalForm gg, hh;
185        DEBOUTLN( cerr, "F = ", F );
186        DEBOUTLN( cerr, "g = ", g );
187        if ( divremt( F, g, gg, hh ) && hh.isZero() ) {
188            ZF.append( CFFactor( g, exp ) );
189            F = gg;
190            theFactors[i] = 1;
191            return true;
192        }
193        else {
194            return liftDegreeFactRec( theFactors, F, recip_lf, f, pk, i+1, d, ZF, exp );
195        }
196    }
197    else  if ( degree( f ) + degree( theFactors[i] ) > d )
198        return false;
199    else {
200        bool ok = liftDegreeFactRec( theFactors, F, recip_lf, pk( recip_lf * f * theFactors[i] ), pk, i+1, d, ZF, exp );
201        if ( ok )
202            theFactors[i] = 1;
203        else
204            ok = liftDegreeFactRec( theFactors, F, recip_lf, f, pk, i+1, d, ZF, exp );
205        return ok;
206    }
207}
208
209static int
210choosePrimes ( int * p, const CanonicalForm & f )
211{
212    int ptr = 0;
213    int i = 0;
214    int maxp = cf_getNumPrimes();
215    int prime;
216
217    while ( ptr < maxp && i < MAX_FP_FAC ) {
218        prime = cf_getPrime( ptr );
219        if ( mod( lc( f ), prime ) != 0 ) {
220            setCharacteristic( prime );
221            if ( isSqrFreeFp( mapinto( f ) ) ) {
222                p[i] = prime;
223                i++;
224            }
225            setCharacteristic( 0 );
226        }
227        ptr++;
228    }
229    return ( i == MAX_FP_FAC );
230}
231
232
233static int
234UnivariateQuadraticLift ( const CanonicalForm &F, const  CanonicalForm & G, const CanonicalForm &H, const modpk & pk, const CanonicalForm & Gamma, CanonicalForm & gk, CanonicalForm & hk )
235{
236    CanonicalForm lf, f, gamma;
237    CanonicalForm a, b, aa, bb, c, g, h, g1, h1, e, modulus, tmp, q, r;
238    int i, j, save;
239    int p = pk.getp(), k = pk.getk();
240    int no_iter = (int)(log(k)/log(2)+2);
241    int * kvals = new int[no_iter];
242
243    DEBOUT( cerr, "quadratic lift called with p = ", p );
244    DEBOUTLN( cerr, "  and k = ", k );
245    for ( j = 0, i = k; i > 1; i = ( i+1 ) / 2, j++ ) kvals[j] = i;
246    kvals[j] = 1;
247
248    save = getCharacteristic(); setCharacteristic( 0 );
249
250    lf = lc( F );
251    f = lf * F;
252    {
253        setCharacteristic( p );
254        g1 = mapinto( lf ) / lc( G ) * G;
255        h1 = mapinto( lf ) / lc( H ) * H;
256        (void)extgcd( g1, h1, a, b );
257        setCharacteristic( 0 );
258    }
259    a = mapinto( a ); b = mapinto( b );
260    g = mapinto( g1 ); h = mapinto( h1 );
261    g = replaceLc( g, lf ); h = replaceLc( h, lf );
262    e = f - g * h;
263    modulus = p;
264    i = 1;
265
266    while ( ! e.isZero() && j > 0 ) {
267        c = div( e, modulus );
268        {
269            j--;
270            setCharacteristic( p, kvals[j+1] );
271            DEBOUT( cerr, "lifting from p^", kvals[j+1] );
272            DEBOUTLN( cerr, " to p^", kvals[j] );
273            c = mapinto( c );
274            DEBOUTLN( cerr, " !!! g = ", mapinto( g ) );
275            g1 = mapinto( lf ) / mapinto( lc( g ) ) * mapinto( g );
276            h1 = mapinto( lf ) / mapinto( lc( h ) ) * mapinto( h );
277//          (void)extgcd( g1, h1, a, b );
278//          DEBOUTLN( cerr, " a = ", aa );
279//          DEBOUTLN( cerr, " b = ", bb );
280            a = mapinto( a ); b = mapinto( b );
281            a += ( ( 1 - a * g1 ) *  a ) % h1;
282            b += ( ( 1 - b * h1 ) *  b ) % g1;
283            DEBOUTLN( cerr, " a = ", a );
284            DEBOUTLN( cerr, " b = ", b );
285            divrem( a * c, h1, q, r );
286            tmp = b * c + q * g1;
287            setCharacteristic( 0 );
288        }
289        a = mapinto( a ); b = mapinto( b );
290        g += mapinto( tmp ) * modulus;
291        h += mapinto( r ) * modulus;
292//      g = replaceLc( g, lf ); h = replaceLc( h, lf );
293        e = f - g * h;
294        modulus = power( CanonicalForm(p), kvals[j] );
295        if ( mod( f - g * h, modulus ) != 0 )
296            DEBOUTLN( cerr, "error at lift stage ", i );
297        i++;
298    }
299    if ( e.isZero() ) {
300        tmp = content( g );
301        gk = g / tmp; hk = h / ( lf / tmp );
302    }
303    else {
304        gk = pk(g); hk = pk(h);
305    }
306    setCharacteristic( save );
307    return e.isZero();
308}
309
310static int
311UnivariateLinearLift ( const CanonicalForm &F, const  CanonicalForm & G, const CanonicalForm &H, const modpk & pk, const CanonicalForm & Gamma, CanonicalForm & gk, CanonicalForm & hk )
312{
313    CanonicalForm lf, f, gamma;
314    CanonicalForm a, b, c, g, h, g1, h1, e, modulus, tmp, q, r;
315    int i, save;
316    int p = pk.getp(), k = pk.getk();
317    save = getCharacteristic(); setCharacteristic( 0 );
318
319    lf = lc( F );
320    f = lf * F;
321    {
322        setCharacteristic( p );
323        g1 = mapinto( lf ) / lc( G ) * G;
324        h1 = mapinto( lf ) / lc( H ) * H;
325        (void)extgcd( g1, h1, a, b );
326        setCharacteristic( 0 );
327    }
328    g = mapinto( g1 ); h = mapinto( h1 );
329    g = replaceLc( g, lf ); h = replaceLc( h, lf );
330    e = f - g * h;
331    modulus = p;
332    i = 1;
333
334    while ( ! e.isZero() && i <= k ) {
335        c = div( e, modulus );
336        {
337            setCharacteristic( p );
338            c = mapinto( c );
339            divrem( a * c, h1, q, r );
340            tmp = b * c + q * g1;
341            setCharacteristic( 0 );
342        }
343        g += mapinto( tmp ) * modulus;
344        h += mapinto( r ) * modulus;
345//      g = replaceLc( g, lf ); h = replaceLc( h, lf );
346        e = f - g * h;
347        modulus *= p;
348        ASSERT( mod( f - g * h, modulus ) == 0, "error at lift stage" );
349        i++;
350    }
351    if ( e.isZero() ) {
352        tmp = content( g );
353        gk = g / tmp; hk = h / ( lf / tmp );
354    }
355    else {
356        gk = pk(g); hk = pk(h);
357    }
358    setCharacteristic( save );
359//    return e.isZero();
360    return (F-gk*hk).isZero();
361}
362
363CFFList
364ZFactorizeUnivariate( const CanonicalForm& ff, bool issqrfree )
365{
366    bool symmsave = isOn( SW_SYMMETRIC_FF );
367    CanonicalForm cont = content( ff );
368    CanonicalForm lf, recip_lf, fp, f, g = ff / cont, dummy1, dummy2;
369    int i, k, exp, n;
370    bool ok;
371    CFFList H, G, F[MAX_FP_FAC];
372    CFFList ZF;
373    int * p = new int [MAX_FP_FAC];
374    int * D = 0;
375    int * Dh = 0;
376    ListIterator<CFFactor> J, I;
377    On( SW_SYMMETRIC_FF );
378    if ( issqrfree )
379        H.append( CFFactor( g, 1 ) );
380    else
381        H = sqrFree( g );
382    for ( J = H; J.hasItem(); ++J ) {
383        f = J.getItem().factor();
384        n = f.degree() / 2 + 1;
385        if ( D != 0 ) {
386            delete [] D;
387            delete [] Dh;
388        }
389        D = new int [n]; D[0] = n;
390        Dh = new int [n]; Dh[0] = n;
391        exp = J.getItem().exp();
392#ifdef TIMING
393        struct tms ts, te;
394        times( &ts );
395#endif
396        ok = choosePrimes( p, f );
397#ifdef TIMING
398        times( &te );
399        cout.setf( ios::fixed, ios::floatfield );
400        cout.precision(2);
401        cout << "time to choose the primes: "
402             << float(te.tms_utime-ts.tms_utime) / CLK_TCK << " sec" << endl;
403#endif
404        if ( ! ok ) {
405            cerr << "factorize warnig: no good prime found" << endl;
406            cerr << f << endl;
407            ZF.append( CFFactor( f, exp ) );
408        }
409        else {
410#ifdef TIMING
411            times( &ts );
412#endif
413            for ( i = 0; i < MAX_FP_FAC; i++ ) {
414                setCharacteristic( p[i] );
415                fp = mapinto( f );
416                F[i] = FpFactorizeUnivariateCZ( fp, true );
417//              if ( p[i] < 23 && fp.degree() < 10 )
418//                  F[i] = FpFactorizeUnivariateB( fp, true );
419//              else
420//                  F[i] = FpFactorizeUnivariateCZ( fp, true );
421                DEBOUT( cerr, "F[i] = ", F[i] );
422                DEBOUTLN( cerr, ", p = ", p[i] );
423            }
424#ifdef TIMING
425        times( &te );
426        cout.setf( ios::fixed, ios::floatfield );
427        cout.precision(2);
428        cout << "time to factorize mod primes: "
429             << float(te.tms_utime-ts.tms_utime) / CLK_TCK << " sec" << endl;
430#endif
431            setCharacteristic( 0 );
432#ifdef DEBUGOUTPUT
433            hprint( D );
434#endif
435            initHG( D, F[0] );
436            hgroup( D );
437#ifdef DEBUGOUTPUT
438            hprint( D );
439#endif
440            for ( i = 1; i < MAX_FP_FAC; i++ ) {
441                initHG( Dh, F[i] );
442                hgroup( Dh );
443#ifdef DEBUGOUTPUT
444                cerr << "Dh = ";
445                hprint( Dh );
446#endif
447                hintersect( D, Dh );
448#ifdef DEBUGOUTPUT
449                cerr << "D = ";
450                hprint( D );
451#endif
452
453            }
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#ifdef TIMING
473            times( &ts );
474#endif
475            while ( I.hasItem() ) {
476                DEBOUTLN( cerr, "factor to lift = ", I.getItem().factor() );
477                if ( isOn( SW_FAC_QUADRATICLIFT ) )
478                    ok = UnivariateQuadraticLift( f, I.getItem().factor(), fp / I.getItem().factor(), pk, lc( f ), dummy1, dummy2 );
479                else
480                    ok = UnivariateLinearLift( f, I.getItem().factor(), fp / I.getItem().factor(), pk, lc( f ), dummy1, dummy2 );
481                if ( ok ) {
482                    // should be done in a more efficient way
483                    DEBOUTLN( cerr, "dummy1 = ", dummy1 );
484                    DEBOUTLN( cerr, "dummy2 = ", dummy2 );
485                    f = dummy2;
486                    fp /= I.getItem().factor();
487                    ZF.append( CFFactor( dummy1, exp ) );
488                    I.remove( 0 ); 
489                    I = F[j];
490                    i = 0;
491                    DEBOUTLN( cerr, "F[j] = ", F[j] );
492                }
493                else {
494                    DEBOUTLN( cerr, "i = ", i );
495                    DEBOUTLN( cerr, "dummy1 = ", dummy1 );
496                    setCharacteristic( 0 );
497//                  theFactors[i] = pk( dummy1 * pk.inverse( lc( dummy1 ) ) );
498                    theFactors[i] = pk( dummy1 );
499                    setCharacteristic( p[j] );
500                    i++;
501                    I++;
502                }
503            }
504#ifdef TIMING
505        times( &te );
506        cout.setf( ios::fixed, ios::floatfield );
507        cout.precision(2);
508        cout << "time to lift the factors: "
509             << float(te.tms_utime-ts.tms_utime) / CLK_TCK << " sec" << endl;
510#endif
511            DEBOUTLN( cerr, "ZF = ", ZF );
512            initHG( Dh, theFactors );
513            hgroup( Dh );
514#ifdef DEBUGOUTPUT
515            cerr << "Dh = ";
516            hprint( Dh );
517#endif
518            hintersect( D, Dh );
519            setCharacteristic( 0 );
520            for ( int l = i; l < F[j].length(); l++ )
521                theFactors[l] = 1;
522            DEBOUTLN( cerr, "theFactors = ", theFactors );
523            DEBOUTLN( cerr, "f = ", f );
524            DEBOUT( cerr, "p = ", pk.getp() );
525            DEBOUTLN( cerr, ", k = ", pk.getk() );
526#ifdef DEBUGOUTPUT
527            hprint( D );
528#endif
529            lf = lc( f );
530            (void)iextgcd( pk.getpk(), lf, dummy1, recip_lf );
531            DEBOUTLN( cerr, "recip_lf = ", recip_lf );
532#ifdef TIMING
533            times( &ts );
534#endif
535            for ( i = 1; i < D[0]; i++ )
536                if ( D[i] != 0 )
537                    while ( liftDegreeFactRec( theFactors, f, recip_lf, lf, pk, 0, i, ZF, exp ) );
538            if ( degree( f ) > 0 )
539                ZF.append( CFFactor( f, exp ) );
540#ifdef TIMING
541        times( &te );
542        cout.setf( ios::fixed, ios::floatfield );
543        cout.precision(2);
544        cout << "time to combine the factors: "
545             << float(te.tms_utime-ts.tms_utime) / CLK_TCK << " sec" << endl;
546#endif
547        }
548    }
549    if ( ZF.getFirst().factor().inCoeffDomain() )
550        ZF.removeFirst();
551    if ( lc( ff ).sign() < 0 )
552        ZF.insert( CFFactor( -cont, 1 ) );
553    else  if ( ! cont.isOne() )
554        ZF.insert( CFFactor( cont, 1 ) );
555    if ( D != 0 ) {
556        delete [] D;
557        delete [] Dh;
558    }
559    if ( ! symmsave )
560        Off( SW_SYMMETRIC_FF );
561    return ZF;
562}
Note: See TracBrowser for help on using the repository browser.