Source file src/pkg/cmd/compile/internal/gc/mpfloat.go

     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		// TODO(gri) why is this needed?
   183		for len(as) > 0 && (as[0] == ' ' || as[0] == '\t') {
   184			as = as[1:]
   185		}
   186	
   187		f, _, err := a.Val.Parse(as, 0)
   188		if err != nil {
   189			yyerror("malformed constant: %s (%v)", as, err)
   190			a.Val.SetFloat64(0)
   191			return
   192		}
   193	
   194		if f.IsInf() {
   195			yyerror("constant too large: %s", as)
   196			a.Val.SetFloat64(0)
   197			return
   198		}
   199	
   200		// -0 becomes 0
   201		if f.Sign() == 0 && f.Signbit() {
   202			a.Val.SetFloat64(0)
   203		}
   204	}
   205	
   206	func (f *Mpflt) String() string {
   207		return f.Val.Text('b', 0)
   208	}
   209	
   210	func (fvp *Mpflt) GoString() string {
   211		// determine sign
   212		sign := ""
   213		f := &fvp.Val
   214		if f.Sign() < 0 {
   215			sign = "-"
   216			f = new(big.Float).Abs(f)
   217		}
   218	
   219		// Don't try to convert infinities (will not terminate).
   220		if f.IsInf() {
   221			return sign + "Inf"
   222		}
   223	
   224		// Use exact fmt formatting if in float64 range (common case):
   225		// proceed if f doesn't underflow to 0 or overflow to inf.
   226		if x, _ := f.Float64(); f.Sign() == 0 == (x == 0) && !math.IsInf(x, 0) {
   227			return fmt.Sprintf("%s%.6g", sign, x)
   228		}
   229	
   230		// Out of float64 range. Do approximate manual to decimal
   231		// conversion to avoid precise but possibly slow Float
   232		// formatting.
   233		// f = mant * 2**exp
   234		var mant big.Float
   235		exp := f.MantExp(&mant) // 0.5 <= mant < 1.0
   236	
   237		// approximate float64 mantissa m and decimal exponent d
   238		// f ~ m * 10**d
   239		m, _ := mant.Float64()                     // 0.5 <= m < 1.0
   240		d := float64(exp) * (math.Ln2 / math.Ln10) // log_10(2)
   241	
   242		// adjust m for truncated (integer) decimal exponent e
   243		e := int64(d)
   244		m *= math.Pow(10, d-float64(e))
   245	
   246		// ensure 1 <= m < 10
   247		switch {
   248		case m < 1-0.5e-6:
   249			// The %.6g format below rounds m to 5 digits after the
   250			// decimal point. Make sure that m*10 < 10 even after
   251			// rounding up: m*10 + 0.5e-5 < 10 => m < 1 - 0.5e6.
   252			m *= 10
   253			e--
   254		case m >= 10:
   255			m /= 10
   256			e++
   257		}
   258	
   259		return fmt.Sprintf("%s%.6ge%+d", sign, m, e)
   260	}
   261	
   262	// complex multiply v *= rv
   263	//	(a, b) * (c, d) = (a*c - b*d, b*c + a*d)
   264	func (v *Mpcplx) Mul(rv *Mpcplx) {
   265		var ac, ad, bc, bd Mpflt
   266	
   267		ac.Set(&v.Real)
   268		ac.Mul(&rv.Real) // ac
   269	
   270		bd.Set(&v.Imag)
   271		bd.Mul(&rv.Imag) // bd
   272	
   273		bc.Set(&v.Imag)
   274		bc.Mul(&rv.Real) // bc
   275	
   276		ad.Set(&v.Real)
   277		ad.Mul(&rv.Imag) // ad
   278	
   279		v.Real.Set(&ac)
   280		v.Real.Sub(&bd) // ac-bd
   281	
   282		v.Imag.Set(&bc)
   283		v.Imag.Add(&ad) // bc+ad
   284	}
   285	
   286	// complex divide v /= rv
   287	//	(a, b) / (c, d) = ((a*c + b*d), (b*c - a*d))/(c*c + d*d)
   288	func (v *Mpcplx) Div(rv *Mpcplx) bool {
   289		if rv.Real.CmpFloat64(0) == 0 && rv.Imag.CmpFloat64(0) == 0 {
   290			return false
   291		}
   292	
   293		var ac, ad, bc, bd, cc_plus_dd Mpflt
   294	
   295		cc_plus_dd.Set(&rv.Real)
   296		cc_plus_dd.Mul(&rv.Real) // cc
   297	
   298		ac.Set(&rv.Imag)
   299		ac.Mul(&rv.Imag)    // dd
   300		cc_plus_dd.Add(&ac) // cc+dd
   301	
   302		// We already checked that c and d are not both zero, but we can't
   303		// assume that c²+d² != 0 follows, because for tiny values of c
   304		// and/or d c²+d² can underflow to zero.  Check that c²+d² is
   305		// nonzero, return if it's not.
   306		if cc_plus_dd.CmpFloat64(0) == 0 {
   307			return false
   308		}
   309	
   310		ac.Set(&v.Real)
   311		ac.Mul(&rv.Real) // ac
   312	
   313		bd.Set(&v.Imag)
   314		bd.Mul(&rv.Imag) // bd
   315	
   316		bc.Set(&v.Imag)
   317		bc.Mul(&rv.Real) // bc
   318	
   319		ad.Set(&v.Real)
   320		ad.Mul(&rv.Imag) // ad
   321	
   322		v.Real.Set(&ac)
   323		v.Real.Add(&bd)         // ac+bd
   324		v.Real.Quo(&cc_plus_dd) // (ac+bd)/(cc+dd)
   325	
   326		v.Imag.Set(&bc)
   327		v.Imag.Sub(&ad)         // bc-ad
   328		v.Imag.Quo(&cc_plus_dd) // (bc+ad)/(cc+dd)
   329	
   330		return true
   331	}
   332	
   333	func (v *Mpcplx) String() string {
   334		return fmt.Sprintf("(%s+%si)", v.Real.String(), v.Imag.String())
   335	}
   336	
   337	func (v *Mpcplx) GoString() string {
   338		var re string
   339		sre := v.Real.CmpFloat64(0)
   340		if sre != 0 {
   341			re = v.Real.GoString()
   342		}
   343	
   344		var im string
   345		sim := v.Imag.CmpFloat64(0)
   346		if sim != 0 {
   347			im = v.Imag.GoString()
   348		}
   349	
   350		switch {
   351		case sre == 0 && sim == 0:
   352			return "0"
   353		case sre == 0:
   354			return im + "i"
   355		case sim == 0:
   356			return re
   357		case sim < 0:
   358			return fmt.Sprintf("(%s%si)", re, im)
   359		default:
   360			return fmt.Sprintf("(%s+%si)", re, im)
   361		}
   362	}
   363	

View as plain text