github.com/gagliardetto/golang-go@v0.0.0-20201020153340-53909ea70814/cmd/compile/internal/gc/mpfloat.go (about)

     1  // Copyright 2009 The Go 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 gc
     6  
     7  import (
     8  	"fmt"
     9  	"math"
    10  	"math/big"
    11  )
    12  
    13  // implements float arithmetic
    14  
    15  const (
    16  	// Maximum size in bits for Mpints before signalling
    17  	// overflow and also mantissa precision for Mpflts.
    18  	Mpprec = 512
    19  	// Turn on for constant arithmetic debugging output.
    20  	Mpdebug = false
    21  )
    22  
    23  // Mpflt represents a floating-point constant.
    24  type Mpflt struct {
    25  	Val big.Float
    26  }
    27  
    28  // Mpcplx represents a complex constant.
    29  type Mpcplx struct {
    30  	Real Mpflt
    31  	Imag Mpflt
    32  }
    33  
    34  // Use newMpflt (not new(Mpflt)!) to get the correct default precision.
    35  func newMpflt() *Mpflt {
    36  	var a Mpflt
    37  	a.Val.SetPrec(Mpprec)
    38  	return &a
    39  }
    40  
    41  // Use newMpcmplx (not new(Mpcplx)!) to get the correct default precision.
    42  func newMpcmplx() *Mpcplx {
    43  	var a Mpcplx
    44  	a.Real = *newMpflt()
    45  	a.Imag = *newMpflt()
    46  	return &a
    47  }
    48  
    49  func (a *Mpflt) SetInt(b *Mpint) {
    50  	if b.checkOverflow(0) {
    51  		// sign doesn't really matter but copy anyway
    52  		a.Val.SetInf(b.Val.Sign() < 0)
    53  		return
    54  	}
    55  	a.Val.SetInt(&b.Val)
    56  }
    57  
    58  func (a *Mpflt) Set(b *Mpflt) {
    59  	a.Val.Set(&b.Val)
    60  }
    61  
    62  func (a *Mpflt) Add(b *Mpflt) {
    63  	if Mpdebug {
    64  		fmt.Printf("\n%v + %v", a, b)
    65  	}
    66  
    67  	a.Val.Add(&a.Val, &b.Val)
    68  
    69  	if Mpdebug {
    70  		fmt.Printf(" = %v\n\n", a)
    71  	}
    72  }
    73  
    74  func (a *Mpflt) AddFloat64(c float64) {
    75  	var b Mpflt
    76  
    77  	b.SetFloat64(c)
    78  	a.Add(&b)
    79  }
    80  
    81  func (a *Mpflt) Sub(b *Mpflt) {
    82  	if Mpdebug {
    83  		fmt.Printf("\n%v - %v", a, b)
    84  	}
    85  
    86  	a.Val.Sub(&a.Val, &b.Val)
    87  
    88  	if Mpdebug {
    89  		fmt.Printf(" = %v\n\n", a)
    90  	}
    91  }
    92  
    93  func (a *Mpflt) Mul(b *Mpflt) {
    94  	if Mpdebug {
    95  		fmt.Printf("%v\n * %v\n", a, b)
    96  	}
    97  
    98  	a.Val.Mul(&a.Val, &b.Val)
    99  
   100  	if Mpdebug {
   101  		fmt.Printf(" = %v\n\n", a)
   102  	}
   103  }
   104  
   105  func (a *Mpflt) MulFloat64(c float64) {
   106  	var b Mpflt
   107  
   108  	b.SetFloat64(c)
   109  	a.Mul(&b)
   110  }
   111  
   112  func (a *Mpflt) Quo(b *Mpflt) {
   113  	if Mpdebug {
   114  		fmt.Printf("%v\n / %v\n", a, b)
   115  	}
   116  
   117  	a.Val.Quo(&a.Val, &b.Val)
   118  
   119  	if Mpdebug {
   120  		fmt.Printf(" = %v\n\n", a)
   121  	}
   122  }
   123  
   124  func (a *Mpflt) Cmp(b *Mpflt) int {
   125  	return a.Val.Cmp(&b.Val)
   126  }
   127  
   128  func (a *Mpflt) CmpFloat64(c float64) int {
   129  	if c == 0 {
   130  		return a.Val.Sign() // common case shortcut
   131  	}
   132  	return a.Val.Cmp(big.NewFloat(c))
   133  }
   134  
   135  func (a *Mpflt) Float64() float64 {
   136  	x, _ := a.Val.Float64()
   137  
   138  	// check for overflow
   139  	if math.IsInf(x, 0) && nsavederrors+nerrors == 0 {
   140  		Fatalf("ovf in Mpflt Float64")
   141  	}
   142  
   143  	return x + 0 // avoid -0 (should not be needed, but be conservative)
   144  }
   145  
   146  func (a *Mpflt) Float32() float64 {
   147  	x32, _ := a.Val.Float32()
   148  	x := float64(x32)
   149  
   150  	// check for overflow
   151  	if math.IsInf(x, 0) && nsavederrors+nerrors == 0 {
   152  		Fatalf("ovf in Mpflt Float32")
   153  	}
   154  
   155  	return x + 0 // avoid -0 (should not be needed, but be conservative)
   156  }
   157  
   158  func (a *Mpflt) SetFloat64(c float64) {
   159  	if Mpdebug {
   160  		fmt.Printf("\nconst %g", c)
   161  	}
   162  
   163  	// convert -0 to 0
   164  	if c == 0 {
   165  		c = 0
   166  	}
   167  	a.Val.SetFloat64(c)
   168  
   169  	if Mpdebug {
   170  		fmt.Printf(" = %v\n", a)
   171  	}
   172  }
   173  
   174  func (a *Mpflt) Neg() {
   175  	// avoid -0
   176  	if a.Val.Sign() != 0 {
   177  		a.Val.Neg(&a.Val)
   178  	}
   179  }
   180  
   181  func (a *Mpflt) SetString(as string) {
   182  	f, _, err := a.Val.Parse(as, 0)
   183  	if err != nil {
   184  		yyerror("malformed constant: %s (%v)", as, err)
   185  		a.Val.SetFloat64(0)
   186  		return
   187  	}
   188  
   189  	if f.IsInf() {
   190  		yyerror("constant too large: %s", as)
   191  		a.Val.SetFloat64(0)
   192  		return
   193  	}
   194  
   195  	// -0 becomes 0
   196  	if f.Sign() == 0 && f.Signbit() {
   197  		a.Val.SetFloat64(0)
   198  	}
   199  }
   200  
   201  func (f *Mpflt) String() string {
   202  	return f.Val.Text('b', 0)
   203  }
   204  
   205  func (fvp *Mpflt) GoString() string {
   206  	// determine sign
   207  	sign := ""
   208  	f := &fvp.Val
   209  	if f.Sign() < 0 {
   210  		sign = "-"
   211  		f = new(big.Float).Abs(f)
   212  	}
   213  
   214  	// Don't try to convert infinities (will not terminate).
   215  	if f.IsInf() {
   216  		return sign + "Inf"
   217  	}
   218  
   219  	// Use exact fmt formatting if in float64 range (common case):
   220  	// proceed if f doesn't underflow to 0 or overflow to inf.
   221  	if x, _ := f.Float64(); f.Sign() == 0 == (x == 0) && !math.IsInf(x, 0) {
   222  		return fmt.Sprintf("%s%.6g", sign, x)
   223  	}
   224  
   225  	// Out of float64 range. Do approximate manual to decimal
   226  	// conversion to avoid precise but possibly slow Float
   227  	// formatting.
   228  	// f = mant * 2**exp
   229  	var mant big.Float
   230  	exp := f.MantExp(&mant) // 0.5 <= mant < 1.0
   231  
   232  	// approximate float64 mantissa m and decimal exponent d
   233  	// f ~ m * 10**d
   234  	m, _ := mant.Float64()                     // 0.5 <= m < 1.0
   235  	d := float64(exp) * (math.Ln2 / math.Ln10) // log_10(2)
   236  
   237  	// adjust m for truncated (integer) decimal exponent e
   238  	e := int64(d)
   239  	m *= math.Pow(10, d-float64(e))
   240  
   241  	// ensure 1 <= m < 10
   242  	switch {
   243  	case m < 1-0.5e-6:
   244  		// The %.6g format below rounds m to 5 digits after the
   245  		// decimal point. Make sure that m*10 < 10 even after
   246  		// rounding up: m*10 + 0.5e-5 < 10 => m < 1 - 0.5e6.
   247  		m *= 10
   248  		e--
   249  	case m >= 10:
   250  		m /= 10
   251  		e++
   252  	}
   253  
   254  	return fmt.Sprintf("%s%.6ge%+d", sign, m, e)
   255  }
   256  
   257  // complex multiply v *= rv
   258  //	(a, b) * (c, d) = (a*c - b*d, b*c + a*d)
   259  func (v *Mpcplx) Mul(rv *Mpcplx) {
   260  	var ac, ad, bc, bd Mpflt
   261  
   262  	ac.Set(&v.Real)
   263  	ac.Mul(&rv.Real) // ac
   264  
   265  	bd.Set(&v.Imag)
   266  	bd.Mul(&rv.Imag) // bd
   267  
   268  	bc.Set(&v.Imag)
   269  	bc.Mul(&rv.Real) // bc
   270  
   271  	ad.Set(&v.Real)
   272  	ad.Mul(&rv.Imag) // ad
   273  
   274  	v.Real.Set(&ac)
   275  	v.Real.Sub(&bd) // ac-bd
   276  
   277  	v.Imag.Set(&bc)
   278  	v.Imag.Add(&ad) // bc+ad
   279  }
   280  
   281  // complex divide v /= rv
   282  //	(a, b) / (c, d) = ((a*c + b*d), (b*c - a*d))/(c*c + d*d)
   283  func (v *Mpcplx) Div(rv *Mpcplx) bool {
   284  	if rv.Real.CmpFloat64(0) == 0 && rv.Imag.CmpFloat64(0) == 0 {
   285  		return false
   286  	}
   287  
   288  	var ac, ad, bc, bd, cc_plus_dd Mpflt
   289  
   290  	cc_plus_dd.Set(&rv.Real)
   291  	cc_plus_dd.Mul(&rv.Real) // cc
   292  
   293  	ac.Set(&rv.Imag)
   294  	ac.Mul(&rv.Imag)    // dd
   295  	cc_plus_dd.Add(&ac) // cc+dd
   296  
   297  	// We already checked that c and d are not both zero, but we can't
   298  	// assume that c²+d² != 0 follows, because for tiny values of c
   299  	// and/or d c²+d² can underflow to zero.  Check that c²+d² is
   300  	// nonzero, return if it's not.
   301  	if cc_plus_dd.CmpFloat64(0) == 0 {
   302  		return false
   303  	}
   304  
   305  	ac.Set(&v.Real)
   306  	ac.Mul(&rv.Real) // ac
   307  
   308  	bd.Set(&v.Imag)
   309  	bd.Mul(&rv.Imag) // bd
   310  
   311  	bc.Set(&v.Imag)
   312  	bc.Mul(&rv.Real) // bc
   313  
   314  	ad.Set(&v.Real)
   315  	ad.Mul(&rv.Imag) // ad
   316  
   317  	v.Real.Set(&ac)
   318  	v.Real.Add(&bd)         // ac+bd
   319  	v.Real.Quo(&cc_plus_dd) // (ac+bd)/(cc+dd)
   320  
   321  	v.Imag.Set(&bc)
   322  	v.Imag.Sub(&ad)         // bc-ad
   323  	v.Imag.Quo(&cc_plus_dd) // (bc+ad)/(cc+dd)
   324  
   325  	return true
   326  }
   327  
   328  func (v *Mpcplx) String() string {
   329  	return fmt.Sprintf("(%s+%si)", v.Real.String(), v.Imag.String())
   330  }
   331  
   332  func (v *Mpcplx) GoString() string {
   333  	var re string
   334  	sre := v.Real.CmpFloat64(0)
   335  	if sre != 0 {
   336  		re = v.Real.GoString()
   337  	}
   338  
   339  	var im string
   340  	sim := v.Imag.CmpFloat64(0)
   341  	if sim != 0 {
   342  		im = v.Imag.GoString()
   343  	}
   344  
   345  	switch {
   346  	case sre == 0 && sim == 0:
   347  		return "0"
   348  	case sre == 0:
   349  		return im + "i"
   350  	case sim == 0:
   351  		return re
   352  	case sim < 0:
   353  		return fmt.Sprintf("(%s%si)", re, im)
   354  	default:
   355  		return fmt.Sprintf("(%s+%si)", re, im)
   356  	}
   357  }