(* ========================================================================= *)
(* Complex transcendental functions.                                         *)
(* ========================================================================= *)

needs "Library/transc.ml";;
needs "Library/floor.ml";;
needs "Complex/complexnumbers.ml";;

unparse_as_infix "exp";;
remove_interface "exp";;

(* ------------------------------------------------------------------------- *)
(* Complex square roots.                                                     *)
(* ------------------------------------------------------------------------- *)

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 COMPLEX_NORM_GE_RE_IM = 
prove (`!z. abs(Re(z)) <= norm(z) /\ abs(Im(z)) <= norm(z)`,
GEN_TAC THEN ONCE_REWRITE_TAC[GSYM POW_2_SQRT_ABS] THEN REWRITE_TAC[complex_norm] THEN CONJ_TAC THEN MATCH_MP_TAC SQRT_MONO_LE THEN ASM_SIMP_TAC[REAL_LE_ADDR; REAL_LE_ADDL; REAL_POW_2; REAL_LE_SQUARE]);;
let CSQRT = 
prove (`!z. csqrt(z) pow 2 = z`,
GEN_TAC THEN REWRITE_TAC[COMPLEX_POW_2; csqrt] THEN COND_CASES_TAC THENL [COND_CASES_TAC THEN ASM_REWRITE_TAC[CX_DEF; complex_mul; RE; IM; REAL_MUL_RZERO; REAL_MUL_LZERO; REAL_SUB_LZERO; REAL_SUB_RZERO; REAL_ADD_LID; COMPLEX_EQ] THEN REWRITE_TAC[REAL_NEG_EQ; GSYM REAL_POW_2] THEN ASM_SIMP_TAC[SQRT_POW_2; REAL_ARITH `~(&0 <= x) ==> &0 <= --x`]; ALL_TAC] THEN REWRITE_TAC[complex_mul; RE; IM] THEN ONCE_REWRITE_TAC[REAL_ARITH `(s * s - (i * s') * (i * s') = s * s - (i * i) * (s' * s')) /\ (s * i * s' + (i * s')* s = &2 * i * s * s')`] THEN REWRITE_TAC[GSYM REAL_POW_2] THEN SUBGOAL_THEN `&0 <= norm(z) + Re(z) /\ &0 <= norm(z) - Re(z)` STRIP_ASSUME_TAC THENL [MP_TAC(SPEC `z:complex` COMPLEX_NORM_GE_RE_IM) THEN REAL_ARITH_TAC; ALL_TAC] THEN ASM_SIMP_TAC[REAL_LE_DIV; REAL_POS; GSYM SQRT_MUL; SQRT_POW_2] THEN REWRITE_TAC[COMPLEX_EQ; RE; IM] THEN CONJ_TAC THENL [ASM_SIMP_TAC[REAL_POW_DIV; REAL_POW2_ABS; REAL_POW_EQ_0; REAL_DIV_REFL] THEN REWRITE_TAC[real_div; REAL_MUL_LID; GSYM REAL_SUB_RDISTRIB] THEN REWRITE_TAC[REAL_ARITH `(m + r) - (m - r) = r * &2`] THEN REWRITE_TAC[GSYM REAL_MUL_ASSOC] THEN CONV_TAC REAL_RAT_REDUCE_CONV THEN REWRITE_TAC[REAL_MUL_RID]; ALL_TAC] THEN REWRITE_TAC[real_div] THEN ONCE_REWRITE_TAC[AC REAL_MUL_AC `(a * b) * a' * b = (a * a') * (b * b:real)`] THEN REWRITE_TAC[REAL_DIFFSQ] THEN REWRITE_TAC[complex_norm; GSYM REAL_POW_2] THEN SIMP_TAC[SQRT_POW_2; REAL_LE_ADD; REWRITE_RULE[GSYM REAL_POW_2] REAL_LE_SQUARE] THEN REWRITE_TAC[REAL_ADD_SUB; GSYM REAL_POW_MUL] THEN REWRITE_TAC[POW_2_SQRT_ABS] THEN REWRITE_TAC[REAL_ABS_MUL; REAL_ABS_INV; REAL_ABS_NUM] THEN ONCE_REWRITE_TAC[AC REAL_MUL_AC `&2 * (i * a') * a * h = i * (&2 * h) * a * a'`] THEN CONV_TAC REAL_RAT_REDUCE_CONV THEN REWRITE_TAC[REAL_MUL_LID; GSYM real_div] THEN ASM_SIMP_TAC[REAL_DIV_REFL; REAL_ABS_ZERO; REAL_MUL_RID]);;
(* ------------------------------------------------------------------------- *) (* Complex exponential. *) (* ------------------------------------------------------------------------- *)
let cexp = new_definition
 `cexp z = Cx(exp(Re z)) * complex(cos(Im z),sin(Im z))`;;
