1signature REAL_ASYMP_DIAG = sig
2
3val pretty_limit : Proof.context -> term -> Pretty.T
4
5val limit_cmd :
6  Proof.context -> (Facts.ref * Token.src list) list list -> string -> string option -> unit
7val limit : Proof.context -> thm list -> term -> term -> Multiseries_Expansion.limit_result
8
9val expansion_cmd :
10   Proof.context -> (Facts.ref * Token.src list) list list -> bool * int ->
11     string -> string option -> unit
12val expansion :
13  Proof.context -> thm list -> bool * int -> term -> term -> term * Asymptotic_Basis.basis
14
15end
16
17structure Real_Asymp_Diag : REAL_ASYMP_DIAG = struct
18
19open Lazy_Eval
20open Multiseries_Expansion
21
22fun pretty_limit _ (Const (\<^const_name>\<open>at_top\<close>, _)) = Pretty.str "\<infinity>"
23  | pretty_limit _ (Const (\<^const_name>\<open>at_bot\<close>, _)) = Pretty.str "-\<infinity>"
24  | pretty_limit _ (Const (\<^const_name>\<open>at_infinity\<close>, _)) = Pretty.str "\<plusminus>\<infinity>"
25  | pretty_limit ctxt (Const (\<^const_name>\<open>at_within\<close>, _) $ c $ 
26        (Const (\<^const_name>\<open>greaterThan\<close>, _) $ _)) = 
27      Pretty.block [Syntax.pretty_term ctxt c, Pretty.str "\<^sup>+"]
28  | pretty_limit ctxt (Const (\<^const_name>\<open>at_within\<close>, _) $ c $ 
29        (Const (\<^const_name>\<open>lessThan\<close>, _) $ _)) = 
30      Pretty.block [Syntax.pretty_term ctxt c, Pretty.str "\<^sup>-"]
31  | pretty_limit ctxt (Const (\<^const_name>\<open>at_within\<close>, _) $ c $ Const ("UNIV", _)) = 
32      Syntax.pretty_term ctxt c
33  | pretty_limit ctxt (Const (\<^const_name>\<open>nhds\<close>, _) $ c) =
34      Syntax.pretty_term ctxt c
35  | pretty_limit _ t = raise TERM ("pretty_limit", [t])
36
37fun reduce_to_at_top flt t = Envir.beta_eta_contract (
38    case flt of
39      \<^term>\<open>at_top :: real filter\<close> => t
40    | \<^term>\<open>at_bot :: real filter\<close> =>
41        Term.betapply (\<^term>\<open>%(f::real\<Rightarrow>real) x. f (-x)\<close>, t)
42    | \<^term>\<open>at_left 0 :: real filter\<close> =>
43        Term.betapply (\<^term>\<open>%(f::real\<Rightarrow>real) x. f (-inverse x)\<close>, t)
44    | \<^term>\<open>at_right 0 :: real filter\<close> =>
45        Term.betapply (\<^term>\<open>%(f::real\<Rightarrow>real) x. f (inverse x)\<close>, t)
46    | \<^term>\<open>at_within :: real => _\<close> $ c $ (\<^term>\<open>greaterThan :: real \<Rightarrow> _\<close> $ c') =>
47        if c aconv c' then
48          Term.betapply (Term.betapply (\<^term>\<open>%(f::real\<Rightarrow>real) c x. f (c + inverse x)\<close>, t), c)
49        else
50          raise TERM ("Unsupported filter for real_limit", [flt])
51    | \<^term>\<open>at_within :: real => _\<close> $ c $ (\<^term>\<open>lessThan :: real \<Rightarrow> _\<close> $ c') =>
52        if c aconv c' then
53          Term.betapply (Term.betapply (\<^term>\<open>%(f::real\<Rightarrow>real) c x. f (c - inverse x)\<close>, t), c)
54        else
55          raise TERM ("Unsupported filter for real_limit", [flt])
56    | _ =>
57        raise TERM ("Unsupported filter for real_limit", [flt]))
58
59fun mk_uminus (\<^term>\<open>uminus :: real => real\<close> $ c) = c
60  | mk_uminus c = Term.betapply (\<^term>\<open>uminus :: real => real\<close>, c)
61
62fun transfer_expansion_from_at_top' flt t = Envir.beta_eta_contract (
63    case flt of
64      \<^term>\<open>at_top :: real filter\<close> => t
65    | \<^term>\<open>at_bot :: real filter\<close> =>
66        Term.betapply (\<^term>\<open>%(f::real\<Rightarrow>real) x. f (-x)\<close>, t)
67    | \<^term>\<open>at_left 0 :: real filter\<close> =>
68        Term.betapply (\<^term>\<open>%(f::real\<Rightarrow>real) x. f (-inverse x)\<close>, t)
69    | \<^term>\<open>at_right 0 :: real filter\<close> =>
70        Term.betapply (\<^term>\<open>%(f::real\<Rightarrow>real) x. f (inverse x)\<close>, t)
71    | \<^term>\<open>at_within :: real => _\<close> $ c $ (\<^term>\<open>greaterThan :: real \<Rightarrow> _\<close> $ c') =>
72        if c aconv c' then
73          Term.betapply (Term.betapply (\<^term>\<open>%(f::real\<Rightarrow>real) c x. f (inverse (x - c))\<close>, t), c)
74        else
75          raise TERM ("Unsupported filter for real_limit", [flt])
76    | \<^term>\<open>at_within :: real => _\<close> $ c $ (\<^term>\<open>lessThan :: real \<Rightarrow> _\<close> $ c') =>
77        if c aconv c' then
78          Term.betapply (Term.betapply (\<^term>\<open>%(f::real\<Rightarrow>real) c x. f (inverse (c - x))\<close>, t), c)
79        else
80          raise TERM ("Unsupported filter for real_limit", [flt])
81    | _ =>
82        raise TERM ("Unsupported filter for real_limit", [flt]))
83
84
85fun transfer_expansion_from_at_top flt =
86  let
87    fun go idx (t as (\<^term>\<open>(powr) :: real => _\<close> $ 
88                 (\<^term>\<open>inverse :: real \<Rightarrow> _\<close> $ Bound n) $ e)) =
89          if n = idx then
90            Envir.beta_eta_contract (\<^term>\<open>(powr) :: real => _\<close> $ Bound n $ mk_uminus e)
91          else
92            t
93      | go idx (t as (\<^term>\<open>(powr) :: real => _\<close> $ (\<^term>\<open>uminus :: real \<Rightarrow> real\<close> $
94              (\<^term>\<open>inverse :: real \<Rightarrow> _\<close> $ Bound n)) $ e)) =
95          if n = idx then
96            Envir.beta_eta_contract (\<^term>\<open>(powr) :: real => _\<close> $
97              (mk_uminus (Bound n)) $ mk_uminus e)
98          else
99            t
100      | go idx (t as (\<^term>\<open>(powr) :: real => _\<close> $ (\<^term>\<open>inverse :: real \<Rightarrow> _\<close> $ 
101              (\<^term>\<open>(-) :: real \<Rightarrow> _\<close> $ Bound n $ c)) $ e)) =
102          if n = idx then
103            Envir.beta_eta_contract (\<^term>\<open>(powr) :: real => _\<close> $
104              (\<^term>\<open>(-) :: real => _\<close> $ Bound n $ c) $ mk_uminus e)
105          else
106            t
107      | go idx (t as (\<^term>\<open>(powr) :: real => _\<close> $ (\<^term>\<open>inverse :: real \<Rightarrow> _\<close> $ 
108              (\<^term>\<open>(-) :: real \<Rightarrow> _\<close> $ c $ Bound n)) $ e)) =
109          if n = idx then
110            Envir.beta_eta_contract (\<^term>\<open>(powr) :: real => _\<close> $
111              (\<^term>\<open>(-) :: real => _\<close> $ c $ Bound n) $ mk_uminus e)
112          else
113            t
114      | go idx (s $ t) = go idx s $ go idx t
115      | go idx (Abs (x, T, t)) = Abs (x, T, go (idx + 1) t)
116      | go _ t = t
117  in
118    transfer_expansion_from_at_top' flt #> go (~1)
119  end
120
121fun gen_limit prep_term prep_flt prep_facts res ctxt facts t flt =
122  let
123    val t = prep_term ctxt t
124    val flt = prep_flt ctxt flt
125    val ctxt = Proof_Context.augment t ctxt
126    val t = reduce_to_at_top flt t
127    val facts = prep_facts ctxt facts
128    val ectxt = mk_eval_ctxt ctxt |> add_facts facts |> set_verbose true
129    val (bnds, basis) = expand_term_bounds ectxt t Asymptotic_Basis.default_basis
130  in
131    res ctxt (limit_of_expansion_bounds ectxt (bnds, basis))
132  end
133
134fun pretty_limit_result ctxt (Exact_Limit lim) = pretty_limit ctxt lim
135  | pretty_limit_result ctxt (Limit_Bounds (l, u)) =
136      let
137        fun pretty s (SOME l) = [Pretty.block [Pretty.str s, pretty_limit ctxt l]]
138          | pretty _ NONE = []
139        val ps = pretty "Lower bound: " l @ pretty "Upper bound: " u
140      in
141        if null ps then Pretty.str "No bounds found." else Pretty.chunks ps
142      end
143
144val limit_cmd =
145  gen_limit 
146    (fn ctxt => 
147      Syntax.parse_term ctxt 
148      #> Type.constraint \<^typ>\<open>real => real\<close> 
149      #> Syntax.check_term ctxt)
150    (fn ctxt => fn flt =>
151        case flt of
152          NONE => \<^term>\<open>at_top :: real filter\<close>
153        | SOME flt =>
154            flt
155            |> Syntax.parse_term ctxt
156            |> Type.constraint \<^typ>\<open>real filter\<close>
157            |> Syntax.check_term ctxt)
158    (fn ctxt => flat o flat o map (map (Proof_Context.get_fact ctxt o fst)))
159    (fn ctxt => pretty_limit_result ctxt #> Pretty.writeln)
160
161val limit = gen_limit Syntax.check_term Syntax.check_term (K I) (K I)
162
163
164fun gen_expansion prep_term prep_flt prep_facts res ctxt facts (n, strict) t flt =
165  let
166    val t = prep_term ctxt t
167    val flt = prep_flt ctxt flt
168    val ctxt = Proof_Context.augment t ctxt
169    val t = reduce_to_at_top flt t
170    val facts = prep_facts ctxt facts
171    val ectxt = mk_eval_ctxt ctxt |> add_facts facts |> set_verbose true
172    val basis = Asymptotic_Basis.default_basis
173    val (thm, basis) = expand_term ectxt t basis
174    val (exp, error) = extract_terms (strict, n) ectxt basis (get_expansion thm)
175    val exp = transfer_expansion_from_at_top flt exp
176    val error =
177      case error of
178        SOME (L $ flt' $ f) => SOME (L $ flt' $ transfer_expansion_from_at_top flt f)
179      | _ => error
180    val t =
181      case error of
182        NONE => exp
183      | SOME err =>
184          Term.betapply (Term.betapply (\<^term>\<open>expansion_with_remainder_term\<close>, exp), err)
185  in
186    res ctxt (t, basis)
187  end
188
189fun print_basis_elem ctxt t =
190  Syntax.pretty_term (Config.put Syntax_Trans.eta_contract false ctxt)
191    (Envir.eta_long [] t)
192
193val expansion_cmd =
194  gen_expansion
195    (fn ctxt => 
196      Syntax.parse_term ctxt 
197      #> Type.constraint \<^typ>\<open>real => real\<close> 
198      #> Syntax.check_term ctxt)
199    (fn ctxt => fn flt =>
200        case flt of
201          NONE => \<^term>\<open>at_top :: real filter\<close>
202        | SOME flt =>
203            flt
204            |> Syntax.parse_term ctxt
205            |> Type.constraint \<^typ>\<open>real filter\<close>
206            |> Syntax.check_term ctxt)
207    (fn ctxt => flat o flat o map (map (Proof_Context.get_fact ctxt o fst)))
208    (fn ctxt => fn (exp, basis) =>
209       Pretty.writeln (Pretty.chunks (
210         [Pretty.str "Expansion:",
211          Pretty.indent 2 (Syntax.pretty_term ctxt exp),
212          Pretty.str "Basis:"] @
213            map (Pretty.indent 2 o print_basis_elem ctxt) (Asymptotic_Basis.get_basis_list basis))))
214
215val expansion = gen_expansion Syntax.check_term (K I) (K I) (K I)
216
217end
218
219local 
220
221fun parse_opts opts dflt =
222  let
223    val p = map (fn (s, p) => Args.$$$ s |-- Args.colon |-- p) opts
224  in
225    Scan.repeat (Scan.first p) >> (fn xs => fold I xs dflt)
226  end
227
228val limit_opts =
229  [("limit", Parse.term >> (fn t => fn {facts, ...} => {limit = SOME t, facts = facts})),
230   ("facts", Parse.and_list1 Parse.thms1 >>
231     (fn thms => fn {limit, facts} => {limit = limit, facts = facts @ thms}))]
232
233val dflt_limit_opts = {limit = NONE, facts = []}
234
235val expansion_opts =
236  [("limit", Parse.term >> (fn t => fn {terms, facts, ...} =>
237     {limit = SOME t, terms = terms, facts = facts})),
238   ("facts", Parse.and_list1 Parse.thms1 >>
239     (fn thms => fn {limit, terms, facts} =>
240       {limit = limit, terms = terms, facts = facts @ thms})),
241   ("terms", Parse.nat -- Scan.optional (Args.parens (Args.$$$ "strict") >> K true) false >>
242     (fn trms => fn {limit, facts, ...} => {limit = limit, terms = trms, facts = facts}))]
243
244val dflt_expansion_opts = {limit = NONE, facts = [], terms = (3, false)}
245
246in
247
248val _ =
249  Outer_Syntax.command \<^command_keyword>\<open>real_limit\<close>
250    "semi-automatically compute limits of real functions"
251    ((Parse.term -- parse_opts limit_opts dflt_limit_opts) >>
252    (fn (t, {limit = flt, facts = thms}) => 
253      (Toplevel.keep (fn state => 
254        Real_Asymp_Diag.limit_cmd (Toplevel.context_of state) thms t flt))))
255
256val _ =
257  Outer_Syntax.command \<^command_keyword>\<open>real_expansion\<close>
258    "semi-automatically compute expansions of real functions"
259    (Parse.term -- parse_opts expansion_opts dflt_expansion_opts >> 
260    (fn (t, {limit = flt, terms = n_strict, facts = thms}) => 
261      (Toplevel.keep (fn state => 
262        Real_Asymp_Diag.expansion_cmd (Toplevel.context_of state) thms (swap n_strict) t flt))))
263
264end