(* ========================================================================== *)
(* FLYSPECK - BOOK FORMALIZATION                                              *)
(*                                                                            *)
(* Lemma: Taylor Series for atn function                                      *)
(* Chapter: Nonlinear Inequalities                                            *)
(* Author: Thomas C. Hales                                                    *)
(* Date: 2010-07-14                                                           *)
(* ========================================================================== *)


(*
This file gives the half-angle identity for atan
    atn (2 x) = 2 atn (...)

It gives a general formula for the nth derivative of catn

It gives the complex Taylor polynomial of catn at Cx(&0).

It gives the real Taylor polynomial of atn at (&0)
*)

module Taylor_atn = struct


let FORCE_EQ = REPEAT (CHANGED_TAC (AP_TERM_TAC ORELSE AP_THM_TAC ORELSE BINOP_TAC)) ;;

let FORCE_MATCH = (MATCH_MP_TAC (TAUT `(a = b) ==> (a ==> b)`)) THEN FORCE_EQ ;;

let FORCE_MATCH_MP_TAC th = 
  MP_TAC th THEN ANTS_TAC THENL[ALL_TAC;FORCE_MATCH
      ];;

(* first we develop the half-angle identity for the atn function *)


let halfatn = new_definition `halfatn x = x / (sqrt(&1 + x pow 2) + &1)`;;


let pos1 = 
prove(
  `!x. &0 < &1 + x pow 2 `,

  MESON_TAC[Collect_geom.pow_g;REAL_ARITH `&0 <= t ==> &0 < &1 + t`];
);;


let ssqrt = new_definition `!x. ssqrt x = (if x < &0 then &0 else sqrt x)`;;


let halfsqrt_ssqrt = 
prove_by_refinement(
  `!x.  sqrt(&1+ x pow 2) = ssqrt(&1 + x pow 2)`,

  (* {{{ proof *)
  [
  REWRITE_TAC[ssqrt;];
  MESON_TAC[pos1;REAL_ARITH `&0 < x ==> ~(x < &0)`];
  ]);;
  (* }}} *)


let pos2 = 
prove (
  `!x. &0 < sqrt(&1 + x pow 2) + &1`,

(* {{{ proof *)
   MESON_TAC[pos1;REAL_ARITH `(&0 <= t ==> &0 < t + &1) /\ (&0 < t ==> &0 <= t)`;SQRT_POS_LE;]
(* }}} *));;


let halfatn_bounds_abs = 
prove_by_refinement(
 `!x.  abs(halfatn x) < &1 `,

(* {{{ proof *)
[
REWRITE_TAC[halfatn;REAL_ABS_DIV];
GEN_TAC;
ASSUME_TAC (ISPEC `x:real` pos2);
ASM_SIMP_TAC[REAL_ARITH `(&0 < x ==> abs(x) = x)/\ (&1 * t = t)`;REAL_LT_LDIV_EQ];
(* *)
REWRITE_TAC[GSYM POW_2_SQRT_ABS];
MATCH_MP_TAC REAL_LET_TRANS;
EXISTS_TAC `sqrt(&1 + x pow 2)`;
CONJ_TAC;
MATCH_MP_TAC SQRT_MONO_LE;
REWRITE_TAC[Collect_geom.pow_g];
ARITH_TAC;
ARITH_TAC;
]
(* }}} *));;


