source: git/factory/fac_multihensel.cc @ 8f5fe5

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