(* ========================================================================= *)
(* The type "real^2" regarded as the complex numbers.                        *)
(*                                                                           *)
(*              (c) Copyright, John Harrison 1998-2008                       *)
(*              (c) Copyright, Valentina Bruno 2010                          *)
(* ========================================================================= *)
needs "Multivariate/integration.ml";;
new_type_abbrev("complex",`:real^2`);;
let prioritize_complex() =
  overload_interface("--",`vector_neg:complex->complex`);
  overload_interface("+",`vector_add:complex->complex->complex`);
  overload_interface("-",`vector_sub:complex->complex->complex`);
  overload_interface("*",`complex_mul:complex->complex->complex`);
  overload_interface("/",`complex_div:complex->complex->complex`);
  overload_interface("pow",`complex_pow:complex->num->complex`);
  overload_interface("inv",`complex_inv:complex->complex`);;
prioritize_complex();;
(* ------------------------------------------------------------------------- *)
(* Real and imaginary parts of a number.                                     *)
(* ------------------------------------------------------------------------- *)
let RE_DEF = new_definition
  `Re(z:complex) = z$1`;;
let IM_DEF = new_definition
  `Im(z:complex) = z$2`;;
let CX_DEF = new_definition
 `Cx(a) = complex(a,&0)`;;
let complex_mul = new_definition
  `w * z = complex(Re(w) * Re(z) - Im(w) * Im(z),
                   Re(w) * Im(z) + Im(w) * Re(z))`;;let complex_inv = new_definition
  `inv(z) = complex(Re(z) / (Re(z) pow 2 + Im(z) pow 2),
                    --(Im(z)) / (Re(z) pow 2 + Im(z) pow 2))`;;let complex_pow = define
  `(x pow 0 = Cx(&1)) /\
   (!n. x pow (SUC n) = x * x pow n)`;;
let COMPLEX_EQ_0 = prove
 (`z = Cx(&0) <=> Re(z) pow 2 + Im(z) pow 2 = &0`,
  REWRITE_TAC[
COMPLEX_EQ; 
CX_DEF; 
RE; 
IM] THEN
  EQ_TAC THEN SIMP_TAC[] THEN CONV_TAC REAL_RAT_REDUCE_CONV THEN
  DISCH_THEN(MP_TAC o MATCH_MP (REAL_ARITH
   `!x y:real. x + y = &0 ==> &0 <= x /\ &0 <= y ==> x = &0 /\ y = &0`)) THEN
  REWRITE_TAC[REAL_POW_2; 
REAL_LE_SQUARE; 
REAL_ENTIRE]);;
 
let COMPLEX_ADD_AC = prove
 (`(m + n = n + m) /\ ((m + n) + p = m + n + p) /\ (m + n + p = n + m + p)`,
  SIMPLE_COMPLEX_ARITH_TAC);;
 
let COMPLEX_MUL_AC = prove
 (`(m * n = n * m) /\ ((m * n) * p = m * n * p) /\ (m * n * p = n * m * p)`,
  SIMPLE_COMPLEX_ARITH_TAC);;
 
let II_NZ = prove
 (`~(ii = Cx(&0))`,
  REWRITE_TAC[ii] THEN SIMPLE_COMPLEX_ARITH_TAC);;
 
let CX_INV = prove
 (`!x. Cx(inv x) = inv(Cx x)`,
  GEN_TAC THEN REWRITE_TAC[
CX_DEF; 
complex_inv; 
RE; 
IM; 
COMPLEX_EQ] THEN
  ASM_CASES_TAC `x = &0` THEN ASM_REWRITE_TAC[] THEN
  CONV_TAC REAL_RAT_REDUCE_CONV THEN
  POP_ASSUM MP_TAC THEN CONV_TAC REAL_FIELD);;
 
let COMPLEX_EXPAND = prove
 (`!z. z = Cx(Re z) + ii * Cx(Im z)`,
  REWRITE_TAC[ii] THEN SIMPLE_COMPLEX_ARITH_TAC);;
 
let COMPLEX_TRAD = prove
 (`!x y. complex(x,y) = Cx(x) + ii * Cx(y)`,
  REWRITE_TAC[ii] THEN SIMPLE_COMPLEX_ARITH_TAC);;
 
let RE_II = prove
 (`Re ii = &0`,
  REWRITE_TAC[ii] THEN SIMPLE_COMPLEX_ARITH_TAC);;
 
let IM_II = prove
 (`Im ii = &1`,
  REWRITE_TAC[ii] THEN SIMPLE_COMPLEX_ARITH_TAC);;
 