let halfatn_bounds = 
prove(
 `!x. -- &1 < halfatn x /\ halfatn x < &1 `,


let halfatn_half = 
prove_by_refinement(
 `!x t. (abs (x) < t ==> abs(halfatn x) < t / &2) `,

(* {{{ proof *)
[
REWRITE_TAC[halfatn;REAL_ABS_DIV];
REPEAT STRIP_TAC;
ASSUME_TAC (ISPEC `x:real` pos2);
ASM_SIMP_TAC[REAL_ARITH `&0 < x ==> &0 < abs(x)`;REAL_LT_LDIV_EQ];
MATCH_MP_TAC REAL_LTE_TRANS;
EXISTS_TAC `t:real`;
ASM_REWRITE_TAC[REAL_ARITH `(t / &2 * x = t * (x / &2)) /\ (t <= t * x/ &2 <=> t * &2 <= t * x)`];
MATCH_MP_TAC REAL_LE_LMUL;
CONJ_TAC;
UNDISCH_TAC `abs x < t`;
REAL_ARITH_TAC;
ASM_SIMP_TAC [REAL_ARITH `&0 < x ==> (abs(x) = x)`];
REWRITE_TAC[REAL_ARITH `(&2 <= x + &1) = (&1 <= x)`];
MATCH_MP_TAC REAL_LE_TRANS;
EXISTS_TAC `sqrt(&1)`;
CONJ_TAC;
REWRITE_TAC[SQRT_1;REAL_ARITH `&1 <= &1`];
MATCH_MP_TAC SQRT_MONO_LE;
CONJ_TAC;
ARITH_TAC;
MATCH_MP_TAC (REAL_ARITH `&0 <= x ==> &1 <= &1 + x`);
REWRITE_TAC[Collect_geom.pow_g];
]
(* }}} *));;


let abs_pass_through = 
prove_by_refinement (
 `(!x f.  (f (-- x) = -- f x) /\ (!y. &0 <= y ==> &0 <= f y)
     ==> (abs (f x) = f (abs x)))`,

(* {{{ proof *)
 [
 REPEAT STRIP_TAC;
 DISJ_CASES_TAC (REAL_ARITH `&0 <= x \/ &0 <= --x`);
   POP_ASSUM MP_TAC;
   POP_ASSUM MP_TAC;
  MESON_TAC[REAL_ARITH `&0 <= x ==> abs(x ) = x`];
 REPEAT (POP_ASSUM MP_TAC);
  MESON_TAC[REAL_ARITH `&0 <= --x ==> abs(x) = -- x`; REAL_ARITH `abs( -- x ) = abs(x)`];
 ]
(* }}} *));;


let atn_abs = 
prove_by_refinement(
 `!x. abs(atn x) = atn (abs x) `,

(* {{{ proof *)
[
GEN_TAC;
MATCH_MP_TAC abs_pass_through;
REWRITE_TAC[ATN_NEG;ATN_POS_LE];
]
(* }}} *));;


let atn_half_range = 
prove_by_refinement (
 `!x. abs(atn (halfatn x)) < pi / &4 `,

(* {{{ proof *)
[
REWRITE_TAC[GSYM ATN_1;atn_abs;ATN_MONO_LT_EQ];
GEN_TAC;
REWRITE_TAC [halfatn_bounds;GSYM REAL_BOUNDS_LT];
]
(* }}} *));;


let tan_one_one = 
prove_by_refinement(
 `!x y. (abs(x) < pi/ &2 /\ (abs y < pi / &2 ) /\ (tan x = tan y) ==> (x = y))`,

(* {{{ proof *)
[
REPEAT STRIP_TAC;
DISJ_CASES_TAC (REAL_ARITH `x < y \/ y < x \/ x = (y:real)`);
REPEAT (POP_ASSUM MP_TAC);
REWRITE_TAC[GSYM REAL_BOUNDS_LT];
MESON_TAC[TAN_MONO_LT_EQ;REAL_ARITH `(x:real<y ) ==> ~(x = y)`];
POP_ASSUM DISJ_CASES_TAC;
REPEAT (POP_ASSUM MP_TAC);
REWRITE_TAC[GSYM REAL_BOUNDS_LT];
MESON_TAC[TAN_MONO_LT_EQ;REAL_ARITH `(x:real<y ) ==> ~(x = y)`];
ASM_REWRITE_TAC[];
]
(* }}} *));;


let abs_lemma = 
prove(
 `!f x. (?n. x = f n) \/ (?n. x = -- f n) <=> (?n. abs(x) = abs(f n))`,

  ASM_MESON_TAC[REAL_ARITH `!x y. abs(x) = abs(y) <=> (x = y)\/ (x = -- y)`]);;


let cos_nz = 
prove_by_refinement (
  `!x. (abs(x) < pi / &2) ==> ~(cos x = &0) `,

(* {{{ proof *)

[
GEN_TAC;
REWRITE_TAC[COS_ZERO_PI;abs_lemma];
ONCE_REWRITE_TAC[TAUT `(a ==> ~b) <=> (b ==> ~a)`];
ONCE_REWRITE_TAC[REAL_ARITH `~(x < y) <=> (y <= x)`];
REPEAT STRIP_TAC;
ASM_REWRITE_TAC[];
MATCH_MP_TAC REAL_LE_TRANS;
EXISTS_TAC`(&n +  &1/ &2) * pi `;
REWRITE_TAC[REAL_ARITH `x <= abs(x)`];
MP_TAC PI_POS;
MP_TAC (REAL_ARITH `&1/ &2 <= (&n + &1/ &2)`);
REWRITE_TAC[REAL_ARITH `pi / &2 = (&1 / &2) * pi`];
ASM_SIMP_TAC[REAL_ARITH `&0 < x ==> &0 <= x`;Real_ext.REAL_LE_RMUL_IMP];
]

(* }}} *));;


let cos_2nz = 
prove_by_refinement(
   `!x. (abs(x) < pi / &4) ==> ~(cos (&2 * x) = &0) `,

(* {{{ proof *)
 [
 STRIP_TAC THEN STRIP_TAC;
 MATCH_MP_TAC cos_nz;
 REWRITE_TAC[REAL_ABS_MUL;REAL_ARITH `abs(&2)= &2 /\ (&2 * x < pi/ &2 <=> x < pi/ &4)`];
 ASM_REWRITE_TAC[];
 ]
(* }}} *));;


let halfatn_double  =
prove_by_refinement(
 `!x. ~(cos (atn (halfatn x)) = &0) /\ ~(cos(&2 * atn (halfatn x)) = &0) `,

(* {{{ proof *)
[
REPLICATE_TAC 2 (STRIP_TAC);
MATCH_MP_TAC cos_nz;
MATCH_MP_TAC REAL_LTE_TRANS;
EXISTS_TAC `pi/ &4`;
REWRITE_TAC[atn_half_range];
MP_TAC PI_POS;
REAL_ARITH_TAC;
MATCH_MP_TAC cos_2nz;
REWRITE_TAC[atn_half_range];
]
(* }}} *));;

let REAL_DIV_MUL2z =  REAL_FIELD 
 `!x y z. (&0 < x)  ==> (y /z = (x pow 2 * y) / (x pow 2* z)) `;;


let atn_half = 
prove_by_refinement (
 `!x. atn x = &2 * atn (halfatn x) `,

(* {{{ proof *)
[
GEN_TAC;
MATCH_MP_TAC tan_one_one;
MATCH_MP_TAC (TAUT `a /\ b /\ (a /\ b ==> c) ==> a /\ b /\ c`);
REPEAT CONJ_TAC;
REWRITE_TAC[GSYM REAL_BOUNDS_LT;ATN_BOUNDS];
(* *)
REWRITE_TAC[REAL_ABS_MUL;REAL_ARITH `abs (&2) = &2`;REAL_ARITH `&2 * x < y / &2 <=> x < y / &4`;atn_half_range];
(* *)
REPEAT STRIP_TAC;
ASSUME_TAC (ISPEC `x:real` halfatn_double);
ASM_SIMP_TAC[TAN_DOUBLE;ATN_TAN];
REWRITE_TAC[halfatn];
ASSUME_TAC (ISPEC `x:real` pos2);
ABBREV_TAC `t = sqrt(&1 + x pow 2) + &1`;
MP_TAC (ISPECL [`t:real`;`&2 * x / t`;`&1 - (x / t) pow 2`] REAL_DIV_MUL2z);
ASM_REWRITE_TAC[];
DISCH_THEN (fun t-> REWRITE_TAC[t]);
ASM_SIMP_TAC[REAL_FIELD `&0 < t ==> t pow 2 * &2 * x / t = t * &2 * x`];
ASM_SIMP_TAC[REAL_FIELD `&0 < t ==> t pow 2 * (&1 - (x / t) pow 2) = t pow 2 - x pow 2`];
EXPAND_TAC "t";
REWRITE_TAC[REAL_FIELD `(a + &1) pow 2 = a pow 2 + &2 * a  + &1`];
ASM_SIMP_TAC[pos1;REAL_ARITH `!x. &0 < x ==> &0 <= x`;SQRT_POW_2];
ASM_REWRITE_TAC[REAL_ARITH `((&1 + v) + &2 * u + &1) - v = (u + &1) * &2`];
UNDISCH_TAC `&0 < t`;
CONV_TAC REAL_FIELD;
]
(* }}} *));;

(* complex taylor for atn *)

prioritize_complex();;

let id1 = COMPLEX_RING `inv (Cx (&1) + z pow 2)   = (inv (Cx (&2))) * ( ( inv (Cx (&1) + z pow 2) * (Cx (&1) - ii *z)) +  (inv (Cx (&1) + z pow 2)) * ( (Cx (&1) + ii * z)))`;;

let id2 = SIMPLE_COMPLEX_ARITH  ` (Cx (&1) + ii * z) * (Cx (&1) - ii * z) = (Cx (&1) - ii * ii * z * z)`;;


let id3 = 
prove_by_refinement (`!u a. a - ii * ii * u = a + u`,

(* {{{ proof *)
[
REWRITE_TAC[ii];
SIMPLE_COMPLEX_ARITH_TAC;
]
(* }}} *));;


let id4 = 
prove_by_refinement (`!z. (Cx (&1) + z pow 2) = (Cx (&1) + ii*z) * (Cx (&1) - ii*z)`,

(* {{{ proof *)
[
REWRITE_TAC[id2;id3;COMPLEX_POW_2];
]
(* }}} *));;


let tactic_list = [SUBGOAL_THEN `(Re (z) = &0) /\ (abs(Im (z)) = &1)` ASSUME_TAC ;
POP_ASSUM MP_TAC ;
ASM_REWRITE_TAC[] ;
REWRITE_TAC[REAL_ARITH `abs(x) = &1 <=> (x = &1 \/ x = -- &1)`;ii] ;
SIMPLE_COMPLEX_ARITH_TAC ;
ASM_MESON_TAC[]];;


let idz = 
prove_by_refinement(
  `!z a.  (Re z = &0 ==> abs(Im z) < &1) /\ ((a = ii \/ a = --ii))  ==>
     ~(Cx (&1) + a * z = Cx (&0))`,

  (* {{{ proof *)
  [
ASSUME_TAC (REAL_ARITH `~(&1 < &1)`);
REPEAT STRIP_TAC;
] @ (tactic_list @ tactic_list));;
  (* }}} *)



let id4a = 
prove_by_refinement (`!z. (Re z = &0 ==> abs(Im z) < &1)
  ==>( inv(Cx (&1) + ii* z) * (Cx (&1)  + ii*z) = Cx (&1))`,

(* {{{ proof *)
[
REPEAT STRIP_TAC;
MATCH_MP_TAC COMPLEX_MUL_LINV;
MATCH_MP_TAC idz;
ASM_REWRITE_TAC[];
]
(* }}} *));;


let id4b = 
prove_by_refinement (`!z. (Re z = &0 ==> abs(Im z) < &1)
  ==>( inv(Cx (&1) -  ii* z) * (Cx (&1)  - ii*z) = Cx (&1))`,

(* {{{ proof *)
[
REPEAT STRIP_TAC;
MATCH_MP_TAC COMPLEX_MUL_LINV;
REWRITE_TAC[SIMPLE_COMPLEX_ARITH `a - ii * z = a + (-- ii) * z`];
MATCH_MP_TAC idz;
ASM_REWRITE_TAC[];
]
(* }}} *));;


let id5 = 
prove_by_refinement (`!z. (Re z = &0 ==> abs(Im z) < &1) ==> ( inv (Cx (&1) + z pow 2) = (inv (Cx (&2))) * ( inv (Cx (&1) + ii * z) + inv (Cx (&1) - ii * z)))`,

(* {{{ proof *)
[
REPEAT STRIP_TAC;
ONCE_REWRITE_TAC[id1];
REWRITE_TAC[id4;COMPLEX_INV_MUL];
REWRITE_TAC[SIMPLE_COMPLEX_ARITH `((a*b)*c + (e*f)*g = (a:complex)*(b*c) + f * (e *g))`];
ASM_SIMP_TAC[id4a;id4b];
REWRITE_TAC[COMPLEX_MUL_RID];
]
(* }}} *));;


let taylor_coeff_catn = new_definition `taylor_coeff_catn n (z:complex) = 
  if (n=0) then catn z else Cx (& (FACT (n-1))) *
   (inv(Cx (&2))) * ( ( (-- ii) pow (n - 1) * ((inv (Cx (&1) + ii * z)) pow n)) +
      ( ii pow (n - 1) * ((inv (Cx (&1) - ii * z)) pow n)))`;;


let taylor_coeff_catn0 = 
prove_by_refinement (
  `taylor_coeff_catn 0 = catn `,

(* {{{ proof *)
 [
 ONCE_REWRITE_TAC[FUN_EQ_THM];
 REWRITE_TAC[taylor_coeff_catn];
 ]
(* }}} *));;


let taylor_coeff_catn1 = 
prove_by_refinement (
   `!z.  (Re z = &0 ==> abs(Im z) < &1) ==> 
       (catn has_complex_derivative (taylor_coeff_catn 1 z)) (at z)`,

(* {{{ proof *)
  [
  REPEAT STRIP_TAC;
    SUBGOAL_THEN `taylor_coeff_catn 1 z = inv (Cx (&1) + z pow 2)` ASSUME_TAC;
  REWRITE_TAC[taylor_coeff_catn;ARITH_RULE `~(1=0) /\ (1-1 = 0) /\ (FACT 0 =1)`;COMPLEX_POW_1;complex_pow;COMPLEX_MUL_LID];
  ASM_SIMP_TAC[id5];
  (* *)
  ASM_SIMP_TAC[HAS_COMPLEX_DERIVATIVE_CATN;];
  ]
(* }}} *));;


let taylor_coeff_catn_pos = 
prove_by_refinement(
  `!n. (n > 0) ==> (taylor_coeff_catn n  = (\z.
		Cx (& (FACT (n-1))) *
   (inv(Cx (&2))) * ( ( (-- ii) pow (n - 1) * ((inv (Cx (&1) + ii * z)) pow n)) +
      ( ii pow (n - 1) * ((inv (Cx (&1) - ii * z)) pow n)))			  ))`,

(* {{{ proof *)
  [
  REPEAT STRIP_TAC;
  ONCE_REWRITE_TAC[FUN_EQ_THM];
  REWRITE_TAC[taylor_coeff_catn];
  ASM_SIMP_TAC[ARITH_RULE `n > 0 ==> ~(n=0)`];
  ]
(* }}} *));;


let taylor_series_inv_pow = 
prove_by_refinement(
  `!n a z. ~(Cx (&1) + a * z = Cx(&0)) ==>
   (((\z. (inv (Cx (&1) + a * z)) pow n) has_complex_derivative 
      (-- Cx(&n) * a * (inv (Cx(&1) + a * z)) pow (n+1))) (at z))`,

  (* {{{ proof *)

  [
  REPEAT STRIP_TAC;
  DISJ_CASES_TAC (ARITH_RULE `(n=0) \/ (n > 0)`);
  ASM_REWRITE_TAC[ARITH_RULE `0+1=1`;complex_pow;COMPLEX_POW_1;SIMPLE_COMPLEX_ARITH `-- Cx (&0) * u = Cx(&0)`;HAS_COMPLEX_DERIVATIVE_CONST];
  ASM_SIMP_TAC[ARITH_RULE `(n>0) ==> (n+1 = 2 + (n-1))`;];
  REWRITE_TAC[COMPLEX_POW_ADD];
  ONCE_REWRITE_TAC[SIMPLE_COMPLEX_ARITH `-- r * s * t * u = r * u * ( -- s * t)`];
  MATCH_MP_TAC HAS_COMPLEX_DERIVATIVE_POW_AT;
  REWRITE_TAC[COMPLEX_POW_INV;GSYM complex_div];
  MATCH_MP_TAC HAS_COMPLEX_DERIVATIVE_INV_AT;
  ASM_REWRITE_TAC[];
  CONV_TAC (PATH_CONV "lr" (ONCE_REWRITE_CONV[SIMPLE_COMPLEX_ARITH `a = Cx(&0) + a * Cx (&1)` ]));
 MATCH_MP_TAC HAS_COMPLEX_DERIVATIVE_ADD;
  REWRITE_TAC[HAS_COMPLEX_DERIVATIVE_CONST];
  MATCH_MP_TAC HAS_COMPLEX_DERIVATIVE_LMUL_AT;
  REWRITE_TAC[HAS_COMPLEX_DERIVATIVE_ID];
  ]);;

  (* }}} *)


let factorial_lemma = 
prove_by_refinement(
  `!n. (n>0) ==> (FACT n = n * FACT (n-1))`,

  (* {{{ proof *)
  [
REPEAT STRIP_TAC;
DISJ_CASES_TAC (ISPEC `n:num` num_CASES);
ASM_MESON_TAC[ARITH_RULE `(n>0) ==> ~(n=0)`];
ASM_MESON_TAC[FACT;ARITH_RULE `SUC n - 1 = n`];
  ]);;
  (* }}} *)


let taylor_coeff_catn_deriv = 
prove_by_refinement(
  `!z n. (Re z = &0 ==> abs(Im z) < &1) ==> 
       ((taylor_coeff_catn n) has_complex_derivative 
              (taylor_coeff_catn (n+1) z)) (at z)`,

(* {{{ proof *)
[
REPEAT STRIP_TAC;
DISJ_CASES_TAC (ARITH_RULE `(n = 0) \/ (n >0)`);
ASM_SIMP_TAC[ taylor_coeff_catn0;taylor_coeff_catn1;ARITH_RULE `0+1=1`];
REWRITE_TAC[taylor_coeff_catn;ARITH_RULE `(n+1)-1 = n`];
ASM_SIMP_TAC[ARITH_RULE `(n>0) ==> ~(n+1=0)`];
ASM_SIMP_TAC[taylor_coeff_catn_pos;factorial_lemma];
(* fact finding *)
SUBGOAL_THEN `!u. Cx (&(n * FACT (n-1))) * u = Cx (&n) * (Cx (&(FACT (n-1))) * u)` MP_TAC;
REWRITE_TAC[CX_MUL;GSYM REAL_OF_NUM_MUL];
SIMPLE_COMPLEX_ARITH_TAC;
DISCH_THEN (fun t-> REWRITE_TAC[t]);
ONCE_REWRITE_TAC[SIMPLE_COMPLEX_ARITH `(a:complex) * b * c * d = b * c * a * d`]; 
ONCE_REWRITE_TAC[SIMPLE_COMPLEX_ARITH `(a:complex) * b * u = (a*b)*u`]; 
MATCH_MP_TAC HAS_COMPLEX_DERIVATIVE_LMUL_AT;
REWRITE_TAC[COMPLEX_ADD_LDISTRIB];
MATCH_MP_TAC HAS_COMPLEX_DERIVATIVE_ADD;
REWRITE_TAC[SIMPLE_COMPLEX_ARITH `a - ii * z = a + (-- ii) * z`];
SUBGOAL_THEN `!a b r. a * b pow n * r = b pow (n-1) * ( -- a * ( -- b) * r)` MP_TAC;
REPEAT STRIP_TAC;
MP_TAC(ARITH_RULE `n>0 ==> n = (n-1) + 1`);
ASM_REWRITE_TAC[];
DISCH_THEN   (fun t -> CONV_TAC (PATH_CONV "lr" (ONCE_REWRITE_CONV[t])));
REWRITE_TAC[COMPLEX_POW_ADD;COMPLEX_POW_1];
SIMPLE_COMPLEX_ARITH_TAC;
DISCH_THEN (fun t->REWRITE_TAC[t]);
CONJ_TAC THEN (MATCH_MP_TAC HAS_COMPLEX_DERIVATIVE_LMUL_AT) ;
(* highly parallel branches *)
REWRITE_TAC[SIMPLE_COMPLEX_ARITH `-- --ii = ii`];
MATCH_MP_TAC  taylor_series_inv_pow;
STRIP_TAC;
ASM_MESON_TAC[idz];
(* *)
MATCH_MP_TAC  taylor_series_inv_pow;
STRIP_TAC;
ASM_MESON_TAC[idz];
]
(* }}} *));;


let ipows2 = 
prove_by_refinement(
  `!n. (--ii ) pow (n + 2) = -- ((-- ii) pow n)`,

  (* {{{ proof *)
  [
  GEN_TAC;
  REWRITE_TAC[COMPLEX_POW_ADD;COMPLEX_POW_2;];
  MATCH_MP_TAC (SIMPLE_COMPLEX_ARITH `(a * a = -- Cx(&1) ) ==> r * a * a = -- r`);
  REWRITE_TAC[ii];
  SIMPLE_COMPLEX_ARITH_TAC;
  ]);;
  (* }}} *)


let ipowsc2 = 
prove_by_refinement(
  `!n. (ii ) pow (n + 2) = -- ((ii) pow n)`,

  (* {{{ proof *)
  [
  GEN_TAC;
  REWRITE_TAC[COMPLEX_POW_ADD;COMPLEX_POW_2;];
  MATCH_MP_TAC (SIMPLE_COMPLEX_ARITH `(a * a = -- Cx(&1) ) ==> r * a * a = -- r`);
  REWRITE_TAC[ii];
  SIMPLE_COMPLEX_ARITH_TAC;
  ]);;
  (* }}} *)



let taylor_coeff0 = 
prove_by_refinement(
  `!n. (taylor_coeff_catn n (Cx (&0)) = if (EVEN n) then (Cx (&0)) else
    Cx (&(FACT (n-1))  * (-- &1) pow ((n - 1) DIV 2)))`,

  (* {{{ proof *)
  [
  GEN_TAC;
  DISJ_CASES_TAC (ISPEC `n:num` EVEN_OR_ODD);
   DISJ_CASES_TAC (ARITH_RULE `(n=0) \/ (n >0)`);
  ASM_REWRITE_TAC[taylor_coeff_catn;GSYM CX_ATN;ATN_0];
  (* *)
  ASM_SIMP_TAC[taylor_coeff_catn_pos];
  REWRITE_TAC[SIMPLE_COMPLEX_ARITH `a + b * Cx(&0) = a /\ a - b * Cx(&0)  =a`];
  REWRITE_TAC[SIMPLE_COMPLEX_ARITH `inv(Cx (&1)) = Cx (&1) /\ a * Cx(&1) = a`;COMPLEX_POW_ONE];
  MATCH_MP_TAC (SIMPLE_COMPLEX_ARITH `c = Cx(&0) ==> a*b*c = Cx(&0)`);
  SUBGOAL_THEN (  `EVEN n  ==> (?k. (n - 1) = 2 * k + 1)`) MP_TAC;
  REWRITE_TAC[EVEN_EXISTS];
   REPEAT STRIP_TAC;
  EXISTS_TAC `m-1`;
  REPEAT (POP_ASSUM MP_TAC);
  ARITH_TAC;
  ASM_REWRITE_TAC[];
  REPEAT STRIP_TAC;
  ASM_REWRITE_TAC[];
  SPEC_TAC (`k:num`,`k:num`);
  INDUCT_TAC;
  REWRITE_TAC[ARITH_RULE `2* 0 + 1 = 1`;COMPLEX_POW_1];
  SIMPLE_COMPLEX_ARITH_TAC;
  ASM_REWRITE_TAC[ipows2;ipowsc2;ARITH_RULE `2* SUC k' + 1 = (2 * k' + 1) + 2`;SIMPLE_COMPLEX_ARITH `--a  + --b = --(a+b) /\ -- Cx(&0) = Cx(&0) `];
  (* ODD *)
  ASM_REWRITE_TAC[GSYM NOT_ODD];
  SUBGOAL_THEN (`ODD n ==> (?k. n = 2 * k + 1)`) MP_TAC;
  REWRITE_TAC[ODD_EXISTS];
  REPEAT STRIP_TAC;
  EXISTS_TAC `m:num`;
  ASM_REWRITE_TAC[];
  ARITH_TAC;
  ASM_REWRITE_TAC[];
   REPEAT STRIP_TAC;
   SUBGOAL_THEN `n > 0` ASSUME_TAC;
  POP_ASSUM MP_TAC;
  ARITH_TAC;
  ASM_SIMP_TAC[taylor_coeff_catn_pos];
  REWRITE_TAC[SIMPLE_COMPLEX_ARITH `a + b * Cx (&0) = a /\ a - b * Cx (&0) = a /\ inv (Cx (&1)) = Cx (&1) /\ a * Cx (&1) = a`;COMPLEX_POW_ONE;ARITH_RULE `(2 * k + 1 ) - 1 = 2 * k`;CX_MUL];
  MATCH_MP_TAC (SIMPLE_COMPLEX_ARITH `(b = c) ==> (a*b = a*c)`);
  SPEC_TAC (`k:num`,`k:num`);
  INDUCT_TAC;
  REWRITE_TAC[ARITH_RULE `2 * 0 =0 /\ 0 DIV 2 = 0`;complex_pow;real_pow];
  SIMPLE_COMPLEX_ARITH_TAC;
  REWRITE_TAC[ARITH_RULE `2 * SUC k' = 2 * k' + 2 /\ (2 * SUC k) DIV 2 = ((2 * k) DIV 2) + 1`;SIMPLE_COMPLEX_ARITH `a * (-- b + -- c) = -- (a * (b+c))`;ipows2;ipowsc2;COMPLEX_POW_ADD;REAL_POW_ADD;CX_MUL;REAL_POW_1];
  ASM_REWRITE_TAC[];
  SIMPLE_COMPLEX_ARITH_TAC;
  ]);;
  (* }}} *)


let term_bound = 
prove_by_refinement(
  `!a n z.  Im (z) = &0 ==>
         norm((Cx(a)* ii) pow n * ((inv (Cx (&1) - Cx(a)* ii * z)) pow (n+1))) 
  <= (abs a) pow n `,

  (* {{{ proof *)
  [
  REPEAT STRIP_TAC;
  REWRITE_TAC[COMPLEX_NORM_MUL;COMPLEX_NORM_II;COMPLEX_NORM_CX;COMPLEX_NORM_POW;COMPLEX_NORM_INV;REAL_ARITH `a * &1 = a`];
  MATCH_MP_TAC (MESON[REAL_LE_LMUL;REAL_ARITH `x = x* &1`]  ( `!x y. &0 <= x /\ y <= &1 ==> x *y <= x`)) ;
  CONJ_TAC;
  MATCH_MP_TAC REAL_POW_LE;
  REWRITE_TAC[REAL_ABS_POS];
  MATCH_MP_TAC (MESON[REAL_POW_LE2;REAL_POW_ONE] `!x. &0 <= x /\ x <= &1 ==> x pow n <= &1`);
  REWRITE_TAC[REAL_LE_INV_EQ;NORM_POS_LE];
  MATCH_MP_TAC (MESON[REAL_LE_INV2;REAL_INV_1;REAL_ARITH `&0 < &1`] `&1 <= x ==> inv (x) <= &1`);
  SUBGOAL_THEN `Im z = &0 ==> Cx(&1) - Cx a * ii * z = complex(&1, -- a * Re(z))` MP_TAC;
  REWRITE_TAC[ii];
  SIMPLE_COMPLEX_ARITH_TAC;
  ASM_REWRITE_TAC[];
  DISCH_THEN (fun t -> REWRITE_TAC[t]);
  (* *)
  MATCH_MP_TAC (MESON[Collect_geom.POW2_COND;REAL_ARITH `&0 <= &1 /\ &1 pow 2 = &1`;NORM_POS_LE] `&1 <= norm x pow 2 ==> &1 <= norm x`);
  REWRITE_TAC[COMPLEX_SQNORM;RE;IM];
  MESON_TAC[Collect_geom.pow_g;REAL_ARITH `&1 pow 2 = &1 /\ (&0 <= t ==> &1 <= &1 + t)`];
  ]);;
  (* }}} *)


let taylor_error_bound = 
prove_by_refinement(
  `!n z. Im(z) = &0 ==> norm(taylor_coeff_catn (n+1) z) <= &(FACT n)`,

  (* {{{ proof *)

  [
  REPEAT STRIP_TAC;
  SIMP_TAC[taylor_coeff_catn_pos;ARITH_RULE `(n+1)>0 /\ ((n+1)-1 = n)`];
  REWRITE_TAC[COMPLEX_NORM_MUL;COMPLEX_NORM_NUM;COMPLEX_NORM_INV];
  REWRITE_TAC[COMPLEX_NORM_NUM];
 MATCH_MP_TAC (MESON[REAL_LE_LMUL;REAL_ARITH `x = x* &1`]  ( `!x y. &0 <= x /\ y <= &1 ==> x *y <= x`)) ;
  CONJ_TAC;
  REWRITE_TAC[REAL_OF_NUM_LE];
  ARITH_TAC;
MATCH_MP_TAC (REAL_ARITH `x <= &1 + &1 ==> inv (&2)* x <= &1`);
MATCH_MP_TAC REAL_LE_TRANS;
EXISTS_TAC `norm (--ii pow n * inv (Cx (&1) + ii * z) pow (n + 1)) + norm ( ii pow n * inv (Cx (&1) - ii * z) pow (n + 1))`;
REWRITE_TAC[NORM_TRIANGLE];
MATCH_MP_TAC REAL_LE_ADD2;
(*  *)
CONJ_TAC;
(* FORCE_MATCH_MP_TAC *)
FORCE_MATCH_MP_TAC (ISPECL [`-- &1`;`n:num`;`z:complex`] term_bound);
ASM_REWRITE_TAC[];
ONCE_REWRITE_TAC[FUN_EQ_THM];
SIMPLE_COMPLEX_ARITH_TAC;
SIMPLE_COMPLEX_ARITH_TAC;
REWRITE_TAC[REAL_ABS_NEG;REAL_ABS_NUM;REAL_POW_ONE];
(* second *)
ONCE_REWRITE_TAC[SIMPLE_COMPLEX_ARITH `ii * z = Cx(&1) * ii * z`];
FORCE_MATCH_MP_TAC (ISPECL [` &1`;`n:num`;`z:complex`] term_bound);
ASM_REWRITE_TAC[];
SIMPLE_COMPLEX_ARITH_TAC;
REWRITE_TAC[REAL_ABS_NEG;REAL_ABS_NUM;REAL_POW_ONE];
  ]);;

  (* }}} *)


let complex_taylor_catn = 
prove_by_refinement(
   ` !n s. (s = {z | Im (z) = &0 }) ==> (!z.
                  (Cx(&0)) IN s /\ z IN s
                  ==> norm
                      (catn z -
                       vsum (0..n)
                       (\i. taylor_coeff_catn i (Cx(&0)) * (z) pow i / Cx (&(FACT i)))) <=
                       norm (z) pow (n + 1) )`,

  (* {{{ proof *)

  [
  GEN_TAC THEN GEN_TAC;
  DISCH_TAC;
  MP_TAC (SPECL[`taylor_coeff_catn`;`n:num`;`s:complex->bool`;`&(FACT n)`] COMPLEX_TAYLOR);
  ASM_REWRITE_TAC[IN_ELIM_THM];
  ANTS_TAC;
  REPEAT (CONJ_TAC THEN (REPEAT STRIP_TAC));
  REWRITE_TAC[convex;IN_ELIM_THM;IM_ADD;IM_CMUL];
  REPEAT STRIP_TAC THEN ASM_REWRITE_TAC[];
  REAL_ARITH_TAC;
  MATCH_MP_TAC HAS_COMPLEX_DERIVATIVE_AT_WITHIN;
  MATCH_MP_TAC taylor_coeff_catn_deriv;
  ASM_REWRITE_TAC[] THEN REAL_ARITH_TAC;
  (*  *)
  MATCH_MP_TAC taylor_error_bound;
  ASM_REWRITE_TAC[];
  DISCH_THEN (MP_TAC o (ISPEC `Cx (&0)`));
  REWRITE_TAC[taylor_coeff_catn0;SIMPLE_COMPLEX_ARITH `z - Cx (&0) = z /\ Im (Cx (&0)) = &0`];
  DISCH_THEN (fun t-> REPEAT STRIP_TAC THEN MP_TAC t) ;
  DISCH_THEN (MP_TAC o (ISPEC `z:complex`));
  DISCH_THEN FORCE_MATCH_MP_TAC;
  ASM_REWRITE_TAC[];
  MATCH_MP_TAC REAL_DIV_LMUL;
  REWRITE_TAC[REAL_OF_NUM_EQ];
  MESON_TAC[ FACT_LT;ARITH_RULE `0 < x ==> ~(x =0)`];
  ]);;

  (* }}} *)


let real_axis = prove_by_refinement( (* not needed *)
  `!z.  (Im z = &0 <=> (?x. z = Cx(x)))`,
  (* {{{ proof *)
  [
  GEN_TAC;
  EQ_TAC;
  DISCH_TAC;
  EXISTS_TAC `(Re z)`;
  POP_ASSUM MP_TAC;
  SIMPLE_COMPLEX_ARITH_TAC;
  REPEAT STRIP_TAC THEN ASM_REWRITE_TAC[IM_CX];
  ]);;
  (* }}} *)


let THREAD_IF = 
prove_by_refinement(
  `!x y z f. (f:A->B) (if x then y else z) = if x then f y else f z`,

  (* {{{ proof *)
  [
  REPEAT GEN_TAC;
  BOOL_CASES_TAC `x:bool` THEN REWRITE_TAC[];
  ]);;
  (* }}} *)



let real_taylor_atn_ver1 = 
prove_by_refinement(
  `!n x.  abs(atn x - sum (0..n) (\i. if EVEN i then &0 else ( -- &1 pow ((i-1) DIV 2)  *      x pow i / &i ))) <= abs(x) pow (n+1)`,

  (* {{{ proof *)
  [
  REPEAT GEN_TAC;
  MP_TAC (ISPECL [`n:num`;`{ z | Im(z) = &0 }`] complex_taylor_catn);
  REWRITE_TAC[];
  DISCH_THEN (fun t -> MP_TAC (ISPEC `Cx (x)` t));
  REWRITE_TAC[VSUM_CX_NUMSEG;COMPLEX_NORM_CX;GSYM CX_SUB;GSYM CX_POW;IN_ELIM_THM;GSYM CX_MUL;GSYM CX_DIV;GSYM CX_ATN;IM_CX;taylor_coeff0;GSYM THREAD_IF];
  FORCE_MATCH;
 ONCE_REWRITE_TAC[FUN_EQ_THM];
  X_GEN_TAC `i:num`;
  BETA_TAC;
  DISJ_CASES_TAC (TAUT `EVEN i \/ ~(EVEN i)`) THEN ASM_REWRITE_TAC[];
  REAL_ARITH_TAC;
  POP_ASSUM MP_TAC;
  REWRITE_TAC[NOT_EVEN;ODD_EXISTS];
  REPEAT STRIP_TAC;
  ASM_REWRITE_TAC[ARITH_RULE `SUC x - 1 = x`;FACT;GSYM REAL_OF_NUM_MUL ];
  ONCE_REWRITE_TAC[REAL_FIELD `((a:real) * b) * c/(d*a) = (b * c/d) * (a/a)`];
  MATCH_MP_TAC (REAL_FIELD `x = &1 ==> (y * x = y)`);
  MATCH_MP_TAC REAL_DIV_REFL;
  REWRITE_TAC[REAL_OF_NUM_EQ];
  MESON_TAC[FACT_LT;ARITH_RULE `0 < x ==> ~(x = 0)`];
  ]);;
  (* }}} *)


let sum_odd  = 
prove_by_refinement(
  `!(g:num->real) n. sum { i | ODD i /\ i IN 0.. 2 * n + 2 } g =   sum (0.. n) (\i. g (2 * i +1))`,

  (* {{{ proof *)
  [
  REPEAT STRIP_TAC;
  FORCE_MATCH_MP_TAC (ISPECL [`\i. (2 * i + 1)`; `g:num->real`;`(0..n)`] SUM_IMAGE);
  BETA_TAC;
  ARITH_TAC;
  REWRITE_TAC[IMAGE;numseg;IN_ELIM_THM;ODD_EXISTS];
  ONCE_REWRITE_TAC[FUN_EQ_THM];
  REWRITE_TAC[IN_ELIM_THM];
  GEN_TAC;
  EQ_TAC;
  REPEAT STRIP_TAC;
  EXISTS_TAC `x':num`;
  POP_ASSUM MP_TAC;
  ARITH_TAC;
  REPEAT (POP_ASSUM MP_TAC);
  ARITH_TAC;
  REPEAT (POP_ASSUM MP_TAC);
  ARITH_TAC;
  REPEAT STRIP_TAC;
  EXISTS_TAC `m:num`;
  REPEAT (POP_ASSUM MP_TAC);
  ARITH_TAC;
  ONCE_REWRITE_TAC[FUN_EQ_THM;];
  REWRITE_TAC[o_THM]]
  (* }}} *)
);;


let sum_even = 
prove_by_refinement(
  `!g n. sum {i | ODD i /\ i IN 0..n } (\i. if EVEN i then &0 else g i) = sum (0..n) (\i. if EVEN i then &0 else g i)`,

  (* {{{ proof *)
  [
  REPEAT GEN_TAC;
  ONCE_REWRITE_TAC [MESON[] `x = y <=> y = x`];
  FORCE_MATCH_MP_TAC (ISPECL [`(\i. if EVEN i then &0 else (g i))`;`{i | ODD i /\ i IN 0..n}`;`(0..n)`] SUM_SUPERSET);
  REPEAT STRIP_TAC;
  SET_TAC[];
  REPEAT (POP_ASSUM MP_TAC);
  REWRITE_TAC[IN_ELIM_THM;numseg];
  MESON_TAC[NOT_ODD];
  REWRITE_TAC[];
  ]);;
  (* }}} *)


let real_taylor_atn = 
prove_by_refinement(
  `!n x. abs(atn x - sum (0..n) (\j. (-- &1 pow j) * x pow (2 * j + 1)/ &(2 * j+ 1))) <= abs(x) pow (2 * n + 3)`,

  (* {{{ proof *)
  [
  REPEAT STRIP_TAC;
  MP_TAC (ISPECL [`2*n + 2`;`x:real`] real_taylor_atn_ver1);
  REWRITE_TAC[ARITH_RULE `(2*n + 2) +1 = 2 *n + 3`];
  REWRITE_TAC[GSYM sum_even];
  REWRITE_TAC[sum_odd];
  SUBGOAL_THEN `!i. ~(EVEN (2 *i + 1))` MP_TAC;
  MP_TAC NOT_EVEN;
  REWRITE_TAC[EVEN_EXISTS;ODD_EXISTS;ARITH_RULE `!m. SUC (m) = m+1`];  
   MESON_TAC[];
  DISCH_TAC;
  ASM_REWRITE_TAC[];
  REWRITE_TAC[ARITH_RULE `((2 * i + 1) - 1) DIV 2 = i`];
  ]);;
  (* }}} *)


let halfatn4 = new_definition `halfatn4 = halfatn o halfatn o halfatn o halfatn`;;


let real_taylor_atn_halfatn4 = 
prove_by_refinement(
  `!n x. abs (atn(halfatn4 x) - sum (0..n) (\j. (-- &1 pow j) * halfatn4 x pow (2 * j + 1)/ &(2 * j+ 1))) <= inv (&8 pow (2 * n + 3))`,

  (* {{{ proof *)
  [
  REPEAT STRIP_TAC;
  MATCH_MP_TAC REAL_LE_TRANS;
  EXISTS_TAC `abs(halfatn4 x) pow (2 * n + 3)`;
  REWRITE_TAC[real_taylor_atn;REAL_INV_POW];
  MATCH_MP_TAC Real_ext.REAL_PROP_LE_POW;
  REWRITE_TAC[REAL_ABS_POS;halfatn4;o_THM];
  MATCH_MP_TAC (REAL_ARITH `x < y ==> x <= y`);
  MP_TAC (ISPEC `x:real` halfatn_bounds_abs);
  REPLICATE_TAC 3(  DISCH_THEN (fun t-> MP_TAC (MATCH_MP halfatn_half t)));
  FORCE_MATCH;
  CONV_TAC REAL_FIELD;
  ]);;
  (* }}} *)


let atn_halfatn4 = 
prove_by_refinement(
  `!x. atn x = &16 * atn(halfatn4 x)`,

  (* {{{ proof *)
  [
  ONCE_REWRITE_TAC[REAL_ARITH `&16 * x = &2 * &2 * &2 * &2 * x`];
  REWRITE_TAC[halfatn4;o_THM;GSYM atn_half];
  ]);;
  (* }}} *)


let real_taylor_atn_halfatn4_a = 
prove_by_refinement(
  `!n x. abs (atn x - &16 * sum (0..n) (\j. (-- &1 pow j) * halfatn4 x pow (2 * j + 1)/ &(2 * j+ 1))) <= inv (&2 pow (6 * n + 5 ))`,

  (* {{{ proof *)
  [
  REPEAT GEN_TAC;
  ONCE_REWRITE_TAC [atn_halfatn4];
  REWRITE_TAC[REAL_ARITH `abs (&16 * x - &16 * y) <=  z <=> abs(x - y) <= z *inv (&2 pow 4)`];
  MP_TAC (ISPECL [`n:num`;`x:real`]  real_taylor_atn_halfatn4);
  FORCE_MATCH;
  REWRITE_TAC[GSYM REAL_INV_MUL;GSYM REAL_POW_ADD;REAL_ARITH `&8 = &2 pow 3`;REAL_POW_POW];
  REPEAT AP_TERM_TAC;
  ARITH_TAC;
  ]);;
  (* }}} *)


let halfatn4_co = new_definition `halfatn4_co x j =  (-- &1 pow j) * halfatn4 x pow (2 * j + 1)/ &(2 * j+ 1)`;;


let atn_bounds_anti = 
prove_by_refinement(
  `!x y. x <= y  ==> abs(atn x - atn y) <= abs(x - y)`,

  (* {{{ proof *)
  [
  REPEAT STRIP_TAC;
  MP_TAC (ISPECL [`atn`;`(\t. inv(&1 + t pow 2))`;`x:real`;`y:real`] REAL_MVT_VERY_SIMPLE);
  REWRITE_TAC[real_interval;IN_ELIM_THM];
  ANTS_TAC;
  ASM_REWRITE_TAC[];
  REPEAT STRIP_TAC;
  MATCH_MP_TAC HAS_REAL_DERIVATIVE_ATREAL_WITHIN;
  REWRITE_TAC[HAS_REAL_DERIVATIVE_ATN];
  REPEAT STRIP_TAC;
  ONCE_REWRITE_TAC[REAL_ARITH `abs(y - x) = abs(x-y)`];
  ASM_REWRITE_TAC[];
  REWRITE_TAC[REAL_ABS_MUL];
  FORCE_MATCH_MP_TAC (ISPECL [`abs(inv (&1 + x' pow 2))`;`&1`;`abs(y-x)`] REAL_LE_RMUL);
  REWRITE_TAC[REAL_ABS_POS;REAL_ABS_INV];
 FORCE_MATCH_MP_TAC (ISPECL [`&1`;`abs(&1 + x' pow 2)`] REAL_LE_INV2);
  CONJ_TAC THEN TRY(REAL_ARITH_TAC);
  MP_TAC (ISPEC `x':real` pos1);
  SIMP_TAC [REAL_ARITH `(&0 < x ==> (abs x = x)) /\  (&1 <= &1 + u <=> &0 <= u)`;Real_ext.REAL_LE_POW_2];
 REAL_ARITH_TAC;
 REAL_ARITH_TAC;
  ]);;
  (* }}} *)

prioritize_real();;


let atn_bounds = 
prove_by_refinement(
  `!x y. abs(atn x - atn y) <= abs(x-y)`,

  (* {{{ proof *)
  [
  REPEAT STRIP_TAC;
  DISJ_CASES_TAC (REAL_ARITH `x<= y \/ y <= x`);
  ASM_SIMP_TAC[atn_bounds_anti];
  ONCE_REWRITE_TAC[REAL_ARITH `abs(x -y) = abs(y-x)`];
  ASM_SIMP_TAC[atn_bounds_anti];
  ]);;
  (* }}} *)


let real_taylor_atn_approx = 
prove_by_refinement(
  `!n x  v eps1 eps2 eps.
    inv (&2 pow (6 * n + 5)) <= eps1 /\
    abs(&16 *sum (0..n) (halfatn4_co x) - v )<= eps2 /\
   (eps1 + eps2 <= eps) ==>
   abs(atn(x) - v) <= eps`,

  (* {{{ proof *)
  [
  REPEAT STRIP_TAC;
    MATCH_MP_TAC REAL_LE_TRANS;
  EXISTS_TAC `( eps1:real) + eps2`;
  ASM_REWRITE_TAC[];
  MATCH_MP_TAC REAL_LE_TRANS;
  ABBREV_TAC `r = &16 * sum(0..n) (halfatn4_co x)`;
  EXISTS_TAC`abs(atn x - r) + abs(r - v)`;
  CONJ_TAC;
  MESON_TAC[REAL_ABS_TRIANGLE;REAL_ARITH `atn x - v = (atn x - r) + r - v`];
 MATCH_MP_TAC (REAL_ARITH `a1 <= b1 /\ a2 <= b2  ==> a1 + a2  <= b1 + b2`);
 ASM_REWRITE_TAC[];
  MATCH_MP_TAC REAL_LE_TRANS;
  EXISTS_TAC `inv(&2 pow (6 * n + 5))`;
  ASM_REWRITE_TAC[];
  EXPAND_TAC "r";
  MP_TAC (ISPECL [`n:num`;`x:real`] real_taylor_atn_halfatn4_a);
  FORCE_MATCH;
  REWRITE_TAC[FUN_EQ_THM;halfatn4_co];
  ]);;
  (* }}} *)


end;;