source: git/factory/fac_multihensel.cc @ 564dd8

spielwiese
Last change on this file since 564dd8 was 564dd8, checked in by Hans Schoenemann <hannes@…>, 4 years ago
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
3
4#include "config.h"
5
6
7#include "cf_assert.h"
8#include "debug.h"
9#include "timing.h"
10
11#include "cf_defs.h"
12#include "cf_eval.h"
13#include "cf_binom.h"
14#include "fac_util.h"
15#include "fac_iterfor.h"
16#include "fac_multihensel.h"
17#include "cf_iter.h"
18
19#ifndef HAVE_NTL
20
21TIMING_DEFINE_PRINT(fac_solve)
22TIMING_DEFINE_PRINT(fac_modpk)
23TIMING_DEFINE_PRINT(fac_corrcoeff)
24TIMING_DEFINE_PRINT(fac_extgcd)
25
26static void
27extgcdrest ( const CanonicalForm & a, const CanonicalForm & b, const CanonicalForm & s, const CanonicalForm & t, const CanonicalForm & c, CanonicalForm & S, CanonicalForm & T, const modpk & /*pk*/ )
28{
29    CanonicalForm sigma = s * c, tau = t * c;
30//    divremainder( sigma, b, T, S, pk );
31//    T = pk( tau + T * a );
32    divrem( sigma, b, T, S );
33    T = tau + T * a;
34}
35
36static void
37solveF ( const CFArray & P, const CFArray & Q, const CFArray & S, const CFArray & T, const CanonicalForm & C, const modpk & pk, int r, CFArray & a )
38{
39    setCharacteristic( pk.getp(), pk.getk() );
40    CanonicalForm g, bb, b = mapinto( C );
41    int j;
42    for ( j = 1; j < r; j++ )
43    {
44        extgcdrest( mapinto( P[j] ), mapinto( Q[j] ), mapinto( S[j] ), mapinto( T[j] ), b, bb, a[j], pk );
45        b = bb;
46    }
47    a[r] = b;
48    setCharacteristic( 0 );
49    for ( j = 1; j <= r; j++ )
50        a[j] = mapinto( a[j] );
51}
52
53static CanonicalForm
54evalF ( const CFArray & P, const CFArray & Q, const CFArray & A, int r )
55{
56    CanonicalForm pprod = 1, sum = 0;
57    for ( int i = 1; i <= r; i++ )
58    {
59        sum += pprod * A[i] * Q[i];
60        pprod *= P[i];
61    }
62    return sum;
63}
64
65static CanonicalForm
66derivAndEval ( const CanonicalForm & f, int n, const Variable & x, const CanonicalForm & a )
67{
68    if ( n == 0 )
69        return f( a, x );
70    else if ( f.degree( x ) < n )
71        return 0;
72    else
73    {
74        CFIterator i;
75        CanonicalForm sum = 0, fact;
76        int min, j;
77        Variable v = Variable( f.level() + 1 );
78        for ( i = swapvar( f, x, v); i.hasTerms() && i.exp() >= n; i++ )
79        {
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
116//static CanonicalForm check_dummy( const CFArray &a, const CFArray & P, const CFArray & Q )
117//{
118//    int i, r = a.size();
119//    CanonicalForm res, prod;
120//    res = 0;
121//    prod = 1;
122//    for ( i = 1; i <= r; i++ )
123//    {
124//        res += prod * a[i] * Q[i];
125//        prod *= P[i];
126//    }
127//    return res;
128//}
129
130static bool check_e( const IteratedFor & e, int k, int m, int * n )
131{
132    int sum = 0;
133    for ( int i = 2; i < k; i++ )
134    {
135        sum += e[i];
136        if ( e[i] > n[i] )
137            return false;
138    }
139    return sum == m+1;
140}
141
142static CanonicalForm modDeltak ( const CanonicalForm & f, const Evaluation & A, int k, int * n )
143{
144    CanonicalForm result = f;
145    for ( int i = 2; i < k; i++ )
146    {
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    {
174        DEBOUTLN( cerr, "there is an error detected, the correction coefficients do not" );
175        DEBOUTLN( cerr, "fulfill equation F(A)" );
176        DEBOUTLN( cerr, "corresponding P " << P );
177        DEBOUTLN( cerr, "              Q " << Q );
178    }
179#endif*/
180    for ( m = 0; m <= h && ( m == 0 || Dm != 0 ); m++ )
181    {
182        Dm = pk( evalF( P, Q, A, r ) - C );
183        if ( Dm != 0 )
184        {
185            if ( k == 2 )
186            {
187                TIMING_START(fac_solve);
188                solveF( P0, Q0, S, T, Dm, pk, r, a );
189                TIMING_END(fac_solve);
190                for ( i = 1; i <= r; i++ )
191                    A[i] -= a[i];
192            }
193            else
194            {
195                IteratedFor e( 2, k-1, m+1 );
196                while ( e.iterations_left() )
197                {
198                    if ( check_e( e, k, m, n ) )
199                    {
200                        g = pk( checkCombination( Dm, I, e, k ) );
201                        if ( ! g.isZero() && ! (g.mvar() > Variable(1)) )
202                        {
203                            prodcomb = prodCombination( I, e, k );
204//                            Dm = Dm - g * prodcomb;
205                            TIMING_START(fac_solve);
206                            solveF( P0, Q0, S, T, g, pk, r, a );
207                            TIMING_END(fac_solve);
208                            for ( i = 1; i <= r; i++ )
209                            {
210//                                A[i] -= a[i] * prodcomb;
211                                A[i] = pk( A[i] - a[i] * prodcomb );
212                            }
213                        }
214                    }
215                    e++;
216                }
217            }
218        }
219        DEBOUTLN( cerr, "the correction coefficients at step " << m );
220        DEBOUTLN( cerr, "are now " << A );
221/*#ifdef DEBUGOUTPUT
222    if ( check_dummy( A, P, Q ) - C != 0 ) {
223        DEBOUTLN( cerr, "there is an error detected, the correction coefficients do not" );
224        DEBOUTLN( cerr, "fulfill equation F(A)" );
225        DEBOUTLN( cerr, "corresponding P " << P );
226        DEBOUTLN( cerr, "              Q " << Q );
227    }
228#endif*/
229    }
230    DEBDECLEVEL( cerr, "findCorrCoeffs" );
231    return A;
232}
233
234
235static bool
236liftStep ( CFArray & P, int k, int r, int t, const modpk & b, const Evaluation & A, const CFArray & lcG, const CanonicalForm & Uk, int * n, int h )
237{
238    DEBINCLEVEL( cerr, "liftStep" );
239    CFArray K( 1, r ), Q( 1, r ), Q0( 1, r ), P0( 1, r ), S( 1, r ), T( 1, r ), alpha( 1, r );
240    CanonicalForm Rm, C, D, xa = Variable(k) - A[k];
241    CanonicalForm xa1 = xa, xa2 = xa*xa;
242    CanonicalForm dummy;
243    int i, m;
244
245    DEBOUTLN( cerr, "we are now performing the liftstep to reach " << Variable(k) );
246    DEBOUTLN( cerr, "the factors so far are " << P );
247    DEBOUTLN( cerr, "modulus p^k= " << b.getpk() << "=" << b.getp() <<"^"<< b.getk() );
248
249    for ( i = 1; i <= r; i++ )
250    {
251        Variable vm = Variable( t + 1 );
252        Variable v1 = Variable(1);
253        K[i] = swapvar( replaceLc( swapvar( P[i], v1, vm ), swapvar( A( lcG[i], k+1, t ), v1, vm ) ), v1, vm );
254        P[i] = A( K[i], k, t );
255    }
256    DEBOUTLN( cerr, "lift K = " << K );
257
258//    d = degree( Uk, Variable( k ) );
259
260    TIMING_START(fac_extgcd);
261    Q[r] = 1;
262    for ( i = r; i > 1; i-- )
263    {
264        Q[i-1] = Q[i] * P[i];
265        P0[i] = A( P[i], 2, k-1 );
266        Q0[i] = A( Q[i], 2, k-1 );
267        extgcd( P0[i], Q0[i], S[i], T[i], b );
268    }
269    P0[1] = A( P[1], 2, k-1 );
270    Q0[1] = A( Q[1], 2, k-1 );
271    extgcd( P0[1], Q0[1], S[1], T[1], b );
272    TIMING_END(fac_extgcd);
273
274    for ( m = 1; m <= n[k]+1; m++ )
275    {
276        TIMING_START(fac_modpk);
277        Rm = modDeltak( prod( K ) - Uk, A, k, n );
278        TIMING_END(fac_modpk);
279#ifdef DEBUGOUTPUT
280        if ( mod( Rm, xa1 ) != 0 )
281        {
282            DEBOUTLN( cerr, "something seems not to be ok with Rm which is " << Rm );
283            DEBOUTLN( cerr, "and should reduce to zero modulo " << xa1 );
284        }
285#endif
286        if ( mod( Rm, xa2 ) != 0 )
287        {
288            C = derivAndEval( Rm, m, Variable( k ), A[k] );
289            D = 1;
290            for ( i = 2; i <= m; i++ ) D *= i;
291            C /= D;
292
293            TIMING_START(fac_corrcoeff);
294            alpha = findCorrCoeffs( P, Q, P0, Q0, S, T, C, A, b, r, k, h, n ); // -> h berechnen
295            TIMING_END(fac_corrcoeff);
296// #ifdef DEBUGOUTPUT
297//             dummy = check_dummy( alpha, P, Q );
298//             if ( b(modDeltak( dummy, A, k, n )) != b(modDeltak( C, A, k, n )) )
299//             {
300//                 DEBOUTLN( cerr, "lift fault" );
301//                 DEBDECLEVEL( cerr, "liftStep" );
302//                 return false;
303//             }
304// #endif
305            for ( i = 1; i <= r; i++ )
306                K[i] = b(K[i] - alpha[i] * xa1);
307            DEBOUTLN( cerr, "the corrected K's are now " << K );
308        }
309        xa1 = xa2;
310        xa2 *= xa;
311    }
312    for ( i = 1; i <= r; i++ )
313        P[i] = K[i];
314    if ( prod( K ) - Uk != 0 )
315    {
316        DEBOUTLN( cerr, "the liftstep produced the wrong result" );
317        DEBOUTLN( cerr, "the product of the factors calculated so far is " << prod(K) );
318        DEBOUTLN( cerr, "and the Uk that should have been reached is " << Uk );
319        DEBDECLEVEL( cerr, "liftStep" );
320        return false;
321    }
322    DEBOUTLN( cerr, "the lift seems ok so far" );
323    DEBDECLEVEL( cerr, "liftStep" );
324    return true; // check for divisibility
325}
326
327bool
328Hensel ( const CanonicalForm & U, CFArray & G, const CFArray & lcG, const Evaluation & A, const modpk & bound, const Variable & /*x*/ )
329{
330    DEBINCLEVEL( cerr, "Hensel" );
331    int k, i, h, t = A.max();
332    bool goodeval = true;
333    CFArray Uk( A.min(), A.max() );
334    int * n = new int[t+1];
335
336    Uk[t] = U;
337    for ( k = t-1; k > 1; k-- )
338    {
339        Uk[k] = Uk[k+1]( A[k+1], Variable( k+1 ) );
340        n[k] = degree( Uk[k], Variable( k ) );
341    }
342    for ( k = A.min(); goodeval && (k <= t); k++ )
343    {
344        h = totaldegree( Uk[k], Variable(A.min()), Variable(k-1) );
345        for ( i = A.min(); i <= k; i++ )
346            n[i] = degree( Uk[k], Variable(i) );
347        goodeval = liftStep( G, k, G.max(), t, bound, A, lcG, Uk[k], n, h );
348        DEBOUTLN( cerr, "Factorization so far: " << G );
349    }
350    DEBDECLEVEL( cerr, "Hensel" );
351    delete[] n;
352    return goodeval;
353}
354#endif
Note: See TracBrowser for help on using the repository browser.