source: git/factory/fac_univar.cc @ 80532d7

spielwiese
Last change on this file since 80532d7 was 80532d7, checked in by Jens Schmidt <schmidt@…>, 26 years ago
* fac_univar.cc (norm): function removed. All references replaced by `euclideanNorm()'. git-svn-id: file:///usr/local/Singular/svn/trunk@1226 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: fac_univar.cc,v 1.17 1998-03-12 14:32:54 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#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
128cmpFactor( const CFFactor & a, const CFFactor & b )
129{
130    CFFactor A( a ), B( b );
131    return degree( A.factor() ) > degree( B.factor() );
132}
133
134static double
135cf2double ( const CanonicalForm & f )
136{
137    CanonicalForm a = f, q, r;
138    double m = 1, res = 0;
139    if ( a.sign() < 0 ) a = -a;
140    while ( ! a.isZero() ) {
141        divrem( a, 10, q, r );
142        res += m * (double)(r.intval());
143        m *= 10;
144        a = q;
145    }
146    if ( f.sign() < 0 ) res = -res;
147    return res;
148}
149
150//{{{ static int kBound ( const CanonicalForm & f, int p )
151//{{{ docu
152//
153// kBound() - return bound of coefficients of factors of f.
154//
155// The bound is returned as an integer k such that p^k is larger
156// than all coefficients of all possible factors of f.  f should
157// be an univariate polynomial over Z.
158//
159// For a discussion of the formula, see the article Mignotte -
160// 'Some Usefull Bounds' in Buchberger, Collins, Loos (eds.) -
161// 'Computer Algebra: Symbolic and Algebraic Computation', 2nd
162// ed.
163//
164// Use by ZFactorizeUnivariate().
165//
166//}}}
167static int
168kBound ( const CanonicalForm & f, int p )
169{
170    return (int)(f.degree() + (double)(ilog2( euclideanNorm(f)+1 ) + 1) / (double)ilog2(p)) + 1;
171}
172//}}}
173
174modpk
175getZFacModulus()
176{
177    return theModulus;
178}
179
180static bool
181liftDegreeFactRec( CFArray & theFactors, CanonicalForm & F, const CanonicalForm & recip_lf, const CanonicalForm & f, const modpk & pk, int i, int d, CFFList & ZF, int exp )
182{
183    if ( i >= theFactors.size() )
184        return false;
185    else  if ( degree( f ) + degree( theFactors[i] ) == d ) {
186        DEBOUTLN( cerr, "ldfr (f) = " << f );
187        DEBOUTLN( cerr, "ldfr (g) = " << theFactors[i] );
188        CanonicalForm g = pp( pk( recip_lf * f * theFactors[i] ) );
189        DEBOUTLN( cerr, "ldfr (pk(f*g)) = " << g );
190        CanonicalForm gg, hh;
191        DEBOUTLN( cerr, "F = " << F );
192        DEBOUTLN( cerr, "g = " << g );
193        if ( divremt( F, g, gg, hh ) && hh.isZero() ) {
194            ZF.append( CFFactor( g, exp ) );
195            F = gg;
196            theFactors[i] = 1;
197            return true;
198        }
199        else {
200            return liftDegreeFactRec( theFactors, F, recip_lf, f, pk, i+1, d, ZF, exp );
201        }
202    }
203    else  if ( degree( f ) + degree( theFactors[i] ) > d )
204        return false;
205    else {
206        bool ok = liftDegreeFactRec( theFactors, F, recip_lf, pk( recip_lf * f * theFactors[i] ), pk, i+1, d, ZF, exp );
207        if ( ok )
208            theFactors[i] = 1;
209        else
210            ok = liftDegreeFactRec( theFactors, F, recip_lf, f, pk, i+1, d, ZF, exp );
211        return ok;
212    }
213}
214
215
216static int
217choosePrimes ( int * p, const CanonicalForm & f )
218{
219    int ptr = 0;
220    int i = 0;
221    int maxp = cf_getNumPrimes();
222    int prime;
223
224    while ( ptr < maxp && i < max_fp_fac ) {
225        prime = cf_getPrime( ptr );
226        if ( mod( lc( f ), prime ) != 0 ) {
227            setCharacteristic( prime );
228            if ( isSqrFreeFp( mapinto( f ) ) ) {
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)+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 = div( 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 = div( 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.