ClosedForm.ml

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

INCLUDE "preamble.ml"

let title _ = <:text<Closed Form>>


(* A wrapper function to call DDMF:-Hypergeometric:-generalTerms *)
(* that does some pre- and post-processing. *)
let call_generalTerms recsys inits u n x lps =
  (* Write recsys in terms of positive shifts (required by generalTerms). *)
  let shift =
    << min(map(a -> op(1, a) - $(n), indets($(recsys), function))) >> in
  let recsys = << subs($(n) = $(n) - $(shift), $(recsys)) >> in
  (* Compute closed forms of the unknown sequences u[i](n). *)
  let p = << indets([$(recsys), $(inits)], name) minus {$(n)} >> in
  let general_terms =
    <<
      DDMF:-Hypergeometric:-generalTerms($(recsys), $(inits), $(n), $(x), $(p))
    >> in
  << [seq($(general_terms)[$(u)[i]], i in $(lps))] >>


(* Compute a closed form expression of a transseries record. The closed form *)
(* is returned as a Maple expression. The record itself is updated, when the *)
(* closed form is computed for the first time. *)
let closed_form_of_transseries ser =
  (* Shorter name for the constituents of the transseries record. *)
  let x = << $(ser):-variable >>
  and ram = << $(ser):-ramification >>
  and u = << $(ser):-sequence_name >>
  and n = << $(ser):-index_name >> in

  (* A local variable to ensure that the final expression is displayed in *)
  (* terms of ser:-variable, i.e., the "transf" given in the sf database, *)
  (* and to make intermediate operations simpler. *)
  let x_loc = << proc() local x_loc; x_loc end proc() >> in

  (* Whether the power factor x^alpha is displayed inside the sum or not. *)
  let pow_factor, x_expon =
    if <:bool< type($(ram) * $(ser):-alpha, integer) >>
    then << 1 >>, << $(ser):-alpha >>
    else << combine($(x) ^ $(ser):-alpha, 'power', 'symbolic') >>, << 0 >> in

  (* Maple expression of an element in the list :-general. *)
  let expr_of_general_term gen =
    let x_part = << $(x_loc) ^ ($(gen):-period / $(ram) + $(x_expon)) >>
    and cf = << $(gen):-cform >> in
    let facs = << `if`(op(0, $(cf)) = `*`, [op($(cf))], [$(cf)]) >> in
    let facs = << [selectremove(has, $(facs), $(n))] >> in
    let out_fac = << `*`(op($(facs)[2])) >>
    and in_fac =
      << map(a -> sort(a, indets(a, name), descending), $(facs)[1]) >> in
    let in_fac = << __DynaMoW_times([op(sort($(in_fac))), $(x_part)]) >> in
    let in_fac = << DDMF:-dynamow_times_to_frac($(in_fac), []) >> in
    << $(out_fac) * Sum($(in_fac), $(n) = 0 .. infinity) >>
  in

  (* Update the transseries record, concerning the closed form sol which *)
  (* is related to index i (which corresponds to u[i](n) * log(x)^i). *)
  let update_index i sol =
    let ln = << ln($(x_loc))^$(i) >> in
    if <:bool< $(sol):-isHyper >> then
      (* If a closed form is found, it is stored in the transseries record. *)
      let gen = CommonTools.symb_list_of_symb << $(sol):-general >> in
      let exprs = List.map expr_of_general_term gen in
      let expr = << `+`(op($(CommonTools.symb_of_symb_list exprs))) >> in
      (* If the solution is finite (no infinite sum involved), it goes to *)
      (* ser:-finite_part, otherwise to ser:-closed_form. *)
      (* In that case, we also update the recurrences and initial values. *)
      if <:bool< $(expr) = 0 >> then
        let rec_inits =
          <<
            proc(xx, u, i)
              local poly, deg, ord, recsys, inits, new_inits, eqns, sol;
              poly := subs(xx = xx^$(ram), $(sol):-poly);
              poly := combine(poly, 'power', 'symbolic');
              deg := degree(poly, xx);
              recsys := $(ser):-recurrences;
              inits := $(ser):-initial_conditions;
              ord := -min(map(a -> op(a) - $(n),
                indets(recsys, specfunc(anything, u))));
              eqns := [seq(op(subs($(n) = ord + j, recsys)), j = 0 .. deg)];
              eqns := subs(
                [seq(u[i](j) = coeff(poly, xx, j), j = 0 .. ord + deg)], eqns);
              new_inits := remove(a -> evalb(eval(op(2, a), u = 0) = 0), eqns);
              eqns := subs(inits, eqns);
              sol := solve(eqns, indets(eqns, specfunc(integer, u)));
              if [sol] = [] then error("Error 666: no solution found.") end if;
              new_inits := select(has, [op(inits), op(sol)],
                indets(map2(op, 1, new_inits), specfunc(integer, u)));
              recsys := select(has, eval(recsys, u[i] = 0), u);
              inits := sort(remove(has, [op(inits), op(new_inits)], u[i]));
              return [recsys, inits];
            end proc($(x_loc), $(u), $(i))
          >> in
        <:unit<
          $(ser):-finite_part := $(ser):-finite_part + $(ln) * $(sol):-poly;
          $(ser):-log_powers := remove(`=`, $(ser):-log_powers, $(i));
          $(ser):-recurrences := $(rec_inits)[1];
          $(ser):-initial_conditions := $(rec_inits)[2];
        >>
      else
        <:unit<
          $(ser):-closed_form :=
            $(ser):-closed_form + $(ln) * ($(sol):-poly + $(expr));
        >>
    else
      (* If no closed form could be found, a symbolic expression in terms *)
      (* of the unknown coefficient sequences u[i](n) is stored in *)
      (* ser:-closed_form. *)
      let x_part = << $(x_loc) ^ ($(n) / $(ram) + $(x_expon)) >> in
      let smnd = << __DynaMoW_times([$(u)[$(i)]($(n)), $(x_part)]) >> in
      <:unit<
        $(ser):-closed_form :=
          $(ser):-closed_form + $(ln) * Sum($(smnd), $(n) = 0 .. infinity);
      >>
  in

  (* If this condition is true, it means that the closed form was never *)
  (* computed before and that we actually have to do something. *)
  if <:bool< $(ser):-closed_form = undefined >>
  then begin
    <:unit< $(ser):-closed_form := 0 >>;
    let lps = << $(ser):-log_powers >> in
    let closed_forms =
      call_generalTerms
        << $(ser):-recurrences >> << $(ser):-initial_conditions >>
        u n << $(x_loc)^(1 / $(ram)) >> lps in
    (* Update the transseries record for all indices i. *)
    List.iter
      (fun i ->
        update_index << $(lps)[$(int: i)] >> << $(closed_forms)[$(int: i)] >>)
      (CommonTools.range <:int< nops($(lps)) >>);
  end;

  (* Now in any case the closed form is available in the record. *)
  <:unit<
    $(ser):-finite_part :=
      DDMF:-replace_and_combine($(ser):-finite_part, $(x_loc), $(x));
    $(ser):-closed_form :=
      DDMF:-replace_and_combine($(ser):-closed_form, $(x_loc), $(x));
  >>;
  let cf = << $(ser):-finite_part + $(ser):-closed_form >> in
  let recsys, inits =
    if <:bool< has($(ser):-closed_form, $(u)) >>
    then << $(ser):-recurrences >>, << $(ser):-initial_conditions >>
    else << [] >>, << [] >> in
  let exp_factor = << combine($(ser):-exp_factor, 'power', 'symbolic') >> in
  << [$(exp_factor) * $(pow_factor) * $(cf), $(recsys), $(inits)] >>


