source: git/factory/cf_resultant.cc @ 75744d

spielwiese
Last change on this file since 75744d was 75744d, checked in by Jens Schmidt <schmidt@…>, 26 years ago
* cf_resultant.cc (subResChain): new function git-svn-id: file:///usr/local/Singular/svn/trunk@648 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 3.2 KB
Line 
1/* emacs edit mode for this file is -*- C++ -*- */
2/* $Id: cf_resultant.cc,v 1.1 1997-09-01 09:01:15 schmidt Exp $ */
3
4//{{{ docu
5//
6// cf_resultant.cc - algorithms for calculating resultants.
7//
8// Header file: cf_algorithm.h
9//
10//}}}
11
12#include <config.h>
13
14#include "assert.h"
15
16#include "canonicalform.h"
17#include "variable.h"
18
19//{{{ CFArray subResChain ( const CanonicalForm & f, const CanonicalForm & g, Variable x )
20//{{{ docu
21//
22// subResChain() - caculate extended subresultant chain.
23//
24// The chain is calculated from f and g with respect to variable
25// x which should not be an algebraic variable.  If f or q equals
26// zero, an array consisting of one zero entry is returned.
27//
28// Note: this is not the standard subresultant chain but the
29// *extended* chain!
30//
31// This algorithm is from the article of R. Loos - 'Generalized
32// Polynomial Remainder Sequences' in B. Buchberger - 'Computer
33// Algebra - Symbolic and Algebraic Computation' with some
34// necessary extensions concerning the calculation of the first
35// step.
36//
37//}}}
38CFArray
39subResChain ( const CanonicalForm & f, const CanonicalForm & g, Variable x )
40{
41    ASSERT( x.level() > 0, "cannot calculate subresultant sequence with respect to algebraic variables" );
42
43    CFArray trivialResult( 0, 0 );
44    CanonicalForm F, G;
45    Variable X;
46
47    // some checks on triviality
48    if ( f.isZero() || g.isZero() ) {
49        trivialResult[0] = 0;
50        return trivialResult;
51    }
52
53    // make x main variable
54    if ( f.mvar() > x || g.mvar() > x ) {
55        if ( f.mvar() > g.mvar() )
56            X = f.mvar();
57        else
58            X = g.mvar();
59        F = swapvar( f, X, x );
60        G = swapvar( g, X, x );
61    }
62    else {
63        X = x;
64        F = f;
65        G = g;
66    }
67    // at this point, we have to calculate the sequence of F and
68    // G in respect to X where X is equal to or greater than the
69    // main variables of F and G
70
71    // initialization of chain
72    int m = degree( F, X );
73    int n = degree( G, X );
74
75    int j = (m <= n) ? n : m-1;
76    int r;
77
78    CFArray S( 0, j+1 );
79    CanonicalForm R;
80    S[j+1] = F; S[j] = G;
81
82    // make sure that S[j+1] is regular and j < n
83    if ( m == n && j > 0 ) {
84        S[j-1] = LC( S[j], X ) * psr( S[j+1], S[j], X );
85        j--;
86    } else if ( m < n ) {
87        S[j-1] = LC( S[j], X ) * LC( S[j], X ) * S[j+1];
88        j--;
89    } else if ( m > n && j > 0 ) {
90        // calculate first step
91        r = degree( S[j], X );
92        R = LC( S[j+1], X );
93
94        // if there was a gap calculate similar polynomial
95        if ( j > r && r >= 0 )
96            S[r] = power( LC( S[j], X ), j - r ) * S[j] * power( R, j - r );
97
98        if ( r > 0 ) {
99            // calculate remainder
100            S[r-1] = psr( S[j+1], S[j], X ) * power( -R, j - r );
101            j = r-1;
102        }
103    }
104
105    while ( j > 0 ) {
106        // at this point, 0 < j < n and S[j+1] is regular
107        r = degree( S[j], X );
108        R = LC( S[j+1], X );
109       
110        // if there was a gap calculate similar polynomial
111        if ( j > r && r >= 0 )
112            S[r] = (power( LC( S[j], X ), j - r ) * S[j]) / power( R, j - r );
113
114        if ( r <= 0 ) break;
115        // calculate remainder
116        S[r-1] = psr( S[j+1], S[j], X ) / power( -R, j - r + 2 );
117
118        j = r-1;
119        // again 0 <= j < r <= jOld and S[j+1] is regular
120    }
121
122    for ( j = 0; j <= S.max(); j++ ) {
123        // reswap variables if necessary
124        if ( X != x ) {
125            S[j] = swapvar( S[j], X, x );
126        }
127    }
128
129    return S;
130}
131//}}}
Note: See TracBrowser for help on using the repository browser.