<< Ïðåäûäóùàÿ ñòð. 30(èç 33 ñòð.)ÎÃËÀÂËÅÍÈÅ Ñëåäóþùàÿ >>
s := {};
for u in fn do
s := s union cf2( u );
od;
s;
end:

coefficient_field := proc( f )
#
# Takes as input either a rational function with
# algebraic-number coefficients -- i.e., a function
# built up from rational numbers, RootOfs and symbols
# or a list or set of such functions.
#
# Returns an ordered list [ s[1], s[2], ..., s[n] ]
# of the RootOfs that appear in f.
# s[j] fails to appear in s[i] if i < j.

131
#
local fcollection, g, gnormal, s;
if type( f, list ) or type( f, set ) then
fcollection := f;
else
fcollection := { f };
fi;
s := {};
for g in fcollection do
gnormal := evala( Normal ( g ) );
s := s union cf2( gnormal );
od;
s := [op(s)];
s := sort( s, proc(a,b) return not evalb( has( a, b ) ) end );
s;
end:

parfrac_summand_decomp := proc( f, x, C_0 )
#
# Takes as input f, a summand of a partial fraction decomposition
# in the indeterminate x over the field C_0. It is assumed that
# f has nontrivial denominator. That is, it is assumed that we can
# write f = N/Q^d, where:
#
# * N and Q are polynomials over C_0
# * degree(N) < degree(Q)
# * Q is monic and irreducible over C_0
# * d is an integer, d >= 1
#
# Returns as output the list [N, Q, d].
#
local N1, D1, L, u, Q1, d, c;
N1 := numer( f );
D1 := denom( f );
L := evala( Factors( D1, C_0 ) );
u := L[1];
if nops( L[2] ) > 1 then
error "Too many factors in denominator of \
partial fractions summand %1" , f ;
fi;
Q1 := L[2][1][1];
d := L[2][1][2];
# f = N1 / (u*Q1^d)
c := lcoeff( Q1, x );
# f = N1 / (u*(c*Q2)^d) = N1 / (u * c^d * Q2^d)
# = (N1 / u / c^d) / Q2^d
return [ N1/u/c^d, Q1/c, d ] ;
end:

132
rational_number_summand_in := proc( a_input )
#
# Takes as input an algebraic number a_input. Suppose
# a_input = a1 + a2, where a1 is a rational number and
# a2 is a Q-linear combination of nontrivial power products
# of RootOfs. Then this procedure returns a1.
#
local a_normal, u;
a_normal := evala( Normal( a_input ) );
if not type( a_normal, algnum ) then
error "%1 is not an algebraic number", a_input;
fi;
if not type( a_normal, â€˜+â€˜ ) then
#
# There is only one summand in a_normal. Replace a_normal
# with { a_normal } so that the subsequent loop through the
# operands of a_normal treats a_normal as a summand.
#
a_normal := { a_normal };
fi;
#
# a_normal is an expanded element of the polynomial algebra
# determined by the RootOfs that appear in a_input. The
# following loop searches the operands that appear in a_normal.
# If one of the summands is a rational number, then the procedure
# returns this number (which is unique); otherwise the procedure
# returns zero.
#
for u in a_normal do
if type( u, rational ) then
return u ;
fi;
od;
return 0 ;
end:

ratfun_relation := proc( s1, s2, x )
#
# Takes the rational functions s1 and s2 with algebraic number
# coefficients as input.
#
# It is assumed that the distinct RootOfs appearing in s1 and s2
# are independent.
#