let RE_MUL_II = prove
 (`!z. Re(z * ii) = --(Im z) /\ Re(ii * z) = --(Im z)`,
  REWRITE_TAC[ii] THEN SIMPLE_COMPLEX_ARITH_TAC);;
 
let IM_MUL_II = prove
 (`!z. Im(z * ii) = Re z /\ Im(ii * z) = Re z`,
  REWRITE_TAC[ii] THEN SIMPLE_COMPLEX_ARITH_TAC);;
 
let RE_MUL_CX = prove
 (`!x z. Re(Cx(x) * z) = x * Re z /\
         Re(z * Cx(x)) = Re z * x`,
  SIMPLE_COMPLEX_ARITH_TAC);;
 
let IM_MUL_CX = prove
 (`!x z. Im(Cx(x) * z) = x * Im z /\
         Im(z * Cx(x)) = Im z * x`,
  SIMPLE_COMPLEX_ARITH_TAC);;
 
let COMPLEX_POLY_CLAUSES = prove
 (`(!x y z. x + (y + z) = (x + y) + z) /\
   (!x y. x + y = y + x) /\
   (!x. Cx(&0) + x = x) /\
   (!x y z. x * (y * z) = (x * y) * z) /\
   (!x y. x * y = y * x) /\
   (!x. Cx(&1) * x = x) /\
   (!x. Cx(&0) * x = Cx(&0)) /\
   (!x y z. x * (y + z) = x * y + x * z) /\
   (!x. x pow 0 = Cx(&1)) /\
   (!x n. x pow (SUC n) = x * x pow n)`,
  REWRITE_TAC[
complex_pow] THEN SIMPLE_COMPLEX_ARITH_TAC)
and COMPLEX_POLY_NEG_CLAUSES = 
prove
 (`(!x. --x = Cx(-- &1) * x) /\
   (!x y. x - y = x + Cx(-- &1) * y)`,
  SIMPLE_COMPLEX_ARITH_TAC);;
 
let COMPLEX_INTEGRAL = prove
   (`(!x. Cx(&0) * x = Cx(&0)) /\
     (!x y z. (x + y = x + z) <=> (y = z)) /\
     (!w x y z. (w * y + x * z = w * z + x * y) <=> (w = x) \/ (y = z))`,
    REWRITE_TAC[
COMPLEX_ENTIRE; SIMPLE_COMPLEX_ARITH
     `(w * y + x * z = w * z + x * y) <=>
      (w - x) * (y - z) = Cx(&0)`] THEN
    SIMPLE_COMPLEX_ARITH_TAC)
  and COMPLEX_RABINOWITSCH = 
prove
   (`!x y:complex. ~(x = y) <=> ?z. (x - y) * z = Cx(&1)`,
    REPEAT GEN_TAC THEN
    GEN_REWRITE_TAC (LAND_CONV o RAND_CONV) [GSYM 
COMPLEX_SUB_0] THEN
    MESON_TAC[
COMPLEX_MUL_RINV; 
COMPLEX_MUL_LZERO;
              SIMPLE_COMPLEX_ARITH `~(Cx(&1) = Cx(&0))`])
  and COMPLEX_IIII = 
prove
   (`ii * ii + Cx(&1) = Cx(&0)`,
    REWRITE_TAC[ii; 
CX_DEF; 
complex_mul; 
complex_add; 
RE; 
IM] THEN
    AP_TERM_TAC THEN BINOP_TAC THEN REAL_ARITH_TAC) in
  let ring,ideal =
    RING_AND_IDEAL_CONV
        (dest_complex_const,mk_complex_const,COMPLEX_RAT_EQ_CONV,
         `(--):complex->complex`,`(+):complex->complex->complex`,
         `(-):complex->complex->complex`,`(inv):complex->complex`,
         `(*):complex->complex->complex`,`(/):complex->complex->complex`,
         `(pow):complex->num->complex`,
         COMPLEX_INTEGRAL,COMPLEX_RABINOWITSCH,COMPLEX_POLY_CONV)
  and ii_tm = `ii` and iiii_tm = concl COMPLEX_IIII in
  (fun tm -> if free_in ii_tm tm then
             MP (ring (mk_imp(iiii_tm,tm))) COMPLEX_IIII
             else ring tm),
  ideal;;
(* ------------------------------------------------------------------------- *)
(* Most basic properties of inverses.                                        *)
(* ------------------------------------------------------------------------- *)
 
