ñòð. 30 |

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

# using addition, subtraction, multiplication and addition --

# 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 |