dual.go

// Copyright ©2018 The Gonum Authors. All rights reserved.
// Use of this source code is governed by a BSD-style
// license that can be found in the LICENSE file.

package dualcmplx

import (
	"fmt"
	"math"
	"math/cmplx"
	"strings"
)

// Number is a float64 precision anti-commutative dual complex number.
type Number struct {
	Real, Dual complex128
}

// Format implements fmt.Formatter.
func (d Number) Format(fs fmt.State, c rune) {
	prec, pOk := fs.Precision()
	if !pOk {
		prec = -1
	}
	width, wOk := fs.Width()
	if !wOk {
		width = -1
	}
	switch c {
	case 'v':
		if fs.Flag('#') {
			fmt.Fprintf(fs, "%T{Real:%#v, Dual:%#v}", d, d.Real, d.Dual)
			return
		}
		if fs.Flag('+') {
			fmt.Fprintf(fs, "{Real:%+v, Dual:%+v}", d.Real, d.Dual)
			return
		}
		c = 'g'
		prec = -1
		fallthrough
	case 'e', 'E', 'f', 'F', 'g', 'G':
		fre := fmtString(fs, c, prec, width, false)
		fim := fmtString(fs, c, prec, width, true)
		fmt.Fprintf(fs, fmt.Sprintf("(%s+%[2]sϵ)", fre, fim), d.Real, d.Dual)
	default:
		fmt.Fprintf(fs, "%%!%c(%T=%[2]v)", c, d)
		return
	}
}

// This is horrible, but it's what we have.
func fmtString(fs fmt.State, c rune, prec, width int, wantPlus bool) string {
	var b strings.Builder
	b.WriteByte('%')
	for _, f := range "0+- " {
		if fs.Flag(int(f)) || (f == '+' && wantPlus) {
			b.WriteByte(byte(f))
		}
	}
	if width >= 0 {
		fmt.Fprint(&b, width)
	}
	if prec >= 0 {
		b.WriteByte('.')
		if prec > 0 {
			fmt.Fprint(&b, prec)
		}
	}
	b.WriteRune(c)
	return b.String()
}

// Add returns the sum of x and y.
func Add(x, y Number) Number {
	return Number{
		Real: x.Real + y.Real,
		Dual: x.Dual + y.Dual,
	}
}

// Sub returns the difference of x and y, x-y.
func Sub(x, y Number) Number {
	return Number{
		Real: x.Real - y.Real,
		Dual: x.Dual - y.Dual,
	}
}

// Mul returns the dual product of x and y, x×y.
func Mul(x, y Number) Number {
	return Number{
		Real: x.Real * y.Real,
		Dual: x.Real*y.Dual + x.Dual*cmplx.Conj(y.Real),
	}
}

// Inv returns the dual inverse of d.
func Inv(d Number) Number {
	return Number{
		Real: 1 / d.Real,
		Dual: -d.Dual / (d.Real * cmplx.Conj(d.Real)),
	}
}

// Conj returns the conjugate of d₁+d₂ϵ, d̅₁+d₂ϵ.
func Conj(d Number) Number {
	return Number{
		Real: cmplx.Conj(d.Real),
		Dual: d.Dual,
	}
}

// Scale returns d scaled by f.
func Scale(f float64, d Number) Number {
	return Number{Real: complex(f, 0) * d.Real, Dual: complex(f, 0) * d.Dual}
}

// Abs returns the absolute value of d.
func Abs(d Number) float64 {
	return cmplx.Abs(d.Real)
}

