source: git/factory/fac_multihensel.cc @ 72dd6e

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