source: git/factory/fac_cantzass.cc @ 9899df

spielwiese
Last change on this file since 9899df was 9899df, checked in by Hans Schoenemann <hannes@…>, 13 years ago
mpz_mypow -> mpz_pow_ui, mpz_mypow_ui ->mpz_ui_pow_ui git-svn-id: file:///usr/local/Singular/svn/trunk@13649 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 8.4 KB
RevLine 
[493c477]1/* emacs edit mode for this file is -*- C++ -*- */
[341696]2/* $Id$ */
[2dd068]3
[b973c0]4#include <config.h>
5
[2dd068]6#include "cf_gmp.h"
7
8#include "assert.h"
[b973c0]9
[2dd068]10#include "cf_defs.h"
11#include "cf_random.h"
12#include "cf_util.h"
13#include "fac_cantzass.h"
14#include "fac_sqrfree.h"
15#include "gmpext.h"
16
17int initializeGMP();
18
19static int initialized = initializeGMP();
20
21static CanonicalForm randomPoly( int n, const Variable & x, const CFRandom & gen );
22
23static CFFList CantorZassenhausFactorFFGF( const CanonicalForm & f, int d, int q, const CFRandom & );
24
25static CFFList CantorZassenhausFactorExt( const CanonicalForm & g, int s, MP_INT * q, const CFRandom & gen );
26
27static CFFList distinctDegreeFactorFFGF ( const CanonicalForm & f, int q );
28
29static CFFList distinctDegreeFactorExt ( const CanonicalForm & f, int p, int n );
30
31// calculates f^m % d
32
33static CanonicalForm powerMod( const CanonicalForm & f, int m, const CanonicalForm & d );
34
35// calculates f^(p^s) % d
36
37static CanonicalForm powerMod( const CanonicalForm & f, int p, int s, const CanonicalForm & d );
38
39// calculates f^((p^s-1)/2) % d
40
41static CanonicalForm powerMod2( const CanonicalForm & f, int p, int s, const CanonicalForm & d );
42
43static CanonicalForm powerMod2( const CanonicalForm & f, MP_INT * q, int s, const CanonicalForm & d );
44
[9c9e2a4]45CFFList FpFactorizeUnivariateCZ( const CanonicalForm& f, bool issqrfree, int numext, const Variable alpha, const Variable beta )
[2dd068]46{
47    CFFList F, G, H, HH;
48    CanonicalForm fac;
49    ListIterator<CFFactor> i, j, k;
50    int d, q, n = 0;
51    bool galoisfield = getGFDegree() > 1;
52    MP_INT qq;
53
54    if ( galoisfield )
[cd86ac]55        q = ipower( getCharacteristic(), getGFDegree() );
[2dd068]56    else
[cd86ac]57        q = getCharacteristic();
[2dd068]58    if ( numext > 0 ) {
[cd86ac]59        if ( numext == 1 )
60            n = getMipo( alpha ).degree();
61        else
62            n = getMipo( alpha ).degree() * getMipo( beta ).degree();
63        mpz_init( &qq );
[9899df]64        mpz_ui_pow_ui ( &qq, q, n );
[2dd068]65    }
66    if ( LC( f ).isOne() )
[cd86ac]67        if ( issqrfree )
68            F.append( CFFactor( f, 1 ) );
69        else
70            F = sqrFreeFp( f );
[2dd068]71    else {
[cd86ac]72        if ( issqrfree )
73            F.append( CFFactor( f / LC( f ), 1 ) );
74        else
75            F = sqrFreeFp( f / LC( f ) );
76        H.append( LC( f ) );
[2dd068]77    }
78    for ( i = F; i.hasItem(); ++i ) {
[cd86ac]79        d = i.getItem().exp();
80        if ( numext > 0 )
81            G = distinctDegreeFactorExt( i.getItem().factor(), q, n );
82        else
83            G = distinctDegreeFactorFFGF( i.getItem().factor(), q );
84        for ( j = G; j.hasItem(); ++j ) {
85            if ( numext > 0 ) {
[0509dde]86             if ( numext == 1 ) {
87                   AlgExtRandomF tmpalpha( alpha );
88                    HH = CantorZassenhausFactorExt( j.getItem().factor(), j.getItem().exp(), &qq, tmpalpha );
89             }
90             else {
91                   AlgExtRandomF tmpalphabeta( alpha, beta );
92                    HH = CantorZassenhausFactorExt( j.getItem().factor(), j.getItem().exp(), &qq, tmpalphabeta );
93             }
[cd86ac]94            }
95            else  if ( galoisfield )
96                HH = CantorZassenhausFactorFFGF( j.getItem().factor(), j.getItem().exp(), q, GFRandom() );
97            else
98                HH = CantorZassenhausFactorFFGF( j.getItem().factor(), j.getItem().exp(), q, FFRandom() );
99            for ( k = HH; k.hasItem(); ++k ) {
100                fac = k.getItem().factor();
101                H.append( CFFactor( fac / LC( fac ), d ) );
102            }
103        }
[2dd068]104    }
105    if ( numext > 0 )
[cd86ac]106        mpz_clear( &qq );
[2dd068]107    return H;
108}
109
110CFFList distinctDegreeFactorFFGF ( const CanonicalForm & f, int q )
111{
112    int i;
113    Variable x = f.mvar();
114    CanonicalForm g = f, h, r = x;
115    CFFList F;
116    i = 1;
117    while ( g.degree(x) > 0 && i <= g.degree(x) ) {
[cd86ac]118        r = powerMod( r, q, g );
119        h = gcd( g, r - x );
120        if ( h.degree(x) > 0 ) {
121            F.append( CFFactor( h, i ) );
122            g /= h;
123        }
124        i++;
[2dd068]125    }
126    ASSERT( g.degree(x) == 0, "fatal fatal" );
127    return F;
128}
129
130CFFList distinctDegreeFactorExt ( const CanonicalForm & f, int p, int n )
131{
132    int i;
133    Variable x = f.mvar();
134    CanonicalForm g = f, h, r = x;
135    CFFList F;
136    i = 1;
137    while ( g.degree(x) > 0 && i <= g.degree(x) ) {
[cd86ac]138        r = powerMod( r, p, n, g );
139        h = gcd( g, r - x );
140        if ( h.degree(x) > 0 ) {
141            F.append( CFFactor( h, i ) );
142            g /= h;
143        }
144        i++;
[2dd068]145    }
146    ASSERT( g.degree(x) == 0, "fatal fatal" );
147    return F;
148}
149
150CFFList CantorZassenhausFactorFFGF( const CanonicalForm & g, int s, int q, const CFRandom & gen )
151{
152    CanonicalForm f = g;
153    CanonicalForm b, f1;
154    int d, d1;
155    Variable x = f.mvar();
[b973c0]156
[2dd068]157    if ( (d=f.degree(x)) == s )
[cd86ac]158        return CFFactor( f, 1 );
[2dd068]159    else while ( 1 ) {
[cd86ac]160        b = randomPoly( d, x, gen );
161        f1 = gcd( b, f );
162        if ( (d1 = f1.degree(x)) > 0 && d1 < d ) {
163            CFFList firstFactor = CantorZassenhausFactorFFGF( f1, s, q, gen );
164            CFFList secondFactor = CantorZassenhausFactorFFGF( f/f1, s, q, gen );
165            return Union( firstFactor, secondFactor );
166        } else {
167            f1 = gcd( f, powerMod2( b, q, s, f ) - 1 );
168            if ( (d1 = f1.degree(x)) > 0 && d1 < d ) {
169                CFFList firstFactor = CantorZassenhausFactorFFGF( f1, s, q, gen );
170                CFFList secondFactor = CantorZassenhausFactorFFGF( f/f1, s, q, gen );
171                return Union( firstFactor, secondFactor );
172            }
173        }
[2dd068]174    }
175}
176
177CFFList CantorZassenhausFactorExt( const CanonicalForm & g, int s, MP_INT * q, const CFRandom & gen )
178{
179    CanonicalForm f = g;
180    CanonicalForm b, f1;
181    int d, d1;
182    Variable x = f.mvar();
[b973c0]183
[2dd068]184    if ( (d=f.degree(x)) == s )
[cd86ac]185        return CFFactor( f, 1 );
[2dd068]186    else while ( 1 ) {
[cd86ac]187        b = randomPoly( d, x, gen );
188        f1 = gcd( b, f );
189        if ( (d1 = f1.degree(x)) > 0 && d1 < d ) {
190            CFFList firstFactor = CantorZassenhausFactorExt( f1, s, q, gen );
191            CFFList secondFactor = CantorZassenhausFactorExt( f/f1, s, q, gen );
192            return Union( firstFactor, secondFactor );
193        } else {
194            f1 = gcd( f, powerMod2( b, q, s, f ) - 1 );
195            if ( (d1 = f1.degree(x)) > 0 && d1 < d ) {
196                CFFList firstFactor = CantorZassenhausFactorExt( f1, s, q, gen );
197                CFFList secondFactor = CantorZassenhausFactorExt( f/f1, s, q, gen );
198                return Union( firstFactor, secondFactor );
199            }
200        }
[2dd068]201    }
202}
203
204CanonicalForm randomPoly( int d, const Variable & x, const CFRandom & g )
205{
206    CanonicalForm result = 0;
207    for ( int i = 0; i < d; i++ )
[cd86ac]208        result += power( x, i ) * g.generate();
[2dd068]209    result += power( x, d );
210    return result;
211}
212
213CanonicalForm powerMod( const CanonicalForm & f, int m, const CanonicalForm & d )
214{
215    CanonicalForm prod = 1;
216    CanonicalForm b = f % d;
217
218    while ( m != 0 ) {
[cd86ac]219        if ( m % 2 != 0 )
220            prod = (prod * b) % d;
221        m /= 2;
222        if ( m != 0 )
223            b = (b * b) % d;
[2dd068]224    }
225    return prod;
226}
227
228CanonicalForm powerMod( const CanonicalForm & f, int p, int s, const CanonicalForm & d )
229{
230    CanonicalForm prod = 1;
231    CanonicalForm b = f % d;
232    int odd;
233
234    MP_INT m;
[b973c0]235
[2dd068]236    mpz_init( &m );
[9899df]237    mpz_ui_pow_ui ( &m, p, s );
[cacfb6]238    while ( mpz_cmp_si( &m, 0 ) != 0 )
239    {
240        odd = mpz_fdiv_q_ui( &m, &m, 2 );
[cd86ac]241        if ( odd != 0 )
242            prod = (prod * b) % d;
243        if ( mpz_cmp_si( &m, 0 ) != 0 )
244            b = (b*b) % d;
[2dd068]245    }
246    mpz_clear( &m );
247    return prod;
248}
249
250CanonicalForm powerMod2( const CanonicalForm & f, int p, int s, const CanonicalForm & d )
251{
252    CanonicalForm prod = 1;
253    CanonicalForm b = f % d;
254    int odd;
255
256    MP_INT m;
[b973c0]257
[2dd068]258    mpz_init( &m );
[9899df]259    mpz_ui_pow_ui ( &m, p, s );
[2dd068]260    mpz_sub_ui( &m, &m, 1 );
[cacfb6]261    mpz_fdiv_q_ui( &m, &m, 2 );
262    while ( mpz_cmp_si( &m, 0 ) != 0 )
263    {
264        odd = mpz_fdiv_q_ui( &m, &m, 2 );
[cd86ac]265        if ( odd != 0 )
266            prod = (prod * b) % d;
267        if ( mpz_cmp_si( &m, 0 ) != 0 )
268            b = (b*b) % d;
[2dd068]269    }
270    mpz_clear( &m );
271    return prod;
272}
273
274CanonicalForm powerMod2( const CanonicalForm & f, MP_INT * q, int s, const CanonicalForm & d )
275{
276    CanonicalForm prod = 1;
277    CanonicalForm b = f % d;
278    int odd;
279
280    MP_INT m;
[b973c0]281
[2dd068]282    mpz_init( &m );
[9899df]283    mpz_pow_ui( &m, q, s );
[2dd068]284    mpz_sub_ui( &m, &m, 1 );
[cacfb6]285    mpz_fdiv_q_ui( &m, &m, 2 );
286    while ( mpz_cmp_si( &m, 0 ) != 0 )
287    {
288        odd = mpz_fdiv_q_ui( &m, &m, 2 );
[cd86ac]289        if ( odd != 0 )
290            prod = (prod * b) % d;
291        if ( mpz_cmp_si( &m, 0 ) != 0 )
292            b = (b*b) % d;
[2dd068]293    }
294    mpz_clear( &m );
295    return prod;
296}
Note: See TracBrowser for help on using the repository browser.