let closed_form_of_linear_comb linear_comb =
  let sers = CommonTools.symb_list_of_symb << map(a -> a[2], $(linear_comb)) >>
  and coeffs = << map(a -> a[1], $(linear_comb)) >> in
  let closed_forms =
    CommonTools.symb_of_symb_list (List.map closed_form_of_transseries sers) in
  let recsys = << map(a -> op(a[2]), $(closed_forms)) >>
  and inits = << map(a -> op(a[3]), $(closed_forms)) >>
  and closed_forms = << map(a -> a[1], $(closed_forms)) >> in
  (* Assemble the linear combination of closed forms. *)
  let cf =
    <<
      `+`(op(zip((a, b) -> `if`(op(0, b) = `+`, map(c -> a * c, b), a * b),
        $(coeffs), $(closed_forms))))
    >> in
  let cf = << `if`(op(0, $(cf)) = `+`, [op($(cf))], [$(cf)]) >> in
  (* Sort the terms according to asymptotic behaviour. *)
  (* This is not perfectly done since for the sums it's difficult to see. *)
  let var = if sers = [] then << x >> else << $(List.hd sers):-variable >> in
  let cf = << selectremove(has, sort($(cf), length), Sum) >> in
  let cf =
    <<
      [op(sort($(cf)[2],
        proc(a, b) local da, db;
          da := DDMF:-my_degree(a, ln($(var)));
          db := DDMF:-my_degree(b, ln($(var)));
          `if`(da = db, [a, b] = sort([a, b]), da > db)
        end proc)), op($(cf)[1])]
    >> in
  (* Simplify coefficients, in particular in front of sums. *)
  let cf =
    <<
      map(a -> `if`(op(0, a) = `*`, ((b, c) -> combine(`*`(op(c))) *
        (`*`(op(b))))(selectremove(has, [op(a)], Sum)), a), $(cf))
    >> in
  (* Rewrite the parts of cf in terms of __DynaMoW_plus, _times and _frac. *)
  let cf =
    <<
      map(a -> `if`(op(0, a) = `*`, DDMF:-dynamow_times_to_frac(__DynaMoW_times(
        map(b -> op(sort(b, (a1, a2) -> `if`(has(a1, Sum) = has(a2, Sum),
          [a1, a2] = sort([a1, a2]), has(a2, Sum)))),
          [selectremove(b -> not(has(b, indets($(var)))), [op(a)])])),
        [Sum, $(var), 1/$(var)]), a), $(cf))
    >> in
  let cf = << `if`(nops($(cf)) > 1, __DynaMoW_plus($(cf)), `+`(op($(cf)))) >> in
  (* If there are unknown sequences (no closed form found), rename them, *)
  (* as to have consecutive indexing. *)
  let seq_names = << map(a -> a[2]:-sequence_name, $(linear_comb)) >> in
  let seq_names =
    << map(a -> `if`(op(0, a) = 'symbol', a, op(0, a)), $(seq_names)) >> in
  if <:bool< nops({op($(seq_names))}) = 1 >> then
    let rename =
      <<
        DDMF:-rename_sequences(
          [$(cf), $(recsys), $(inits)], $(seq_names)[1])[1]
      >> in
    << $(rename)[1] >>, << $(rename)[2] >>, << $(rename)[3] >>
  else cf, recsys, inits


let_service ClosedForm
  (linear_comb : any maple)
  (notation : string) :

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

  let linear_comb = << DDMF:-decode_linear_combination($(linear_comb)) >> in
  let cf, recsys, inits = closed_form_of_linear_comb linear_comb in
  let cf = << DDMF:-maple_to_dynamow($(cf)) >> in

  <:par<Taylor coefficients: <:dmath< $(str: notation) = $(symb: cf). >>>>,
  (cf, << DDMF:-encode_linear_combination($(linear_comb)) >>)

Generated by GNU Enscript 1.6.5.90.