source: git/factory/fac_util.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: 6.2 KB
Line 
1/* emacs edit mode for this file is -*- C++ -*- */
2
3#include "config.h"
4
5#include "cf_assert.h"
6
7#include "cf_defs.h"
8#include "canonicalform.h"
9#include "cf_iter.h"
10#include "fac_util.h"
11
12static CanonicalForm PK, PKHALF;
13
14static CanonicalForm mappk ( const CanonicalForm& );
15
16static CanonicalForm mappksymmetric ( const CanonicalForm& );
17
18
19modpk::modpk()
20{
21    p = 0;
22    k = 0;
23    pk = 1;
24    pkhalf = 0;
25}
26
27modpk::modpk( int q, int l )
28{
29    p = q;
30    k = l;
31    pk = power( CanonicalForm( p ), k );
32    pkhalf = pk / 2;
33}
34
35modpk::modpk( const modpk & m )
36{
37    p = m.p;
38    k = m.k;
39    pk = m.pk;
40    pkhalf = m.pkhalf;
41}
42
43modpk&
44modpk::operator= ( const modpk & m )
45{
46    if ( this != &m ) {
47        p = m.p;
48        k = m.k;
49        pk = m.pk;
50        pkhalf = m.pkhalf;
51    }
52    return *this;
53}
54
55CanonicalForm
56modpk::inverse( const CanonicalForm & f, bool symmetric ) const
57{
58    CanonicalForm u, r0 = this->operator()( f, false ), r1 = pk, q0 = 1, q1 = 0;
59    while ( ( r0 > 0 ) && ( r1 > 0 ) ) {
60        u = r0 / r1;
61        r0 = r0 % r1;
62        q0 = u*q1 + q0;
63        if ( r0 > 0 ) {
64            u = r1 / r0;
65            r1 = r1 % r0;
66            q1 = u*q0 + q1;
67        }
68    }
69    if ( r0 == 0 )
70        return this->operator()( pk-q1, symmetric );
71    else
72        return this->operator()( q0, symmetric );
73}
74
75CanonicalForm
76modpk::operator() ( const CanonicalForm & f, bool symmetric ) const
77{
78    PKHALF = pkhalf;
79    PK = pk;
80    if ( symmetric )
81        return mapdomain( f, mappksymmetric );
82    else
83        return mapdomain( f, mappk );
84}
85
86CanonicalForm
87replaceLc( const CanonicalForm & f, const CanonicalForm & c )
88{
89    if ( f.inCoeffDomain() )
90        return c;
91    else
92        return f + ( c - LC( f ) ) * power( f.mvar(), degree( f ) );
93}
94
95CanonicalForm
96remainder( const CanonicalForm & f, const CanonicalForm & g, const modpk & pk )
97{
98    ASSERT( (f.inCoeffDomain() || f.isUnivariate()) && (g.inCoeffDomain() || g.isUnivariate()) && (f.inCoeffDomain() || g.inCoeffDomain() || f.mvar() == g.mvar()), "can not build remainder" );
99    if ( f.inCoeffDomain() )
100        if ( g.inCoeffDomain() )
101            return pk( f % g );
102        else
103            return pk( f );
104    else {
105        Variable x = f.mvar();
106        CanonicalForm result = f;
107        int degg = g.degree();
108        CanonicalForm invlcg = pk.inverse( g.lc() );
109        CanonicalForm gg = pk( g*invlcg );
110        if((gg.lc().isOne()))
111        {
112          while ( result.degree() >= degg )
113          {
114            result -= pk(lc( result ) * gg) * power( x, result.degree() - degg );
115            result=pk(result);
116          }
117        }
118        else
119        // no inverse found
120        {
121          CanonicalForm ic=icontent(g);
122          if (!ic.isOne())
123          {
124            gg=g/ic;
125            return remainder(f,gg,pk);
126          }
127          while ( result.degree() >= degg )
128          {
129            if (gg.lc().isZero()) return result;
130            CanonicalForm lcgf = result.lc() / gg.lc();
131            if (lcgf.inZ())
132              gg = pk( g*lcgf );
133            else
134            {
135              //printf("!\n\n");
136              return result;
137            }
138            result -=  gg * power( x, result.degree() - degg );
139            result=pk(result);
140          }
141        }
142        return result;
143    }
144}
145
146void
147divremainder( const CanonicalForm & f, const CanonicalForm & g, CanonicalForm & quot, CanonicalForm & rem, const modpk & pk )
148{
149    ASSERT( (f.inCoeffDomain() || f.isUnivariate()) && (g.inCoeffDomain() || g.isUnivariate()) && (f.inCoeffDomain() || g.inCoeffDomain() || f.mvar() == g.mvar()), "can not build remainder" );
150    if ( f.inCoeffDomain() )
151        if ( g.inCoeffDomain() ) {
152            divrem( f, g, quot, rem );
153            quot = pk( quot );
154            rem = pk( rem );
155        }
156        else {
157            quot = 0;
158            rem = pk( f );
159        }
160    else {
161        Variable x = f.mvar();
162        CanonicalForm invlcg = pk.inverse( g.lc() );
163        rem = f;
164        int degg = g.degree();
165        while ( rem.degree() >= degg ) {
166            quot += pk( lc( rem ) * invlcg ) * power( x, rem.degree() - degg );
167            rem = pk( rem - lc( rem ) * invlcg * g * power( x, rem.degree() - degg ) );
168        }
169    }
170}
171
172CanonicalForm
173mappksymmetric ( const CanonicalForm & f )
174{
175    CanonicalForm result = mod( f, PK );
176    if ( result > PKHALF )
177        return result - PK;
178    else
179        return result;
180}
181
182CanonicalForm
183mappk ( const CanonicalForm & f )
184{
185    return mod( f, PK );
186}
187
188void
189extgcd ( const CanonicalForm & a, const CanonicalForm & b, CanonicalForm & S, CanonicalForm & T, const modpk & pk )
190{
191    int p = pk.getp(), k = pk.getk(), j;
192    CanonicalForm amodp, bmodp, smodp, tmodp, s, t, sigma, tau, e;
193    CanonicalForm modulus = p, sigmat, taut, q;
194
195    setCharacteristic( p );
196    {
197        amodp = mapinto( a ); bmodp = mapinto( b );
198        (void)extgcd( amodp, bmodp, smodp, tmodp );
199    }
200    setCharacteristic( 0 );
201    s = mapinto( smodp ); t = mapinto( tmodp );
202
203    for ( j = 1; j < k; j++ ) {
204        e = ( 1 - s * a - t * b ) / modulus;
205        setCharacteristic( p );
206        {
207            e = mapinto( e );
208            sigmat = smodp * e;
209            taut = tmodp * e;
210            divrem( sigmat, bmodp, q, sigma );
211            tau = taut + q * amodp;
212        }
213        setCharacteristic( 0 );
214        s += mapinto( sigma ) * modulus;
215        t += mapinto( tau ) * modulus;
216        modulus *= p;
217    }
218    S = s; T = t;
219}
220
221CanonicalForm
222sum ( const CFArray & a, int f, int l )
223{
224    if ( f < a.min() ) f = a.min();
225    if ( l > a.max() ) l = a.max();
226    CanonicalForm s = 0;
227    for ( int i = f; i <= l; i++ )
228        s += a[i];
229    return s;
230}
231
232CanonicalForm
233prod ( const CFArray & a, int f, int l )
234{
235    if ( f < a.min() ) f = a.min();
236    if ( l > a.max() ) l = a.max();
237    CanonicalForm p = 1;
238    for ( int i = f; i <= l; i++ )
239        p *= a[i];
240    return p;
241}
242
243CanonicalForm
244sum ( const CFArray & a )
245{
246    return sum( a, a.min(), a.max() );
247}
248
249CanonicalForm
250prod ( const CFArray & a )
251{
252    return prod( a, a.min(), a.max() );
253}
254
255CanonicalForm
256crossprod ( const CFArray & a, const CFArray & b )
257{
258    ASSERT( a.size() == b.size(), "array size mismatch" );
259    CanonicalForm s = 0;
260    int fa = a.min();
261    int fb = b.min();
262    int n = a.max();
263    for ( ; fa <= n; fa++, fb++ )
264        s += a[fa] * b[fb];
265    return s;
266}
Note: See TracBrowser for help on using the repository browser.