let COMPLEX_INV_INV = prove
 (`!x:complex. inv(inv x) = x`,
  GEN_TAC THEN ASM_CASES_TAC `x = Cx(&0)` THEN
  ASM_REWRITE_TAC[
COMPLEX_INV_0] THEN
  POP_ASSUM MP_TAC THEN
  MAP_EVERY (fun t -> MP_TAC(SPEC t 
COMPLEX_MUL_RINV))
   [`x:complex`; `inv(x):complex`] THEN
  CONV_TAC COMPLEX_RING);;
 
let COMPLEX_INV_EQ_0 = prove
 (`!x. inv x = Cx(&0) <=> x = Cx(&0)`,
  GEN_TAC THEN ASM_CASES_TAC `x = Cx(&0)` THEN
  ASM_REWRITE_TAC[
COMPLEX_INV_0] THEN POP_ASSUM MP_TAC THEN
  CONV_TAC COMPLEX_FIELD);;
 
let COMPLEX_INV_EQ_1 = prove
 (`!x. inv x = Cx(&1) <=> x = Cx(&1)`,
  GEN_TAC THEN ASM_CASES_TAC `x = Cx(&0)` THEN
  ASM_REWRITE_TAC[
COMPLEX_INV_0] THEN POP_ASSUM MP_TAC THEN
  CONV_TAC COMPLEX_FIELD);;
 
