source: git/factory/fac_multihensel.cc @ 4a7a45

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