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

spielwiese
Last change on this file since 564dd8 was 564dd8, checked in by Hans Schoenemann <hannes@…>, 3 years ago
factorize in Z[x,..] w/o NTL
  • Property mode set to 100644
File size: 11.1 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_algorithm.h"
13#include "fac_multivar.h"
14#include "fac_univar.h"
15#include "cf_reval.h"
16#include "cf_map.h"
17#include "fac_util.h"
18#include "cf_binom.h"
19#include "cf_iter.h"
20#include "cf_primes.h"
21#include "fac_distrib.h"
22#include "fac_multihensel.h"
23
24#ifndef HAVE_NTL
25void out_cf(const char *s1,const CanonicalForm &f,const char *s2);
26void out_cff(CFFList &L);
27
28TIMING_DEFINE_PRINT(fac_content)
29TIMING_DEFINE_PRINT(fac_findeval)
30TIMING_DEFINE_PRINT(fac_distrib)
31TIMING_DEFINE_PRINT(fac_hensel)
32
33static CFArray
34conv_to_factor_array( const CFFList & L )
35{
36    int n;
37    CFFListIterator I = L;
38    bool negate = false;
39
40    if ( ! I.hasItem() )
41        n = 0;
42    else  if ( I.getItem().factor().inBaseDomain() ) {
43        negate = I.getItem().factor().sign() < 0;
44        I++;
45        n = L.length();
46    }
47    else
48        n = L.length() + 1;
49    CFFListIterator J = I;
50    while ( J.hasItem() ) {
51        n += J.getItem().exp() - 1;
52        J++;
53    }
54    CFArray result( 1, n-1 );
55    int i, j, k;
56    i = 1;
57    while ( I.hasItem() ) {
58        k = I.getItem().exp();
59        for ( j = 1; j <= k; j++ ) {
60            result[i] = I.getItem().factor();
61            i++;
62        }
63        I++;
64    }
65    if ( negate )
66        result[1] = -result[1];
67    return result;
68}
69
70static modpk
71coeffBound ( const CanonicalForm & f, int p )
72{
73    int * degs = degrees( f );
74    int M = 0, i, k = f.level();
75    for ( i = 1; i <= k; i++ )
76        M += degs[i];
77    CanonicalForm b = 2 * maxNorm( f ) * power( CanonicalForm( 3 ), M );
78    CanonicalForm B = p;
79    k = 1;
80    while ( B < b ) {
81        B *= p;
82        k++;
83    }
84    return modpk( p, k );
85}
86
87// static bool
88// nonDivisors ( CanonicalForm omega, CanonicalForm delta, const CFArray & F, CFArray & d )
89// {
90//     DEBOUTLN( cerr, "nondivisors omega = " << omega );
91//     DEBOUTLN( cerr, "nondivisors delta = " << delta );
92//     DEBOUTLN( cerr, "nondivisors F = " << F );
93//     CanonicalForm q, r;
94//     int k = F.size();
95//     d = CFArray( 0, k );
96//     d[0] = delta * omega;
97//     for ( int i = 1; i <= k; i++ ) {
98//         q = abs(F[i]);
99//         for ( int j = i-1; j >= 0; j-- ) {
100//             r = d[j];
101//             do {
102//                 r = gcd( r, q );
103//                 q = q / r;
104//             } while ( r != 1 );
105//             if ( q == 1 )
106//                 return false;
107//         }
108//         d[i] = q;
109//     }
110//     return true;
111// }
112
113static void
114findEvaluation ( const CanonicalForm & U, const CanonicalForm & V, const CanonicalForm & omega, const CFFList & F, Evaluation & A, CanonicalForm & U0, CanonicalForm & delta, CFArray & D, int r )
115{
116    DEBINCLEVEL( cerr, "findEvaluation" );
117    CanonicalForm Vn;
118    CFFListIterator I;
119    int j;
120    bool found = false;
121    CFArray FF = CFArray( 1, F.length() );
122    if ( r > 0 )
123        A.nextpoint();
124    while ( ! found )
125    {
126        Vn = A( V );
127        if ( Vn != 0 )
128        {
129            U0 = A( U );
130            DEBOUTLN( cerr, "U0 = " << U0 );
131            if ( isSqrFree( U0 ) )
132            {
133                delta = content( U0 );
134                DEBOUTLN( cerr, "content( U0 ) = " << delta );
135                for ( I = F, j = 1; I.hasItem(); I++, j++ )
136                    FF[j] = A( I.getItem().factor() );
137                found = nonDivisors( omega, delta, FF, D );
138            }
139            else
140            {
141                DEBOUTLN( cerr, "not sqrfree : " << sqrFree( U0 ) );
142            }
143        }
144        if ( ! found )
145            A.nextpoint();
146    }
147    DEBDECLEVEL( cerr, "findEvaluation" );
148}
149
150#ifdef HAVE_NTL
151VAR int prime_number=0;
152void find_good_prime(const CanonicalForm &f, int &start)
153{
154  if (! f.inBaseDomain() )
155  {
156    CFIterator i = f;
157    for(;;)
158    {
159      if  ( i.hasTerms() )
160      {
161        find_good_prime(i.coeff(),start);
162        if (0==cf_getSmallPrime(start)) return;
163        if((i.exp()!=0) && ((i.exp() % cf_getSmallPrime(start))==0))
164        {
165          start++;
166          i=f;
167        }
168        else  i++;
169      }
170      else break;
171    }
172  }
173  else
174  {
175    if (f.inZ())
176    {
177      if (0==cf_getSmallPrime(start)) return;
178      while((!f.isZero()) && (mod(f,cf_getSmallPrime(start))==0))
179      {
180        start++;
181        if (0==cf_getSmallPrime(start)) return;
182      }
183    }
184/* should not happen!
185    else if (f.inQ())
186    {
187      while((f.den()!=0) && (mod(f.den(),cf_getSmallPrime(start))==0))
188      {
189        start++;
190      }
191      while((f.num()!=0) && (mod(f.num(),cf_getSmallPrime(start))==0))
192      {
193        start++;
194      }
195    }
196    else
197          cout <<"??"<< f <<"\n";
198*/
199  }
200}
201#endif
202
203static CFArray ZFactorizeMulti ( const CanonicalForm & arg )
204{
205    DEBINCLEVEL( cerr, "ZFactorizeMulti" );
206    CFMap M;
207    CanonicalForm UU, U = compress( arg, M );
208    CanonicalForm delta, omega, V = LC( U, 1 );
209    int t = U.level();
210    CFFList F = factorize( V );
211    CFFListIterator I, J;
212    CFArray G, lcG, D;
213    int i, j, r, maxdeg;
214    REvaluation A( 2, t, IntRandom( 50 ) );
215    CanonicalForm U0;
216    CanonicalForm ft, ut, gt, d;
217    modpk b;
218    bool negate = false;
219
220    DEBOUTLN( cerr, "-----------------------------------------------------" );
221    DEBOUTLN( cerr, "trying to factorize U = " << U );
222    DEBOUTLN( cerr, "U is a polynomial of level = " << arg.level() );
223    DEBOUTLN( cerr, "U will be factorized with respect to variable " << Variable(1) );
224    DEBOUTLN( cerr, "the leading coefficient of U with respect to that variable is " << V );
225    DEBOUTLN( cerr, "which is factorized as " << F );
226
227    maxdeg = 0;
228    for ( i = 2; i <= t; i++ )
229    {
230        j = U.degree( Variable( i ) );
231        if ( j > maxdeg ) maxdeg = j;
232    }
233
234    if ( F.getFirst().factor().inCoeffDomain() )
235    {
236        omega = F.getFirst().factor();
237        F.removeFirst();
238        if ( omega < 0 )
239        {
240            negate = true;
241            omega = -omega;
242            U = -U;
243        }
244    }
245    else
246        omega = 1;
247
248    bool goodeval = false;
249    r = 0;
250
251//    for ( i = 0; i < 10*t; i++ )
252//        A.nextpoint();
253
254    while ( ! goodeval )
255    {
256        TIMING_START(fac_findeval);
257        findEvaluation( U, V, omega, F, A, U0, delta, D, r );
258        TIMING_END(fac_findeval);
259        DEBOUTLN( cerr, "the evaluation point to reduce to an univariate problem is " << A );
260        DEBOUTLN( cerr, "corresponding delta = " << delta );
261        DEBOUTLN( cerr, "              omega = " << omega );
262        DEBOUTLN( cerr, "              D     = " << D );
263        DEBOUTLN( cerr, "now factorize the univariate polynomial " << U0 );
264        G = conv_to_factor_array( factorize( U0, false ) );
265        DEBOUTLN( cerr, "which factorizes into " << G );
266        #ifdef HAVE_NTL
267        {
268          int i=prime_number;
269          find_good_prime(arg,i);
270          find_good_prime(U0,i);
271          find_good_prime(U,i);
272          int p=cf_getSmallPrime(i);
273          //printf("found:p=%d (%d)\n",p,i);
274          if (p==0)
275          {
276            return conv_to_factor_array(CFFactor(arg,1));
277            //printf("out of primes - switch ot non-NTL\n");
278          }
279          else if (((i==0)||(i!=prime_number)))
280          {
281            b = coeffBound( U, p );
282            prime_number=i;
283          }
284          // p!=0:
285          modpk bb=coeffBound(U0,p);
286          if (bb.getk() > b.getk() ) b=bb;
287          bb=coeffBound(arg,p);
288          if (bb.getk() > b.getk() ) b=bb;
289        }
290        #else
291        b = coeffBound( U, getZFacModulus().getp() );
292        if ( getZFacModulus().getpk() > b.getpk() )
293            b = getZFacModulus();
294        #endif
295        //printf("p=%d, k=%d\n",b.getp(),b.getk());
296        DEBOUTLN( cerr, "the coefficient bound of the factors of U is " << b.getpk() );
297
298        r = G.size();
299        lcG = CFArray( 1, r );
300        UU = U;
301        DEBOUTLN( cerr, "now trying to distribute the leading coefficients ..." );
302        TIMING_START(fac_distrib);
303        goodeval = distributeLeadingCoeffs( UU, G, lcG, F, D, delta, omega, A, r );
304        TIMING_END(fac_distrib);
305#ifdef DEBUGOUTPUT
306        if ( goodeval )
307        {
308            DEBOUTLN( cerr, "the univariate factors after distribution are " << G );
309            DEBOUTLN( cerr, "the distributed leading coeffs are " << lcG );
310            DEBOUTLN( cerr, "U may have changed and is now " << UU );
311            DEBOUTLN( cerr, "which has leading coefficient " << LC( UU, Variable(1) ) );
312
313            if ( LC( UU, Variable(1) ) != prod( lcG ) || A(UU) != prod( G ) )
314            {
315                DEBOUTLN( cerr, "!!! distribution was not correct !!!" );
316                DEBOUTLN( cerr, "product of leading coeffs is " << prod( lcG ) );
317                DEBOUTLN( cerr, "product of univariate factors is " << prod( G ) );
318                DEBOUTLN( cerr, "the new U is evaluated as " << A(UU) );
319            }
320            else
321                DEBOUTLN( cerr, "leading coeffs correct" );
322        }
323        else
324        {
325            DEBOUTLN( cerr, "we have found a bad evaluation point" );
326        }
327#endif
328        if ( goodeval )
329        {
330            TIMING_START(fac_hensel);
331            goodeval = Hensel( UU, G, lcG, A, b, Variable(1) );
332            TIMING_END(fac_hensel);
333        }
334    }
335    for ( i = 1; i <= r; i++ )
336    {
337        G[i] /= icontent( G[i] );
338        G[i] = M(G[i]);
339    }
340    // negate noch beachten !
341    if ( negate )
342        G[1] = -G[1];
343    DEBDECLEVEL( cerr, "ZFactorMulti" );
344    return G;
345}
346
347CFFList ZFactorizeMultivariate ( const CanonicalForm & f, bool issqrfree )
348{
349    CFFList G, F, R;
350    CFArray GG;
351    CFFListIterator i, j;
352    CFMap M;
353    CanonicalForm g, cont;
354    Variable v1, vm;
355    int k, m, n;
356
357    v1 = Variable(1);
358    if ( issqrfree )
359        F = CFFactor( f, 1 );
360    else
361        F = sqrFree( f );
362
363    for ( i = F; i.hasItem(); i++ )
364    {
365        if ( i.getItem().factor().inCoeffDomain() )
366            R.append( CFFactor( i.getItem().factor(), i.getItem().exp() ) );
367        else
368        {
369            TIMING_START(fac_content);
370            g = compress( i.getItem().factor(), M );
371            // now after compress g contains Variable(1)
372            vm = g.mvar();
373            g = swapvar( g, v1, vm );
374            cont = content( g );
375            g = swapvar( g / cont, v1, vm );
376            cont = swapvar( cont, v1, vm );
377            n = i.getItem().exp();
378            TIMING_END(fac_content);
379            DEBOUTLN( cerr, "now after content ..." );
380            if ( g.isUnivariate() )
381            {
382                G = factorize( g, true );
383                for ( j = G; j.hasItem(); j++ )
384                    if ( ! j.getItem().factor().isOne() )
385                        R.append( CFFactor( M( j.getItem().factor() ), n ) );
386            }
387            else
388            {
389                GG = ZFactorizeMulti( g );
390                m = GG.max();
391                for ( k = GG.min(); k <= m; k++ )
392                    if ( ! GG[k].isOne() )
393                        R.append( CFFactor( M( GG[k] ), n ) );
394            }
395            G = factorize( cont, true );
396            for ( j = G; j.hasItem(); j++ )
397                if ( ! j.getItem().factor().isOne() )
398                    R.append( CFFactor( M( j.getItem().factor() ), n ) );
399        }
400    }
401    return R;
402}
403#endif
Note: See TracBrowser for help on using the repository browser.