// PowReal returns d**p, the base-d exponential of p.
//
// Special cases are (in order):
//
//	PowReal(NaN+xϵ, ±0) = 1+NaNϵ for any x
//	Pow(0+xϵ, y) = 0+Infϵ for all y < 1.
//	Pow(0+xϵ, y) = 0 for all y > 1.
//	PowReal(x, ±0) = 1 for any x
//	PowReal(1+xϵ, y) = 1+xyϵ for any y
//	Pow(Inf, y) = +Inf+NaNϵ for y > 0
//	Pow(Inf, y) = +0+NaNϵ for y < 0
//	PowReal(x, 1) = x for any x
//	PowReal(NaN+xϵ, y) = NaN+NaNϵ
//	PowReal(x, NaN) = NaN+NaNϵ
//	PowReal(-1, ±Inf) = 1
//	PowReal(x+0ϵ, +Inf) = +Inf+NaNϵ for |x| > 1
//	PowReal(x+yϵ, +Inf) = +Inf for |x| > 1
//	PowReal(x, -Inf) = +0+NaNϵ for |x| > 1
//	PowReal(x, +Inf) = +0+NaNϵ for |x| < 1
//	PowReal(x+0ϵ, -Inf) = +Inf+NaNϵ for |x| < 1
//	PowReal(x, -Inf) = +Inf-Infϵ for |x| < 1
//	PowReal(+Inf, y) = +Inf for y > 0
//	PowReal(+Inf, y) = +0 for y < 0
//	PowReal(-Inf, y) = Pow(-0, -y)
func PowReal(d Number, p float64) Number {
	switch {
	case p == 0:
		switch {
		case cmplx.IsNaN(d.Real):
			return Number{Real: 1, Dual: cmplx.NaN()}
		case d.Real == 0, cmplx.IsInf(d.Real):
			return Number{Real: 1}
		}
	case p == 1:
		if cmplx.IsInf(d.Real) {
			d.Dual = cmplx.NaN()
		}
		return d
	case math.IsInf(p, 1):
		if d.Real == -1 {
			return Number{Real: 1, Dual: cmplx.NaN()}
		}
		if Abs(d) > 1 {
			if d.Dual == 0 {
				return Number{Real: cmplx.Inf(), Dual: cmplx.NaN()}
			}
			return Number{Real: cmplx.Inf(), Dual: cmplx.Inf()}
		}
		return Number{Real: 0, Dual: cmplx.NaN()}
	case math.IsInf(p, -1):
		if d.Real == -1 {
			return Number{Real: 1, Dual: cmplx.NaN()}
		}
		if Abs(d) > 1 {
			return Number{Real: 0, Dual: cmplx.NaN()}
		}
		if d.Dual == 0 {
			return Number{Real: cmplx.Inf(), Dual: cmplx.NaN()}
		}
		return Number{Real: cmplx.Inf(), Dual: cmplx.Inf()}
	case math.IsNaN(p):
		return Number{Real: cmplx.NaN(), Dual: cmplx.NaN()}
	case d.Real == 0:
		if p < 1 {
			return Number{Real: d.Real, Dual: cmplx.Inf()}
		}
		return Number{Real: d.Real}
	case cmplx.IsInf(d.Real):
		if p < 0 {
			return Number{Real: 0, Dual: cmplx.NaN()}
		}
		return Number{Real: cmplx.Inf(), Dual: cmplx.NaN()}
	}
	return Pow(d, Number{Real: complex(p, 0)})
}

// Pow returns d**p, the base-d exponential of p.
func Pow(d, p Number) Number {
	return Exp(Mul(p, Log(d)))
}

// Sqrt returns the square root of d.
//
// Special cases are:
//
//	Sqrt(+Inf) = +Inf
//	Sqrt(±0) = (±0+Infϵ)
//	Sqrt(x < 0) = NaN
//	Sqrt(NaN) = NaN
func Sqrt(d Number) Number {
	return PowReal(d, 0.5)
}

// Exp returns e**q, the base-e exponential of d.
//
// Special cases are:
//
//	Exp(+Inf) = +Inf
//	Exp(NaN) = NaN
//
// Very large values overflow to 0 or +Inf.
// Very small values underflow to 1.
func Exp(d Number) Number {
	fn := cmplx.Exp(d.Real)
	if imag(d.Real) == 0 {
		return Number{Real: fn, Dual: fn * d.Dual}
	}
	conj := cmplx.Conj(d.Real)
	return Number{
		Real: fn,
		Dual: ((fn - cmplx.Exp(conj)) / (d.Real - conj)) * d.Dual,
	}
}

// Log returns the natural logarithm of d.
//
// Special cases are:
//
//	Log(+Inf) = (+Inf+0ϵ)
//	Log(0) = (-Inf±Infϵ)
//	Log(x < 0) = NaN
//	Log(NaN) = NaN
func Log(d Number) Number {
	fn := cmplx.Log(d.Real)
	switch {
	case d.Real == 0:
		return Number{
			Real: fn,
			Dual: complex(math.Copysign(math.Inf(1), real(d.Real)), math.NaN()),
		}
	case imag(d.Real) == 0:
		return Number{
			Real: fn,
			Dual: d.Dual / d.Real,
		}
	case cmplx.IsInf(d.Real):
		return Number{
			Real: fn,
			Dual: 0,
		}
	}
	conj := cmplx.Conj(d.Real)
	return Number{
		Real: fn,
		Dual: ((fn - cmplx.Log(conj)) / (d.Real - conj)) * d.Dual,
	}
}