source: git/factory/cf_ops.cc @ e4eea0e

spielwiese
Last change on this file since e4eea0e was e4eea0e, checked in by Jens Schmidt <schmidt@…>, 27 years ago
* cf_ops.cc (resultant): assertion added (resultant): new variable flipFactor (resultant): handling of trivial cases fixed (resultant): some other minor bug fixes * cf_ops.cc (resultant): doc fix git-svn-id: file:///usr/local/Singular/svn/trunk@537 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 8.1 KB
Line 
1/* emacs edit mode for this file is -*- C++ -*- */
2/* $Id: cf_ops.cc,v 1.3 1997-07-19 08:09:28 schmidt Exp $ */
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
12CanonicalForm
13psr( const CanonicalForm &f, const CanonicalForm &g, const Variable & x )
14{
15    return ( power( LC(g,x), degree(f,x)-degree(g,x)+1 ) * f ) % g;
16}
17
18CanonicalForm
19psq( const CanonicalForm &f, const CanonicalForm &g, const Variable & x )
20{
21    return ( power( LC(g,x), degree(f,x)-degree(g,x)+1 ) * f ) / g;
22}
23
24void
25psqr( const CanonicalForm &f, const CanonicalForm &g, CanonicalForm &q, CanonicalForm &r, const Variable& x )
26{
27    divrem( power( LC(g,x), degree(f,x)-degree(g,x)+1 ) * f, g, q, r );
28}
29
30static void swapvar_rec ( const CanonicalForm &f, CanonicalForm &result, const CanonicalForm &term );
31
32static void swapvar_between ( const CanonicalForm &f, CanonicalForm &result, const CanonicalForm &term, int expx2 );
33
34static Variable sv_x1, sv_x2;
35
36CanonicalForm
37swapvar ( const CanonicalForm &f, const Variable &x1, const Variable &x2 )
38{
39    ASSERT( x1.level() > 0 && x2.level() > 0, "cannot swap algebraic Variables" );
40    if ( f.inCoeffDomain() || x1 == x2 || ( x1 > f.mvar() && x2 > f.mvar() ) )
41        return f;
42    else {
43        CanonicalForm result = 0;
44        if ( x1 > x2 ) {
45            sv_x1 = x2; sv_x2 = x1;
46        }
47        else {
48            sv_x1 = x1; sv_x2 = x2;
49        }
50        if ( f.mvar() < sv_x2 )
51            swapvar_between( f, result, 1, 0 );
52        else
53            swapvar_rec( f, result, 1 );
54        return result;
55    }
56}
57
58void
59swapvar_rec ( const CanonicalForm &f, CanonicalForm &result, const CanonicalForm &term )
60{
61    if ( f.inCoeffDomain() || f.mvar() < sv_x1 )
62        result += term * f;
63    else  if ( f.mvar() == sv_x2 )
64        for ( CFIterator i = f; i.hasTerms(); i++ )
65            swapvar_between( i.coeff(), result, term, i.exp() );
66    else  if ( f.mvar() < sv_x2 )
67        swapvar_between( f, result, term, 0 );
68    else
69        for ( CFIterator i = f; i.hasTerms(); i++ )
70            swapvar_rec( i.coeff(), result, term*power( f.mvar(), i.exp() ) );
71}
72
73void
74swapvar_between ( const CanonicalForm &f, CanonicalForm &result, const CanonicalForm &term, int expx2 )
75{
76    if ( f.inCoeffDomain() || f.mvar() < sv_x1 )
77        result += term * power( sv_x1, expx2 ) * f;
78    else  if ( f.mvar() == sv_x1 )
79        for ( CFIterator i = f; i.hasTerms(); i++ )
80            result += power( sv_x2, i.exp() ) * term * i.coeff() * power( sv_x1, expx2 );
81    else
82        for ( CFIterator i = f; i.hasTerms(); i++ )
83            swapvar_between( i.coeff(), result, term*power( f.mvar(), i.exp() ), expx2 );
84}
85
86CanonicalForm
87apply ( const CanonicalForm & f, void (*mf)( CanonicalForm &, int & ) )
88{
89    if ( f.inCoeffDomain() ) {
90        int exp = 0;
91        CanonicalForm result;
92        mf( result, exp );
93        ASSERT( exp == 0, "illegal result, do not know what variable to use" );
94        return result;
95    }
96    else {
97        CanonicalForm result, coeff;
98        CFIterator i;
99        int exp;
100        Variable x = f.mvar();
101        for ( i = f; i.hasTerms(); i++ ) {
102            coeff = i.coeff();
103            exp = i.exp();
104            mf( coeff, exp );
105            if ( ! coeff.isZero() )
106                result += power( x, exp ) * coeff;
107        }
108        return result;
109    }
110}
111
112void
113degreesRec ( const CanonicalForm & f, int * degs )
114{
115    if ( ! f.inCoeffDomain() ) {
116        int level = f.level();
117        int deg = f.degree();
118        if ( degs[level] < deg )
119            degs[level] = f.degree();
120        for ( CFIterator i = f; i.hasTerms(); i++ )
121            degreesRec( i.coeff(), degs );
122    }
123}
124
125int *
126degrees ( const CanonicalForm & f, int * degs )
127{
128    if ( f.inCoeffDomain() )
129        return 0;
130    else {
131        int level = f.level();
132        if ( degs == 0 )
133            degs = new int[level+1];
134        for ( int i = 0; i <= level; i++ )
135            degs[i] = 0;
136        degreesRec( f, degs );
137        return degs;
138    }
139}
140
141CanonicalForm
142mapdomain ( const CanonicalForm & f, CanonicalForm (*mf)( const CanonicalForm & ) )
143{
144    if ( f.inCoeffDomain() )
145        return mf( f );
146    else {
147        CanonicalForm result = 0;
148        CFIterator i;
149        Variable x = f.mvar();
150        for ( i = f; i.hasTerms(); i++ )
151            result += power( x, i.exp() ) * mapdomain( i.coeff(), mf );
152        return result;
153    }
154}
155
156int
157totaldegree ( const CanonicalForm & f )
158{
159    if ( f.isZero() )
160        return -1;
161    else if ( f.inCoeffDomain() )
162        return 0;
163    else {
164        CFIterator i;
165        int cdeg = 0, dummy;
166        for ( i = f; i.hasTerms(); i++ )
167            if ( (dummy = totaldegree( i.coeff() ) + i.exp()) > cdeg )
168                cdeg = dummy;
169        return cdeg;
170    }
171}
172
173
174int
175totaldegree ( const CanonicalForm & f, const Variable & v1, const Variable & v2 )
176{
177    if ( f.isZero() )
178        return -1;
179    else if ( v1 > v2 )
180        return 0;
181    else if ( f.inCoeffDomain() )
182        return 0;
183    else if ( f.mvar() < v1 )
184        return 0;
185    else if ( f.mvar() == v1 )
186        return f.degree();
187    else if ( f.mvar() > v2 ) {
188        CFIterator i;
189        int cdeg = 0, dummy;
190        for ( i = f; i.hasTerms(); i++ )
191            if ( (dummy = totaldegree( i.coeff(), v1, v2 )) > cdeg )
192                cdeg = dummy;
193        return cdeg;
194    }
195    else {
196        CFIterator i;
197        int cdeg = 0, dummy;
198        for ( i = f; i.hasTerms(); i++ )
199            if ( (dummy = totaldegree( i.coeff(), v1, v2 ) + i.exp()) > cdeg )
200                cdeg = dummy;
201        return cdeg;
202    }
203}
204
205//{{{ CanonicalForm resultant( const CanonicalForm & f, const CanonicalForm & g, const Variable & x )
206//{{{ docu
207//
208// resultant() - return resultant of f and g with respect to x.
209//
210// We calculate the resultant using a subresultant PSR.
211//
212// flipFactor: Res(f, g) = flipFactor * Res(g, f)
213// F, G: f and g with x as main variable
214// pi, pi1, pi2: used to compute PSR
215// delta:
216// bi, Hi:
217//
218//}}}
219CanonicalForm
220resultant( const CanonicalForm & f, const CanonicalForm & g, const Variable & x )
221{
222    CanonicalForm Hi, bi, pi, pi1, pi2, F, G;
223    int delta, flipFactor;
224    Variable v;
225
226    ASSERT( x.level() > 0, "cannot calculate resultant in respect to algebraic variables" );
227
228    // some checks on triviality.  We will not use degree( v )
229    // here because this may involve variable swapping.
230    if ( f.isZero() || g.isZero() )
231        return 0;
232    if ( f.mvar() < x )
233        return power( f, g.degree( x ) );
234    if ( g.mvar() < x )
235        return power( g, f.degree( x ) );
236
237    // make x main variale
238    if ( f.mvar() > x || g.mvar() > x ) {
239        if ( f.mvar() > g.mvar() )
240            v = f.mvar();
241        else
242            v = g.mvar();
243        F = swapvar( f, v, x );
244        G = swapvar( g, v, x );
245    }
246    else {
247        v = x;
248        F = f;
249        G = g;
250    }
251    // at this point, we have to calculate resultant( F, G, v )
252    // where v is equal to or greater than the main variables
253    // of F and G
254
255    // trivial case: F or G in R.  Swapping will not occur
256    // when calling degree( v ).
257    if ( F.degree( v ) < 1 )
258        return power( f, G.degree( v ) );
259    if ( G.degree( v ) < 1 )
260        return power( g, F.degree( v ) );
261
262    // start the pseudo remainder sequence
263    if ( F.degree( v ) >= G.degree( v ) ) {
264        pi = F; pi1 = G;
265        flipFactor = 1;
266    }
267    else {
268        if ( (F.degree( v ) * G.degree( v )) % 2 )
269            flipFactor = -1;
270        else
271            flipFactor = 1;
272        pi = G; pi1 = F;
273    }
274
275    delta = degree( pi, v ) - degree( pi1, v );
276    Hi = power( LC( pi1, v ), delta );
277
278    // Ist hier nicht if und else zweig vertauscht ???
279    if ( (delta+1) % 2 )
280        bi = 1;
281    else
282        bi = -1;
283
284    // Ist pi1.isZero vielleich schneller ???
285    while ( degree( pi1, v ) >= 0 ) {
286        pi2 = psr( pi, pi1, v );
287        pi2 = pi2 / bi;
288        pi = pi1; pi1 = pi2;
289        if ( degree( pi1, v ) >= 0 ) {
290            delta = degree( pi, v ) - degree( pi1, v );
291
292            // Ist hier nicht if und else zweig vertauscht ???
293            if ( (delta+1) % 2 )
294                bi = LC( pi, v ) * power( Hi, delta );
295            else
296                bi = -LC( pi, v ) * power( Hi, delta );
297
298            // Was ist f"ur delta == 0 ???
299            Hi = power( LC( pi1, v ), delta ) / power( Hi, delta-1 );
300        }
301    }
302
303    // f and g have non-trivial common divisor
304    // if ( degree( pi, v ) > 0 )
305    // return 0;
306
307    // undo variable swap
308    if ( v == x )
309        // Gibt man hier nicht den letzten Rest der PSR zur"uck
310        // und nicht den Korrekturterm Hi ???
311        return Hi * flipFactor;
312    else
313        return swapvar( Hi, v, x ) * flipFactor;
314}
315//}}}
316
317static CanonicalForm
318cden ( const CanonicalForm & f )
319{
320    if ( f.inCoeffDomain() )
321        return f.den();
322    else {
323        CFIterator i;
324        CanonicalForm cd = 1;
325        for ( i = f; i.hasTerms(); i++ )
326            cd = lcm( cd, cden( i.coeff() ) );
327        return cd;
328    }
329}
330
331CanonicalForm
332common_den ( const CanonicalForm & f )
333{
334    if ( getCharacteristic() == 0 && isOn( SW_RATIONAL ) ) {
335        Off( SW_RATIONAL );
336        CanonicalForm cd = cden( f );
337        On( SW_RATIONAL );
338        return cd;
339    }
340    else
341        return 1;
342}
Note: See TracBrowser for help on using the repository browser.