133
# Returns as output a list [j1,j2] of integers such that
#
# j1*s1 + j2*s2 = hâ€™/h
#
# for some rational function h.
#
# The list will be nonzero if possible; i.e., if there exist
# integers n1, n2, not both zero, such that n1*s1 + n2*s2 = gâ€™/g
# for some rational function g, then the output of
# ratfun_relation is UNEQUAL to [0,0].
#
# If j1,j2 are not both zero, then [j1,j2] will be
# in lowest possible terms. That is, if [j1,j2] = [c*k1,c*k2]
# for some integers c, k1, k2 such that k1*s1 + k2*s2 = fâ€™/f
# for some rational function f, then c is equal to
# either 1 or -1.
#
local C_0, sp, i, u, s_fld, s_cancel, psd, A, Q, Qprime,
alpha, B, alpha1, alpha2, scq, fld_denoms, n;
if s1 = 0 then return [1,0] ; fi;
if s2 = 0 then return [0,1] ; fi;
C_0 := coefficient_field( [ s1, s2 ] );
sp := [0,0];
s_fld := [0,0];
s_cancel := [0,0];
fld_denoms := [{},{}];
#
#
# sp will be a two-element list consisting of the
# partial-fraction decompositions of s1 and s2 over C_0.
# s_fld (for Fraction of a Logarithmic Derivative) will be
# a two-element list such that s_fld[i] is the "part" of
# sp[i] which is the log-derivative of an nth root of a
# rational function over C_0 for some n. s_cancel will be a
# two-element list such that
#
# sp[i] = s_fld[i] + s_cancel[i];
#
# we will want s_cancel[1] and s_cancel[2] to "cancel" in
# the sense that some Z-linear combination of these two
# rational functions should be zero; i.e., that they differ
# multiplicatively by a rational number.
#
#
sp[1] := convert( s1, parfrac, x, C_0 );
sp[2] := convert( s2, parfrac, x, C_0 );
for i from 1 to 2 do
if not type( sp[i], â€˜+â€˜ ) then
#

134
# There is only one summand in the partial-fraction
# decomposition of s_i. Replace s_i with { s_i } so that
# the subsequent loop through the operands of s_i treats
# s_i as a summand.
#
sp[i] := { sp[i] };
fi;
for u in sp[i] do
if type( u, polynom ) then
#
# Summand is of the form c*x^d with d >= 0;
# place in s_cancel[i]
#
s_cancel[i] := s_cancel[i] + u;
next;
fi;
psd := parfrac_summand_decomp( u, x, C_0 );
if psd[3] > 1 then
#
# Summand is of the form A/Q^d with d > 1;
# place in s_cancel[i]
#
s_cancel[i] := s_cancel[i] + u;
else
#
# Summand is of the form A/Q, Q irreducible over C_0;
# write A = alpha*diff(Q,x) + B, where deg B < deg Qâ€™.
# Then write alpha = alpha1 + alpha2, where
# alpha1 is a rational number and alpha2 is a
# Q-linear combination of nontrivial power products of
# RootOfs. Then,
#
# A/Q = ((alpha1+alpha2)*Qâ€™ + B)/Q
# = (alpha1 * Qâ€™ / Q) + (alpha2*Qâ€™+B)/Q
# = f1 + f2.
#
# Add f1 to s_fld[i], f2 to s_cancel[i].
#
#
#
A := psd[1];
Q := psd[2];
Qprime := diff( Q, x );
alpha := quo( A, Qprime, x, â€™Bâ€™ );
if degree( alpha, x ) > 0 then
error "Unexpected high-degree numerator in \
partial fractions summand %1", A/Q;
fi;
alpha1 := rational_number_summand_in( alpha );

135
fld_denoms[i] := fld_denoms[i] union { denom( alpha1 ) };
alpha2 := alpha - alpha1;
s_fld[i] := s_fld[i] + alpha1 * Qprime / Q;
s_cancel[i] := s_cancel[i] + (alpha2 * Qprime + B)/Q;
fi;
od;
od;
if s_cancel[1] = 0 then
return [ lcm(op(fld_denoms[1])), 0 ] ;
fi;
if s_cancel[2] = 0 then
return [0, lcm(op(fld_denoms[2]))] ;
fi;
scq := evala( Normal( s_cancel[1]/s_cancel[2] ) );
if type( scq, rational ) then
#
# The s_cancel parts differ by a multiplicative rational number, scq.
# Say scq = n1/n2; then s_cancel[1]/s_cancel[2] = n1/n2, so that
# n2*s_cancel[1] - n1*s_cancel[2] = 0 and therefore (by hypotheses)
# n2*s1 - n1*s2 is the log-derivative of the nth root of
# a rational function, so that
# n*n2*s1 - n*n1*s2 is the log-derivative of a rational function.
# Here, n is the LCM of the denominators of the rational numbers
# alpha1 found above.
#
n := lcm( op(fld_denoms[1]), op(fld_denoms[2]) );
return [ n * denom( scq ), (-1) * n * numer( scq ) ] ;
else
#
# The s_cancel parts fail to cancel
#
return [0, 0] ;
fi;
end:

compute_torus := proc( r1, r2, x )
#
# This procedure assumes that r1 and r2 are rational functions
# in the indeterminate x satisfying the following property:
# Let
 << Ïðåäûäóùàÿ ñòð. 30(èç 33 ñòð.)ÎÃËÀÂËÅÍÈÅ Ñëåäóþùàÿ >>