source: git/factory/fac_util.cc @ ec970e

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