Changeset a2c0302 in git


Ignore:
Timestamp:
Apr 22, 1997, 5:40:32 PM (27 years ago)
Author:
Jens Schmidt <schmidt@…>
Branches:
(u'spielwiese', '5b153614cbc72bfa198d75b1e9e33dab2645d9fe')
Children:
56f7c59395cec760c41b7311172e225a5cc099eb
Parents:
eef3aacee99ae0cf329b82837ba229a249b16451
Message:
o some spurious preprocessr directives removed
o #define MAX_FP_FAC changed to a const max_fp_fac
o calls to hprint() changed to calls of a macro named
  DEBOUTHPRINT()
o initHG(): warning added because of some strange error in
  the internals of factory.  we will wait for some example
  to search for the error.
o kBound(), UnivariatequadraticLift(): added cast to double
  in the calculation of the bound.  at least one compiler
  subjects an ambiguity at this point.
o ZFactorizeUnivariate():
  - spurious variable CFFList G removed
  - sequence 'if ( D != 0 ) delete D;' changed to 'delete
    D;'
  - warning added because of some strange error in the
    internals of factory.  we will wait for some example
    to search for the error.
o calls to macro DEBOUTLN changed to new calling syntax


git-svn-id: file:///usr/local/Singular/svn/trunk@186 2c84dea3-7e68-4137-9b89-c4e89433aadc
File:
1 edited

Legend:

