1 (* ========================================================================= *)
3 (* ========================================================================= *)
5 needs "Complex/complex_transc.ml";;
9 (* ------------------------------------------------------------------------- *)
10 (* Define cube roots (it doesn't matter which one we choose here) *)
11 (* ------------------------------------------------------------------------- *)
13 let ccbrt = new_definition
14 `ccbrt(z) = if z = Cx(&0) then Cx(&0) else cexp(clog(z) / Cx(&3))`;;
17 (`!z. ccbrt(z) pow 3 = z`,
18 GEN_TAC THEN REWRITE_TAC[ccbrt] THEN COND_CASES_TAC THEN
19 ASM_REWRITE_TAC[] THENL [CONV_TAC COMPLEX_RING; ALL_TAC] THEN
20 REWRITE_TAC[COMPLEX_FIELD `z pow 3 = z * z * z:complex`; GSYM CEXP_ADD] THEN
21 REWRITE_TAC[COMPLEX_FIELD `z / Cx(&3) + z / Cx(&3) + z / Cx(&3) = z`] THEN
22 ASM_SIMP_TAC[CEXP_CLOG]);;
24 (* ------------------------------------------------------------------------- *)
25 (* The reduction to a simpler form. *)
26 (* ------------------------------------------------------------------------- *)
28 let CUBIC_REDUCTION = COMPLEX_FIELD
30 x = y - b / (Cx(&3) * a) /\
31 p = (Cx(&3) * a * c - b pow 2) / (Cx(&9) * a pow 2) /\
32 q = (Cx(&9) * a * b * c - Cx(&2) * b pow 3 - Cx(&27) * a pow 2 * d) /
34 ==> (a * x pow 3 + b * x pow 2 + c * x + d = Cx(&0) <=>
35 y pow 3 + Cx(&3) * p * y - Cx(&2) * q = Cx(&0))`;;
37 (* ------------------------------------------------------------------------- *)
38 (* The solutions of the special form. *)
39 (* ------------------------------------------------------------------------- *)
41 let CUBIC_BASIC = COMPLEX_FIELD
43 s pow 2 = q pow 2 + p pow 3 /\
44 (s1 pow 3 = if p = Cx(&0) then Cx(&2) * q else q + s) /\
45 s2 = --s1 * (Cx(&1) + i * t) / Cx(&2) /\
46 s3 = --s1 * (Cx(&1) - i * t) / Cx(&2) /\
47 i pow 2 + Cx(&1) = Cx(&0) /\
49 ==> if p = Cx(&0) then
50 (y pow 3 + Cx(&3) * p * y - Cx(&2) * q = Cx(&0) <=>
51 y = s1 \/ y = s2 \/ y = s3)
54 (y pow 3 + Cx(&3) * p * y - Cx(&2) * q = Cx(&0) <=>
55 (y = s1 - p / s1 \/ y = s2 - p / s2 \/ y = s3 - p / s3))`;;
57 (* ------------------------------------------------------------------------- *)
58 (* Explicit formula for the roots (doesn't matter which square/cube roots). *)
59 (* ------------------------------------------------------------------------- *)
63 ==> let p = (Cx(&3) * a * c - b pow 2) / (Cx(&9) * a pow 2)
64 and q = (Cx(&9) * a * b * c - Cx(&2) * b pow 3 - Cx(&27) * a pow 2 * d) /
65 (Cx(&54) * a pow 3) in
66 let s = csqrt(q pow 2 + p pow 3) in
67 let s1 = if p = Cx(&0) then ccbrt(Cx(&2) * q) else ccbrt(q + s) in
68 let s2 = --s1 * (Cx(&1) + ii * csqrt(Cx(&3))) / Cx(&2)
69 and s3 = --s1 * (Cx(&1) - ii * csqrt(Cx(&3))) / Cx(&2) in
71 a * x pow 3 + b * x pow 2 + c * x + d = Cx(&0) <=>
72 x = s1 - b / (Cx(&3) * a) \/
73 x = s2 - b / (Cx(&3) * a) \/
74 x = s3 - b / (Cx(&3) * a)
77 (a * x pow 3 + b * x pow 2 + c * x + d = Cx(&0) <=>
78 x = s1 - p / s1 - b / (Cx(&3) * a) \/
79 x = s2 - p / s2 - b / (Cx(&3) * a) \/
80 x = s3 - p / s3 - b / (Cx(&3) * a))`,
81 DISCH_TAC THEN REPEAT LET_TAC THEN
82 ABBREV_TAC `y = x + b / (Cx(&3) * a)` THEN
84 `a * x pow 3 + b * x pow 2 + c * x + d = Cx(&0) <=>
85 y pow 3 + Cx(&3) * p * y - Cx(&2) * q = Cx(&0)`
87 [MATCH_MP_TAC CUBIC_REDUCTION THEN ASM_REWRITE_TAC[] THEN
88 EXPAND_TAC "y" THEN CONV_TAC COMPLEX_RING;
90 ONCE_REWRITE_TAC[COMPLEX_RING `x = a - b <=> x + b = (a:complex)`] THEN
91 ASM_REWRITE_TAC[] THEN MATCH_MP_TAC CUBIC_BASIC THEN
93 [`ii`; `csqrt(Cx(&3))`; `csqrt (q pow 2 + p pow 3)`] THEN
94 ASM_REWRITE_TAC[] THEN REPEAT CONJ_TAC THENL
95 [ASM_MESON_TAC[CSQRT];
97 MP_TAC COMPLEX_POW_II_2 THEN CONV_TAC COMPLEX_RING;
98 ASM_MESON_TAC[CSQRT]]);;