source: git/factory/fac_multivar.cc @ ba5e9e

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