source: git/factory/fac_multihensel.cc @ 79592ac

spielwiese
Last change on this file since 79592ac was 362fc67, checked in by Martin Lee <martinlee84@…>, 12 years ago
chg: remove $Id$
  • Property mode set to 100644
File size: 10.7 KB
Line 
1/* emacs edit mode for this file is -*- C++ -*- */
2
3#include "config.h"
4
5#include "cf_assert.h"
6#include "debug.h"
7#include "timing.h"
8
9#include "cf_defs.h"
10#include "cf_eval.h"
11#include "cf_binom.h"
12#include "fac_util.h"
13#include "fac_iterfor.h"
14#include "cf_iter.h"
15
16
17TIMING_DEFINE_PRINT(fac_solve)
18TIMING_DEFINE_PRINT(fac_modpk)
19TIMING_DEFINE_PRINT(fac_corrcoeff)
20TIMING_DEFINE_PRINT(fac_extgcd)
21
22static void
23extgcdrest ( const CanonicalForm & a, const CanonicalForm & b, const CanonicalForm & s, const CanonicalForm & t, const CanonicalForm & c, CanonicalForm & S, CanonicalForm & T, const modpk & /*pk*/ )
24{
25    CanonicalForm sigma = s * c, tau = t * c;
26//    divremainder( sigma, b, T, S, pk );
27//    T = pk( tau + T * a );
28    divrem( sigma, b, T, S );
29    T = tau + T * a;
30}
31
32static void
33solveF ( const CFArray & P, const CFArray & Q, const CFArray & S, const CFArray & T, const CanonicalForm & C, const modpk & pk, int r, CFArray & a )
34{
35    setCharacteristic( pk.getp(), pk.getk() );
36    CanonicalForm g, bb, b = mapinto( C );
37    int j;
38    for ( j = 1; j < r; j++ )
39    {
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    {
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    {
70        CFIterator i;
71        CanonicalForm sum = 0, fact;
72        int min, j;
73        Variable v = Variable( f.level() + 1 );
74        for ( i = swapvar( f, x, v); i.hasTerms() && i.exp() >= n; i++ )
75        {
76            fact = 1;
77            min = i.exp() - n;
78            for ( j = i.exp(); j > min; j-- )
79                fact *= j;
80            sum += fact * i.coeff() * power( v, min );
81        }
82        return sum( a, v );
83    }
84}
85
86static CanonicalForm
87checkCombination ( const CanonicalForm & W, const Evaluation & a, const IteratedFor & e, int k )
88{
89    CanonicalForm dW = W;
90    int i, j;
91    for ( i = k-1; i >= 2 && ! dW.isZero(); i-- )
92        dW = derivAndEval( dW, e[i], Variable( i ), a[i] );
93    if ( ! dW.isZero() ) {
94        CanonicalForm fact = 1;
95        for ( i = 2; i < k; i++ )
96            for ( j = 2; j <= e[i]; j++ )
97                fact *= j;
98        dW /= fact;
99    }
100    return dW;
101}
102
103static CanonicalForm
104prodCombination ( const Evaluation & a, const IteratedFor & e, int k )
105{
106    CanonicalForm p = 1;
107    for ( int i = k-1; i > 1; i-- )
108        p *= binomialpower( Variable(i), -a[i], e[i] );
109    return p;
110}
111
112//static CanonicalForm check_dummy( const CFArray &a, const CFArray & P, const CFArray & Q )
113//{
114//    int i, r = a.size();
115//    CanonicalForm res, prod;
116//    res = 0;
117//    prod = 1;
118//    for ( i = 1; i <= r; i++ )
119//    {
120//        res += prod * a[i] * Q[i];
121//        prod *= P[i];
122//    }
123//    return res;
124//}
125
126static bool check_e( const IteratedFor & e, int k, int m, int * n )
127{
128    int sum = 0;
129    for ( int i = 2; i < k; i++ )
130    {
131        sum += e[i];
132        if ( e[i] > n[i] )
133            return false;
134    }
135    return sum == m+1;
136}
137
138static CanonicalForm modDeltak ( const CanonicalForm & f, const Evaluation & A, int k, int * n )
139{
140    CanonicalForm result = f;
141    for ( int i = 2; i < k; i++ )
142    {
143        result.mod( binomialpower( Variable(i), -A[i], n[i]+1 ) );
144    }
145    return result;
146}
147
148static CFArray
149findCorrCoeffs ( 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 )
150{
151    DEBINCLEVEL( cerr, "findCorrCoeffs" );
152    int i, m;
153    CFArray A(1,r), a(1,r);
154    CanonicalForm C0, Dm, g, prodcomb;
155
156    TIMING_START(fac_solve);
157    C0 = pk( I( C, 2, k-1 ), true );
158    solveF( P0, Q0, S, T, 1, pk, r, a );
159    TIMING_END(fac_solve);
160
161    DEBOUTLN( cerr, "trying to find correction coefficients for " << C );
162    DEBOUTLN( cerr, "which evaluates to " << C0 );
163
164    for ( i = 1; i <= r; i++ )
165        A[i] = remainder( pk( a[i] * C0 ), P0[i], pk );
166    DEBOUTLN( cerr, "the first approximation of the correction coefficients is " << A );
167/*#ifdef DEBUGOUTPUT
168    if ( check_dummy( A, P, Q ) - C != 0 )
169    {
170        DEBOUTLN( cerr, "there is an error detected, the correction coefficients do not" );
171        DEBOUTLN( cerr, "fulfill equation F(A)" );
172        DEBOUTLN( cerr, "corresponding P " << P );
173        DEBOUTLN( cerr, "              Q " << Q );
174    }
175#endif*/
176    for ( m = 0; m <= h && ( m == 0 || Dm != 0 ); m++ )
177    {
178        Dm = pk( evalF( P, Q, A, r ) - C );
179        if ( Dm != 0 )
180        {
181            if ( k == 2 )
182            {
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            {
191                IteratedFor e( 2, k-1, m+1 );
192                while ( e.iterations_left() )
193                {
194                    if ( check_e( e, k, m, n ) )
195                    {
196                        g = pk( checkCombination( Dm, I, e, k ) );
197                        if ( ! g.isZero() && ! (g.mvar() > Variable(1)) )
198                        {
199                            prodcomb = prodCombination( I, e, k );
200//                            Dm = Dm - g * prodcomb;
201                            TIMING_START(fac_solve);
202                            solveF( P0, Q0, S, T, g, pk, r, a );
203                            TIMING_END(fac_solve);
204                            for ( i = 1; i <= r; i++ )
205                            {
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    DEBOUTLN( cerr, "modulus p^k= " << b.getpk() << "=" << b.getp() <<"^"<< b.getk() );
244
245    for ( i = 1; i <= r; i++ )
246    {
247        Variable vm = Variable( t + 1 );
248        Variable v1 = Variable(1);
249        K[i] = swapvar( replaceLc( swapvar( P[i], v1, vm ), swapvar( A( lcG[i], k+1, t ), v1, vm ) ), v1, vm );
250        P[i] = A( K[i], k, t );
251    }
252    DEBOUTLN( cerr, "lift K = " << K );
253
254//    d = degree( Uk, Variable( k ) );
255
256    TIMING_START(fac_extgcd);
257    Q[r] = 1;
258    for ( i = r; i > 1; i-- )
259    {
260        Q[i-1] = Q[i] * P[i];
261        P0[i] = A( P[i], 2, k-1 );
262        Q0[i] = A( Q[i], 2, k-1 );
263        extgcd( P0[i], Q0[i], S[i], T[i], b );
264    }
265    P0[1] = A( P[1], 2, k-1 );
266    Q0[1] = A( Q[1], 2, k-1 );
267    extgcd( P0[1], Q0[1], S[1], T[1], b );
268    TIMING_END(fac_extgcd);
269
270    for ( m = 1; m <= n[k]+1; m++ )
271    {
272        TIMING_START(fac_modpk);
273        Rm = modDeltak( prod( K ) - Uk, A, k, n );
274        TIMING_END(fac_modpk);
275#ifdef DEBUGOUTPUT
276        if ( mod( Rm, xa1 ) != 0 )
277        {
278            DEBOUTLN( cerr, "something seems not to be ok with Rm which is " << Rm );
279            DEBOUTLN( cerr, "and should reduce to zero modulo " << xa1 );
280        }
281#endif
282        if ( mod( Rm, xa2 ) != 0 )
283        {
284            C = derivAndEval( Rm, m, Variable( k ), A[k] );
285            D = 1;
286            for ( i = 2; i <= m; i++ ) D *= i;
287            C /= D;
288
289            TIMING_START(fac_corrcoeff);
290            alpha = findCorrCoeffs( P, Q, P0, Q0, S, T, C, A, b, r, k, h, n ); // -> h berechnen
291            TIMING_END(fac_corrcoeff);
292// #ifdef DEBUGOUTPUT
293//             dummy = check_dummy( alpha, P, Q );
294//             if ( b(modDeltak( dummy, A, k, n )) != b(modDeltak( C, A, k, n )) )
295//             {
296//                 DEBOUTLN( cerr, "lift fault" );
297//                 DEBDECLEVEL( cerr, "liftStep" );
298//                 return false;
299//             }
300// #endif
301            for ( i = 1; i <= r; i++ )
302                K[i] = b(K[i] - alpha[i] * xa1);
303            DEBOUTLN( cerr, "the corrected K's are now " << K );
304        }
305        xa1 = xa2;
306        xa2 *= xa;
307    }
308    for ( i = 1; i <= r; i++ )
309        P[i] = K[i];
310    if ( prod( K ) - Uk != 0 )
311    {
312        DEBOUTLN( cerr, "the liftstep produced the wrong result" );
313        DEBOUTLN( cerr, "the product of the factors calculated so far is " << prod(K) );
314        DEBOUTLN( cerr, "and the Uk that should have been reached is " << Uk );
315        DEBDECLEVEL( cerr, "liftStep" );
316        return false;
317    }
318    DEBOUTLN( cerr, "the lift seems ok so far" );
319    DEBDECLEVEL( cerr, "liftStep" );
320    return true; // check for divisibility
321}
322
323bool
324Hensel ( const CanonicalForm & U, CFArray & G, const CFArray & lcG, const Evaluation & A, const modpk & bound, const Variable & /*x*/ )
325{
326    DEBINCLEVEL( cerr, "Hensel" );
327    int k, i, h, t = A.max();
328    bool goodeval = true;
329    CFArray Uk( A.min(), A.max() );
330    int * n = new int[t+1];
331
332    Uk[t] = U;
333    for ( k = t-1; k > 1; k-- )
334    {
335        Uk[k] = Uk[k+1]( A[k+1], Variable( k+1 ) );
336        n[k] = degree( Uk[k], Variable( k ) );
337    }
338    for ( k = A.min(); goodeval && (k <= t); k++ )
339    {
340        h = totaldegree( Uk[k], Variable(A.min()), Variable(k-1) );
341        for ( i = A.min(); i <= k; i++ )
342            n[i] = degree( Uk[k], Variable(i) );
343        goodeval = liftStep( G, k, G.max(), t, bound, A, lcG, Uk[k], n, h );
344        DEBOUTLN( cerr, "Factorization so far: " << G );
345    }
346    DEBDECLEVEL( cerr, "Hensel" );
347    delete[] n;
348    return goodeval;
349}
Note: See TracBrowser for help on using the repository browser.