Merge 613237d into 65e7c60

airy/internal/amos/amos.go

@@ -12,7 +12,6 @@ import (
 
 /*
 The AMOS functions are included in SLATEC, and the SLATEC guide (http://www.netlib.org/slatec/guide) explicitly states:
-
 "The Library is in the public domain and distributed by the Energy
 Science and Technology Software Center."
 Mention of AMOS's inclusion in SLATEC goes back at least to this 1985 technical report from Sandia National Labs: http://infoserve.sandia.gov/sand_doc/1985/851018.pdf
@@ -554,6 +553,8 @@ func Zbknu(ZR, ZI, FNU float64, KODE, N int, YR, YI []float64, NZ int, TOL, ELIM
 
 	var I, IFLAG, INU, K, KFLAG, KK, KMAX, KODED, IDUM, J, IC, INUB, NW int
 
+	var sinh, cosh complex128
+
 	var tmp complex128
 	var CSSR, CSRR, BRY [4]float64
 	var CYR, CYI [3]float64
@@ -618,7 +619,11 @@ func Zbknu(ZR, ZI, FNU float64, KODE, N int, YR, YI []float64, NZ int, TOL, ELIM
 	SMUI = imag(tmp)
 	FMUR = SMUR * DNU
 	FMUI = SMUI * DNU
-	FMUR, FMUI, CSHR, CSHI, CCHR, CCHI = Zshch(FMUR, FMUI, CSHR, CSHI, CCHR, CCHI)
+	sinh, cosh = Zshch(complex(FMUR, FMUI))
+	CSHR = real(sinh)
+	CSHI = imag(sinh)
+	CCHR = real(cosh)
+	CCHI = imag(cosh)
 	if DNU == 0.0E0 {
 		goto Ten
 	}
@@ -2200,20 +2205,9 @@ Ten:
 	return ZRR, ZRI, S1R, S1I, S2R, S2I, NZ, ASCLE, ALIM, IUF
 }
 
-// ZSHCH COMPUTES THE COMPLEX HYPERBOLIC FUNCTIONS CSH=SINH(X+iY) AND
-// CCH=COSH(X+I*Y), WHERE I**2=-1.
-// TODO(btracey): use cmplx.Sinh and cmplx.Cosh.
-func Zshch(ZR, ZI, CSHR, CSHI, CCHR, CCHI float64) (ZRout, ZIout, CSHRout, CSHIout, CCHRout, CCHIout float64) {
-	var CH, CN, SH, SN float64
-	SH = math.Sinh(ZR)
-	CH = math.Cosh(ZR)
-	SN = dsin(ZI)
-	CN = dcos(ZI)
-	CSHR = SH * CN
-	CSHI = CH * SN
-	CCHR = CH * CN
-	CCHI = SH * SN
-	return ZR, ZI, CSHR, CSHI, CCHR, CCHI
+// Zshch computes the hyperbolic sin and cosine of the input z.
+func Zshch(z complex128) (sinh, cosh complex128) {
+	return cmplx.Sinh(z), cmplx.Cosh(z)
 }
 
 func dmax(a, b float64) float64 {

airy/internal/amos/amos_test.go

@@ -9,6 +9,8 @@ import (
 	"math/rand"
 	"strconv"
 	"testing"
+
+	"github.com/gonum/floats"
 )
 
 type input struct {
@@ -370,8 +372,8 @@ func zbknutest(t *testing.T, x []float64, is []int, tol float64, n int, yr, yi [
 	sameF64(t, "zbknu elim", ELIMfort, ELIMamos)
 	sameF64(t, "zbknu alim", ALIMfort, ALIMamos)
 
-	sameF64S(t, "zbknu yr", YRfort, YRamos)
-	sameF64S(t, "zbknu yi", YIfort, YIamos)
+	sameF64SApprox(t, "zbknu yr", YRfort, YRamos, 1e-12)
+	sameF64SApprox(t, "zbknu yi", YIfort, YIamos, 1e-12)
 }
 
 func zairytest(t *testing.T, x []float64, kode, id int) {
@@ -383,8 +385,8 @@ func zairytest(t *testing.T, x []float64, kode, id int) {
 	AIRfort, AIIfort, NZfort := zairyOrig(ZR, ZI, ID, KODE)
 	AIRamos, AIIamos, NZamos := Zairy(ZR, ZI, ID, KODE)
 
-	sameF64(t, "zairy air", AIRfort, AIRamos)
-	sameF64(t, "zairy aii", AIIfort, AIIamos)
+	sameF64Approx(t, "zairy air", AIRfort, AIRamos, 1e-12)
+	sameF64Approx(t, "zairy aii", AIIfort, AIIamos, 1e-12)
 	sameInt(t, "zairy nz", NZfort, NZamos)
 }
 
@@ -425,8 +427,8 @@ func zacaitest(t *testing.T, x []float64, is []int, tol float64, n int, yr, yi [
 	sameF64(t, "zacai elim", ELIMfort, ELIMamos)
 	sameF64(t, "zacai elim", ALIMfort, ALIMamos)
 
-	sameF64S(t, "zacai yr", YRfort, YRamos)
-	sameF64S(t, "zacai yi", YIfort, YIamos)
+	sameF64SApprox(t, "zacai yr", YRfort, YRamos, 1e-12)
+	sameF64SApprox(t, "zacai yi", YIfort, YIamos, 1e-12)
 }
 
 func sameF64(t *testing.T, str string, c, native float64) {
@@ -441,6 +443,24 @@ func sameF64(t *testing.T, str string, c, native float64) {
 	t.Errorf("Case %s: Float64 mismatch. c = %v, native = %v\n cb: %v, nb: %v\n", str, c, native, cb, nb)
 }
 
+func sameF64Approx(t *testing.T, str string, c, native, tol float64) {
+	if math.IsNaN(c) && math.IsNaN(native) {
+		return
+	}
+	if floats.EqualWithinAbsOrRel(c, native, tol, tol) {
+		return
+	}
+	// Have a much looser tolerance for correctness when the values are large.
+	// Floating point noise makes the relative tolerance difference greater for
+	// higher values.
+	if c > 1e200 && floats.EqualWithinAbsOrRel(c, native, 10, 10) {
+		return
+	}
+	cb := math.Float64bits(c)
+	nb := math.Float64bits(native)
+	t.Errorf("Case %s: Float64 mismatch. c = %v, native = %v\n cb: %v, nb: %v\n", str, c, native, cb, nb)
+}
+
 func sameInt(t *testing.T, str string, c, native int) {
 	if c != native {
 		t.Errorf("Case %s: Int mismatch. c = %v, native = %v.", str, c, native)
@@ -455,3 +475,12 @@ func sameF64S(t *testing.T, str string, c, native []float64) {
 		sameF64(t, str+"_idx_"+strconv.Itoa(i), v, native[i])
 	}
 }
+
+func sameF64SApprox(t *testing.T, str string, c, native []float64, tol float64) {
+	if len(c) != len(native) {
+		panic(str)
+	}
+	for i, v := range c {
+		sameF64Approx(t, str+"_idx_"+strconv.Itoa(i), v, native[i], tol)
+	}
+}