let COMPLEX_DIV_POW = prove
 (`!x:complex n k:num.
      ~(x= Cx(&0)) /\ k <= n /\ ~(k = 0)
      ==> x pow (n - k) = x pow n / x pow k`,
 
let COMPLEX_DIV_POW2 = prove
 (`!z m n. ~(z = Cx(&0))
           ==> z pow m / z pow n =
               if n <= m then z pow (m - n) else inv(z pow (n - m))`,
  REPEAT STRIP_TAC THEN COND_CASES_TAC THEN
  ASM_SIMP_TAC[
COMPLEX_POW_EQ_0; COMPLEX_FIELD
    `~(b = Cx(&0)) /\ ~(c = Cx(&0))
     ==> (a / b = inv c <=> a * c = b)`] THEN
  ASM_SIMP_TAC[
COMPLEX_POW_EQ_0; COMPLEX_FIELD
   `~(b = Cx(&0)) ==> (a / b = c <=> b * c = a)`] THEN
  REWRITE_TAC[GSYM 
COMPLEX_POW_ADD] THEN AP_TERM_TAC THEN ASM_ARITH_TAC);;
 
let csqrt = new_definition
  `csqrt(z) = if Im(z) = &0 then
                if &0 <= Re(z) then complex(sqrt(Re(z)),&0)
                else complex(&0,sqrt(--Re(z)))
              else complex(sqrt((norm(z) + Re(z)) / &2),
                           (Im(z) / abs(Im(z))) *
                           sqrt((norm(z) - Re(z)) / &2))`;;let CSQRT_UNIQUE = prove
 (`!s z. s pow 2 = z /\ (&0 < Re s \/ Re s = &0 /\ &0 <= Im s)
         ==> csqrt z = s`,
  REPEAT GEN_TAC THEN DISCH_THEN(CONJUNCTS_THEN ASSUME_TAC) THEN
  FIRST_X_ASSUM(SUBST_ALL_TAC o SYM) THEN
  MP_TAC(SPEC `(s:complex) pow 2` 
CSQRT) THEN
  SIMP_TAC[COMPLEX_RING `a pow 2 = b pow 2 <=> a = b \/ a = --b:complex`] THEN
  STRIP_TAC THEN ASM_REWRITE_TAC[COMPLEX_RING `--z = z <=> z = Cx(&0)`] THEN
  FIRST_ASSUM(MP_TAC o AP_TERM `Re`) THEN
  FIRST_X_ASSUM(MP_TAC o AP_TERM `Im`) THEN
  REWRITE_TAC[
RE_NEG; 
IM_NEG; 
COMPLEX_EQ; 
RE_CX; 
IM_CX] THEN
  MP_TAC(SPEC `(s:complex) pow 2` 
CSQRT_PRINCIPAL) THEN
  POP_ASSUM MP_TAC THEN REAL_ARITH_TAC);;
 
let CSQRT_EQ_0 = prove
 (`!z. csqrt z = Cx(&0) <=> z = Cx(&0)`,
  GEN_TAC THEN MP_TAC (SPEC `z:complex` 
CSQRT) THEN CONV_TAC COMPLEX_RING);;
 
let COMPLEX_SUB_POW = prove
 (`!x y n.
        1 <= n ==> x pow n - y pow n =
                   (x - y) * vsum(0..n-1) (\i. x pow i * y pow (n - 1 - i))`,
 
let IN_SEGMENT_CX = prove
 (`!a b x. Cx(x) 
IN segment[Cx(a),Cx(b)] <=>
                a <= x /\ x <= b \/ b <= x /\ x <= a`,
  REPEAT STRIP_TAC THEN REWRITE_TAC[segment; 
IN_ELIM_THM] THEN
  REWRITE_TAC[
COMPLEX_CMUL; GSYM 
CX_ADD; 
CX_INJ; GSYM 
CX_MUL] THEN
  ASM_CASES_TAC `a:real = b` THENL
   [ASM_REWRITE_TAC[REAL_ARITH `(&1 - u) * b + u * b = b`] THEN
    ASM_CASES_TAC `x:real = b` THEN ASM_REWRITE_TAC[REAL_LE_ANTISYM] THEN
    EXISTS_TAC `&0` THEN REWRITE_TAC[
REAL_POS];
    ALL_TAC] THEN
  EQ_TAC THENL
   [DISCH_THEN(X_CHOOSE_THEN `u:real`
     (CONJUNCTS_THEN2 STRIP_ASSUME_TAC SUBST1_TAC)) THEN
    REWRITE_TAC[REAL_ARITH `a <= (&1 - u) * a + u * b <=> &0 <= u * (b - a)`;
      REAL_ARITH `b <= (&1 - u) * a + u * b <=> &0 <= (&1 - u) * (a - b)`;
      REAL_ARITH `(&1 - u) * a + u * b <= a <=> &0 <= u * (a - b)`;
      REAL_ARITH `(&1 - u) * a + u * b <= b <=> &0 <= (&1 - u) * (b - a)`] THEN
    DISJ_CASES_TAC(REAL_ARITH `a <= b \/ b <= a`) THENL
     [DISJ1_TAC; DISJ2_TAC] THEN
    CONJ_TAC THEN MATCH_MP_TAC 
REAL_LE_MUL THEN
    ASM_REAL_ARITH_TAC;
    ALL_TAC] THEN
  STRIP_TAC THENL
   [SUBGOAL_THEN `&0 < b - a` ASSUME_TAC THENL
     [ASM_REAL_ARITH_TAC;
      EXISTS_TAC `(x - a:real) / (b - a)`];
    SUBGOAL_THEN `&0 < a - b` ASSUME_TAC THENL
     [ASM_REAL_ARITH_TAC;
      EXISTS_TAC `(a - x:real) / (a - b)`]] THEN
  (CONJ_TAC THENL
    [ALL_TAC; UNDISCH_TAC `~(a:real = b)` THEN CONV_TAC REAL_FIELD]) THEN
  ASM_SIMP_TAC[
REAL_LE_LDIV_EQ; 
REAL_LE_RDIV_EQ] THEN
  ASM_REAL_ARITH_TAC);;
 
let IN_SEGMENT_CX_GEN = prove
 (`!a b x.
        x 
IN segment[Cx a,Cx b] <=>
        Im(x) = &0 /\ (a <= Re x /\ Re x <= b \/ b <= Re x /\ Re x <= a)`,
  REPEAT STRIP_TAC THEN REWRITE_TAC[GSYM real] THEN
  ASM_CASES_TAC `real x` THENL
   [FIRST_X_ASSUM(SUBST1_TAC o SYM o REWRITE_RULE[
REAL]) THEN
    REWRITE_TAC[
IN_SEGMENT_CX; 
REAL_CX; 
RE_CX] THEN REAL_ARITH_TAC;
    ASM_MESON_TAC[
REAL_SEGMENT; 
REAL_CX]]);;
 
let VSUM_GP = prove
 (`!x m n.
        vsum(m..n) (\i. x pow i) =
                if n < m then Cx(&0)
                else if x = Cx(&1) then Cx(&((n + 1) - m))
                else (x pow m - x pow (SUC n)) / (Cx(&1) - x)`,
 
let VSUM_GP_OFFSET = prove
 (`!x m n. vsum(m..m+n) (\i. x pow i) =
                if x = Cx(&1) then Cx(&n) + Cx(&1)
                else x pow m * (Cx(&1) - x pow (SUC n)) / (Cx(&1) - x)`,
  REPEAT GEN_TAC THEN REWRITE_TAC[
VSUM_GP; ARITH_RULE `~(m + n < m:num)`] THEN
  COND_CASES_TAC THEN ASM_REWRITE_TAC[] THENL
   [REWRITE_TAC[REAL_OF_NUM_ADD; GSYM 
CX_ADD] THEN
    AP_TERM_TAC THEN AP_TERM_TAC THEN ARITH_TAC;
    REWRITE_TAC[
complex_div; 
complex_pow; 
COMPLEX_POW_ADD] THEN
    SIMPLE_COMPLEX_ARITH_TAC]);;
 
let COMPLEX_SUB_POLYFUN = prove
 (`!a x y n.
   1 <= n
   ==> vsum(0..n) (\i. a i * x pow i) - vsum(0..n) (\i. a i * y pow i) =
       (x - y) *
       vsum(0..n-1) (\j. vsum(j+1..n) (\i. a i * y pow (i - j - 1)) * x pow j)`,
 
let COMPLEX_SUB_POLYFUN_ALT = prove
 (`!a x y n.
    1 <= n
    ==> vsum(0..n) (\i. a i * x pow i) - vsum(0..n) (\i. a i * y pow i) =
        (x - y) *
        vsum(0..n-1) (\j. vsum(0..n-j-1) (\k. a(j+k+1) * y pow k) * x pow j)`,
  REPEAT STRIP_TAC THEN ASM_SIMP_TAC[
COMPLEX_SUB_POLYFUN] THEN AP_TERM_TAC THEN
  MATCH_MP_TAC 
VSUM_EQ_NUMSEG THEN X_GEN_TAC `j:num` THEN REPEAT STRIP_TAC THEN
  REWRITE_TAC[] THEN AP_THM_TAC THEN AP_TERM_TAC THEN
  MATCH_MP_TAC 
VSUM_EQ_GENERAL_INVERSES THEN
  MAP_EVERY EXISTS_TAC
   [`\i. i - (j + 1)`; `\k. j + k + 1`] THEN
  REWRITE_TAC[
IN_NUMSEG] THEN REPEAT STRIP_TAC THEN
  TRY(BINOP_TAC THEN AP_TERM_TAC) THEN ASM_ARITH_TAC);;
 
let COMPLEX_POLYFUN_EQ_CONST = prove
 (`!n c k. (!z. vsum(0..n) (\i. c i * z pow i) = k) <=>
           c 0 = k /\ (!i. i 
IN 1..n ==> c i = Cx(&0))`,
  REPEAT GEN_TAC THEN MATCH_MP_TAC 
EQ_TRANS THEN EXISTS_TAC
   `!x. vsum(0..n) (\i. (if i = 0 then c 0 - k else c i) * x pow i) =
        Cx(&0)` THEN
  CONJ_TAC THENL
   [SIMP_TAC[
VSUM_CLAUSES_LEFT; 
LE_0; 
complex_pow; 
COMPLEX_MUL_RID] THEN
    REWRITE_TAC[COMPLEX_RING `(c - k) + s = Cx(&0) <=> c + s = k`] THEN
    AP_TERM_TAC THEN ABS_TAC THEN AP_THM_TAC THEN AP_TERM_TAC THEN
    AP_TERM_TAC THEN MATCH_MP_TAC 
VSUM_EQ THEN GEN_TAC THEN
    REWRITE_TAC[
IN_NUMSEG] THEN
    COND_CASES_TAC THEN ASM_REWRITE_TAC[ARITH];
    REWRITE_TAC[
COMPLEX_POLYFUN_EQ_0; 
IN_NUMSEG; 
LE_0] THEN
    GEN_REWRITE_TAC LAND_CONV [MESON[]
     `(!n. P n) <=> P 0 /\ (!n. ~(n = 0) ==> P n)`] THEN
    SIMP_TAC[
LE_0; 
COMPLEX_SUB_0] THEN MESON_TAC[
LE_1]]);;
 
let CPRODUCT_CLAUSES_NUMSEG = prove
 (`(!m. cproduct(m..0) f = if m = 0 then f(0) else Cx(&1)) /\
   (!m n. cproduct(m..SUC n) f = if m <= SUC n then cproduct(m..n) f * f(SUC n)
                                 else cproduct(m..n) f)`,
 
let th = prove
 (`(!f g s.   (!x. x 
IN s ==> f(x) = g(x))
              ==> cproduct s (\i. f(i)) = cproduct s g) /\
   (!f g a b. (!i. a <= i /\ i <= b ==> f(i) = g(i))
              ==> cproduct(a..b) (\i. f(i)) = cproduct(a..b) g) /\
   (!f g p.   (!x. p x ==> f x = g x)
              ==> cproduct {y | p y} (\i. f(i)) = cproduct {y | p y} g)`,