1// Copyright ©2018 The Gonum Authors. All rights reserved.
2// Use of this source code is governed by a BSD-style
3// license that can be found in the LICENSE file.
4
5package dualcmplx
6
7import (
8	"fmt"
9	"math"
10	"math/cmplx"
11	"strings"
12)
13
14// Number is a float64 precision anti-commutative dual complex number.
15type Number struct {
16	Real, Dual complex128
17}
18
19// Format implements fmt.Formatter.
20func (d Number) Format(fs fmt.State, c rune) {
21	prec, pOk := fs.Precision()
22	if !pOk {
23		prec = -1
24	}
25	width, wOk := fs.Width()
26	if !wOk {
27		width = -1
28	}
29	switch c {
30	case 'v':
31		if fs.Flag('#') {
32			fmt.Fprintf(fs, "%T{Real:%#v, Dual:%#v}", d, d.Real, d.Dual)
33			return
34		}
35		if fs.Flag('+') {
36			fmt.Fprintf(fs, "{Real:%+v, Dual:%+v}", d.Real, d.Dual)
37			return
38		}
39		c = 'g'
40		prec = -1
41		fallthrough
42	case 'e', 'E', 'f', 'F', 'g', 'G':
43		fre := fmtString(fs, c, prec, width, false)
44		fim := fmtString(fs, c, prec, width, true)
45		fmt.Fprintf(fs, fmt.Sprintf("(%s+%[2]sϵ)", fre, fim), d.Real, d.Dual)
46	default:
47		fmt.Fprintf(fs, "%%!%c(%T=%[2]v)", c, d)
48		return
49	}
50}
51
52// This is horrible, but it's what we have.
53func fmtString(fs fmt.State, c rune, prec, width int, wantPlus bool) string {
54	var b strings.Builder
55	b.WriteByte('%')
56	for _, f := range "0+- " {
57		if fs.Flag(int(f)) || (f == '+' && wantPlus) {
58			b.WriteByte(byte(f))
59		}
60	}
61	if width >= 0 {
62		fmt.Fprint(&b, width)
63	}
64	if prec >= 0 {
65		b.WriteByte('.')
66		if prec > 0 {
67			fmt.Fprint(&b, prec)
68		}
69	}
70	b.WriteRune(c)
71	return b.String()
72}
73
74// Add returns the sum of x and y.
75func Add(x, y Number) Number {
76	return Number{
77		Real: x.Real + y.Real,
78		Dual: x.Dual + y.Dual,
79	}
80}
81
82// Sub returns the difference of x and y, x-y.
83func Sub(x, y Number) Number {
84	return Number{
85		Real: x.Real - y.Real,
86		Dual: x.Dual - y.Dual,
87	}
88}
89
90// Mul returns the dual product of x and y, x×y.
91func Mul(x, y Number) Number {
92	return Number{
93		Real: x.Real * y.Real,
94		Dual: x.Real*y.Dual + x.Dual*cmplx.Conj(y.Real),
95	}
96}
97
98// Inv returns the dual inverse of d.
99func Inv(d Number) Number {
100	return Number{
101		Real: 1 / d.Real,
102		Dual: -d.Dual / (d.Real * cmplx.Conj(d.Real)),
103	}
104}
105
106// Conj returns the conjugate of d₁+d₂ϵ, d̅₁+d₂ϵ.
107func Conj(d Number) Number {
108	return Number{
109		Real: cmplx.Conj(d.Real),
110		Dual: d.Dual,
111	}
112}
113
114// Scale returns d scaled by f.
115func Scale(f float64, d Number) Number {
116	return Number{Real: complex(f, 0) * d.Real, Dual: complex(f, 0) * d.Dual}
117}
118
119// Abs returns the absolute value of d.
120func Abs(d Number) float64 {
121	return cmplx.Abs(d.Real)
122}
123
124// PowReal returns d**p, the base-d exponential of p.
125//
126// Special cases are (in order):
127//	PowReal(NaN+xϵ, ±0) = 1+NaNϵ for any x
128//	Pow(0+xϵ, y) = 0+Infϵ for all y < 1.
129//	Pow(0+xϵ, y) = 0 for all y > 1.
130//	PowReal(x, ±0) = 1 for any x
131//	PowReal(1+xϵ, y) = 1+xyϵ for any y
132//	Pow(Inf, y) = +Inf+NaNϵ for y > 0
133//	Pow(Inf, y) = +0+NaNϵ for y < 0
134//	PowReal(x, 1) = x for any x
135//	PowReal(NaN+xϵ, y) = NaN+NaNϵ
136//	PowReal(x, NaN) = NaN+NaNϵ
137//	PowReal(-1, ±Inf) = 1
138//	PowReal(x+0ϵ, +Inf) = +Inf+NaNϵ for |x| > 1
139//	PowReal(x+yϵ, +Inf) = +Inf for |x| > 1
140//	PowReal(x, -Inf) = +0+NaNϵ for |x| > 1
141//	PowReal(x, +Inf) = +0+NaNϵ for |x| < 1
142//	PowReal(x+0ϵ, -Inf) = +Inf+NaNϵ for |x| < 1
143//	PowReal(x, -Inf) = +Inf-Infϵ for |x| < 1
144//	PowReal(+Inf, y) = +Inf for y > 0
145//	PowReal(+Inf, y) = +0 for y < 0
146//	PowReal(-Inf, y) = Pow(-0, -y)
147func PowReal(d Number, p float64) Number {
148	switch {
149	case p == 0:
150		switch {
151		case cmplx.IsNaN(d.Real):
152			return Number{Real: 1, Dual: cmplx.NaN()}
153		case d.Real == 0, cmplx.IsInf(d.Real):
154			return Number{Real: 1}
155		}
156	case p == 1:
157		if cmplx.IsInf(d.Real) {
158			d.Dual = cmplx.NaN()
159		}
160		return d
161	case math.IsInf(p, 1):
162		if d.Real == -1 {
163			return Number{Real: 1, Dual: cmplx.NaN()}
164		}
165		if Abs(d) > 1 {
166			if d.Dual == 0 {
167				return Number{Real: cmplx.Inf(), Dual: cmplx.NaN()}
168			}
169			return Number{Real: cmplx.Inf(), Dual: cmplx.Inf()}
170		}
171		return Number{Real: 0, Dual: cmplx.NaN()}
172	case math.IsInf(p, -1):
173		if d.Real == -1 {
174			return Number{Real: 1, Dual: cmplx.NaN()}
175		}
176		if Abs(d) > 1 {
177			return Number{Real: 0, Dual: cmplx.NaN()}
178		}
179		if d.Dual == 0 {
180			return Number{Real: cmplx.Inf(), Dual: cmplx.NaN()}
181		}
182		return Number{Real: cmplx.Inf(), Dual: cmplx.Inf()}
183	case math.IsNaN(p):
184		return Number{Real: cmplx.NaN(), Dual: cmplx.NaN()}
185	case d.Real == 0:
186		if p < 1 {
187			return Number{Real: d.Real, Dual: cmplx.Inf()}
188		}
189		return Number{Real: d.Real}
190	case cmplx.IsInf(d.Real):
191		if p < 0 {
192			return Number{Real: 0, Dual: cmplx.NaN()}
193		}
194		return Number{Real: cmplx.Inf(), Dual: cmplx.NaN()}
195	}
196	return Pow(d, Number{Real: complex(p, 0)})
197}
198
199// Pow returns d**p, the base-d exponential of p.
200func Pow(d, p Number) Number {
201	return Exp(Mul(p, Log(d)))
202}
203
204// Sqrt returns the square root of d.
205//
206// Special cases are:
207//	Sqrt(+Inf) = +Inf
208//	Sqrt(±0) = (±0+Infϵ)
209//	Sqrt(x < 0) = NaN
210//	Sqrt(NaN) = NaN
211func Sqrt(d Number) Number {
212	return PowReal(d, 0.5)
213}
214
215// Exp returns e**q, the base-e exponential of d.
216//
217// Special cases are:
218//	Exp(+Inf) = +Inf
219//	Exp(NaN) = NaN
220// Very large values overflow to 0 or +Inf.
221// Very small values underflow to 1.
222func Exp(d Number) Number {
223	fn := cmplx.Exp(d.Real)
224	if imag(d.Real) == 0 {
225		return Number{Real: fn, Dual: fn * d.Dual}
226	}
227	conj := cmplx.Conj(d.Real)
228	return Number{
229		Real: fn,
230		Dual: ((fn - cmplx.Exp(conj)) / (d.Real - conj)) * d.Dual,
231	}
232}
233
234// Log returns the natural logarithm of d.
235//
236// Special cases are:
237//	Log(+Inf) = (+Inf+0ϵ)
238//	Log(0) = (-Inf±Infϵ)
239//	Log(x < 0) = NaN
240//	Log(NaN) = NaN
241func Log(d Number) Number {
242	fn := cmplx.Log(d.Real)
243	switch {
244	case d.Real == 0:
245		return Number{
246			Real: fn,
247			Dual: complex(math.Copysign(math.Inf(1), real(d.Real)), math.NaN()),
248		}
249	case imag(d.Real) == 0:
250		return Number{
251			Real: fn,
252			Dual: d.Dual / d.Real,
253		}
254	case cmplx.IsInf(d.Real):
255		return Number{
256			Real: fn,
257			Dual: 0,
258		}
259	}
260	conj := cmplx.Conj(d.Real)
261	return Number{
262		Real: fn,
263		Dual: ((fn - cmplx.Log(conj)) / (d.Real - conj)) * d.Dual,
264	}
265}
266