source: git/factory/fac_multivar.cc @ e89e56

spielwiese
Last change on this file since e89e56 was e89e56, checked in by Hans Schönemann <hannes@…>, 16 years ago
*hannes: fact.tst git-svn-id: file:///usr/local/Singular/svn/trunk@10621 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 15.5 KB
Line 
1/* emacs edit mode for this file is -*- C++ -*- */
2/* $Id: fac_multivar.cc,v 1.14 2008-03-17 17:44:04 Singular Exp $ */
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
22void out_cf(char *s1,const CanonicalForm &f,char *s2);
23void out_cff(CFFList &L);
24
25TIMING_DEFINE_PRINT(fac_content);
26TIMING_DEFINE_PRINT(fac_findeval);
27TIMING_DEFINE_PRINT(fac_distrib);
28TIMING_DEFINE_PRINT(fac_hensel);
29
30static CFArray
31conv_to_factor_array( const CFFList & L )
32{
33    int n;
34    CFFListIterator I = L;
35    bool negate = false;
36
37    if ( ! I.hasItem() )
38        n = 0;
39    else  if ( I.getItem().factor().inBaseDomain() ) {
40        negate = I.getItem().factor().sign() < 0;
41        I++;
42        n = L.length();
43    }
44    else
45        n = L.length() + 1;
46    CFFListIterator J = I;
47    while ( J.hasItem() ) {
48        n += J.getItem().exp() - 1;
49        J++;
50    }
51    CFArray result( 1, n-1 );
52    int i, j, k;
53    i = 1;
54    while ( I.hasItem() ) {
55        k = I.getItem().exp();
56        for ( j = 1; j <= k; j++ ) {
57            result[i] = I.getItem().factor();
58            i++;
59        }
60        I++;
61    }
62    if ( negate )
63        result[1] = -result[1];
64    return result;
65}
66
67static modpk
68coeffBound ( const CanonicalForm & f, int p )
69{
70    int * degs = degrees( f );
71    int M = 0, i, k = f.level();
72    for ( i = 1; i <= k; i++ )
73        M += degs[i];
74    CanonicalForm b = 2 * maxNorm( f ) * power( CanonicalForm( 3 ), M );
75    CanonicalForm B = p;
76    k = 1;
77    while ( B < b ) {
78        B *= p;
79        k++;
80    }
81    return modpk( p, k );
82}
83
84// static bool
85// nonDivisors ( CanonicalForm omega, CanonicalForm delta, const CFArray & F, CFArray & d )
86// {
87//     DEBOUTLN( cerr, "nondivisors omega = " << omega );
88//     DEBOUTLN( cerr, "nondivisors delta = " << delta );
89//     DEBOUTLN( cerr, "nondivisors F = " << F );
90//     CanonicalForm q, r;
91//     int k = F.size();
92//     d = CFArray( 0, k );
93//     d[0] = delta * omega;
94//     for ( int i = 1; i <= k; i++ ) {
95//         q = abs(F[i]);
96//         for ( int j = i-1; j >= 0; j-- ) {
97//             r = d[j];
98//             do {
99//                 r = gcd( r, q );
100//                 q = q / r;
101//             } while ( r != 1 );
102//             if ( q == 1 )
103//                 return false;
104//         }
105//         d[i] = q;
106//     }
107//     return true;
108// }
109
110static void
111findEvaluation ( const CanonicalForm & U, const CanonicalForm & V, const CanonicalForm & omega, const CFFList & F, Evaluation & A, CanonicalForm & U0, CanonicalForm & delta, CFArray & D, int r )
112{
113    DEBINCLEVEL( cerr, "findEvaluation" );
114    CanonicalForm Vn;
115    CFFListIterator I;
116    int j;
117    bool found = false;
118    CFArray FF = CFArray( 1, F.length() );
119    if ( r > 0 )
120        A.nextpoint();
121    while ( ! found )
122    {
123        Vn = A( V );
124        if ( Vn != 0 )
125        {
126            U0 = A( U );
127            DEBOUTLN( cerr, "U0 = " << U0 );
128            if ( isSqrFree( U0 ) )
129            {
130                delta = content( U0 );
131                DEBOUTLN( cerr, "content( U0 ) = " << delta );
132                for ( I = F, j = 1; I.hasItem(); I++, j++ )
133                    FF[j] = A( I.getItem().factor() );
134                found = nonDivisors( omega, delta, FF, D );
135            }
136            else
137            {
138                DEBOUTLN( cerr, "not sqrfree : " << sqrFree( U0 ) );
139            }
140        }
141        if ( ! found )
142            A.nextpoint();
143    }
144    DEBDECLEVEL( cerr, "findEvaluation" );
145}
146
147#ifdef HAVE_NTL
148int prime_number=0;
149void find_good_prime(const CanonicalForm &f, int &start)
150{
151  if (! f.inBaseDomain() )
152  {
153    int l = f.level();
154    CFIterator i = f;
155    for(;;)
156    {
157      if  ( i.hasTerms() )
158      {
159        find_good_prime(i.coeff(),start);
160        if((i.exp()!=0) && ((i.exp() % cf_getSmallPrime(start))==0))
161        {
162          start++;
163          i=f;
164        }
165        else  i++;
166      }
167      else break;
168    }
169  }
170  else
171  {
172    if (f.inZ())
173    {
174      while((f!=0) && (mod(f,cf_getSmallPrime(start))==0))
175      {
176        start++;
177      }
178    }
179/* should not happen!
180    else if (f.inQ())
181    {
182      while((f.den()!=0) && (mod(f.den(),cf_getSmallPrime(start))==0))
183      {
184        start++;
185      }
186      while((f.num()!=0) && (mod(f.num(),cf_getSmallPrime(start))==0))
187      {
188        start++;
189      }
190    }
191    else
192          cout <<"??"<< f <<"\n";
193*/
194  }
195}
196#endif
197
198static CFArray ZFactorizeMulti ( const CanonicalForm & arg )
199{
200    DEBINCLEVEL( cerr, "ZFactorizeMulti" );
201    CFMap M;
202    CanonicalForm UU, U = compress( arg, M );
203    CanonicalForm delta, omega, V = LC( U, 1 );
204    int t = U.level();
205    CFFList F = factorize( V );
206    CFFListIterator I, J;
207    CFArray G, lcG, D;
208    int i, j, k, m, r, maxdeg, h;
209    REvaluation A( 2, t, IntRandom( 50 ) );
210    CanonicalForm U0;
211    CanonicalForm ft, ut, gt, d;
212    modpk b;
213    bool negate = false;
214
215    DEBOUTLN( cerr, "-----------------------------------------------------" );
216    DEBOUTLN( cerr, "trying to factorize U = " << U );
217    DEBOUTLN( cerr, "U is a polynomial of level = " << arg.level() );
218    DEBOUTLN( cerr, "U will be factorized with respect to variable " << Variable(1) );
219    DEBOUTLN( cerr, "the leading coefficient of U with respect to that variable is " << V );
220    DEBOUTLN( cerr, "which is factorized as " << F );
221
222    maxdeg = 0;
223    for ( i = 2; i <= t; i++ )
224    {
225        j = U.degree( Variable( i ) );
226        if ( j > maxdeg ) maxdeg = j;
227    }
228
229    if ( F.getFirst().factor().inCoeffDomain() )
230    {
231        omega = F.getFirst().factor();
232        F.removeFirst();
233        if ( omega < 0 )
234        {
235            negate = true;
236            omega = -omega;
237            U = -U;
238        }
239    }
240    else
241        omega = 1;
242
243    bool goodeval = false;
244    r = 0;
245
246//    for ( i = 0; i < 10*t; i++ )
247//        A.nextpoint();
248
249    while ( ! goodeval )
250    {
251        TIMING_START(fac_findeval);
252        findEvaluation( U, V, omega, F, A, U0, delta, D, r );
253        TIMING_END(fac_findeval);
254        DEBOUTLN( cerr, "the evaluation point to reduce to an univariate problem is " << A );
255        DEBOUTLN( cerr, "corresponding delta = " << delta );
256        DEBOUTLN( cerr, "              omega = " << omega );
257        DEBOUTLN( cerr, "              D     = " << D );
258        DEBOUTLN( cerr, "now factorize the univariate polynomial " << U0 );
259        G = conv_to_factor_array( factorize( U0, false ) );
260        DEBOUTLN( cerr, "which factorizes into " << G );
261        #ifdef HAVE_NTL
262        {
263          int i=prime_number;
264          find_good_prime(arg,i);
265          find_good_prime(U0,i);
266          find_good_prime(U,i);
267          int p=cf_getSmallPrime(i);
268          //printf("found:p=%d (%d)\n",p,i);
269          if ((i==0)||(i!=prime_number))
270          {
271            b = coeffBound( U, p );
272            prime_number=i;
273          }
274          modpk bb=coeffBound(U0,p);
275          if (bb.getk() > b.getk() ) b=bb;
276          bb=coeffBound(arg,p);
277          if (bb.getk() > b.getk() ) b=bb;
278        }
279        #else
280        b = coeffBound( U, getZFacModulus().getp() );
281        if ( getZFacModulus().getpk() > b.getpk() )
282            b = getZFacModulus();
283        #endif
284        //printf("p=%d, k=%d\n",b.getp(),b.getk());
285        DEBOUTLN( cerr, "the coefficient bound of the factors of U is " << b.getpk() );
286
287        r = G.size();
288        lcG = CFArray( 1, r );
289        UU = U;
290        DEBOUTLN( cerr, "now trying to distribute the leading coefficients ..." );
291        TIMING_START(fac_distrib);
292        goodeval = distributeLeadingCoeffs( UU, G, lcG, F, D, delta, omega, A, r );
293        TIMING_END(fac_distrib);
294#ifdef DEBUGOUTPUT
295        if ( goodeval )
296        {
297            DEBOUTLN( cerr, "the univariate factors after distribution are " << G );
298            DEBOUTLN( cerr, "the distributed leading coeffs are " << lcG );
299            DEBOUTLN( cerr, "U may have changed and is now " << UU );
300            DEBOUTLN( cerr, "which has leading coefficient " << LC( UU, Variable(1) ) );
301
302            if ( LC( UU, Variable(1) ) != prod( lcG ) || A(UU) != prod( G ) )
303            {
304                DEBOUTLN( cerr, "!!! distribution was not correct !!!" );
305                DEBOUTLN( cerr, "product of leading coeffs is " << prod( lcG ) );
306                DEBOUTLN( cerr, "product of univariate factors is " << prod( G ) );
307                DEBOUTLN( cerr, "the new U is evaluated as " << A(UU) );
308            }
309            else
310                DEBOUTLN( cerr, "leading coeffs correct" );
311        }
312        else
313        {
314            DEBOUTLN( cerr, "we have found a bad evaluation point" );
315        }
316#endif
317        if ( goodeval )
318        {
319            TIMING_START(fac_hensel);
320            goodeval = Hensel( UU, G, lcG, A, b, Variable(1) );
321            TIMING_END(fac_hensel);
322        }
323    }
324    for ( i = 1; i <= r; i++ )
325    {
326        G[i] /= icontent( G[i] );
327        G[i] = M(G[i]);
328    }
329    // negate noch beachten !
330    if ( negate )
331        G[1] = -G[1];
332    DEBDECLEVEL( cerr, "ZFactorMulti" );
333    return G;
334}
335
336CFFList ZFactorizeMultivariate ( const CanonicalForm & f, bool issqrfree )
337{
338    CFFList G, F, R;
339    CFArray GG;
340    CFFListIterator i, j;
341    CFMap M;
342    CanonicalForm g, cont;
343    Variable v1, vm;
344    int k, m, n;
345
346    v1 = Variable(1);
347    if ( issqrfree )
348        F = CFFactor( f, 1 );
349    else
350        F = sqrFree( f );
351
352    for ( i = F; i.hasItem(); i++ )
353    {
354        if ( i.getItem().factor().inCoeffDomain() )
355        {
356            if ( ! i.getItem().factor().isOne() )
357                R.append( CFFactor( i.getItem().factor(), i.getItem().exp() ) );
358        }
359        else
360        {
361            TIMING_START(fac_content);
362            g = compress( i.getItem().factor(), M );
363            // now after compress g contains Variable(1)
364            vm = g.mvar();
365            g = swapvar( g, v1, vm );
366            cont = content( g );
367            g = swapvar( g / cont, v1, vm );
368            cont = swapvar( cont, v1, vm );
369            n = i.getItem().exp();
370            TIMING_END(fac_content);
371            DEBOUTLN( cerr, "now after content ..." );
372            if ( g.isUnivariate() )
373            {
374                G = factorize( g, true );
375                for ( j = G; j.hasItem(); j++ )
376                    if ( ! j.getItem().factor().isOne() )
377                        R.append( CFFactor( M( j.getItem().factor() ), n ) );
378            }
379            else
380            {
381                GG = ZFactorizeMulti( g );
382                m = GG.max();
383                for ( k = GG.min(); k <= m; k++ )
384                    if ( ! GG[k].isOne() )
385                        R.append( CFFactor( M( GG[k] ), n ) );
386            }
387            G = factorize( cont, true );
388            for ( j = G; j.hasItem(); j++ )
389                if ( ! j.getItem().factor().isOne() )
390                    R.append( CFFactor( M( j.getItem().factor() ), n ) );
391        }
392    }
393    return R;
394}
395
396static CFArray FpFactorizeMulti ( const CanonicalForm & arg )
397{
398    out_cf("FpFactorizeMulti:",arg,"\n");
399    DEBINCLEVEL( cerr, "FpFactorizeMulti" );
400    CFMap M;
401    CanonicalForm UU, U = compress( arg, M );
402    CanonicalForm delta, omega, V = LC( U, 1 );
403    int t = U.level();
404    CFFList F = factorize( V );
405    CFFListIterator I, J;
406    CFArray G, lcG, D;
407    int i, j, k, m, r, maxdeg, h;
408    REvaluation A( 2, t, FFRandom() );
409    CanonicalForm U0;
410    CanonicalForm ft, ut, gt, d;
411    modpk b;
412    bool negate = false;
413
414    DEBOUTLN( cerr, "-----------------------------------------------------" );
415    DEBOUTLN( cerr, "trying to factorize U = " << U );
416    DEBOUTLN( cerr, "U is a polynomial of level = " << arg.level() );
417    DEBOUTLN( cerr, "U will be factorized with respect to variable " << Variable(1) );
418    DEBOUTLN( cerr, "the leading coefficient of U with respect to that variable is " << V );
419    DEBOUTLN( cerr, "which is factorized as " << F );
420
421    maxdeg = 0;
422    out_cf("try:",U,"\n");
423    for ( i = 2; i <= t; i++ )
424    {
425        j = U.degree( Variable( i ) );
426        if ( j > maxdeg ) maxdeg = j;
427    }
428
429    if ( F.getFirst().factor().inCoeffDomain() )
430    {
431        omega = F.getFirst().factor();
432        F.removeFirst();
433        if ( omega < 0 )
434        {
435            negate = true;
436            omega = -omega;
437            U = -U;
438        }
439    }
440    else
441        omega = 1;
442
443    bool goodeval = false;
444    r = 0;
445
446//    for ( i = 0; i < 10*t; i++ )
447//        A.nextpoint();
448
449    while ( ! goodeval )
450    {
451        TIMING_START(fac_findeval);
452        findEvaluation( U, V, omega, F, A, U0, delta, D, r );
453        out_cf("univ.U0:",U0,"\n");
454        TIMING_END(fac_findeval);
455        DEBOUTLN( cerr, "the evaluation point to reduce to an univariate problem is " << A );
456        DEBOUTLN( cerr, "corresponding delta = " << delta );
457        DEBOUTLN( cerr, "              omega = " << omega );
458        DEBOUTLN( cerr, "              D     = " << D );
459        DEBOUTLN( cerr, "now factorize the univariate polynomial " << U0 );
460        G = conv_to_factor_array( factorize( U0, false ) );
461        printf("conv_to_factor_array\n");
462        DEBOUTLN( cerr, "which factorizes into " << G );
463
464        r = G.size();
465        lcG = CFArray( 1, r );
466        UU = U;
467        //if ( goodeval )
468        {
469            printf("start hensel\n");
470            TIMING_START(fac_hensel);
471            goodeval = Hensel( UU, G, lcG, A, b, Variable(1) );
472            TIMING_END(fac_hensel);
473        }
474    }
475    for ( i = 1; i <= r; i++ )
476    {
477        G[i] /= icontent( G[i] );
478        G[i] = M(G[i]);
479    }
480    // negate noch beachten !
481    if ( negate )
482        G[1] = -G[1];
483    DEBDECLEVEL( cerr, "ZFactorMulti" );
484    return G;
485}
486CFFList FpFactorizeMultivariate ( const CanonicalForm & f, bool issqrfree )
487{
488    CFFList G, F, R;
489    CFArray GG;
490    CFFListIterator i, j;
491    CFMap M;
492    CanonicalForm g, cont;
493    Variable v1, vm;
494    int k, m, n;
495
496    v1 = Variable(1);
497    if ( issqrfree )
498        F = CFFactor( f, 1 );
499    else
500        F = sqrFree( f );
501
502    printf("nach sqrFree:\n");
503    out_cff(F);
504    for ( i = F; i.hasItem(); i++ )
505    {
506        out_cf("consider:",i.getItem().factor(),"\n");
507        if ( i.getItem().factor().inCoeffDomain() )
508        {
509            if ( ! i.getItem().factor().isOne() )
510                R.append( CFFactor( i.getItem().factor(), i.getItem().exp() ) );
511        }
512        else
513        {
514            TIMING_START(fac_content);
515            g = compress( i.getItem().factor(), M );
516            // now after compress g contains Variable(1)
517            vm = g.mvar();
518            g = swapvar( g, v1, vm );
519            cont = content( g );
520            g = swapvar( g / cont, v1, vm );
521            cont = swapvar( cont, v1, vm );
522            n = i.getItem().exp();
523            TIMING_END(fac_content);
524            DEBOUTLN( cerr, "now after content ..." );
525            if ( g.isUnivariate() )
526            {
527                G = factorize( g, true );
528                for ( j = G; j.hasItem(); j++ )
529                    if ( ! j.getItem().factor().isOne() )
530                        R.append( CFFactor( M( j.getItem().factor() ), n ) );
531            }
532            else
533            {
534                GG = FpFactorizeMulti( g );
535                m = GG.max();
536                for ( k = GG.min(); k <= m; k++ )
537                    if ( ! GG[k].isOne() )
538                        R.append( CFFactor( M( GG[k] ), n ) );
539            }
540            out_cf("try cont:",cont,"\n");
541            G = factorize( cont, true );
542            out_cff(G);
543            for ( j = G; j.hasItem(); j++ )
544                if ( ! j.getItem().factor().isOne() )
545                    R.append( CFFactor( M( j.getItem().factor() ), n ) );
546        }
547    }
548    return R;
549}
Note: See TracBrowser for help on using the repository browser.