vendor/gonum.org/v1/gonum/internal/math32/math_test.go

   1 // Copyright ©2015 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 math32
   6
   7 import (
   8         "math"
   9         "testing"
  10         "testing/quick"
  11
  12         "gonum.org/v1/gonum/floats"
  13 )
  14
  15 const tol = 1e-7
  16
  17 func TestAbs(t *testing.T) {
  18         f := func(x float32) bool {
  19                 y := Abs(x)
  20                 return y == float32(math.Abs(float64(x)))
  21         }
  22         if err := quick.Check(f, nil); err != nil {
  23                 t.Error(err)
  24         }
  25 }
  26
  27 func TestCopySign(t *testing.T) {
  28         f := func(x struct{ X, Y float32 }) bool {
  29                 y := Copysign(x.X, x.Y)
  30                 return y == float32(math.Copysign(float64(x.X), float64(x.Y)))
  31         }
  32         if err := quick.Check(f, nil); err != nil {
  33                 t.Error(err)
  34         }
  35 }
  36
  37 func TestHypot(t *testing.T) {
  38         // tol is increased for Hypot to avoid failures
  39         // related to https://github.com/gonum/gonum/issues/110.
  40         const tol = 1e-6
  41         f := func(x struct{ X, Y float32 }) bool {
  42                 y := Hypot(x.X, x.Y)
  43                 if math.Hypot(float64(x.X), float64(x.Y)) > math.MaxFloat32 {
  44                         return true
  45                 }
  46                 return floats.EqualWithinRel(float64(y), math.Hypot(float64(x.X), float64(x.Y)), tol)
  47         }
  48         if err := quick.Check(f, nil); err != nil {
  49                 t.Error(err)
  50         }
  51 }
  52
  53 func TestInf(t *testing.T) {
  54         if float64(Inf(1)) != math.Inf(1) || float64(Inf(-1)) != math.Inf(-1) {
  55                 t.Error("float32(inf) not infinite")
  56         }
  57 }
  58
  59 func TestIsInf(t *testing.T) {
  60         posInf := float32(math.Inf(1))
  61         negInf := float32(math.Inf(-1))
  62         if !IsInf(posInf, 0) || !IsInf(negInf, 0) || !IsInf(posInf, 1) || !IsInf(negInf, -1) || IsInf(posInf, -1) || IsInf(negInf, 1) {
  63                 t.Error("unexpected isInf value")
  64         }
  65         f := func(x struct {
  66                 F    float32
  67                 Sign int
  68         }) bool {
  69                 y := IsInf(x.F, x.Sign)
  70                 return y == math.IsInf(float64(x.F), x.Sign)
  71         }
  72         if err := quick.Check(f, nil); err != nil {
  73                 t.Error(err)
  74         }
  75 }
  76
  77 func TestIsNaN(t *testing.T) {
  78         f := func(x float32) bool {
  79                 y := IsNaN(x)
  80                 return y == math.IsNaN(float64(x))
  81         }
  82         if err := quick.Check(f, nil); err != nil {
  83                 t.Error(err)
  84         }
  85 }
  86
  87 func TestNaN(t *testing.T) {
  88         if !math.IsNaN(float64(NaN())) {
  89                 t.Errorf("float32(nan) is a number: %f", NaN())
  90         }
  91 }
  92
  93 func TestSignbit(t *testing.T) {
  94         f := func(x float32) bool {
  95                 return Signbit(x) == math.Signbit(float64(x))
  96         }
  97         if err := quick.Check(f, nil); err != nil {
  98                 t.Error(err)
  99         }
 100 }
 101
 102 func TestSqrt(t *testing.T) {
 103         f := func(x float32) bool {
 104                 y := Sqrt(x)
 105                 if IsNaN(y) && IsNaN(sqrt(x)) {
 106                         return true
 107                 }
 108                 return floats.EqualWithinRel(float64(y), float64(sqrt(x)), tol)
 109         }
 110         if err := quick.Check(f, nil); err != nil {
 111                 t.Error(err)
 112         }
 113 }
 114
 115 // Copyright 2009 The Go Authors. All rights reserved.
 116 // Use of this source code is governed by a BSD-style
 117 // license that can be found in the LICENSE file.
 118
 119 // The original C code and the long comment below are
 120 // from FreeBSD's /usr/src/lib/msun/src/e_sqrt.c and
 121 // came with this notice.  The go code is a simplified
 122 // version of the original C.
 123 //
 124 // ====================================================
 125 // Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
 126 //
 127 // Developed at SunPro, a Sun Microsystems, Inc. business.
 128 // Permission to use, copy, modify, and distribute this
 129 // software is freely granted, provided that this notice
 130 // is preserved.
 131 // ====================================================
 132 //
 133 // __ieee754_sqrt(x)
 134 // Return correctly rounded sqrt.
 135 //           -----------------------------------------
 136 //           | Use the hardware sqrt if you have one |
 137 //           -----------------------------------------
 138 // Method:
 139 //   Bit by bit method using integer arithmetic. (Slow, but portable)
 140 //   1. Normalization
 141 //      Scale x to y in [1,4) with even powers of 2:
 142 //      find an integer k such that  1 <= (y=x*2**(2k)) < 4, then
 143 //              sqrt(x) = 2**k * sqrt(y)
 144 //   2. Bit by bit computation
 145 //      Let q  = sqrt(y) truncated to i bit after binary point (q = 1),
 146 //           i                                                   0
 147 //                                     i+1         2
 148 //          s  = 2*q , and      y  =  2   * ( y - q  ).          (1)
 149 //           i      i            i                 i
 150 //
 151 //      To compute q    from q , one checks whether
 152 //                  i+1       i
 153 //
 154 //                            -(i+1) 2
 155 //                      (q + 2      )  <= y.                     (2)
 156 //                        i
 157 //                                                            -(i+1)
 158 //      If (2) is false, then q   = q ; otherwise q   = q  + 2      .
 159 //                             i+1   i             i+1   i
 160 //
 161 //      With some algebraic manipulation, it is not difficult to see
 162 //      that (2) is equivalent to
 163 //                             -(i+1)
 164 //                      s  +  2       <= y                       (3)
 165 //                       i                i
 166 //
 167 //      The advantage of (3) is that s  and y  can be computed by
 168 //                                    i      i
 169 //      the following recurrence formula:
 170 //          if (3) is false
 171 //
 172 //          s     =  s  ,       y    = y   ;                     (4)
 173 //           i+1      i          i+1    i
 174 //
 175 //      otherwise,
 176 //                         -i                      -(i+1)
 177 //          s     =  s  + 2  ,  y    = y  -  s  - 2              (5)
 178 //           i+1      i          i+1    i     i
 179 //
 180 //      One may easily use induction to prove (4) and (5).
 181 //      Note. Since the left hand side of (3) contain only i+2 bits,
 182 //            it does not necessary to do a full (53-bit) comparison
 183 //            in (3).
 184 //   3. Final rounding
 185 //      After generating the 53 bits result, we compute one more bit.
 186 //      Together with the remainder, we can decide whether the
 187 //      result is exact, bigger than 1/2ulp, or less than 1/2ulp
 188 //      (it will never equal to 1/2ulp).
 189 //      The rounding mode can be detected by checking whether
 190 //      huge + tiny is equal to huge, and whether huge - tiny is
 191 //      equal to huge for some floating point number "huge" and "tiny".
 192 //
 193 func sqrt(x float32) float32 {
 194         // special cases
 195         switch {
 196         case x == 0 || IsNaN(x) || IsInf(x, 1):
 197                 return x
 198         case x < 0:
 199                 return NaN()
 200         }
 201         ix := math.Float32bits(x)
 202         // normalize x
 203         exp := int((ix >> shift) & mask)
 204         if exp == 0 { // subnormal x
 205                 for ix&1<<shift == 0 {
 206                         ix <<= 1
 207                         exp--
 208                 }
 209                 exp++
 210         }
 211         exp -= bias // unbias exponent
 212         ix &^= mask << shift
 213         ix |= 1 << shift
 214         if exp&1 == 1 { // odd exp, double x to make it even
 215                 ix <<= 1
 216         }
 217         exp >>= 1 // exp = exp/2, exponent of square root
 218         // generate sqrt(x) bit by bit
 219         ix <<= 1
 220         var q, s uint32               // q = sqrt(x)
 221         r := uint32(1 << (shift + 1)) // r = moving bit from MSB to LSB
 222         for r != 0 {
 223                 t := s + r
 224                 if t <= ix {
 225                         s = t + r
 226                         ix -= t
 227                         q += r
 228                 }
 229                 ix <<= 1
 230                 r >>= 1
 231         }
 232         // final rounding
 233         if ix != 0 { // remainder, result not exact
 234                 q += q & 1 // round according to extra bit
 235         }
 236         ix = q>>1 + uint32(exp-1+bias)<<shift // significand + biased exponent
 237         return math.Float32frombits(ix)
 238 }