source: git/factory/fac_multihensel.cc @ f3a82f4

jengelh-datetimespielwiese
Last change on this file since f3a82f4 was f3a82f4, checked in by Martin Lee <martinlee84@…>, 11 years ago
fix: compilation errors with --enable-debugoutput
  • Property mode set to 100644
File size: 10.7 KB
Line 
1/* emacs edit mode for this file is -*- C++ -*- */
2/* $Id$ */
3
4#include "config.h"
5
6#include "cf_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    {
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    {
56        sum += pprod * A[i] * Q[i];
57        pprod *= P[i];
58    }
59    return sum;
60}
61
62static CanonicalForm
63derivAndEval ( const CanonicalForm & f, int n, const Variable & x, const CanonicalForm & a )
64{
65    if ( n == 0 )
66        return f( a, x );
67    else if ( f.degree( x ) < n )
68        return 0;
69    else
70    {
71        CFIterator i;
72        CanonicalForm sum = 0, fact;
73        int min, j;
74        Variable v = Variable( f.level() + 1 );
75        for ( i = swapvar( f, x, v); i.hasTerms() && i.exp() >= n; i++ )
76        {
77            fact = 1;
78            min = i.exp() - n;
79            for ( j = i.exp(); j > min; j-- )
80                fact *= j;
81            sum += fact * i.coeff() * power( v, min );
82        }
83        return sum( a, v );
84    }
85}
86
87static CanonicalForm
88checkCombination ( const CanonicalForm & W, const Evaluation & a, const IteratedFor & e, int k )
89{
90    CanonicalForm dW = W;
91    int i, j;
92    for ( i = k-1; i >= 2 && ! dW.isZero(); i-- )
93        dW = derivAndEval( dW, e[i], Variable( i ), a[i] );
94    if ( ! dW.isZero() ) {
95        CanonicalForm fact = 1;
96        for ( i = 2; i < k; i++ )
97            for ( j = 2; j <= e[i]; j++ )
98                fact *= j;
99        dW /= fact;
100    }
101    return dW;
102}
103
104static CanonicalForm
105prodCombination ( const Evaluation & a, const IteratedFor & e, int k )
106{
107    CanonicalForm p = 1;
108    for ( int i = k-1; i > 1; i-- )
109        p *= binomialpower( Variable(i), -a[i], e[i] );
110    return p;
111}
112
113//static CanonicalForm check_dummy( const CFArray &a, const CFArray & P, const CFArray & Q )
114//{
115//    int i, r = a.size();
116//    CanonicalForm res, prod;
117//    res = 0;
118//    prod = 1;
119//    for ( i = 1; i <= r; i++ )
120//    {
121//        res += prod * a[i] * Q[i];
122//        prod *= P[i];
123//    }
124//    return res;
125//}
126
127static bool check_e( const IteratedFor & e, int k, int m, int * n )
128{
129    int sum = 0;
130    for ( int i = 2; i < k; i++ )
131    {
132        sum += e[i];
133        if ( e[i] > n[i] )
134            return false;
135    }
136    return sum == m+1;
137}
138
139static CanonicalForm modDeltak ( const CanonicalForm & f, const Evaluation & A, int k, int * n )
140{
141    CanonicalForm result = f;
142    for ( int i = 2; i < k; i++ )
143    {
144        result.mod( binomialpower( Variable(i), -A[i], n[i]+1 ) );
145    }
146    return result;
147}
148
149static CFArray
150findCorrCoeffs ( 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 )
151{
152    DEBINCLEVEL( cerr, "findCorrCoeffs" );
153    int i, m;
154    CFArray A(1,r), a(1,r);
155    CanonicalForm C0, Dm, g, prodcomb;
156
157    TIMING_START(fac_solve);
158    C0 = pk( I( C, 2, k-1 ), true );
159    solveF( P0, Q0, S, T, 1, pk, r, a );
160    TIMING_END(fac_solve);
161
162    DEBOUTLN( cerr, "trying to find correction coefficients for " << C );
163    DEBOUTLN( cerr, "which evaluates to " << C0 );
164
165    for ( i = 1; i <= r; i++ )
166        A[i] = remainder( pk( a[i] * C0 ), P0[i], pk );
167    DEBOUTLN( cerr, "the first approximation of the correction coefficients is " << A );
168/*#ifdef DEBUGOUTPUT
169    if ( check_dummy( A, P, Q ) - C != 0 )
170    {
171        DEBOUTLN( cerr, "there is an error detected, the correction coefficients do not" );
172        DEBOUTLN( cerr, "fulfill equation F(A)" );
173        DEBOUTLN( cerr, "corresponding P " << P );
174        DEBOUTLN( cerr, "              Q " << Q );
175    }
176#endif*/
177    for ( m = 0; m <= h && ( m == 0 || Dm != 0 ); m++ )
178    {
179        Dm = pk( evalF( P, Q, A, r ) - C );
180        if ( Dm != 0 )
181        {
182            if ( k == 2 )
183            {
184                TIMING_START(fac_solve);
185                solveF( P0, Q0, S, T, Dm, pk, r, a );
186                TIMING_END(fac_solve);
187                for ( i = 1; i <= r; i++ )
188                    A[i] -= a[i];
189            }
190            else
191            {
192                IteratedFor e( 2, k-1, m+1 );
193                while ( e.iterations_left() )
194                {
195                    if ( check_e( e, k, m, n ) )
196                    {
197                        g = pk( checkCombination( Dm, I, e, k ) );
198                        if ( ! g.isZero() && ! (g.mvar() > Variable(1)) )
199                        {
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                            {
207//                                A[i] -= a[i] * prodcomb;
208                                A[i] = pk( A[i] - a[i] * prodcomb );
209                            }
210                        }
211                    }
212                    e++;
213                }
214            }
215        }
216        DEBOUTLN( cerr, "the correction coefficients at step " << m );
217        DEBOUTLN( cerr, "are now " << A );
218/*#ifdef DEBUGOUTPUT
219    if ( check_dummy( A, P, Q ) - C != 0 ) {
220        DEBOUTLN( cerr, "there is an error detected, the correction coefficients do not" );
221        DEBOUTLN( cerr, "fulfill equation F(A)" );
222        DEBOUTLN( cerr, "corresponding P " << P );
223        DEBOUTLN( cerr, "              Q " << Q );
224    }
225#endif*/
226    }
227    DEBDECLEVEL( cerr, "findCorrCoeffs" );
228    return A;
229}
230
231
232static bool
233liftStep ( CFArray & P, int k, int r, int t, const modpk & b, const Evaluation & A, const CFArray & lcG, const CanonicalForm & Uk, int * n, int h )
234{
235    DEBINCLEVEL( cerr, "liftStep" );
236    CFArray K( 1, r ), Q( 1, r ), Q0( 1, r ), P0( 1, r ), S( 1, r ), T( 1, r ), alpha( 1, r );
237    CanonicalForm Rm, C, D, xa = Variable(k) - A[k];
238    CanonicalForm xa1 = xa, xa2 = xa*xa;
239    CanonicalForm dummy;
240    int i, m;
241
242    DEBOUTLN( cerr, "we are now performing the liftstep to reach " << Variable(k) );
243    DEBOUTLN( cerr, "the factors so far are " << P );
244    DEBOUTLN( cerr, "modulus p^k= " << b.getpk() << "=" << b.getp() <<"^"<< b.getk() );
245
246    for ( i = 1; i <= r; i++ )
247    {
248        Variable vm = Variable( t + 1 );
249        Variable v1 = Variable(1);
250        K[i] = swapvar( replaceLc( swapvar( P[i], v1, vm ), swapvar( A( lcG[i], k+1, t ), v1, vm ) ), v1, vm );
251        P[i] = A( K[i], k, t );
252    }
253    DEBOUTLN( cerr, "lift K = " << K );
254
255//    d = degree( Uk, Variable( k ) );
256
257    TIMING_START(fac_extgcd);
258    Q[r] = 1;
259    for ( i = r; i > 1; i-- )
260    {
261        Q[i-1] = Q[i] * P[i];
262        P0[i] = A( P[i], 2, k-1 );
263        Q0[i] = A( Q[i], 2, k-1 );
264        extgcd( P0[i], Q0[i], S[i], T[i], b );
265    }
266    P0[1] = A( P[1], 2, k-1 );
267    Q0[1] = A( Q[1], 2, k-1 );
268    extgcd( P0[1], Q0[1], S[1], T[1], b );
269    TIMING_END(fac_extgcd);
270
271    for ( m = 1; m <= n[k]+1; m++ )
272    {
273        TIMING_START(fac_modpk);
274        Rm = modDeltak( prod( K ) - Uk, A, k, n );
275        TIMING_END(fac_modpk);
276#ifdef DEBUGOUTPUT
277        if ( mod( Rm, xa1 ) != 0 )
278        {
279            DEBOUTLN( cerr, "something seems not to be ok with Rm which is " << Rm );
280            DEBOUTLN( cerr, "and should reduce to zero modulo " << xa1 );
281        }
282#endif
283        if ( mod( Rm, xa2 ) != 0 )
284        {
285            C = derivAndEval( Rm, m, Variable( k ), A[k] );
286            D = 1;
287            for ( i = 2; i <= m; i++ ) D *= i;
288            C /= D;
289
290            TIMING_START(fac_corrcoeff);
291            alpha = findCorrCoeffs( P, Q, P0, Q0, S, T, C, A, b, r, k, h, n ); // -> h berechnen
292            TIMING_END(fac_corrcoeff);
293// #ifdef DEBUGOUTPUT
294//             dummy = check_dummy( alpha, P, Q );
295//             if ( b(modDeltak( dummy, A, k, n )) != b(modDeltak( C, A, k, n )) )
296//             {
297//                 DEBOUTLN( cerr, "lift fault" );
298//                 DEBDECLEVEL( cerr, "liftStep" );
299//                 return false;
300//             }
301// #endif
302            for ( i = 1; i <= r; i++ )
303                K[i] = b(K[i] - alpha[i] * xa1);
304            DEBOUTLN( cerr, "the corrected K's are now " << K );
305        }
306        xa1 = xa2;
307        xa2 *= xa;
308    }
309    for ( i = 1; i <= r; i++ )
310        P[i] = K[i];
311    if ( prod( K ) - Uk != 0 )
312    {
313        DEBOUTLN( cerr, "the liftstep produced the wrong result" );
314        DEBOUTLN( cerr, "the product of the factors calculated so far is " << prod(K) );
315        DEBOUTLN( cerr, "and the Uk that should have been reached is " << Uk );
316        DEBDECLEVEL( cerr, "liftStep" );
317        return false;
318    }
319    DEBOUTLN( cerr, "the lift seems ok so far" );
320    DEBDECLEVEL( cerr, "liftStep" );
321    return true; // check for divisibility
322}
323
324bool
325Hensel ( const CanonicalForm & U, CFArray & G, const CFArray & lcG, const Evaluation & A, const modpk & bound, const Variable & /*x*/ )
326{
327    DEBINCLEVEL( cerr, "Hensel" );
328    int k, i, h, t = A.max();
329    bool goodeval = true;
330    CFArray Uk( A.min(), A.max() );
331    int * n = new int[t+1];
332
333    Uk[t] = U;
334    for ( k = t-1; k > 1; k-- )
335    {
336        Uk[k] = Uk[k+1]( A[k+1], Variable( k+1 ) );
337        n[k] = degree( Uk[k], Variable( k ) );
338    }
339    for ( k = A.min(); goodeval && (k <= t); k++ )
340    {
341        h = totaldegree( Uk[k], Variable(A.min()), Variable(k-1) );
342        for ( i = A.min(); i <= k; i++ )
343            n[i] = degree( Uk[k], Variable(i) );
344        goodeval = liftStep( G, k, G.max(), t, bound, A, lcG, Uk[k], n, h );
345        DEBOUTLN( cerr, "Factorization so far: " << G );
346    }
347    DEBDECLEVEL( cerr, "Hensel" );
348    delete[] n;
349    return goodeval;
350}
Note: See TracBrowser for help on using the repository browser.