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