github.com/jingcheng-WU/gonum@v0.9.1-0.20210323123734-f1a2a11a8f7b/num/hyperdual/hyperdual.go (about)

     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  
     5  package hyperdual
     6  
     7  import (
     8  	"fmt"
     9  	"math"
    10  	"strings"
    11  )
    12  
    13  // Number is a float64 precision hyperdual number.
    14  type Number struct {
    15  	Real, E1mag, E2mag, E1E2mag float64
    16  }
    17  
    18  var negZero = math.Float64frombits(1 << 63)
    19  
    20  // Format implements fmt.Formatter.
    21  func (d Number) Format(fs fmt.State, c rune) {
    22  	prec, pOk := fs.Precision()
    23  	if !pOk {
    24  		prec = -1
    25  	}
    26  	width, wOk := fs.Width()
    27  	if !wOk {
    28  		width = -1
    29  	}
    30  	switch c {
    31  	case 'v':
    32  		if fs.Flag('#') {
    33  			fmt.Fprintf(fs, "%T{Real:%#v, E1mag:%#v, E2mag:%#v, E1E2mag:%#v}", d, d.Real, d.E1mag, d.E2mag, d.E1E2mag)
    34  			return
    35  		}
    36  		if fs.Flag('+') {
    37  			fmt.Fprintf(fs, "{Real:%+v, E1mag:%+v, E2mag:%+v, E1E2mag:%+v}", d.Real, d.E1mag, d.E2mag, d.E1E2mag)
    38  			return
    39  		}
    40  		c = 'g'
    41  		prec = -1
    42  		fallthrough
    43  	case 'e', 'E', 'f', 'F', 'g', 'G':
    44  		fre := fmtString(fs, c, prec, width, false)
    45  		fim := fmtString(fs, c, prec, width, true)
    46  		fmt.Fprintf(fs, fmt.Sprintf("(%s%[2]sϵ₁%[2]sϵ₂%[2]sϵ₁ϵ₂)", fre, fim), d.Real, d.E1mag, d.E2mag, d.E1E2mag)
    47  	default:
    48  		fmt.Fprintf(fs, "%%!%c(%T=%[2]v)", c, d)
    49  		return
    50  	}
    51  }
    52  
    53  // This is horrible, but it's what we have.
    54  func fmtString(fs fmt.State, c rune, prec, width int, wantPlus bool) string {
    55  	var b strings.Builder
    56  	b.WriteByte('%')
    57  	for _, f := range "0+- " {
    58  		if fs.Flag(int(f)) || (f == '+' && wantPlus) {
    59  			b.WriteByte(byte(f))
    60  		}
    61  	}
    62  	if width >= 0 {
    63  		fmt.Fprint(&b, width)
    64  	}
    65  	if prec >= 0 {
    66  		b.WriteByte('.')
    67  		if prec > 0 {
    68  			fmt.Fprint(&b, prec)
    69  		}
    70  	}
    71  	b.WriteRune(c)
    72  	return b.String()
    73  }
    74  
    75  // Add returns the sum of x and y.
    76  func Add(x, y Number) Number {
    77  	return Number{
    78  		Real:    x.Real + y.Real,
    79  		E1mag:   x.E1mag + y.E1mag,
    80  		E2mag:   x.E2mag + y.E2mag,
    81  		E1E2mag: x.E1E2mag + y.E1E2mag,
    82  	}
    83  }
    84  
    85  // Sub returns the difference of x and y, x-y.
    86  func Sub(x, y Number) Number {
    87  	return Number{
    88  		Real:    x.Real - y.Real,
    89  		E1mag:   x.E1mag - y.E1mag,
    90  		E2mag:   x.E2mag - y.E2mag,
    91  		E1E2mag: x.E1E2mag - y.E1E2mag,
    92  	}
    93  }
    94  
    95  // Mul returns the hyperdual product of x and y.
    96  func Mul(x, y Number) Number {
    97  	return Number{
    98  		Real:    x.Real * y.Real,
    99  		E1mag:   x.Real*y.E1mag + x.E1mag*y.Real,
   100  		E2mag:   x.Real*y.E2mag + x.E2mag*y.Real,
   101  		E1E2mag: x.Real*y.E1E2mag + x.E1mag*y.E2mag + x.E2mag*y.E1mag + x.E1E2mag*y.Real,
   102  	}
   103  }
   104  
   105  // Inv returns the hyperdual inverse of d.
   106  //
   107  // Special cases are:
   108  //	Inv(±Inf) = ±0-0ϵ₁-0ϵ₂±0ϵ₁ϵ₂
   109  //	Inv(±0) = ±Inf-Infϵ₁-Infϵ₂±Infϵ₁ϵ₂
   110  func Inv(d Number) Number {
   111  	if d.Real == 0 {
   112  		return Number{
   113  			Real:    1 / d.Real,
   114  			E1mag:   math.Inf(-1),
   115  			E2mag:   math.Inf(-1),
   116  			E1E2mag: 1 / d.Real, // Return a signed inf from a signed zero.
   117  		}
   118  	}
   119  	d2 := d.Real * d.Real
   120  	return Number{
   121  		Real:    1 / d.Real,
   122  		E1mag:   -d.E1mag / d2,
   123  		E2mag:   -d.E2mag / d2,
   124  		E1E2mag: -d.E1E2mag/d2 + 2*d.E1mag*d.E2mag/(d2*d.Real),
   125  	}
   126  }
   127  
   128  // Scale returns d scaled by f.
   129  func Scale(f float64, d Number) Number {
   130  	return Number{Real: f * d.Real, E1mag: f * d.E1mag, E2mag: f * d.E2mag, E1E2mag: f * d.E1E2mag}
   131  }
   132  
   133  // Abs returns the absolute value of d.
   134  func Abs(d Number) Number {
   135  	if math.Float64bits(d.Real)&(1<<63) == 0 {
   136  		return d
   137  	}
   138  	return Scale(-1, d)
   139  }