# TruncatedSeries.ml

```(* Copyright INRIA and Microsoft Corporation, 2008-2013. *)

INCLUDE "preamble.ml"

let title _ = <:text<Truncated Series>>

(* Compute the truncated version (without prefactor) of a transseries record, *)
(* up to the order ord and w.r.t. the variable var. *)
(* We don't use tser:-variable on purpose, for having more flexibility. *)
(* This is exploited below (see the use of my_var). *)
let truncate_transseries ord tser var =
<<
proc(s, ord, x, \$)
local i, j, u, n, r, lp, fin, lnx, smnd, trunc,
recs, shifts, m, rec_ord, sys, sol;
u := s:-sequence_name;
n := s:-index_name;
r := s:-ramification;
lp := s:-log_powers;

# Deal with the finite part and truncate it.
fin := DDMF:-transform_back(s:-finite_part, s:-variable, x);
fin := expand(subs(ln(x^(1/r)) = lnx, fin));
fin := `if`(op(0, fin) = `+`, [op(fin)], [fin]);
fin := select(a -> DDMF:-my_degree(a, x) < ord, fin);
fin := subs(lnx = ln(x^(1/r)), `+`(op(fin)));
if lp = [] then return x ^ s:-alpha * fin end if;

# Construct the truncated expression with undetermined coefficients.
smnd := add(u[lp[j]](n) * ln(x^(1/r))^lp[j], j = 1 .. nops(lp)) * x^(n/r);
trunc := add(subs(n = i, smnd), i = 0 .. r*ord - 1);
trunc := subs(s:-initial_conditions, trunc);

# The maximal order of the recurrences.
recs := s:-recurrences;
shifts := map(a -> op(a) - n, indets(recs, specfunc(anything, u)));
m := max(shifts);
if m <> 0 then recs := subs(n = n - m, recs) end if;
rec_ord := m - min(shifts);

# Solve for the undetermined coefficients and plug in.
# TODO Refine! Don't solve this big system in one stroke.
sys := [seq(op(subs(n = i, recs)), i = rec_ord .. r*ord - 1)];
sys := select(has, subs(s:-initial_conditions, sys), u);
sol := solve(sys, [seq(seq(u[lp[j]](i), j = 1 .. nops(lp)),
i = rec_ord .. r*ord - 1)]);
if nops(sol) <> 1 then error "Error in truncate_transseries." end if;
trunc := subs(sol[1], trunc);

return x ^ s:-alpha * (fin + trunc);
end proc(\$(tser), \$(int: ord), \$(var))
>>

(* Deal with a linear combination of equivalent (infinite) transseries, *)
(* i.e., sharing the same prefactor modulo an integer power of x, *)
(* and compute its truncation of order ord. *)
let deal_with_eqv ord eqv =
let alpha1 = << (\$(eqv)[1,2]):-alpha >> in
let min_alpha =
<< \$(alpha1) + min(map(a -> a[2]:-alpha - \$(alpha1), \$(eqv))) >> in

(* If alpha is an integer, it must be taken into account for the O-term. *)
let corr_alpha =
if <:bool< type(\$(min_alpha), integer) >>
then << 0 >> else << \$(min_alpha) >> in

(* Introduce a local variable, this makes simplifications more convenient. *)
let my_var = << proc() local my_var; my_var end proc() >> in

(* Compute the truncation for each transseries with the appropriate order. *)
let truncs =
List.map
(fun a ->
truncate_transseries
(ord - <:int< \$(a):-alpha - \$(corr_alpha) >>) a my_var)
(CommonTools.symb_list_of_symb << map2(op, 2, \$(eqv)) >>) in
let truncs = CommonTools.symb_of_symb_list truncs in

(* Combine the truncations, do some syntactic simplifications that *)
(* are supposed to yield nicer-looking formulas, and add O-term. *)
let trunc = << add(\$(eqv)[i,1] * \$(truncs)[i], i = 1 .. nops(\$(eqv))) >>
and rams = << map(a -> a[2]:-ramification, \$(eqv)) >> in
(* Test whether all transseries have the same ramification. *)
if <:bool< nops({op(\$(rams))}) <> 1 >> then raise BadRamification;
let ram = << \$(rams)[1] >> in
let trunc =
<<
proc(x)
local expr, lnx, cont, terms, par, norm_sum, norm_fac, norm_term,
ln_power;
expr := subs(ln(x ^ (1 / \$(ram))) = lnx, \$(trunc) / x ^ \$(corr_alpha));
expr := collect(expr, x, factor);

# Extract content
if op(0, expr) = `+` then
cont := map(a -> `if`(op(0, a) = `*`, {op(a)}, {a}), [op(expr)]);
cont := `*`(op(remove(has, `intersect`(op(cont)), [x, lnx])));
expr := map(a -> a / cont, expr);
terms := [op(expr)];
else cont := 1; terms := [expr] end if;
terms := map(a -> `if`(DDMF:-my_degree(a, x) = 0 and degree(a, lnx) > 0,
expand(a), a), terms);
terms := map(a -> `if`(op(0, a) = `+`, op(a), a), terms);

# Sort the expression according to asymptotic behaviour.
terms := sort(terms,
proc(a, b) local da, db;
da := DDMF:-my_degree(a, x);
db := DDMF:-my_degree(b, x);
if da = db then da := -degree(a, lnx); db := -degree(b, lnx) end if;
`if`(da = db, [a, b] = sort([a, b]), da < db)
end proc);
par := [lnx, op(sort(remove(has, indets(terms, name), [x, lnx])))];
norm_sum := proc(s)
__DynaMoW_plus(map(a -> `if`(op(0, a) = `*`,
DDMF:-dynamow_times_to_frac(__DynaMoW_times([op(a)]), lnx),
`if`(type(a, complex(rational)), DDMF:-complex_to_frac(a), a)),
[op(sort(s, par, descending))]))
end proc;
norm_fac := proc(f)
if op(0, f) = `+` then return norm_sum(f) end if;
if op(0, f) = `^` and op(0, op(1, f)) = `+` then
return norm_sum(op(1, f)) ^ op(2, f);
end if;
return f;
end proc;
norm_term := proc(t)
sort(map(norm_fac, [op(t)]),
proc(a, b) local t1, t2;
t1 := `if`(has(a, x), 2, `if`(has(a, lnx), 1, 0));
t2 := `if`(has(b, x), 2, `if`(has(b, lnx), 1, 0));
`if`(t1 = t2, [a, b] = sort([a, b]), t1 < t2)
end proc);
end proc;
terms := map(a ->
`if`(op(0, a) = `*`, DDMF:-dynamow_times_to_frac(
__DynaMoW_times(norm_term(a)), [x, lnx]),
`if`(type(a, complex(rational)), DDMF:-complex_to_frac(a), a)),
terms);

# Put everything together and add O-term if needed.
cont := factor(cont * x ^ \$(corr_alpha));
if map(a -> op(a[2]:-recurrences), \$(eqv)) <> [] then
terms := remove(`=`, terms, 0);
ln_power := max(0, map(a -> op(a[2]:-log_powers), \$(eqv)));
terms := [op(terms), O(lnx ^ ln_power * x ^ \$(int: ord))];
end if;
if nops(terms) = 1 then
expr := cont * terms[1];
expr := __DynaMoW_times(`if`(op(0, expr) = `*`, [op(expr)], [expr]));
else
cont := `if`(op(0, cont) = `*`, [op(cont)], [cont]);
terms := map(a -> `if`(op(0, a) = `+`, op(a), a), terms);
terms := map(a -> `if`(op(0, a) = `*`, DDMF:-dynamow_times_to_frac(
__DynaMoW_times(sort([op(a)])), [x, lnx]), a), terms);
expr := __DynaMoW_times([op(cont), __DynaMoW_plus(terms)]);
end if;
lnx := ln(x ^ (1 / \$(ram)));
return expr;
end proc(\$(my_var))
>> in

(* Replace local variable by var and assemble the final expression. *)
let var = << \$(eqv)[1,2]:-variable >> in
let trunc = << DDMF:-replace_and_combine(\$(trunc), \$(my_var), \$(var)) >> in
let ef = << combine(\$(eqv)[1,2]:-exp_factor, 'power', 'symbolic') >> in
let ef = << `if`(op(0, \$(ef)) = `*`, [op(\$(ef))], [\$(ef)]) >> in
let trunc = << remove(`=`, [op(\$(ef)), op(op(1, \$(trunc)))], 1) >> in
let trunc =
<<
sort(\$(trunc), (a, b) ->
`if`(has(a, O) = has(b, O), [a, b] = sort([a, b]), has(b, O)))
>> in
<<
`if`(nops(\$(trunc)) > 1, DDMF:-dynamow_times_to_frac(
__DynaMoW_times(\$(trunc)), [O, \$(var), 1/\$(var)]), `*`(op(\$(trunc))))
>>

let_service TruncatedSeries
(linear_comb : any maple)
(notation : string)
(order : int = Constants.default_order) :

DC.sec_entities * any maple with { title = title } =

let linear_comb = << DDMF:-decode_linear_combination(\$(linear_comb)) >> in

(* Sort the transseries into equivalence classes such that all series *)
(* in a class have the same prefactor (= exp_factor * x^alpha) modulo *)
(* an integer power of x. *)
let eqvs =
<<
DDMF:-equivalence_classes(\$(linear_comb),
(a, b) -> evalb(a[2]:-exp_factor = b[2]:-exp_factor and
type(a[2]:-alpha - b[2]:-alpha, integer)))
>> in

(* Compute the truncated series for each class separately and add them up. *)
(* We don't expect that any simplification is possible, except of the type *)
(* exp(I*x)*(...) + exp(-I*x)*(...). *)
(* TODO: Deal with such special cases and rewrite with sin and cos? *)
let truncs =
CommonTools.symb_of_symb_list
(List.map (deal_with_eqv order) (CommonTools.symb_list_of_symb eqvs)) in
let first_terms =
<<
`if`(nops(\$(truncs)) > 1,
__DynaMoW_plus(sort(\$(truncs), length)), `+`(op(\$(truncs))))
>> in
let first_terms = << DDMF:-maple_to_dynamow(\$(first_terms)) >> in

<:par<First terms:<:dmath< \$(str: notation) = \$(symb: first_terms). >>>>,
first_terms
```

Generated by GNU Enscript 1.6.5.90.