source: git/factory/fac_multivar.cc @ 464b18

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