/
bigcomplex.go
90 lines (66 loc) · 1.93 KB
/
bigcomplex.go
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
package fractal
import "math/big"
// BigComplex
type BigComplex struct {
r, i big.Float
}
/* Multiplying two complex numbers looks like a quadratic equation
*
* (x + yi)(w + zi) == w * w + xzi + wyi + yzi^2
* And since i^2 = -1, this becomes ww-yz + (xz+wy)i
*/
func (z *BigComplex) Mul(x, y *BigComplex) *BigComplex {
rSquared := new(big.Float).Copy(&x.r)
rSquared.Mul(rSquared, &y.r)
iSquared := new(big.Float).Copy(&x.i)
iSquared.Mul(iSquared, &y.i)
r := new(big.Float).Copy(rSquared)
r.Sub(r, iSquared)
i := new(big.Float).Copy(&x.r)
i.Mul(i, &y.i)
tmp := new(big.Float).Copy(&x.i)
tmp.Mul(tmp, &y.r)
i.Add(i, tmp)
z.r = *r
z.i = *i
return z
}
func (z *BigComplex) Add(x, y *BigComplex) *BigComplex {
z.r.Add(&x.r, &y.r)
z.i.Add(&x.i, &y.i)
return z
}
func (z *BigComplex) Sub(x, y *BigComplex) *BigComplex {
z.r.Sub(&x.r, &y.r)
z.i.Sub(&x.i, &y.i)
return z
}
/*
x + yi / u + vi = (xu + yv) + (-xv + yu)i / u^2 + v^2
*/
func (z *BigComplex) Quo(x, y *BigComplex) *BigComplex {
divisor := new(big.Float).Set(&y.r)
divisor.Mul(divisor, divisor) // u^2
tmp := new(big.Float).Set(&y.i)
tmp.Mul(tmp, tmp) // v^2
divisor.Add(divisor, tmp) // u^2 + v^2
xu := new(big.Float).SetFloat64(0.0).Mul(&x.r, &y.r)
yv := new(big.Float).SetFloat64(0.0).Mul(&x.i, &y.i)
negXv := new(big.Float).SetFloat64(0.0).Mul(&x.r, &y.i)
negXv.Neg(negXv)
yu := new(big.Float).SetFloat64(0.0).Mul(&x.i, &y.r)
dividend_real := new(big.Float).Add(xu, yv)
dividend_imag := new(big.Float).Add(negXv, yu)
z.r = *new(big.Float).Quo(dividend_real, divisor)
z.i = *new(big.Float).Quo(dividend_imag, divisor)
return z
}
// sqrt(x.real^2 + x.Imag^2)
func Abs(x *BigComplex) *big.Float {
realsquared := new(big.Float).Set(&x.r)
realsquared.Mul(realsquared, realsquared)
imagsquared := new(big.Float).Set(&x.i)
imagsquared.Mul(imagsquared, imagsquared)
realsquared.Add(realsquared, imagsquared)
return realsquared.Sqrt(realsquared)
}