source: git/factory/fac_multihensel.cc @ d754b7

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