source: git/factory/fac_multivar.cc @ 96ce32

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