TruncatedSeries.ml

(* Copyright INRIA and Microsoft Corporation, 2008-2013. *)
(* DDMF is distributed under CeCILL-B license. *)

INCLUDE "preamble.ml"

exception BadRamification

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.