let CX_CEXP = 
prove (`!x. Cx(exp x) = cexp(Cx x)`,
REWRITE_TAC[cexp; CX_DEF; RE; IM; SIN_0; COS_0] THEN REWRITE_TAC[GSYM CX_DEF; GSYM CX_MUL; REAL_MUL_RID]);;
let CEXP_0 = 
prove (`cexp(Cx(&0)) = Cx(&1)`,
REWRITE_TAC[GSYM CX_CEXP; REAL_EXP_0]);;
let CEXP_ADD = 
prove (`!w z. cexp(w + z) = cexp(w) * cexp(z)`,
REWRITE_TAC[COMPLEX_EQ; cexp; complex_mul; complex_add; RE; IM; CX_DEF] THEN REWRITE_TAC[REAL_EXP_ADD; SIN_ADD; COS_ADD] THEN CONV_TAC REAL_RING);;
let CEXP_MUL = 
prove (`!n z. cexp(Cx(&n) * z) = cexp(z) pow n`,
INDUCT_TAC THEN REWRITE_TAC[complex_pow] THEN REWRITE_TAC[COMPLEX_MUL_LZERO; CEXP_0] THEN REWRITE_TAC[GSYM REAL_OF_NUM_SUC; COMPLEX_ADD_RDISTRIB; CX_ADD] THEN ASM_REWRITE_TAC[CEXP_ADD; COMPLEX_MUL_LID] THEN REWRITE_TAC[COMPLEX_MUL_AC]);;
let CEXP_NONZERO = 
prove (`!z. ~(cexp z = Cx(&0))`,
GEN_TAC THEN REWRITE_TAC[cexp; COMPLEX_ENTIRE; CX_INJ; REAL_EXP_NZ] THEN REWRITE_TAC[CX_DEF; RE; IM; COMPLEX_EQ] THEN MP_TAC(SPEC `Im z` SIN_CIRCLE) THEN CONV_TAC REAL_RING);;
let CEXP_NEG_LMUL = 
prove (`!z. cexp(--z) * cexp(z) = Cx(&1)`,
REWRITE_TAC[GSYM CEXP_ADD; COMPLEX_ADD_LINV; CEXP_0]);;
let CEXP_NEG_RMUL = 
prove (`!z. cexp(z) * cexp(--z) = Cx(&1)`,
REWRITE_TAC[GSYM CEXP_ADD; COMPLEX_ADD_RINV; CEXP_0]);;
let CEXP_NEG = 
prove (`!z. cexp(--z) = inv(cexp z)`,
let CEXP_SUB = 
prove (`!w z. cexp(w - z) = cexp(w) / cexp(z)`,
(* ------------------------------------------------------------------------- *) (* Complex trig functions. *) (* ------------------------------------------------------------------------- *)
let ccos = new_definition
  `ccos z = (cexp(ii * z) + cexp(--ii * z)) / Cx(&2)`;;
let csin = new_definition
  `csin z = (cexp(ii * z) - cexp(--ii * z)) / (Cx(&2) * ii)`;;
let CX_CSIN,CX_CCOS = (CONJ_PAIR o prove) (`(!x. Cx(sin x) = csin(Cx x)) /\ (!x. Cx(cos x) = ccos(Cx x))`, REWRITE_TAC[csin; ccos; cexp; CX_DEF; ii; RE; IM; complex_mul; complex_add; REAL_MUL_RZERO; REAL_MUL_LZERO; REAL_SUB_RZERO; REAL_MUL_LID; complex_neg; REAL_EXP_0; REAL_ADD_LID; REAL_MUL_LNEG; REAL_NEG_0; REAL_ADD_RID; complex_sub; SIN_NEG; COS_NEG; GSYM REAL_MUL_2; GSYM real_sub; complex_div; REAL_SUB_REFL; complex_inv; REAL_SUB_RNEG] THEN CONJ_TAC THEN GEN_TAC THEN CONV_TAC REAL_RAT_REDUCE_CONV THEN REWRITE_TAC[REAL_MUL_RZERO] THEN AP_TERM_TAC THEN AP_THM_TAC THEN AP_TERM_TAC THEN CONV_TAC REAL_RING);;
let CSIN_0 = 
prove (`csin(Cx(&0)) = Cx(&0)`,
REWRITE_TAC[GSYM CX_CSIN; SIN_0]);;
let CCOS_0 = 
prove (`ccos(Cx(&0)) = Cx(&1)`,
REWRITE_TAC[GSYM CX_CCOS; COS_0]);;
let CSIN_CIRCLE = 
prove (`!z. csin(z) pow 2 + ccos(z) pow 2 = Cx(&1)`,
GEN_TAC THEN REWRITE_TAC[csin; ccos] THEN MP_TAC(SPEC `ii * z` CEXP_NEG_LMUL) THEN MP_TAC COMPLEX_POW_II_2 THEN REWRITE_TAC[COMPLEX_MUL_LNEG] THEN CONV_TAC COMPLEX_FIELD);;
let CSIN_ADD = 
prove (`!w z. csin(w + z) = csin(w) * ccos(z) + ccos(w) * csin(z)`,
REPEAT GEN_TAC THEN REWRITE_TAC[csin; ccos; COMPLEX_ADD_LDISTRIB; CEXP_ADD] THEN MP_TAC COMPLEX_POW_II_2 THEN CONV_TAC COMPLEX_FIELD);;
let CCOS_ADD = 
prove (`!w z. ccos(w + z) = ccos(w) * ccos(z) - csin(w) * csin(z)`,
REPEAT GEN_TAC THEN REWRITE_TAC[csin; ccos; COMPLEX_ADD_LDISTRIB; CEXP_ADD] THEN MP_TAC COMPLEX_POW_II_2 THEN CONV_TAC COMPLEX_FIELD);;
let CSIN_NEG = 
prove (`!z. csin(--z) = --(csin(z))`,
REWRITE_TAC[csin; COMPLEX_MUL_LNEG; COMPLEX_MUL_RNEG; COMPLEX_NEG_NEG] THEN GEN_TAC THEN MP_TAC COMPLEX_POW_II_2 THEN CONV_TAC COMPLEX_FIELD);;
let CCOS_NEG = 
prove (`!z. ccos(--z) = ccos(z)`,
REWRITE_TAC[ccos; COMPLEX_MUL_LNEG; COMPLEX_MUL_RNEG; COMPLEX_NEG_NEG] THEN GEN_TAC THEN MP_TAC COMPLEX_POW_II_2 THEN CONV_TAC COMPLEX_FIELD);;
let CSIN_DOUBLE = 
prove (`!z. csin(Cx(&2) * z) = Cx(&2) * csin(z) * ccos(z)`,
REWRITE_TAC[COMPLEX_RING `Cx(&2) * x = x + x`; CSIN_ADD] THEN CONV_TAC COMPLEX_RING);;
let CCOS_DOUBLE = 
prove (`!z. ccos(Cx(&2) * z) = (ccos(z) pow 2) - (csin(z) pow 2)`,
REWRITE_TAC[COMPLEX_RING `Cx(&2) * x = x + x`; CCOS_ADD] THEN CONV_TAC COMPLEX_RING);;
(* ------------------------------------------------------------------------- *) (* Euler and de Moivre formulas. *) (* ------------------------------------------------------------------------- *)
let CEXP_EULER = 
prove (`!z. cexp(ii * z) = ccos(z) + ii * csin(z)`,
REWRITE_TAC[ccos; csin] THEN MP_TAC COMPLEX_POW_II_2 THEN CONV_TAC COMPLEX_FIELD);;
let DEMOIVRE = 
prove (`!z n. (ccos z + ii * csin z) pow n = ccos(Cx(&n) * z) + ii * csin(Cx(&n) * z)`,
REWRITE_TAC[GSYM CEXP_EULER; GSYM CEXP_MUL] THEN REWRITE_TAC[COMPLEX_MUL_AC]);;
(* ------------------------------------------------------------------------- *) (* Some lemmas. *) (* ------------------------------------------------------------------------- *)
let EXISTS_COMPLEX = 
prove (`!P. (?z. P (Re z) (Im z)) <=> ?x y. P x y`,
MESON_TAC[RE; IM; COMPLEX]);;
let COMPLEX_UNIMODULAR_POLAR = 
prove (`!z. (norm z = &1) ==> ?x. z = complex(cos(x),sin(x))`,
GEN_TAC THEN DISCH_THEN(MP_TAC o C AP_THM `2` o AP_TERM `(pow):real->num->real`) THEN REWRITE_TAC[complex_norm] THEN SIMP_TAC[REAL_POW_2; REWRITE_RULE[REAL_POW_2] SQRT_POW_2; REAL_LE_SQUARE; REAL_LE_ADD] THEN REWRITE_TAC[GSYM REAL_POW_2; REAL_MUL_LID] THEN DISCH_THEN(X_CHOOSE_TAC `t:real` o MATCH_MP CIRCLE_SINCOS) THEN EXISTS_TAC `t:real` THEN ASM_REWRITE_TAC[COMPLEX_EQ; RE; IM]);;
let SIN_INTEGER_2PI = 
prove (`!n. integer n ==> sin((&2 * pi) * n) = &0`,
REWRITE_TAC[integer; REAL_ARITH `abs(x) = &n <=> x = &n \/ x = -- &n`] THEN REPEAT STRIP_TAC THEN ASM_REWRITE_TAC[REAL_MUL_RNEG; SIN_NEG] THEN REWRITE_TAC[GSYM REAL_MUL_ASSOC; SIN_DOUBLE] THEN REWRITE_TAC[REAL_ARITH `pi * &n = &n * pi`; SIN_NPI] THEN REWRITE_TAC[REAL_MUL_LZERO; REAL_MUL_RZERO; REAL_NEG_0]);;
let COS_INTEGER_2PI = 
prove (`!n. integer n ==> cos((&2 * pi) * n) = &1`,
REWRITE_TAC[integer; REAL_ARITH `abs(x) = &n <=> x = &n \/ x = -- &n`] THEN REPEAT STRIP_TAC THEN ASM_REWRITE_TAC[REAL_MUL_RNEG; COS_NEG] THEN REWRITE_TAC[GSYM REAL_MUL_ASSOC; COS_DOUBLE] THEN REWRITE_TAC[REAL_ARITH `pi * &n = &n * pi`; SIN_NPI; COS_NPI] THEN REWRITE_TAC[REAL_POW_POW] THEN ONCE_REWRITE_TAC[MULT_SYM] THEN REWRITE_TAC[GSYM REAL_POW_POW; REAL_POW_2] THEN CONV_TAC REAL_RAT_REDUCE_CONV THEN REWRITE_TAC[REAL_POW_ONE; REAL_SUB_RZERO]);;
let SINCOS_PRINCIPAL_VALUE = 
prove (`!x. ?y. (--pi < y /\ y <= pi) /\ (sin(y) = sin(x) /\ cos(y) = cos(x))`,
GEN_TAC THEN EXISTS_TAC `pi - (&2 * pi) * frac((pi - x) / (&2 * pi))` THEN CONJ_TAC THENL [SIMP_TAC[REAL_ARITH `--p < p - x <=> x < (&2 * p) * &1`; REAL_ARITH `p - x <= p <=> (&2 * p) * &0 <= x`; REAL_LT_LMUL_EQ; REAL_LE_LMUL_EQ; REAL_LT_MUL; PI_POS; REAL_OF_NUM_LT; ARITH; FLOOR_FRAC]; REWRITE_TAC[FRAC_FLOOR; REAL_SUB_LDISTRIB] THEN SIMP_TAC[REAL_DIV_LMUL; REAL_ENTIRE; REAL_OF_NUM_EQ; ARITH; REAL_LT_IMP_NZ; PI_POS; REAL_ARITH `a - (a - b - c):real = b + c`; SIN_ADD; COS_ADD] THEN SIMP_TAC[FLOOR_FRAC; SIN_INTEGER_2PI; COS_INTEGER_2PI] THEN CONV_TAC REAL_RING]);;
(* ------------------------------------------------------------------------- *) (* Complex logarithms (the conventional principal value). *) (* ------------------------------------------------------------------------- *)
let clog = new_definition
 `clog z = @w. cexp(w) = z /\ --pi < Im(w) /\ Im(w) <= pi`;;
let CLOG_WORKS = 
prove (`!z. ~(z = Cx(&0)) ==> cexp(clog z) = z /\ --pi < Im(clog z) /\ Im(clog z) <= pi`,
GEN_TAC THEN DISCH_TAC THEN REWRITE_TAC[clog] THEN CONV_TAC SELECT_CONV THEN REWRITE_TAC[cexp; EXISTS_COMPLEX] THEN EXISTS_TAC `ln(norm(z:complex))` THEN SUBGOAL_THEN `exp(ln(norm(z:complex))) = norm(z)` SUBST1_TAC THENL [ASM_MESON_TAC[REAL_EXP_LN; COMPLEX_NORM_NZ]; ALL_TAC] THEN MP_TAC(SPEC `z / Cx(norm z)` COMPLEX_UNIMODULAR_POLAR) THEN ANTS_TAC THENL [ASM_SIMP_TAC[COMPLEX_NORM_DIV; COMPLEX_NORM_CX] THEN ASM_SIMP_TAC[COMPLEX_ABS_NORM; REAL_DIV_REFL; COMPLEX_NORM_ZERO]; ALL_TAC] THEN DISCH_THEN(X_CHOOSE_THEN `x:real` STRIP_ASSUME_TAC) THEN MP_TAC(SPEC `x:real` SINCOS_PRINCIPAL_VALUE) THEN MATCH_MP_TAC MONO_EXISTS THEN X_GEN_TAC `y:real` THEN STRIP_TAC THEN ASM_REWRITE_TAC[] THEN ASM_MESON_TAC[CX_INJ; COMPLEX_DIV_LMUL; COMPLEX_NORM_ZERO]);;
let CEXP_CLOG = 
prove (`!z. ~(z = Cx(&0)) ==> cexp(clog z) = z`,
SIMP_TAC[CLOG_WORKS]);;
(* ------------------------------------------------------------------------- *) (* Unwinding number. *) (* ------------------------------------------------------------------------- *)
let unwinding = new_definition
 `unwinding(z) = (z - clog(cexp z)) / (Cx(&2 * pi) * ii)`;;
let COMPLEX_II_NZ = 
prove (`~(ii = Cx(&0))`,
MP_TAC COMPLEX_POW_II_2 THEN CONV_TAC COMPLEX_RING);;
let UNWINDING_2PI = 
prove (`Cx(&2 * pi) * ii * unwinding(z) = z - clog(cexp z)`,
REWRITE_TAC[unwinding; COMPLEX_MUL_ASSOC] THEN MATCH_MP_TAC COMPLEX_DIV_LMUL THEN REWRITE_TAC[COMPLEX_ENTIRE; CX_INJ; COMPLEX_II_NZ] THEN MP_TAC PI_POS THEN REAL_ARITH_TAC);;
(* ------------------------------------------------------------------------- *) (* An example of how to get nice identities with unwinding number. *) (* ------------------------------------------------------------------------- *)
let CLOG_MUL = 
prove (`!w z. ~(w = Cx(&0)) /\ ~(z = Cx(&0)) ==> clog(w * z) = clog(w) + clog(z) - Cx(&2 * pi) * ii * unwinding(clog w + clog z)`,
REWRITE_TAC[UNWINDING_2PI; COMPLEX_RING `w + z - ((w + z) - c) = c:complex`] THEN ASM_SIMP_TAC[CEXP_ADD; CEXP_CLOG]);;