source: git/factory/fac_multihensel.cc @ 173e86

spielwiese
Last change on this file since 173e86 was 173e86, checked in by Jens Schmidt <schmidt@…>, 27 years ago
#include <config.h> added git-svn-id: file:///usr/local/Singular/svn/trunk@138 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 9.5 KB
Line 
1// emacs edit mode for this file is -*- C++ -*-
2// $Id: fac_multihensel.cc,v 1.3 1997-04-08 10:29:10 schmidt Exp $
3
4/*
5 * $Log: not supported by cvs2svn $
6 * Revision 1.2  1997/03/27 09:45:49  schmidt
7 * debug output rewritten
8 * some spurious preprocessor directives removed
9 *
10 * Revision 1.1  1996/12/06 14:46:22  schmidt
11 * Initial revision
12 *
13 */
14
15#include <config.h>
16
17#include "assert.h"
18#include "debug.h"
19#include "timing.h"
20
21#include "cf_defs.h"
22#include "cf_eval.h"
23#include "cf_map.h"
24#include "cf_binom.h"
25#include "fac_util.h"
26#include "fac_iterfor.h"
27#include "cf_iter.h"
28
29
30TIMING_DEFINE_PRINT(fac_solve);
31TIMING_DEFINE_PRINT(fac_modpk);
32TIMING_DEFINE_PRINT(fac_corrcoeff);
33TIMING_DEFINE_PRINT(fac_extgcd);
34
35static void
36extgcdrest ( const CanonicalForm & a, const CanonicalForm & b, const CanonicalForm & s, const CanonicalForm & t, const CanonicalForm & c, CanonicalForm & S, CanonicalForm & T, const modpk & pk )
37{
38    CanonicalForm sigma = s * c, tau = t * c;
39//    divremainder( sigma, b, T, S, pk );
40//    T = pk( tau + T * a );
41    divrem( sigma, b, T, S );
42    T = tau + T * a;
43}
44
45static void
46solveF ( const CFArray & P, const CFArray & Q, const CFArray & S, const CFArray & T, const CanonicalForm & C, const modpk & pk, int r, CFArray & a )
47{
48    setCharacteristic( pk.getp(), pk.getk() );
49    CanonicalForm g, bb, b = mapinto( C );
50    int j;
51    for ( j = 1; j < r; j++ ) {
52        extgcdrest( mapinto( P[j] ), mapinto( Q[j] ), mapinto( S[j] ), mapinto( T[j] ), b, bb, a[j], pk );
53        b = bb;
54    }
55    a[r] = b;
56    setCharacteristic( 0 );
57    for ( j = 1; j <= r; j++ )
58        a[j] = mapinto( a[j] );
59}
60
61static CanonicalForm
62evalF ( const CFArray & P, const CFArray & Q, const CFArray & A, int r )
63{
64    CanonicalForm pprod = 1, sum = 0;
65    for ( int i = 1; i <= r; i++ ) {
66        sum += pprod * A[i] * Q[i];
67        pprod *= P[i];
68    }
69    return sum;
70}
71
72static CanonicalForm
73derivAndEval ( const CanonicalForm & f, int n, const Variable & x, const CanonicalForm & a )
74{
75    if ( n == 0 )
76        return f( a, x );
77    else if ( f.degree( x ) < n )
78        return 0;
79    else {
80        CFIterator i;
81        CanonicalForm sum = 0, fact;
82        int min, j;
83        Variable v = Variable( f.level() + 1 );
84        for ( i = swapvar( f, x, v); i.hasTerms() && i.exp() >= n; i++ ) {
85            fact = 1;
86            min = i.exp() - n;
87            for ( j = i.exp(); j > min; j-- )
88                fact *= j;
89            sum += fact * i.coeff() * power( v, min );
90        }
91        return sum( a, v );
92    }
93}
94
95static CanonicalForm
96checkCombination ( const CanonicalForm & W, const Evaluation & a, const IteratedFor & e, int k )
97{
98    CanonicalForm dW = W;
99    int i, j;
100    for ( i = k-1; i >= 2 && ! dW.isZero(); i-- )
101        dW = derivAndEval( dW, e[i], Variable( i ), a[i] );
102    if ( ! dW.isZero() ) {
103        CanonicalForm fact = 1;
104        for ( i = 2; i < k; i++ )
105            for ( j = 2; j <= e[i]; j++ )
106                fact *= j;
107        dW /= fact;
108    }
109    return dW;
110}
111
112static CanonicalForm
113prodCombination ( const Evaluation & a, const IteratedFor & e, int k )
114{
115    CanonicalForm p = 1;
116    for ( int i = k-1; i > 1; i-- )
117        p *= binomialpower( Variable(i), -a[i], e[i] );
118    return p;
119}
120
121static CanonicalForm
122check_dummy( const CFArray &a, const CFArray & P, const CFArray & Q )
123{
124    int i, r = a.size();
125    CanonicalForm res, prod;
126    res = 0;
127    prod = 1;
128    for ( i = 1; i <= r; i++ ) {
129        res += prod * a[i] * Q[i];
130        prod *= P[i];
131    }
132    return res;
133}
134
135static bool
136check_e( const IteratedFor & e, int k, int m, int * n )
137{
138    int sum = 0;
139    for ( int i = 2; i < k; i++ ) {
140        sum += e[i];
141        if ( e[i] > n[i] )
142            return false;
143    }
144    return sum == m+1;
145}
146
147static CanonicalForm
148modDeltak ( const CanonicalForm & f, const Evaluation & A, int k, int * n )
149{
150    CanonicalForm result = f;
151    for ( int i = 2; i < k; i++ ) {
152        result.mod( binomialpower( Variable(i), -A[i], n[i]+1 ) );
153    }
154    return result;
155}
156
157static CFArray
158findCorrCoeffs ( const CFArray & P, const CFArray & Q, const CFArray & P0, const CFArray & Q0, const CFArray & S, const CFArray & T, const CanonicalForm & C, const Evaluation & I, const modpk & pk, int r, int k, int h, int * n )
159{
160    DEBINCLEVEL( cerr, "findCorrCoeffs" );
161    int i, m;
162    CFArray A(1,r), a(1,r);
163    CanonicalForm C0, Dm, g, prodcomb;
164
165    TIMING_START(fac_solve);
166    C0 = pk( I( C, 2, k-1 ), true );
167    solveF( P0, Q0, S, T, 1, pk, r, a );
168    TIMING_END(fac_solve);
169
170    DEBOUTLN( cerr, "trying to find correction coefficients for ", C );
171    DEBOUTLN( cerr, "which evaluates to ", C0 );
172
173    for ( i = 1; i <= r; i++ )
174        A[i] = remainder( pk( a[i] * C0 ), P0[i], pk );
175    DEBOUTLN( cerr, "the first approximation of the correction coefficients is ", A );
176#ifdef DEBUGOUTPUT
177    if ( check_dummy( A, P, Q ) - C != 0 ) {
178        DEBOUTLN( cerr, "there is an error detected, the correction coefficients do not", ' ' );
179        DEBOUTLN( cerr, "fulfill equation F(A)", ' ' );
180        DEBOUTLN( cerr, "corresponding P ", P );
181        DEBOUTLN( cerr, "              Q ", Q );
182    }
183#endif
184    for ( m = 0; m <= h && ( m == 0 || Dm != 0 ); m++ ) {
185        Dm = pk( evalF( P, Q, A, r ) - C );
186        if ( Dm != 0 ) {
187            if ( k == 2 ) {
188                TIMING_START(fac_solve);
189                solveF( P0, Q0, S, T, Dm, pk, r, a );
190                TIMING_END(fac_solve);
191                for ( i = 1; i <= r; i++ )
192                    A[i] -= a[i];
193            }
194            else {
195                IteratedFor e( 2, k-1, m+1 );
196                while ( e.iterations_left() ) {
197                    if ( check_e( e, k, m, n ) ) {
198                        g = pk( checkCombination( Dm, I, e, k ) );
199                        if ( ! g.isZero() && ! (g.mvar() > Variable(1)) ) {
200                            prodcomb = prodCombination( I, e, k );
201//                          Dm = Dm - g * prodcomb;
202                            TIMING_START(fac_solve);
203                            solveF( P0, Q0, S, T, g, pk, r, a );
204                            TIMING_END(fac_solve);
205                            for ( i = 1; i <= r; i++ ) {
206//                              A[i] -= a[i] * prodcomb;
207                                A[i] = pk( A[i] - a[i] * prodcomb );
208                            }
209                        }
210                    }
211                    e++;
212                }
213            }
214        }
215        DEBOUTLN( cerr, "the correction coefficients at step ", m );
216        DEBOUTLN( cerr, "are now ", A );
217#ifdef DEBUGOUTPUT
218    if ( check_dummy( A, P, Q ) - C != 0 ) {
219        DEBOUTLN( cerr, "there is an error detected, the correction coefficients do not", ' ' );
220        DEBOUTLN( cerr, "fulfill equation F(A)", ' ' );
221        DEBOUTLN( cerr, "corresponding P ", P );
222        DEBOUTLN( cerr, "              Q ", Q );
223    }
224#endif
225    }
226    DEBDECLEVEL( cerr, "findCorrCoeffs" );
227    return A;
228}
229
230
231static bool
232liftStep ( CFArray & P, int k, int r, int t, const modpk & b, const Evaluation & A, const CFArray & lcG, const CanonicalForm & Uk, int * n, int h )
233{
234    DEBINCLEVEL( cerr, "liftStep" );
235    CFArray K( 1, r ), Q( 1, r ), Q0( 1, r ), P0( 1, r ), S( 1, r ), T( 1, r ), alpha( 1, r );
236    CanonicalForm Rm, C, D, xa = Variable(k) - A[k];
237    CanonicalForm xa1 = xa, xa2 = xa*xa;
238    CanonicalForm dummy;
239    int i, m;
240
241    DEBOUTLN( cerr, "we are now performing the liftstep to reach ", Variable(k) );
242    DEBOUTLN( cerr, "the factors so far are ", P );
243
244    for ( i = 1; i <= r; i++ ) {
245        Variable vm = Variable( t + 1 );
246        Variable v1 = Variable(1);
247        K[i] = swapvar( replaceLc( swapvar( P[i], v1, vm ), swapvar( A( lcG[i], k+1, t ), v1, vm ) ), v1, vm );
248        P[i] = A( K[i], k, t );
249    }
250    DEBOUTLN( cerr, "lift K = ", K );
251
252//    d = degree( Uk, Variable( k ) );
253
254    TIMING_START(fac_extgcd);
255    Q[r] = 1;
256    for ( i = r; i > 1; i-- ) {
257        Q[i-1] = Q[i] * P[i];
258        P0[i] = A( P[i], 2, k-1 );
259        Q0[i] = A( Q[i], 2, k-1 );
260        extgcd( P0[i], Q0[i], S[i], T[i], b );
261    }
262    P0[1] = A( P[1], 2, k-1 );
263    Q0[1] = A( Q[1], 2, k-1 );
264    extgcd( P0[1], Q0[1], S[1], T[1], b );
265    TIMING_END(fac_extgcd);
266
267    for ( m = 1; m <= n[k]+1; m++ ) {
268        TIMING_START(fac_modpk);
269        Rm = modDeltak( prod( K ) - Uk, A, k, n );
270        TIMING_END(fac_modpk);
271#ifdef DEBUGOUTPUT
272        if ( mod( Rm, xa1 ) != 0 ) {
273            DEBOUTLN( cerr, "something seems not to be ok with Rm which is ", Rm );
274            DEBOUTLN( cerr, "and should reduce to zero modulo ", xa1 );
275        }
276#endif
277        if ( mod( Rm, xa2 ) != 0 ) {
278            C = derivAndEval( Rm, m, Variable( k ), A[k] );
279            D = 1;
280            for ( i = 2; i <= m; i++ ) D *= i;
281            C /= D;
282
283            TIMING_START(fac_corrcoeff);
284            alpha = findCorrCoeffs( P, Q, P0, Q0, S, T, C, A, b, r, k, h, n ); // -> h berechnen
285            TIMING_END(fac_corrcoeff);
286// #ifdef DEBUGOUTPUT
287//          dummy = check_dummy( alpha, P, Q );
288//          if ( b(modDeltak( dummy, A, k, n )) != b(modDeltak( C, A, k, n )) ) {
289//              DEBOUTLN( cerr, "lift fault", ' ' );
290//              DEBDECLEVEL( cerr, "liftStep" );
291//              return false;
292//          }
293// #endif
294            for ( i = 1; i <= r; i++ )
295                K[i] = b(K[i] - alpha[i] * xa1);
296            DEBOUTLN( cerr, "the corrected K's are now ", K );
297        }
298        xa1 = xa2;
299        xa2 *= xa;
300    }
301    for ( i = 1; i <= r; i++ )
302        P[i] = K[i];
303    if ( prod( K ) - Uk != 0 ) {
304        DEBOUTLN( cerr, "the liftstep produced the wrong result", ' ' );
305        DEBOUTLN( cerr, "the product of the factors calculated so far is ", prod(K) );
306        DEBOUTLN( cerr, "and the Uk that should have been reached is ", Uk );
307        DEBDECLEVEL( cerr, "liftStep" );
308        return false;
309    }
310    DEBOUTLN( cerr, "the lift seems ok so far", ' ' );
311    DEBDECLEVEL( cerr, "liftStep" );
312    return true; // check for divisibility
313}
314
315bool
316Hensel ( const CanonicalForm & U, CFArray & G, const CFArray & lcG, const Evaluation & A, const modpk & bound, const Variable & x )
317{
318    DEBINCLEVEL( cerr, "Hensel" );
319    int k, i, h, t = A.max();
320    bool goodeval = true;
321    CFArray Uk( A.min(), A.max() );
322    int * n = new int[t+1];
323
324    Uk[t] = U;
325    for ( k = t-1; k > 1; k-- ) {
326        Uk[k] = Uk[k+1]( A[k+1], Variable( k+1 ) );
327        n[k] = degree( Uk[k], Variable( k ) );
328    }
329    for ( k = A.min(); goodeval && (k <= t); k++ ) {
330        h = totaldegree( Uk[k], Variable(A.min()), Variable(k-1) );
331        for ( i = A.min(); i <= k; i++ )
332            n[i] = degree( Uk[k], Variable(i) );
333        goodeval = liftStep( G, k, G.max(), t, bound, A, lcG, Uk[k], n, h );
334        DEBOUTLN( cerr, "Factorization so far: ", G );
335    }
336    DEBDECLEVEL( cerr, "Hensel" );
337    return goodeval;
338}
Note: See TracBrowser for help on using the repository browser.