Unmodified
Added
Removed
  • factory/fac_univar.cc

    reef3aa ra2c0302  
    11// emacs edit mode for this file is -*- C++ -*-
    2 // $Id: fac_univar.cc,v 1.6 1997-04-08 10:33:19 schmidt Exp $
     2// $Id: fac_univar.cc,v 1.7 1997-04-22 15:40:32 schmidt Exp $
    33
    44/*
    55$Log: not supported by cvs2svn $
     6Revision 1.6  1997/04/08 10:33:19  schmidt
     7#include <config.h> added
     8
    69Revision 1.5  1997/03/27 09:54:41  schmidt
    710timing output changed to TIMING
     
    3033
    3134*/
     35
    3236
    3337#include <config.h>
     
    5357TIMING_DEFINE_PRINT(fac_combineFactors);
    5458
    55 #define MAX_FP_FAC 3
     59
     60const int max_fp_fac = 3;
    5661
    5762static modpk theModulus;
    5863
     64// !!! this should be placed in cf_gcd.h
    5965CanonicalForm
    6066iextgcd ( const CanonicalForm & f, const CanonicalForm & g, CanonicalForm & a, CanonicalForm & b );
    6167
     68
    6269#ifdef DEBUGOUTPUT
     70#define DEBOUTHPRINT(stream, msg, hg) \
     71{stream << deb_level_msg << msg, stream.flush(); hprint( hg ); stream << endl;}
     72
    6373static void
    6474hprint ( int * a )
     
    7282        i++;
    7383    }
    74     cerr << ")" << endl << endl;
    75 }
    76 #endif
     84    cerr << ")";
     85}
     86#else /* DEBUGOUTPUT */
     87#define DEBOUTHPRINT(stream, msg, hg)
     88#endif /* DEBUGOUTPUT */
    7789
    7890static void
     
    123135    for ( i = F; i.hasItem(); i++ )
    124136        if ( (k = i.getItem().factor().degree()) < n )
    125             if ( k != 0 )
     137            if ( k == -1 ) {
     138                WARN( k == -1, "there occured an error.  factory was not able to factorize\n"
     139                      "correctly mod p.  Please send the example which caused\n"
     140                      "this error to the authors.  Nonetheless we will go on with the\n"
     141                      "calculations hoping the result will be correct.  Thank you." );
     142            }
     143            else if ( k != 0 )
    126144                a[k] = 1;
    127145}
     
    134152    for ( i = 1; i < m; i++ )
    135153        if ( (k = F[i].degree()) < n )
    136             if ( k != 0 )
     154            if ( k == -1 ) {
     155                WARN( k == -1, "there occured an error.  factory was not able to factorize\n"
     156                      "correctly mod p.  Please send the example which caused\n"
     157                      "this error to the authors.  Nonetheless we will go on with the\n"
     158                      "calculations hoping the result will be correct.  Thank you." );
     159            }
     160            else if ( k != 0 )
    137161                a[k] = 1;
    138162}
     
    167191    CanonicalForm sum = 0;
    168192    for ( i = f; i.hasTerms(); i++ ) sum += i.coeff() * i.coeff();
    169     DEBOUTLN( cerr, "sum = ", sum );
     193    DEBOUTLN( cerr, "sum = " << sum );
    170194    return sqrt( cf2double( sum ) );
    171195}
     
    174198kBound ( const CanonicalForm & f, int p )
    175199{
    176     DEBOUTLN( cerr, "lc(f) = ", lc(f) );
    177     return (int)((f.degree()*log(2)+log( fabs(cf2double(lc(f))) )+log( dnorm( f ) )) / log( p ) + 0.5) + 1;
     200    DEBOUTLN( cerr, "lc(f) = " << lc(f) );
     201    return (int)((f.degree()*log(2)+log( fabs(cf2double(lc(f))) )+log( dnorm( f ) )) / log( (double)p ) + 0.5) + 1;
    178202}
    179203
     
    190214        return false;
    191215    else  if ( degree( f ) + degree( theFactors[i] ) == d ) {
    192         DEBOUTLN( cerr, "ldfr (f) = ", f );
    193         DEBOUTLN( cerr, "ldfr (g) = ", theFactors[i] );
     216        DEBOUTLN( cerr, "ldfr (f) = " << f );
     217        DEBOUTLN( cerr, "ldfr (g) = " << theFactors[i] );
    194218        CanonicalForm g = pp( pk( recip_lf * f * theFactors[i] ) );
    195         DEBOUTLN( cerr, "ldfr (pk(f*g)) = ", g );
     219        DEBOUTLN( cerr, "ldfr (pk(f*g)) = " << g );
    196220        CanonicalForm gg, hh;
    197         DEBOUTLN( cerr, "F = ", F );
    198         DEBOUTLN( cerr, "g = ", g );
     221        DEBOUTLN( cerr, "F = " << F );
     222        DEBOUTLN( cerr, "g = " << g );
    199223        if ( divremt( F, g, gg, hh ) && hh.isZero() ) {
    200224            ZF.append( CFFactor( g, exp ) );
     
    219243}
    220244
     245
    221246static int
    222247choosePrimes ( int * p, const CanonicalForm & f )
     
    227252    int prime;
    228253
    229     while ( ptr < maxp && i < MAX_FP_FAC ) {
     254    while ( ptr < maxp && i < max_fp_fac ) {
    230255        prime = cf_getPrime( ptr );
    231256        if ( mod( lc( f ), prime ) != 0 ) {
     
    239264        ptr++;
    240265    }
    241     return ( i == MAX_FP_FAC );
     266    return ( i == max_fp_fac );
    242267}
    243268
     
    250275    int i, j, save;
    251276    int p = pk.getp(), k = pk.getk();
    252     int no_iter = (int)(log(k)/log(2)+2);
     277    int no_iter = (int)(log( (double)k )/log(2)+2);
    253278    int * kvals = new int[no_iter];
    254279
    255     DEBOUTSL( cerr );
    256     DEBOUT( cerr, "quadratic lift called with p = " << p << "  and k = " << k );
    257     DEBOUTENDL( cerr );
     280    DEBOUTLN( cerr, "quadratic lift called with p = " << p << "  and k = " << k );
    258281    for ( j = 0, i = k; i > 1; i = ( i+1 ) / 2, j++ ) kvals[j] = i;
    259282    kvals[j] = 1;
     
    282305            j--;
    283306            setCharacteristic( p, kvals[j+1] );
    284             DEBOUTSL( cerr );
    285             DEBOUT( cerr, "lifting from p^" << kvals[j+1] << " to p^" << kvals[j] );
    286             DEBOUTENDL( cerr );
     307            DEBOUTLN( cerr, "lifting from p^" << kvals[j+1] << " to p^" << kvals[j] );
    287308            c = mapinto( c );
    288             DEBOUTLN( cerr, " !!! g = ", mapinto( g ) );
     309            DEBOUTLN( cerr, " !!! g = " << mapinto( g ) );
    289310            g1 = mapinto( lf ) / mapinto( lc( g ) ) * mapinto( g );
    290311            h1 = mapinto( lf ) / mapinto( lc( h ) ) * mapinto( h );
    291312//          (void)extgcd( g1, h1, a, b );
    292 //          DEBOUTLN( cerr, " a = ", aa );
    293 //          DEBOUTLN( cerr, " b = ", bb );
     313//          DEBOUTLN( cerr, " a = " << aa );
     314//          DEBOUTLN( cerr, " b = " << bb );
    294315            a = mapinto( a ); b = mapinto( b );
    295316            a += ( ( 1 - a * g1 ) *  a ) % h1;
    296317            b += ( ( 1 - b * h1 ) *  b ) % g1;
    297             DEBOUTLN( cerr, " a = ", a );
    298             DEBOUTLN( cerr, " b = ", b );
     318            DEBOUTLN( cerr, " a = " << a );
     319            DEBOUTLN( cerr, " b = " << b );
    299320            divrem( a * c, h1, q, r );
    300321            tmp = b * c + q * g1;
     
    308329        modulus = power( CanonicalForm(p), kvals[j] );
    309330        if ( mod( f - g * h, modulus ) != 0 )
    310             DEBOUTLN( cerr, "error at lift stage ", i );
     331            DEBOUTLN( cerr, "error at lift stage " << i );
    311332        i++;
    312333    }
     
    375396}
    376397
     398
    377399CFFList
    378400ZFactorizeUnivariate( const CanonicalForm& ff, bool issqrfree )
     
    383405    int i, k, exp, n;
    384406    bool ok;
    385     CFFList H, G, F[MAX_FP_FAC];
     407    CFFList H, F[max_fp_fac];
    386408    CFFList ZF;
    387     int * p = new int [MAX_FP_FAC];
     409    int * p = new int [max_fp_fac];
    388410    int * D = 0;
    389411    int * Dh = 0;
    390412    ListIterator<CFFactor> J, I;
     413
     414    DEBINCLEVEL( cerr, "ZFactorizeUnivariate" );
    391415    On( SW_SYMMETRIC_FF );
     416
     417    // get squarefree decomposition of f
    392418    if ( issqrfree )
    393419        H.append( CFFactor( g, 1 ) );
    394420    else
    395421        H = sqrFree( g );
     422
     423    DEBOUTLN( cerr, "H = " << H );
     424
     425    // cycle through squarefree factors of f
    396426    for ( J = H; J.hasItem(); ++J ) {
    397427        f = J.getItem().factor();
    398428        if ( f.inCoeffDomain() ) continue;
    399429        n = f.degree() / 2 + 1;
    400         if ( D != 0 ) {
    401             delete [] D;
    402             delete [] Dh;
    403         }
     430        delete [] D;
     431        delete [] Dh;
    404432        D = new int [n]; D[0] = n;
    405433        Dh = new int [n]; Dh[0] = n;
    406434        exp = J.getItem().exp();
     435
     436        // choose primes to factor f
    407437        TIMING_START(fac_choosePrimes);
    408438        ok = choosePrimes( p, f );
    409439        TIMING_END_AND_PRINT(fac_choosePrimes, "time to choose the primes: ");
    410440        if ( ! ok ) {
    411             DEBOUTLN( cerr, "error: no good prime found to factorize ", f );
    412             ASSERT( 0, "error: no good prime found" );
     441            DEBOUTLN( cerr, "warning: no good prime found to factorize " << f );
     442            WARN( ok, "there occured an error.  We went out of primes p\n"
     443                  "to factorize mod p.  Please send the example which caused\n"
     444                  "this error to the authors.  Nonetheless we will go on with the\n"
     445                  "calculations hoping the result will be correct.  Thank you.");
    413446            ZF.append( CFFactor( f, exp ) );
    414         }
    415         else {
    416             TIMING_START(fac_facModPrimes);
    417             for ( i = 0; i < MAX_FP_FAC; i++ ) {
    418                 setCharacteristic( p[i] );
    419                 fp = mapinto( f );
    420                 F[i] = FpFactorizeUnivariateCZ( fp, true );
     447            continue;
     448        }
     449
     450        // factorize f modulo certain primes
     451        TIMING_START(fac_facModPrimes);
     452        for ( i = 0; i < max_fp_fac; i++ ) {
     453            setCharacteristic( p[i] );
     454            fp = mapinto( f );
     455            F[i] = FpFactorizeUnivariateCZ( fp, true );
    421456//              if ( p[i] < 23 && fp.degree() < 10 )
    422457//                  F[i] = FpFactorizeUnivariateB( fp, true );
    423458//              else
    424459//                  F[i] = FpFactorizeUnivariateCZ( fp, true );
    425                 DEBOUTSL( cerr );
    426                 DEBOUT( cerr, "F[i] = " << F[i] << ", p = " << p[i] );
    427                 DEBOUTENDL( cerr );
     460            DEBOUTLN( cerr, "F[i] = " << F[i] << ", p = " << p[i] );
     461        }
     462        TIMING_END_AND_PRINT(fac_facModPrimes, "time to factorize mod primes: ");
     463        setCharacteristic( 0 );
     464
     465        // do some strange things with the D's
     466        initHG( D, F[0] );
     467        hgroup( D );
     468        DEBOUTHPRINT( cerr, "D = ", D );
     469        for ( i = 1; i < max_fp_fac; i++ ) {
     470            initHG( Dh, F[i] );
     471            hgroup( Dh );
     472            DEBOUTHPRINT( cerr, "Dh = ", Dh );
     473            hintersect( D, Dh );
     474            DEBOUTHPRINT( cerr, "D = ", D );
     475        }
     476
     477        // look which p gives the shortest factorization of f mod p
     478        // j: index of that p in p[]
     479        int min, j;
     480        min = F[0].length(), j = 0;
     481        for ( i = 1; i < max_fp_fac; i++ ) {
     482            if ( min >= F[i].length() ) {
     483                j = i; min = F[i].length();
    428484            }
    429             TIMING_END_AND_PRINT(fac_facModPrimes, "time to factorize mod primes: ");
    430             setCharacteristic( 0 );
    431 #ifdef DEBUGOUTPUT
    432             DEBOUTLN( cerr, "D = ", ' ' );
    433             hprint( D );
    434 #endif
    435             initHG( D, F[0] );
    436             hgroup( D );
    437 #ifdef DEBUGOUTPUT
    438             DEBOUTLN( cerr, "D = ", ' ' );
    439             hprint( D );
    440 #endif
    441             for ( i = 1; i < MAX_FP_FAC; i++ ) {
    442                 initHG( Dh, F[i] );
    443                 hgroup( Dh );
    444 #ifdef DEBUGOUTPUT
    445                 DEBOUTLN( cerr, "Dh = ", ' ' );
    446                 hprint( Dh );
    447 #endif
    448                 hintersect( D, Dh );
    449 #ifdef DEBUGOUTPUT
    450                 DEBOUTLN( cerr, "D = ", ' ' );
    451                 hprint( D );
    452 #endif
    453 
    454             }
    455             int min, j;
    456             min = F[0].length(), j = 0;
    457             for ( i = 1; i < MAX_FP_FAC; i++ ) {
    458                 if ( min >= F[i].length() ) {
    459                     j = i; min = F[i].length();
    460                 }
    461             }
    462             k = kBound( f, p[j] );
    463             CFArray theFactors( F[j].length() );
     485        }
     486        k = kBound( f, p[j] );
     487        CFArray theFactors( F[j].length() );
    464488//          pk = power( CanonicalForm( p[j] ), k );
    465489//          pkhalf = pk / 2;
    466             modpk pk( p[j], k );
    467             DEBOUTLN( cerr, "coeff bound = ", pk.getpk() );
    468             theModulus = pk;
    469             setCharacteristic( p[j] );
    470             fp = mapinto( f );
    471             F[j].sort( cmpFactor );
    472             I = F[j]; i = 0;
    473             TIMING_START(fac_liftFactors);
    474             while ( I.hasItem() ) {
    475                 DEBOUTLN( cerr, "factor to lift = ", I.getItem().factor() );
    476                 if ( isOn( SW_FAC_QUADRATICLIFT ) )
    477                     ok = UnivariateQuadraticLift( f, I.getItem().factor(), fp / I.getItem().factor(), pk, lc( f ), dummy1, dummy2 );
    478                 else
    479                     ok = UnivariateLinearLift( f, I.getItem().factor(), fp / I.getItem().factor(), pk, lc( f ), dummy1, dummy2 );
    480                 if ( ok ) {
    481                     // should be done in a more efficient way
    482                     DEBOUTLN( cerr, "dummy1 = ", dummy1 );
    483                     DEBOUTLN( cerr, "dummy2 = ", dummy2 );
    484                     f = dummy2;
    485                     fp /= I.getItem().factor();
    486                     ZF.append( CFFactor( dummy1, exp ) );
    487                     I.remove( 0 );
    488                     I = F[j];
    489                     i = 0;
    490                     DEBOUTLN( cerr, "F[j] = ", F[j] );
    491                 }
    492                 else {
    493                     DEBOUTLN( cerr, "i = ", i );
    494                     DEBOUTLN( cerr, "dummy1 = ", dummy1 );
    495                     setCharacteristic( 0 );
     490        modpk pk( p[j], k );
     491        DEBOUTLN( cerr, "coeff bound = " << pk.getpk() );
     492        theModulus = pk;
     493        setCharacteristic( p[j] );
     494        fp = mapinto( f );
     495        F[j].sort( cmpFactor );
     496        I = F[j]; i = 0;
     497        TIMING_START(fac_liftFactors);
     498        while ( I.hasItem() ) {
     499            DEBOUTLN( cerr, "factor to lift = " << I.getItem().factor() );
     500            if ( isOn( SW_FAC_QUADRATICLIFT ) )
     501                ok = UnivariateQuadraticLift( f, I.getItem().factor(), fp / I.getItem().factor(), pk, lc( f ), dummy1, dummy2 );
     502            else
     503                ok = UnivariateLinearLift( f, I.getItem().factor(), fp / I.getItem().factor(), pk, lc( f ), dummy1, dummy2 );
     504            if ( ok ) {
     505                // should be done in a more efficient way
     506                DEBOUTLN( cerr, "dummy1 = " << dummy1 );
     507                DEBOUTLN( cerr, "dummy2 = " << dummy2 );
     508                f = dummy2;
     509                fp /= I.getItem().factor();
     510                ZF.append( CFFactor( dummy1, exp ) );
     511                I.remove( 0 );
     512                I = F[j];
     513                i = 0;
     514                DEBOUTLN( cerr, "F[j] = " << F[j] );
     515            }
     516            else {
     517                DEBOUTLN( cerr, "i = " << i );
     518                DEBOUTLN( cerr, "dummy1 = " << dummy1 );
     519                setCharacteristic( 0 );
    496520//                  theFactors[i] = pk( dummy1 * pk.inverse( lc( dummy1 ) ) );
    497                     theFactors[i] = pk( dummy1 );
    498                     setCharacteristic( p[j] );
    499                     i++;
    500                     I++;
    501                 }
     521                theFactors[i] = pk( dummy1 );
     522                setCharacteristic( p[j] );
     523                i++;
     524                I++;
    502525            }
    503             TIMING_END_AND_PRINT(fac_liftFactors, "time to lift the factors: ");
    504             DEBOUTLN( cerr, "ZF = ", ZF );
    505             initHG( Dh, theFactors );
    506             hgroup( Dh );
    507 #ifdef DEBUGOUTPUT
    508             DEBOUTLN( cerr, "Dh = ", ' ' );
    509             hprint( Dh );
    510 #endif
    511             hintersect( D, Dh );
    512             setCharacteristic( 0 );
    513             for ( int l = i; l < F[j].length(); l++ )
    514                 theFactors[l] = 1;
    515             DEBOUTLN( cerr, "theFactors = ", theFactors );
    516             DEBOUTLN( cerr, "f = ", f );
    517             DEBOUTSL( cerr );
    518             DEBOUT( cerr, "p = " << pk.getp() << ", k = " << pk.getk() );
    519             DEBOUTENDL( cerr );
    520 #ifdef DEBUGOUTPUT
    521             DEBOUTLN( cerr, "D = ", ' ' );
    522             hprint( D );
    523 #endif
    524             lf = lc( f );
    525             (void)iextgcd( pk.getpk(), lf, dummy1, recip_lf );
    526             DEBOUTLN( cerr, "recip_lf = ", recip_lf );
    527             TIMING_START(fac_combineFactors);
    528             for ( i = 1; i < D[0]; i++ )
    529                 if ( D[i] != 0 )
    530                     while ( liftDegreeFactRec( theFactors, f, recip_lf, lf, pk, 0, i, ZF, exp ) );
    531             if ( degree( f ) > 0 )
    532                 ZF.append( CFFactor( f, exp ) );
    533             TIMING_END_AND_PRINT(fac_combineFactors, "time to combine the factors: ");
    534         }
    535     }
     526        }
     527        TIMING_END_AND_PRINT(fac_liftFactors, "time to lift the factors: ");
     528        DEBOUTLN( cerr, "ZF = " << ZF );
     529        initHG( Dh, theFactors );
     530        hgroup( Dh );
     531        DEBOUTHPRINT( cerr, "Dh = ", Dh );
     532        hintersect( D, Dh );
     533        setCharacteristic( 0 );
     534        for ( int l = i; l < F[j].length(); l++ )
     535            theFactors[l] = 1;
     536        DEBOUTLN( cerr, "theFactors = " << theFactors );
     537        DEBOUTLN( cerr, "f = " << f );
     538        DEBOUTLN( cerr, "p = " << pk.getp() << ", k = " << pk.getk() );
     539        DEBOUTHPRINT( cerr, "D = ", D );
     540        lf = lc( f );
     541        (void)iextgcd( pk.getpk(), lf, dummy1, recip_lf );
     542        DEBOUTLN( cerr, "recip_lf = " << recip_lf );
     543        TIMING_START(fac_combineFactors);
     544        for ( i = 1; i < D[0]; i++ )
     545            if ( D[i] != 0 )
     546                while ( liftDegreeFactRec( theFactors, f, recip_lf, lf, pk, 0, i, ZF, exp ) );
     547        if ( degree( f ) > 0 )
     548            ZF.append( CFFactor( f, exp ) );
     549        TIMING_END_AND_PRINT(fac_combineFactors, "time to combine the factors: ");
     550    }
     551
     552    // brush up our result
    536553    if ( ZF.getFirst().factor().inCoeffDomain() )
    537554        ZF.removeFirst();
     
    540557    else
    541558        ZF.insert( CFFactor( cont, 1 ) );
    542     if ( D != 0 ) {
    543         delete [] D;
    544         delete [] Dh;
    545     }
     559    delete [] D;
     560    delete [] Dh;
    546561    if ( ! symmsave )
    547562        Off( SW_SYMMETRIC_FF );
     563
     564    DEBDECLEVEL( cerr, "ZFactorizeUnivariate" );
    548565    return ZF;
    549566}
Note: See TracChangeset for help on using the changeset viewer.