Changeset 125195 in git for factory


Ignore:
Timestamp:
Jul 23, 1996, 11:18:35 AM (28 years ago)
Author:
Rüdiger Stobbe <stobbe@…>
Branches:
(u'spielwiese', '82fc009ea2b0098c1a4896c841bb70860976bdfc')
Children:
aef6ab2c9227b84444bc58f9c4dd89ca55f754f2
Parents:
70161e0eb6f97b4a4008f8a9dee8629adf6bb83c
Message:
Version in work.


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

Legend:

Unmodified
Added
Removed
  • factory/fac_multivar.cc

    r70161e r125195  
    11// emacs edit mode for this file is -*- C++ -*-
    2 // $Id: fac_multivar.cc,v 1.0 1996-05-17 10:59:45 stobbe Exp $
     2// $Id: fac_multivar.cc,v 1.1 1996-07-23 09:18:35 stobbe Exp $
    33
    44/*
    55$Log: not supported by cvs2svn $
     6Revision 1.0  1996/05/17 10:59:45  stobbe
     7Initial revision
     8
    69*/
    710
    8 //#define TIMING
    9 
    10 #ifndef NDEBUG
    11 //#define DEBUGOUTPUT
    12 #endif
    13 
    14 #include <time.h>
    15 #include <sys/times.h>
     11#define TIMING
     12#undef DEBUGOUTPUT
     13
     14#include "timing.h"
    1615
    1716#include "assert.h"
     
    2322#include "fac_util.h"
    2423#include "cf_binom.h"
    25 #include "fac_util.h"
    26 #include "fac_diophant.h"
    27 #include "fac_iterfor.h"
    2824#include "cf_iter.h"
    2925
    3026
    31 static CanonicalForm
    32 check_dummy( const CFArray &a, const CFArray & P, const CFArray & Q )
    33 {
    34     int i, r = a.size();
    35     CanonicalForm res, prod;
    36     res = 0;
    37     prod = 1;
    38     for ( i = 1; i <= r; i++ ) {
    39         res += prod * a[i] * Q[i];
    40         prod *= P[i];
    41     }
    42     return res;
    43 }
    44 
    45 static bool
    46 check_e( const IteratedFor & e, int k, int m, int * n )
    47 {
    48     int sum = 0;
    49     for ( int i = 2; i < k; i++ ) {
    50         sum += e[i];
    51         if ( e[i] > n[i] )
    52             return false;
    53     }
    54     return sum == m+1;
    55 }
     27TIMING_DEFINE_PRINT(fac_content);
     28TIMING_DEFINE_PRINT(fac_findeval);
     29TIMING_DEFINE_PRINT(fac_distrib);
     30
    5631
    5732static CFArray
     
    6035    int n;
    6136    CFFListIterator I = L;
     37    bool negate = false;
     38
    6239    if ( ! I.hasItem() )
    6340        n = 0;
    6441    else  if ( I.getItem().factor().inBaseDomain() ) {
     42        negate = I.getItem().factor().sign() < 0;
    6543        I++;
    6644        n = L.length();
     
    8462        I++;
    8563    }
     64    if ( negate )
     65        result[1] = -result[1];
    8666    return result;
    8767}
    8868
    89 CanonicalForm
    90 maxCoeff( const CanonicalForm & f )
    91 {
    92     if ( f.inCoeffDomain() )
    93         return abs( f );
    94     else {
    95         CanonicalForm M = 0, m;
    96         for ( CFIterator i = f; i.hasTerms(); i++ )
    97             if ( (m = maxCoeff( i.coeff() )) > M )
    98                 M = m;
    99         return M;
    100     }
    101 }
    102 
    103 modpk
     69static modpk
    10470coeffBound ( const CanonicalForm & f, int p )
    10571{
     
    12187nonDivisors ( CanonicalForm omega, CanonicalForm delta, const CFArray & F, CFArray & d )
    12288{
    123     DEBOUTLN( cerr, "omega = ", omega );
    124     DEBOUTLN( cerr, "delta = ", delta );
    125     DEBOUTLN( cerr, "F = ", F );
     89    DEBOUTLN( cerr, "nondivisors omega = ", omega );
     90    DEBOUTLN( cerr, "nondivisors delta = ", delta );
     91    DEBOUTLN( cerr, "nondivisors F = ", F );
    12692    CanonicalForm q, r;
    12793    int k = F.size();
     
    145111
    146112static void
    147 findEvaluation ( const CanonicalForm & U, const CanonicalForm & V, const CanonicalForm & omega, const CFFList & F, Evaluation & A, CanonicalForm & U0, CanonicalForm & delta, int r )
     113findEvaluation ( const CanonicalForm & U, const CanonicalForm & V, const CanonicalForm & omega, const CFFList & F, Evaluation & A, CanonicalForm & U0, CanonicalForm & delta, CFArray & D, int r )
    148114{
    149115    CanonicalForm Vn;
     
    151117    int j;
    152118    bool found = false;
    153     CFArray FF = CFArray( 1, F.length() ), D;
     119    CFArray FF = CFArray( 1, F.length() );
    154120    if ( r > 0 )
    155121        A.nextpoint();
     
    158124        if ( Vn != 0 ) {
    159125            U0 = A( U );
    160             DEBOUTLN( cerr, "U0 = ", U0 );
     126            DEBOUTLN( cerr, "findev U0 = ", U0 );
    161127            if ( isSqrFree( U0 ) ) {
    162128                delta = content( U0 );
    163                 DEBOUTLN( cerr, "content( U0 ) = ", delta );
     129                DEBOUTLN( cerr, "findev content( U0 ) = ", delta );
    164130                for ( I = F, j = 1; I.hasItem(); I++, j++ )
    165131                    FF[j] = A( I.getItem().factor() );
     
    167133            }
    168134            else {
    169                 DEBOUTLN( cerr, "not sqrfree :", sqrFree( U0 ) );
     135                DEBOUTLN( cerr, "findev not sqrfree :", sqrFree( U0 ) );
    170136            }
    171137        }
     
    175141}
    176142
    177 static CanonicalForm
    178 modDeltak ( const CanonicalForm & f, const Evaluation & A, int k, int * n )
    179 {
    180     CanonicalForm result = f;
    181     for ( int i = 2; i < k; i++ ) {
    182         result.mod( binomialpower( Variable(i), -A[i], n[i]+1 ) );
    183     }
    184     return result;
    185 }
    186 
    187 CFArray
    188 findCorrCoeffs ( const CFArray & P, const CFArray & Q, const CanonicalForm & C, const Evaluation & I, const modpk & pk, int r, int k, int h, int * n )
    189 {
    190 #ifdef TIMING
    191     static float totalsolve = 0.0;
    192     struct tms ts, te;
    193 #endif
    194     int i, m;
    195     CFArray A(1,r), a(1,r), P0(1,r), Q0(1,r);
    196     CanonicalForm C0, Dm, g, prodcomb;
    197 
    198 #ifdef TIMING
    199     times( &ts );
    200 #endif
    201     for ( i = 1; i <= r; i++ ) {
    202         P0[i] = I( P[i], 2, k-1 );
    203         Q0[i] = I( Q[i], 2, k-1 );
    204     }
    205     C0 = pk( I( C, 2, k-1 ), true );
    206     DEBOUTLN( cerr, "C  = ", C );
    207     DEBOUTLN( cerr, "C0 = ", C0 );
    208     DEBOUTLN( cerr, "P0 = ", P0 );
    209     DEBOUTLN( cerr, "Q0 = ", Q0 );
    210     solveF( P0, Q0, 1, pk, r, a );
    211 //    solveF( P0, Q0, C0, pk, r, a );
    212     DEBOUTLN( cerr, "a  = ", a );
    213 #ifdef TIMING
    214     times( &te );
    215     cerr.setf( ios::fixed, ios::floatfield );
    216     cerr.precision(2);
    217     cerr << "time used for solve: "
    218          << float(te.tms_utime-ts.tms_utime) / CLK_TCK << " sec, total time = ";
    219     totalsolve += float(te.tms_utime-ts.tms_utime) / CLK_TCK;
    220     cerr << totalsolve << endl;
    221 #endif
    222     for ( i = 1; i <= r; i++ )
    223         A[i] = remainder( pk( a[i] * C0 ), P0[i], pk );
    224     for ( m = 0; m <= h && ( m == 0 || Dm != 0 ); m++ ) {
    225         Dm = pk( evalF( P, Q, A, r ) - C );
    226         if ( Dm != 0 ) {
    227             if ( k == 2 ) {
    228                 solveF( P0, Q0, Dm, pk, r, a );
    229                 for ( i = 1; i <= r; i++ )
    230                     A[i] -= a[i];
    231             }
    232             else {
    233                 IteratedFor e( k-1, m+1 );
    234                 while ( e.iterations_left() ) {
    235                     if ( check_e( e, k, m, n ) ) {
    236                         g = pk( checkCombination( Dm, I, e, k ) );
    237                         if ( ! g.isZero() && ! (g.mvar() > Variable(1)) ) {
    238                             prodcomb = prodCombination( I, e, k );
    239 //                          Dm = Dm - g * prodcomb;
    240                             solveF( P0, Q0, g, pk, r, a );
    241                             for ( i = 1; i <= r; i++ ) {
    242 //                              A[i] -= a[i] * prodcomb;
    243                                 A[i] = pk( A[i] - a[i] * prodcomb );
    244                             }
    245                         }
    246                     }
    247                     e++;
    248                 }
    249             }
    250         }
    251     }
    252     return A;
    253 }
    254                
    255 
    256 bool
    257 liftStep ( CFArray & P, int k, int r, int t, const modpk & b, const Evaluation & A, const CFArray & lcG, const CanonicalForm & Uk, int * n, int d, int h )
    258 {
    259     CFArray K( 1, r ), Q( 1, r ), alpha( 1, r );
    260     CanonicalForm Rm, C, D, xa = Variable(k) - A[k];
    261     CanonicalForm xa1 = xa, xa2 = xa*xa;
    262     CanonicalForm dummy;
    263     int i, m;
    264 
    265     DEBOUTLN( cerr, " P = ", P );
    266     DEBOUTLN( cerr, " lc= ", lcG );
    267     for ( i = 1; i <= r; i++ ) {
    268         Variable vm = Variable( t + 1 );
    269         Variable v1 = Variable(1);
    270         K[i] = swapvar( replaceLc( swapvar( P[i], v1, vm ), A( lcG[i], k+1, t )), v1, vm );
    271     }
    272     DEBOUTLN( cerr, " K = ", K );
    273 
    274 //    d = degree( Uk, Variable( k ) );
    275 
    276     Q[r] = 1;
    277     for ( i = r; i > 1; i-- )
    278         Q[i-1] = Q[i] * P[i];
    279     for ( m = 1; m <= n[k]+1; m++ ) {
    280         Rm = modDeltak( prod( K ) - Uk, A, k, n );
    281 #ifdef 0
    282 #ifdef DEBUGOUTPUT
    283         cerr << "p(K)-Uk = " << prod( K ) - Uk << endl;
    284         cerr << "p(K)-Uk mod deltak = " << Rm << endl;
    285         cerr << "xa1 = " << xa1 << endl;
    286         cerr << "Rm % xa1 = " << Rm % xa1 << endl;
    287         if ( mod( Rm, xa1 ) != 0 )
    288             cerr << "Rm not ok : " << Rm << endl;
    289 #endif
    290 #endif
    291         if ( mod( Rm, xa2 ) != 0 ) {
    292             C = derivAndEval( Rm, m, Variable( k ), A[k] );
    293             D = 1;
    294             for ( i = 2; i <= m; i++ ) D *= i;
    295             C /= D;
    296 
    297             alpha = findCorrCoeffs( P, Q, C, A, b, r, k, h, n ); // -> h berechnen
    298             dummy = check_dummy( alpha, P, Q );
    299             if ( b(modDeltak( dummy, A, k, n )) != b(modDeltak( C, A, k, n )) ) {
    300 #ifdef DEBUGOUTPUT
    301                 cerr << "fault in multi diophant" << endl
    302                      << " alpha = " << alpha << endl
    303                      << " dummy = " << modDeltak( dummy, A, k, n ) << endl
    304                      << " C = " << modDeltak( C, A, k, n ) << endl;
    305 #endif
    306                 cerr << "fault" << endl;
    307                 return false;
    308             }
    309             for ( i = 1; i <= r; i++ )
    310                 K[i] = b(K[i] - alpha[i] * xa1);
    311         }
    312         xa1 = xa2;
    313         xa2 *= xa;
    314     }
    315     for ( i = 1; i <= r; i++ )
    316         P[i] = K[i];
    317     DEBOUTLN( cerr, "K = ", K );
    318     if ( prod( K ) - Uk != 0 )
    319         return false;
    320     return true; // check for divisibility
    321 }
    322 
    323 CFArray
     143static CFArray
    324144ZFactorizeMulti ( const CanonicalForm & arg )
    325145{
     146    DEBOUTLN( cerr, "MULTIFACTOR START ----------------------------------- level = ", arg.level() );
    326147    CFMap M;
    327148    CanonicalForm UU, U = compress( arg, M );
     
    330151    CFFList F = factorize( V );
    331152    CFFListIterator I, J;
    332     CFArray G, lcG;
     153    CFArray G, lcG, D;
    333154    int i, j, k, m, r, maxdeg, h;
    334155    REvaluation A( 2, t, IntRandom( 100 ) );
     
    339160
    340161#ifdef DEBUGOUTPUT
    341     cerr << "U = " << U << endl;
     162    cerr << "fac U = " << U << endl;
    342163    for ( i = 1; i <= level( U ); i++ )
    343         cerr << " deg(U," << Variable(i) << ") = "
     164        cerr << "fac deg(U," << Variable(i) << ") = "
    344165             << degree( U, Variable(i) ) << endl
    345              << " lc(U," << Variable(i) << ") = "
     166             << "fac lc(U," << Variable(i) << ") = "
    346167             << LC( U, Variable(i) ) << endl;
    347168#endif
     
    370191//      A.nextpoint();
    371192
     193    DEBOUTLN( cerr, "---------------------------------------", ' ' );
     194
    372195    while ( ! goodeval ) {
    373         findEvaluation( U, V, omega, F, A, U0, delta, r );
    374         cerr << "A = " << A << endl;
    375         DEBOUTLN( cerr, "evaluation = ", A );
    376         DEBOUTLN( cerr, "fac(U0) = ", factorize( U0, false ) );
    377         G = conv_to_factor_array( factorize( U0, true ) );
     196        TIMING_START(fac_findeval);
     197        findEvaluation( U, V, omega, F, A, U0, delta, D, r );
     198        TIMING_END(fac_findeval);
     199        DEBOUTLN( cerr, "fac evaluation = ", A );
     200        G = conv_to_factor_array( factorize( U0, false ) );
     201        DEBOUTLN( cerr, "fac fac(U0) = ", G );
    378202        b = coeffBound( U, getZFacModulus().getp() );
    379203#ifdef DEBUGOUTPUT
    380         cout << "(fac: U  = " << U << endl
     204        cerr << "p^k(" << U.level() << ") = " << b.getp() << "^" << b.getk() << endl;
     205        cerr << "(fac: U  = " << U << endl
    381206             << "      U0 = " << U0 << endl
    382207             << "      V  = " << V << endl
     
    386211             << "      b  = " << b.getp() << "^" << b.getk() << endl
    387212             << "      ub = " << b.getp() << "^" << getZFacModulus().getk() << endl
    388              << "      G  = " << G << " )" << endl;
     213             << "      G  = " << G << endl
     214             << "      D  = " << D << " )" << endl;
    389215#endif
    390216        r = G.size();
    391         DEBOUTLN( cerr, "SIZE OF UNIFAC = ", r );
     217        DEBOUTLN( cerr, "fac SIZE OF UNIFAC = ", r );
    392218        lcG = CFArray( 1, r );
    393219        for ( j = 1; j <= r; j ++ )
    394220            lcG[j] = 1;
    395221
     222        TIMING_START(fac_distrib);
    396223        goodeval = true;
    397224        for ( I = F; goodeval && I.hasItem(); I++ ) {
     
    409236            goodeval = (m == 0);
    410237        }
     238        TIMING_END(fac_distrib);
    411239        if ( goodeval ) {
    412             DEBOUTLN( cerr, "good evaluation, lcG = ", lcG );
    413             DEBOUTLN( cerr, "                   G = ", G );
    414             DEBOUTLN( cerr, "delta = ", delta );
    415240            if ( delta != 1 ) {
    416241                for ( j = 1; j <= r; j++ ) {
     
    422247                    delta /= gt;
    423248                }
    424                 DEBOUTLN( cerr, "delta = ", delta );
     249                DEBOUTLN( cerr, "fac delta = ", delta );
    425250                if ( delta != 1 ) {
    426251                    for ( j = 1; j <= r; j++ ) {
     
    435260            else
    436261                UU = U;
    437             DEBOUTLN( cerr, "good evaluation, lcG = ", lcG );
    438             DEBOUTLN( cerr, "                   G = ", G );
    439             CFArray Uk( 1, t );
    440             Uk[t] = UU;
    441             int * n = new int[t+1];
    442             for ( k = t-1; k > 1; k-- ) {
    443                 Uk[k] = Uk[k+1]( A[k+1], Variable( k+1 ) );
    444                 n[k] = degree( Uk[k], Variable( k ) );
    445             }
    446             for ( k = 2; goodeval && (k <= t); k++ ) {
    447                 DEBOUTLN( cerr, "begin lifting -------> k  = ", k );
    448                 DEBOUTLN( cerr, "                       Uk = ", Uk[k] ); 
    449                 h = totaldegree( Uk[k], Variable(2), Variable(k-1) );
    450                 DEBOUTLN( cerr, "                       h  = ", h ); 
    451                 for ( i = 2; i <= k; i++ )
    452                     n[i] = degree( Uk[k], Variable(i) );
    453                 goodeval = liftStep( G, k, r, t, b, A, lcG, Uk[k], n, maxdeg, h );
    454             }
     262            DEBOUTLN( cerr, "fac good evaluation, lcG = ", lcG );
     263            DEBOUTLN( cerr, "fac                    G = ", G );
     264            DEBOUTLN( cerr, "fac delta = ", delta );
     265            DEBOUTLN( cerr, "fac omega = ", omega );
     266            for ( j = 1; j <= r; j++ ) {
     267                gt = A( lcG[j] );
     268                if ( gt != lc( G[j] ) ) {
     269                    gt = lc( G[j] ) / gt;
     270                    lcG[j] *= gt;
     271                    omega /= gt;
     272                }
     273            }
     274            DEBOUTLN( cerr, "fac good evaluation, lcG = ", lcG );
     275            DEBOUTLN( cerr, "fac                    G = ", G );
     276            DEBOUTLN( cerr, "fac delta = ", delta );
     277            DEBOUTLN( cerr, "fac omega = ", omega );
     278            goodeval = Hensel( UU, G, lcG, A, b, Variable(1) );
    455279        }
    456280    }
     
    460284    if ( negate )
    461285        G[1] = -G[1];
     286    DEBOUTLN( cerr, "MULTIFACTOR END   ----------------------------------- level = ", arg.level() );
    462287    return G;
    463288}
     
    486311        }
    487312        else {
     313            TIMING_START(fac_content);
    488314            g = compress( i.getItem().factor(), M );
    489315            // now after compress g contains Variable(1)
     
    494320            cont = swapvar( cont, v1, vm );
    495321            n = i.getItem().exp();
     322            TIMING_END(fac_content);
     323            DEBOUTLN( cerr, "now after content ...", ' ' );
    496324            if ( g.isUnivariate() ) {
    497325                G = factorize( g, true );
Note: See TracChangeset for help on using the changeset viewer.