math: add functions; update tests and special cases
Added special cases to comments for asin.go and fabs.go.
Added Trunc() to floor.go and floor_386.s. Fixed formatting
error in hypot_386.s Added new functions Acosh, Asinh,
Atanh, Copysign, Erf, Erfc, Expm1, and Log1p. Added
386 FPU version of Fmod. Added tests, benchmarks, and
precision to expected results in all_test.go. Edited
makefile so it all compiles.
R=rsc
CC=golang-dev
https://golang.org/cl/195052
diff --git a/src/pkg/math/Makefile b/src/pkg/math/Makefile
index b10df65..f228084 100644
--- a/src/pkg/math/Makefile
+++ b/src/pkg/math/Makefile
@@ -15,6 +15,7 @@
exp_386.$O\
fabs_386.$O\
floor_386.$O\
+ fmod_386.$O\
hypot_386.$O\
log_386.$O\
sin_386.$O\
@@ -25,17 +26,24 @@
$(OFILES_$(GOARCH))
ALLGOFILES=\
+ acosh.go\
asin.go\
+ asinh.go\
atan.go\
+ atanh.go\
atan2.go\
bits.go\
const.go\
+ copysign.go\
+ erf.go\
exp.go\
+ expm1.go\
fabs.go\
floor.go\
fmod.go\
hypot.go\
log.go\
+ log1p.go\
pow.go\
pow10.go\
sin.go\
diff --git a/src/pkg/math/acosh.go b/src/pkg/math/acosh.go
new file mode 100644
index 0000000..1f0d3f3
--- /dev/null
+++ b/src/pkg/math/acosh.go
@@ -0,0 +1,59 @@
+// Copyright 2010 The Go Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+package math
+
+
+// The original C code, the long comment, and the constants
+// below are from FreeBSD's /usr/src/lib/msun/src/e_acosh.c
+// and came with this notice. The go code is a simplified
+// version of the original C.
+//
+// ====================================================
+// Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+//
+// Developed at SunPro, a Sun Microsystems, Inc. business.
+// Permission to use, copy, modify, and distribute this
+// software is freely granted, provided that this notice
+// is preserved.
+// ====================================================
+//
+//
+// __ieee754_acosh(x)
+// Method :
+// Based on
+// acosh(x) = log [ x + sqrt(x*x-1) ]
+// we have
+// acosh(x) := log(x)+ln2, if x is large; else
+// acosh(x) := log(2x-1/(sqrt(x*x-1)+x)) if x>2; else
+// acosh(x) := log1p(t+sqrt(2.0*t+t*t)); where t=x-1.
+//
+// Special cases:
+// acosh(x) is NaN with signal if x<1.
+// acosh(NaN) is NaN without signal.
+//
+
+// Acosh(x) calculates the inverse hyperbolic cosine of x.
+//
+// Special cases are:
+// Acosh(x) = NaN if x < 1
+// Acosh(NaN) = NaN
+func Acosh(x float64) float64 {
+ const (
+ Ln2 = 6.93147180559945286227e-01 // 0x3FE62E42FEFA39EF
+ Large = 1 << 28 // 2^28
+ )
+ switch {
+ case x < 1 || IsNaN(x):
+ return NaN()
+ case x == 1:
+ return 0
+ case x >= Large:
+ return Log(x) + Ln2 // x > 2^28
+ case x > 2:
+ return Log(2*x - 1/(x+Sqrt(x*x-1))) // 2^28 > x > 2
+ }
+ t := x - 1
+ return Log1p(t + Sqrt(2*t+t*t)) // 2 >= x > 1
+}
diff --git a/src/pkg/math/all_test.go b/src/pkg/math/all_test.go
index 9490e06..59fdd9e 100644
--- a/src/pkg/math/all_test.go
+++ b/src/pkg/math/all_test.go
@@ -22,41 +22,83 @@
1.8253080916808550e+00,
-8.6859247685756013e+00,
}
+// The expected results below were computed by the high precision calculators
+// at http://keisan.casio.com/. More exact input values (array vf[], above)
+// were obtained by printing them with "%.26f". The answers were calculated
+// to 26 digits (by using the "Digit number" drop-down control of each
+// calculator). Twenty-six digits were chosen so that the answers would be
+// accurate even for a float128 type.
var acos = []float64{
- 1.0496193546107222e+00,
- 6.858401281366443e-01,
- 1.598487871457716e+00,
- 2.095619936147586e+00,
- 2.7053008467824158e-01,
- 1.2738121680361776e+00,
- 1.0205369421140630e+00,
- 1.2945003481781246e+00,
- 1.3872364345374451e+00,
- 2.6231510803970464e+00,
+ 1.0496193546107222142571536e+00,
+ 6.8584012813664425171660692e-01,
+ 1.5984878714577160325521819e+00,
+ 2.0956199361475859327461799e+00,
+ 2.7053008467824138592616927e-01,
+ 1.2738121680361776018155625e+00,
+ 1.0205369421140629186287407e+00,
+ 1.2945003481781246062157835e+00,
+ 1.3872364345374451433846657e+00,
+ 2.6231510803970463967294145e+00,
+}
+var acosh = []float64{
+ 2.4743347004159012494457618e+00,
+ 2.8576385344292769649802701e+00,
+ 7.2796961502981066190593175e-01,
+ 2.4796794418831451156471977e+00,
+ 3.0552020742306061857212962e+00,
+ 2.044238592688586588942468e+00,
+ 2.5158701513104513595766636e+00,
+ 1.99050839282411638174299e+00,
+ 1.6988625798424034227205445e+00,
+ 2.9611454842470387925531875e+00,
}
var asin = []float64{
- 5.2117697218417440e-01,
- 8.8495619865825236e-01,
- -2.7691544662819413e-02,
- -5.2482360935268932e-01,
- 1.3002662421166553e+00,
- 2.9698415875871901e-01,
- 5.5025938468083364e-01,
- 2.7629597861677200e-01,
- 1.8355989225745148e-01,
- -1.0523547536021498e+00,
+ 5.2117697218417440497416805e-01,
+ 8.8495619865825236751471477e-01,
+ -02.769154466281941332086016e-02,
+ -5.2482360935268931351485822e-01,
+ 1.3002662421166552333051524e+00,
+ 2.9698415875871901741575922e-01,
+ 5.5025938468083370060258102e-01,
+ 2.7629597861677201301553823e-01,
+ 1.83559892257451475846656e-01,
+ -1.0523547536021497774980928e+00,
+}
+var asinh = []float64{
+ 2.3083139124923523427628243e+00,
+ 2.743551594301593620039021e+00,
+ -2.7345908534880091229413487e-01,
+ -2.3145157644718338650499085e+00,
+ 2.9613652154015058521951083e+00,
+ 1.7949041616585821933067568e+00,
+ 2.3564032905983506405561554e+00,
+ 1.7287118790768438878045346e+00,
+ 1.3626658083714826013073193e+00,
+ -2.8581483626513914445234004e+00,
}
var atan = []float64{
- 1.3725902621296217e+00,
- 1.4422906096452980e+00,
- -2.7011324359471755e-01,
- -1.3738077684543379e+00,
- 1.4673921193587666e+00,
- 1.2415173565870167e+00,
- 1.3818396865615167e+00,
- 1.2194305844639670e+00,
- 1.0696031952318783e+00,
- -1.4561721938838085e+00,
+ 1.372590262129621651920085e+00,
+ 1.442290609645298083020664e+00,
+ -2.7011324359471758245192595e-01,
+ -1.3738077684543379452781531e+00,
+ 1.4673921193587666049154681e+00,
+ 1.2415173565870168649117764e+00,
+ 1.3818396865615168979966498e+00,
+ 1.2194305844639670701091426e+00,
+ 1.0696031952318783760193244e+00,
+ -1.4561721938838084990898679e+00,
+}
+var atanh = []float64{
+ 5.4651163712251938116878204e-01,
+ 1.0299474112843111224914709e+00,
+ -2.7695084420740135145234906e-02,
+ -5.5072096119207195480202529e-01,
+ 1.9943940993171843235906642e+00,
+ 3.01448604578089708203017e-01,
+ 5.8033427206942188834370595e-01,
+ 2.7987997499441511013958297e-01,
+ 1.8459947964298794318714228e-01,
+ -1.3273186910532645867272502e+00,
}
var ceil = []float64{
5.0000000000000000e+00,
@@ -70,17 +112,89 @@
2.0000000000000000e+00,
-8.0000000000000000e+00,
}
+var copysign = []float64{
+ -4.9790119248836735e+00,
+ -7.7388724745781045e+00,
+ -2.7688005719200159e-01,
+ -5.0106036182710749e+00,
+ -9.6362937071984173e+00,
+ -2.9263772392439646e+00,
+ -5.2290834314593066e+00,
+ -2.7279399104360102e+00,
+ -1.8253080916808550e+00,
+ -8.6859247685756013e+00,
+}
+var cos = []float64{
+ 2.634752140995199110787593e-01,
+ 1.148551260848219865642039e-01,
+ 9.6191297325640768154550453e-01,
+ 2.938141150061714816890637e-01,
+ -9.777138189897924126294461e-01,
+ -9.7693041344303219127199518e-01,
+ 4.940088096948647263961162e-01,
+ -9.1565869021018925545016502e-01,
+ -2.517729313893103197176091e-01,
+ -7.39241351595676573201918e-01,
+}
+var cosh = []float64{
+ 7.2668796942212842775517446e+01,
+ 1.1479413465659254502011135e+03,
+ 1.0385767908766418550935495e+00,
+ 7.5000957789658051428857788e+01,
+ 7.655246669605357888468613e+03,
+ 9.3567491758321272072888257e+00,
+ 9.331351599270605471131735e+01,
+ 7.6833430994624643209296404e+00,
+ 3.1829371625150718153881164e+00,
+ 2.9595059261916188501640911e+03,
+}
+var erf = []float64{
+ 5.1865354817738701906913566e-01,
+ 7.2623875834137295116929844e-01,
+ -3.123458688281309990629839e-02,
+ -5.2143121110253302920437013e-01,
+ 8.2704742671312902508629582e-01,
+ 3.2101767558376376743993945e-01,
+ 5.403990312223245516066252e-01,
+ 3.0034702916738588551174831e-01,
+ 2.0369924417882241241559589e-01,
+ -7.8069386968009226729944677e-01,
+}
+var erfc = []float64{
+ 4.8134645182261298093086434e-01,
+ 2.7376124165862704883070156e-01,
+ 1.0312345868828130999062984e+00,
+ 1.5214312111025330292043701e+00,
+ 1.7295257328687097491370418e-01,
+ 6.7898232441623623256006055e-01,
+ 4.596009687776754483933748e-01,
+ 6.9965297083261411448825169e-01,
+ 7.9630075582117758758440411e-01,
+ 1.7806938696800922672994468e+00,
+}
var exp = []float64{
- 1.4533071302642137e+02,
- 2.2958822575694450e+03,
- 7.5814542574851666e-01,
- 6.6668778421791010e-03,
- 1.5310493273896035e+04,
- 1.8659907517999329e+01,
- 1.8662167355098713e+02,
- 1.5301332413189379e+01,
- 6.2047063430646876e+00,
- 1.6894712385826522e-04,
+ 1.4533071302642137507696589e+02,
+ 2.2958822575694449002537581e+03,
+ 7.5814542574851666582042306e-01,
+ 6.6668778421791005061482264e-03,
+ 1.5310493273896033740861206e+04,
+ 1.8659907517999328638667732e+01,
+ 1.8662167355098714543942057e+02,
+ 1.5301332413189378961665788e+01,
+ 6.2047063430646876349125085e+00,
+ 1.6894712385826521111610438e-04,
+}
+var expm1 = []float64{
+ 5.105047796122957327384770212e-02,
+ 8.046199708567344080562675439e-02,
+ -2.764970978891639815187418703e-03,
+ -4.8871434888875355394330300273e-02,
+ 1.0115864277221467777117227494e-01,
+ 2.969616407795910726014621657e-02,
+ 5.368214487944892300914037972e-02,
+ 2.765488851131274068067445335e-02,
+ 1.842068661871398836913874273e-02,
+ -8.3193870863553801814961137573e-02,
}
var floor = []float64{
4.0000000000000000e+00,
@@ -95,115 +209,150 @@
-9.0000000000000000e+00,
}
var fmod = []float64{
- 4.1976150232653000e-02,
- 2.2611275254218955e+00,
- 3.2317941087942760e-02,
- 4.9893963817289251e+00,
- 3.6370629280158270e-01,
- 1.2208682822681062e+00,
- 4.7709165685406934e+00,
- 1.8161802686919694e+00,
- 8.7345954159572500e-01,
- 1.3140752314243987e+00,
+ 4.197615023265299782906368e-02,
+ 2.261127525421895434476482e+00,
+ 3.231794108794261433104108e-02,
+ 4.989396381728925078391512e+00,
+ 3.637062928015826201999516e-01,
+ 1.220868282268106064236690e+00,
+ 4.770916568540693347699744e+00,
+ 1.816180268691969246219742e+00,
+ 8.734595415957246977711748e-01,
+ 1.314075231424398637614104e+00,
}
var log = []float64{
- 1.6052314626930630e+00,
- 2.0462560018708768e+00,
- -1.2841708730962657e+00,
- 1.6115563905281544e+00,
- 2.2655365644872018e+00,
- 1.0737652208918380e+00,
- 1.6542360106073545e+00,
- 1.0035467127723465e+00,
- 6.0174879014578053e-01,
- 2.1617038728473527e+00,
+ 1.605231462693062999102599e+00,
+ 2.0462560018708770653153909e+00,
+ -1.2841708730962657801275038e+00,
+ 1.6115563905281545116286206e+00,
+ 2.2655365644872016636317461e+00,
+ 1.0737652208918379856272735e+00,
+ 1.6542360106073546632707956e+00,
+ 1.0035467127723465801264487e+00,
+ 6.0174879014578057187016475e-01,
+ 2.161703872847352815363655e+00,
}
var log10 = []float64{
- 6.9714316642508291e-01,
- 8.8867769017393205e-01,
- -5.5770832400658930e-01,
- 6.9989004768229943e-01,
- 9.8391002850684232e-01,
- 4.6633031029295153e-01,
- 7.1842557117242328e-01,
- 4.3583479968917772e-01,
- 2.6133617905227037e-01,
- 9.3881606348649405e-01,
+ 6.9714316642508290997617083e-01,
+ 8.886776901739320576279124e-01,
+ -5.5770832400658929815908236e-01,
+ 6.998900476822994346229723e-01,
+ 9.8391002850684232013281033e-01,
+ 4.6633031029295153334285302e-01,
+ 7.1842557117242328821552533e-01,
+ 4.3583479968917773161304553e-01,
+ 2.6133617905227038228626834e-01,
+ 9.3881606348649405716214241e-01,
+}
+var log1p = []float64{
+ 4.8590257759797794104158205e-02,
+ 7.4540265965225865330849141e-02,
+ -2.7726407903942672823234024e-03,
+ -5.1404917651627649094953380e-02,
+ 9.1998280672258624681335010e-02,
+ 2.8843762576593352865894824e-02,
+ 5.0969534581863707268992645e-02,
+ 2.6913947602193238458458594e-02,
+ 1.8088493239630770262045333e-02,
+ -9.0865245631588989681559268e-02,
}
var pow = []float64{
- 9.5282232631648415e+04,
- 5.4811599352999900e+07,
- 5.2859121715894400e-01,
- 9.7587991957286472e-06,
- 4.3280643293460450e+09,
- 8.4406761805034551e+02,
- 1.6946633276191194e+05,
- 5.3449040147551940e+02,
- 6.6881821384514159e+01,
- 2.0609869004248744e-09,
+ 9.5282232631648411840742957e+04,
+ 5.4811599352999901232411871e+07,
+ 5.2859121715894396531132279e-01,
+ 9.7587991957286474464259698e-06,
+ 4.328064329346044846740467e+09,
+ 8.4406761805034547437659092e+02,
+ 1.6946633276191194947742146e+05,
+ 5.3449040147551939075312879e+02,
+ 6.688182138451414936380374e+01,
+ 2.0609869004248742886827439e-09,
}
var sin = []float64{
- -9.6466616586009283e-01,
- 9.9338225271646543e-01,
- -2.7335587039794395e-01,
- 9.5586257685042800e-01,
- -2.0994210667799692e-01,
- 2.1355787807998605e-01,
- -8.6945689711673619e-01,
- 4.0195666811555783e-01,
- 9.6778633541688000e-01,
- -6.7344058690503452e-01,
+ -9.6466616586009283766724726e-01,
+ 9.9338225271646545763467022e-01,
+ -2.7335587039794393342449301e-01,
+ 9.5586257685042792878173752e-01,
+ -2.099421066779969164496634e-01,
+ 2.135578780799860532750616e-01,
+ -8.694568971167362743327708e-01,
+ 4.019566681155577786649878e-01,
+ 9.6778633541687993721617774e-01,
+ -6.734405869050344734943028e-01,
}
var sinh = []float64{
- 7.2661916084208533e+01,
- 1.1479409110035194e+03,
- -2.8043136512812520e-01,
- -7.4994290911815868e+01,
- 7.6552466042906761e+03,
- 9.3031583421672010e+00,
- 9.3308157558281088e+01,
- 7.6179893137269143e+00,
- 3.0217691805496156e+00,
- -2.9595057572444951e+03,
+ 7.2661916084208532301448439e+01,
+ 1.1479409110035194500526446e+03,
+ -2.8043136512812518927312641e-01,
+ -7.499429091181587232835164e+01,
+ 7.6552466042906758523925934e+03,
+ 9.3031583421672014313789064e+00,
+ 9.330815755828109072810322e+01,
+ 7.6179893137269146407361477e+00,
+ 3.021769180549615819524392e+00,
+ -2.95950575724449499189888e+03,
}
var sqrt = []float64{
- 2.2313699659365484e+00,
- 2.7818829009464263e+00,
- 5.2619393496314792e-01,
- 2.2384377628763938e+00,
- 3.1042380236055380e+00,
- 1.7106657298385224e+00,
- 2.2867189227054791e+00,
- 1.6516476350711160e+00,
- 1.3510396336454586e+00,
- 2.9471892997524950e+00,
+ 2.2313699659365484748756904e+00,
+ 2.7818829009464263511285458e+00,
+ 5.2619393496314796848143251e-01,
+ 2.2384377628763938724244104e+00,
+ 3.1042380236055381099288487e+00,
+ 1.7106657298385224403917771e+00,
+ 2.286718922705479046148059e+00,
+ 1.6516476350711159636222979e+00,
+ 1.3510396336454586262419247e+00,
+ 2.9471892997524949215723329e+00,
}
var tan = []float64{
- -3.6613165650402277e+00,
- 8.6490023264859754e+00,
- -2.8417941955033615e-01,
- 3.2532901859747287e+00,
- 2.1472756403802937e-01,
- -2.1860091071106700e-01,
- -1.7600028178723679e+00,
- -4.3898089147528178e-01,
- -3.8438855602011305e+00,
- 9.1098879337768517e-01,
+ -3.661316565040227801781974e+00,
+ 8.64900232648597589369854e+00,
+ -2.8417941955033612725238097e-01,
+ 3.253290185974728640827156e+00,
+ 2.147275640380293804770778e-01,
+ -2.18600910711067004921551e-01,
+ -1.760002817872367935518928e+00,
+ -4.389808914752818126249079e-01,
+ -3.843885560201130679995041e+00,
+ 9.10988793377685105753416e-01,
}
var tanh = []float64{
- 9.9990531206936328e-01,
- 9.9999962057085307e-01,
- -2.7001505097318680e-01,
- -9.9991110943061700e-01,
- 9.9999999146798441e-01,
- 9.9427249436125233e-01,
- 9.9994257600983156e-01,
- 9.9149409509772863e-01,
- 9.4936501296239700e-01,
- -9.9999994291374019e-01,
+ 9.9990531206936338549262119e-01,
+ 9.9999962057085294197613294e-01,
+ -2.7001505097318677233756845e-01,
+ -9.9991110943061718603541401e-01,
+ 9.9999999146798465745022007e-01,
+ 9.9427249436125236705001048e-01,
+ 9.9994257600983138572705076e-01,
+ 9.9149409509772875982054701e-01,
+ 9.4936501296239685514466577e-01,
+ -9.9999994291374030946055701e-01,
+}
+var trunc = []float64{
+ 4.0000000000000000e+00,
+ 7.0000000000000000e+00,
+ -0.0000000000000000e+00,
+ -5.0000000000000000e+00,
+ 9.0000000000000000e+00,
+ 2.0000000000000000e+00,
+ 5.0000000000000000e+00,
+ 2.0000000000000000e+00,
+ 1.0000000000000000e+00,
+ -8.0000000000000000e+00,
}
// arguments and expected results for special cases
+var vfacoshSC = []float64{
+ Inf(-1),
+ 0.5,
+ NaN(),
+}
+var acoshSC = []float64{
+ NaN(),
+ NaN(),
+ NaN(),
+}
+
var vfasinSC = []float64{
NaN(),
-Pi,
@@ -215,6 +364,17 @@
NaN(),
}
+var vfasinhSC = []float64{
+ Inf(-1),
+ Inf(1),
+ NaN(),
+}
+var asinhSC = []float64{
+ Inf(-1),
+ Inf(1),
+ NaN(),
+}
+
var vfatanSC = []float64{
NaN(),
}
@@ -222,6 +382,21 @@
NaN(),
}
+var vfatanhSC = []float64{
+ -Pi,
+ -1,
+ 1,
+ Pi,
+ NaN(),
+}
+var atanhSC = []float64{
+ NaN(),
+ Inf(-1),
+ Inf(1),
+ NaN(),
+ NaN(),
+}
+
var vfceilSC = []float64{
Inf(-1),
Inf(1),
@@ -233,6 +408,33 @@
NaN(),
}
+var vfcopysignSC = []float64{
+ Inf(-1),
+ Inf(1),
+ NaN(),
+}
+var copysignSC = []float64{
+ Inf(-1),
+ Inf(-1),
+ NaN(),
+}
+
+var vferfSC = []float64{
+ Inf(-1),
+ Inf(1),
+ NaN(),
+}
+var erfSC = []float64{
+ -1,
+ 1,
+ NaN(),
+}
+var erfcSC = []float64{
+ 2,
+ 0,
+ NaN(),
+}
+
var vfexpSC = []float64{
Inf(-1),
Inf(1),
@@ -243,6 +445,11 @@
Inf(1),
NaN(),
}
+var expm1SC = []float64{
+ -1,
+ Inf(1),
+ NaN(),
+}
var vffmodSC = [][2]float64{
[2]float64{Inf(-1), Inf(-1)},
@@ -359,6 +566,21 @@
NaN(),
}
+var vflog1pSC = []float64{
+ Inf(-1),
+ -Pi,
+ -1,
+ Inf(1),
+ NaN(),
+}
+var log1pSC = []float64{
+ NaN(),
+ NaN(),
+ Inf(-1),
+ Inf(1),
+ NaN(),
+}
+
var vfpowSC = [][2]float64{
[2]float64{-Pi, Pi},
[2]float64{-Pi, -Pi},
@@ -494,8 +716,9 @@
func TestAcos(t *testing.T) {
for i := 0; i < len(vf); i++ {
- if f := Acos(vf[i] / 10); !close(acos[i], f) {
- t.Errorf("Acos(%g) = %g, want %g\n", vf[i]/10, f, acos[i])
+ a := vf[i] / 10
+ if f := Acos(a); !close(acos[i], f) {
+ t.Errorf("Acos(%g) = %g, want %g\n", a, f, acos[i])
}
}
for i := 0; i < len(vfasinSC); i++ {
@@ -505,10 +728,25 @@
}
}
+func TestAcosh(t *testing.T) {
+ for i := 0; i < len(vf); i++ {
+ a := 1 + Fabs(vf[i])
+ if f := Acosh(a); !veryclose(acosh[i], f) {
+ t.Errorf("Acosh(%g) = %g, want %g\n", a, f, acosh[i])
+ }
+ }
+ for i := 0; i < len(vfacoshSC); i++ {
+ if f := Acosh(vfacoshSC[i]); !alike(acoshSC[i], f) {
+ t.Errorf("Acosh(%g) = %g, want %g\n", vfacoshSC[i], f, acoshSC[i])
+ }
+ }
+}
+
func TestAsin(t *testing.T) {
for i := 0; i < len(vf); i++ {
- if f := Asin(vf[i] / 10); !veryclose(asin[i], f) {
- t.Errorf("Asin(%g) = %g, want %g\n", vf[i]/10, f, asin[i])
+ a := vf[i] / 10
+ if f := Asin(a); !veryclose(asin[i], f) {
+ t.Errorf("Asin(%g) = %g, want %g\n", a, f, asin[i])
}
}
for i := 0; i < len(vfasinSC); i++ {
@@ -518,6 +756,19 @@
}
}
+func TestAsinh(t *testing.T) {
+ for i := 0; i < len(vf); i++ {
+ if f := Asinh(vf[i]); !veryclose(asinh[i], f) {
+ t.Errorf("Asinh(%g) = %g, want %g\n", vf[i], f, asinh[i])
+ }
+ }
+ for i := 0; i < len(vfasinhSC); i++ {
+ if f := Asinh(vfasinhSC[i]); !alike(asinhSC[i], f) {
+ t.Errorf("Asinh(%g) = %g, want %g\n", vfasinhSC[i], f, asinhSC[i])
+ }
+ }
+}
+
func TestAtan(t *testing.T) {
for i := 0; i < len(vf); i++ {
if f := Atan(vf[i]); !veryclose(atan[i], f) {
@@ -531,6 +782,20 @@
}
}
+func TestAtanh(t *testing.T) {
+ for i := 0; i < len(vf); i++ {
+ a := vf[i] / 10
+ if f := Atanh(a); !veryclose(atanh[i], f) {
+ t.Errorf("Atanh(%g) = %g, want %g\n", a, f, atanh[i])
+ }
+ }
+ for i := 0; i < len(vfatanhSC); i++ {
+ if f := Atanh(vfatanhSC[i]); !alike(atanhSC[i], f) {
+ t.Errorf("Atanh(%g) = %g, want %g\n", vfatanhSC[i], f, atanhSC[i])
+ }
+ }
+}
+
func TestCeil(t *testing.T) {
for i := 0; i < len(vf); i++ {
if f := Ceil(vf[i]); ceil[i] != f {
@@ -544,6 +809,63 @@
}
}
+func TestCopysign(t *testing.T) {
+ for i := 0; i < len(vf); i++ {
+ if f := Copysign(vf[i], -1); copysign[i] != f {
+ t.Errorf("Copysign(%g, -1) = %g, want %g\n", vf[i], f, copysign[i])
+ }
+ }
+ for i := 0; i < len(vfcopysignSC); i++ {
+ if f := Copysign(vfcopysignSC[i], -1); !alike(copysignSC[i], f) {
+ t.Errorf("Copysign(%g, -1) = %g, want %g\n", vfcopysignSC[i], f, copysignSC[i])
+ }
+ }
+}
+
+func TestCos(t *testing.T) {
+ for i := 0; i < len(vf); i++ {
+ if f := Cos(vf[i]); !close(cos[i], f) {
+ t.Errorf("Cos(%g) = %g, want %g\n", vf[i], f, cos[i])
+ }
+ }
+}
+
+func TestCosh(t *testing.T) {
+ for i := 0; i < len(vf); i++ {
+ if f := Cosh(vf[i]); !veryclose(cosh[i], f) {
+ t.Errorf("Cosh(%g) = %g, want %g\n", vf[i], f, cosh[i])
+ }
+ }
+}
+
+func TestErf(t *testing.T) {
+ for i := 0; i < len(vf); i++ {
+ a := vf[i] / 10
+ if f := Erf(a); !veryclose(erf[i], f) {
+ t.Errorf("Erf(%g) = %g, want %g\n", a, f, erf[i])
+ }
+ }
+ for i := 0; i < len(vferfSC); i++ {
+ if f := Erf(vferfSC[i]); !alike(erfSC[i], f) {
+ t.Errorf("Erf(%g) = %g, want %g\n", vferfSC[i], f, erfSC[i])
+ }
+ }
+}
+
+func TestErfc(t *testing.T) {
+ for i := 0; i < len(vf); i++ {
+ a := vf[i] / 10
+ if f := Erfc(a); !veryclose(erfc[i], f) {
+ t.Errorf("Erfc(%g) = %g, want %g\n", a, f, erfc[i])
+ }
+ }
+ for i := 0; i < len(vferfSC); i++ {
+ if f := Erfc(vferfSC[i]); !alike(erfcSC[i], f) {
+ t.Errorf("Erfc(%g) = %g, want %g\n", vferfSC[i], f, erfcSC[i])
+ }
+ }
+}
+
func TestExp(t *testing.T) {
for i := 0; i < len(vf); i++ {
if f := Exp(vf[i]); !close(exp[i], f) {
@@ -557,6 +879,20 @@
}
}
+func TestExpm1(t *testing.T) {
+ for i := 0; i < len(vf); i++ {
+ a := vf[i] / 100
+ if f := Expm1(a); !veryclose(expm1[i], f) {
+ t.Errorf("Expm1(%.26fg) = %.26fg, want %.26fg\n", a, f, expm1[i])
+ }
+ }
+ for i := 0; i < len(vfexpSC); i++ {
+ if f := Expm1(vfexpSC[i]); !alike(expm1SC[i], f) {
+ t.Errorf("Expm1(%g) = %g, want %g\n", vfexpSC[i], f, expm1SC[i])
+ }
+ }
+}
+
func TestFloor(t *testing.T) {
for i := 0; i < len(vf); i++ {
if f := Floor(vf[i]); floor[i] != f {
@@ -572,8 +908,8 @@
func TestFmod(t *testing.T) {
for i := 0; i < len(vf); i++ {
- if f := Fmod(10, vf[i]); !close(fmod[i], f) {
- t.Errorf("Fmod(10, %g) = %g, want %g\n", vf[i], f, fmod[i])
+ if f := Fmod(10, vf[i]); fmod[i] != f { /*!close(fmod[i], f)*/
+ t.Errorf("Fmod(10, %g) = %g, want %g\n", vf[i], f, fmod[i])
}
}
for i := 0; i < len(vffmodSC); i++ {
@@ -631,6 +967,24 @@
}
}
+func TestLog1p(t *testing.T) {
+ for i := 0; i < len(vf); i++ {
+ a := vf[i] / 100
+ if f := Log1p(a); !veryclose(log1p[i], f) {
+ t.Errorf("Log1p(%g) = %g, want %g\n", a, f, log1p[i])
+ }
+ }
+ a := float64(9)
+ if f := Log1p(a); f != Ln10 {
+ t.Errorf("Log1p(%g) = %g, want %g\n", a, f, Ln10)
+ }
+ for i := 0; i < len(vflogSC); i++ {
+ if f := Log1p(vflog1pSC[i]); !alike(log1pSC[i], f) {
+ t.Errorf("Log1p(%g) = %g, want %g\n", vflog1pSC[i], f, log1pSC[i])
+ }
+ }
+}
+
func TestPow(t *testing.T) {
for i := 0; i < len(vf); i++ {
if f := Pow(10, vf[i]); !close(pow[i], f) {
@@ -694,6 +1048,19 @@
}
}
+func TestTrunc(t *testing.T) {
+ for i := 0; i < len(vf); i++ {
+ if f := Trunc(vf[i]); trunc[i] != f {
+ t.Errorf("Trunc(%g) = %g, want %g\n", vf[i], f, trunc[i])
+ }
+ }
+ for i := 0; i < len(vfceilSC); i++ {
+ if f := Trunc(vfceilSC[i]); !alike(ceilSC[i], f) {
+ t.Errorf("Trunc(%g) = %g, want %g\n", vfceilSC[i], f, ceilSC[i])
+ }
+ }
+}
+
// Check that math functions of high angle values
// return similar results to low angle values
func TestLargeSin(t *testing.T) {
@@ -718,7 +1085,6 @@
}
}
-
func TestLargeTan(t *testing.T) {
large := float64(100000 * Pi)
for i := 0; i < len(vf); i++ {
@@ -758,9 +1124,15 @@
// Benchmarks
-func BenchmarkAtan(b *testing.B) {
+func BenchmarkAcos(b *testing.B) {
for i := 0; i < b.N; i++ {
- Atan(.5)
+ Acos(.5)
+ }
+}
+
+func BenchmarkAcosh(b *testing.B) {
+ for i := 0; i < b.N; i++ {
+ Acosh(.5)
}
}
@@ -770,9 +1142,105 @@
}
}
-func BenchmarkAcos(b *testing.B) {
+func BenchmarkAsinh(b *testing.B) {
for i := 0; i < b.N; i++ {
- Acos(.5)
+ Asinh(.5)
+ }
+}
+
+func BenchmarkAtan(b *testing.B) {
+ for i := 0; i < b.N; i++ {
+ Atan(.5)
+ }
+}
+
+func BenchmarkAtanh(b *testing.B) {
+ for i := 0; i < b.N; i++ {
+ Atanh(.5)
+ }
+}
+
+func BenchmarkCeil(b *testing.B) {
+ for i := 0; i < b.N; i++ {
+ Ceil(.5)
+ }
+}
+
+func BenchmarkCopysign(b *testing.B) {
+ for i := 0; i < b.N; i++ {
+ Copysign(.5, -1)
+ }
+}
+
+func BenchmarkCos(b *testing.B) {
+ for i := 0; i < b.N; i++ {
+ Cos(.5)
+ }
+}
+
+func BenchmarkCosh(b *testing.B) {
+ for i := 0; i < b.N; i++ {
+ Cosh(2.5)
+ }
+}
+
+func BenchmarkErf(b *testing.B) {
+ for i := 0; i < b.N; i++ {
+ Erf(.5)
+ }
+}
+
+func BenchmarkErfc(b *testing.B) {
+ for i := 0; i < b.N; i++ {
+ Erfc(.5)
+ }
+}
+
+func BenchmarkExp(b *testing.B) {
+ for i := 0; i < b.N; i++ {
+ Exp(.5)
+ }
+}
+
+func BenchmarkExpm1(b *testing.B) {
+ for i := 0; i < b.N; i++ {
+ Expm1(.5)
+ }
+}
+
+func BenchmarkFloor(b *testing.B) {
+ for i := 0; i < b.N; i++ {
+ Floor(.5)
+ }
+}
+
+func BenchmarkFmod(b *testing.B) {
+ for i := 0; i < b.N; i++ {
+ Fmod(10, 3)
+ }
+}
+
+func BenchmarkHypot(b *testing.B) {
+ for i := 0; i < b.N; i++ {
+ Hypot(3, 4)
+ }
+}
+
+func BenchmarkLog(b *testing.B) {
+ for i := 0; i < b.N; i++ {
+ Log(.5)
+ }
+}
+
+func BenchmarkLog10(b *testing.B) {
+ for i := 0; i < b.N; i++ {
+ Log10(.5)
+ }
+}
+
+func BenchmarkLog1p(b *testing.B) {
+ for i := 0; i < b.N; i++ {
+ Log1p(.5)
}
}
@@ -788,8 +1256,43 @@
}
}
+func BenchmarkSin(b *testing.B) {
+ for i := 0; i < b.N; i++ {
+ Sin(.5)
+ }
+}
+
+func BenchmarkSinh(b *testing.B) {
+ for i := 0; i < b.N; i++ {
+ Sinh(2.5)
+ }
+}
+
func BenchmarkSqrt(b *testing.B) {
for i := 0; i < b.N; i++ {
Sqrt(10)
}
}
+
+func BenchmarkSqrtGo(b *testing.B) {
+ for i := 0; i < b.N; i++ {
+ SqrtGo(10)
+ }
+}
+
+func BenchmarkTan(b *testing.B) {
+ for i := 0; i < b.N; i++ {
+ Tan(.5)
+ }
+}
+
+func BenchmarkTanh(b *testing.B) {
+ for i := 0; i < b.N; i++ {
+ Tanh(2.5)
+ }
+}
+func BenchmarkTrunc(b *testing.B) {
+ for i := 0; i < b.N; i++ {
+ Trunc(.5)
+ }
+}
diff --git a/src/pkg/math/asin.go b/src/pkg/math/asin.go
index a9df663..63cac1c 100644
--- a/src/pkg/math/asin.go
+++ b/src/pkg/math/asin.go
@@ -6,13 +6,16 @@
/*
- Floating-point sine and cosine.
+ Floating-point arcsine and arccosine.
They are implemented by computing the arctangent
after appropriate range reduction.
*/
// Asin returns the arcsine of x.
+//
+// Special case is:
+// Asin(x) = NaN if x < -1 or x > 1
func Asin(x float64) float64 {
sign := false
if x < 0 {
@@ -37,4 +40,7 @@
}
// Acos returns the arccosine of x.
+//
+// Special case is:
+// Acos(x) = NaN if x < -1 or x > 1
func Acos(x float64) float64 { return Pi/2 - Asin(x) }
diff --git a/src/pkg/math/asinh.go b/src/pkg/math/asinh.go
new file mode 100644
index 0000000..b38bbd7
--- /dev/null
+++ b/src/pkg/math/asinh.go
@@ -0,0 +1,72 @@
+// Copyright 2010 The Go Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+package math
+
+
+// The original C code, the long comment, and the constants
+// below are from FreeBSD's /usr/src/lib/msun/src/s_asinh.c
+// and came with this notice. The go code is a simplified
+// version of the original C.
+//
+// ====================================================
+// Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+//
+// Developed at SunPro, a Sun Microsystems, Inc. business.
+// Permission to use, copy, modify, and distribute this
+// software is freely granted, provided that this notice
+// is preserved.
+// ====================================================
+//
+//
+// asinh(x)
+// Method :
+// Based on
+// asinh(x) = sign(x) * log [ |x| + sqrt(x*x+1) ]
+// we have
+// asinh(x) := x if 1+x*x=1,
+// := sign(x)*(log(x)+ln2)) for large |x|, else
+// := sign(x)*log(2|x|+1/(|x|+sqrt(x*x+1))) if|x|>2, else
+// := sign(x)*log1p(|x| + x^2/(1 + sqrt(1+x^2)))
+//
+
+// Asinh(x) calculates the inverse hyperbolic sine of x.
+//
+// Special cases are:
+// Asinh(+Inf) = +Inf
+// Asinh(-Inf) = -Inf
+// Asinh(NaN) = NaN
+func Asinh(x float64) float64 {
+ const (
+ Ln2 = 6.93147180559945286227e-01 // 0x3FE62E42FEFA39EF
+ NearZero = 1.0 / (1 << 28) // 2^-28
+ Large = 1 << 28 // 2^28
+ )
+ // TODO(rsc): Remove manual inlining of IsNaN, IsInf
+ // when compiler does it for us
+ // special cases
+ if x != x || x > MaxFloat64 || x < -MaxFloat64 { // IsNaN(x) || IsInf(x, 0)
+ return x
+ }
+ sign := false
+ if x < 0 {
+ x = -x
+ sign = true
+ }
+ var temp float64
+ switch {
+ case x > Large:
+ temp = Log(x) + Ln2 // |x| > 2^28
+ case x > 2:
+ temp = Log(2*x + 1/(Sqrt(x*x+1)+x)) // 2^28 > |x| > 2.0
+ case x < NearZero:
+ temp = x // |x| < 2^-28
+ default:
+ temp = Log1p(x + x*x/(1+Sqrt(1+x*x))) // 2.0 > |x| > 2^-28
+ }
+ if sign {
+ temp = -temp
+ }
+ return temp
+}
diff --git a/src/pkg/math/atanh.go b/src/pkg/math/atanh.go
new file mode 100644
index 0000000..72ae2a6
--- /dev/null
+++ b/src/pkg/math/atanh.go
@@ -0,0 +1,79 @@
+// Copyright 2010 The Go Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+package math
+
+
+// The original C code, the long comment, and the constants
+// below are from FreeBSD's /usr/src/lib/msun/src/e_atanh.c
+// and came with this notice. The go code is a simplified
+// version of the original C.
+//
+// ====================================================
+// Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+//
+// Developed at SunPro, a Sun Microsystems, Inc. business.
+// Permission to use, copy, modify, and distribute this
+// software is freely granted, provided that this notice
+// is preserved.
+// ====================================================
+//
+//
+// __ieee754_atanh(x)
+// Method :
+// 1. Reduce x to positive by atanh(-x) = -atanh(x)
+// 2. For x>=0.5
+// 1 2x x
+// atanh(x) = --- * log(1 + -------) = 0.5 * log1p(2 * --------)
+// 2 1 - x 1 - x
+//
+// For x<0.5
+// atanh(x) = 0.5*log1p(2x+2x*x/(1-x))
+//
+// Special cases:
+// atanh(x) is NaN if |x| > 1 with signal;
+// atanh(NaN) is that NaN with no signal;
+// atanh(+-1) is +-INF with signal.
+//
+
+// Atanh(x) calculates the inverse hyperbolic tangent of x.
+//
+// Special cases are:
+// Atanh(x) = NaN if x < -1 or x > 1
+// Atanh(1) = +Inf
+// Atanh(-1) = -Inf
+// Atanh(NaN) = NaN
+func Atanh(x float64) float64 {
+ const NearZero = 1.0 / (1 << 28) // 2^-28
+ // TODO(rsc): Remove manual inlining of IsNaN
+ // when compiler does it for us
+ // special cases
+ switch {
+ case x < -1 || x > 1 || x != x: // x < -1 || x > 1 || IsNaN(x):
+ return NaN()
+ case x == 1:
+ return Inf(1)
+ case x == -1:
+ return Inf(-1)
+ }
+ sign := false
+ if x < 0 {
+ x = -x
+ sign = true
+ }
+ var temp float64
+ switch {
+ case x < NearZero:
+ temp = x
+ case x < 0.5:
+ temp = x + x
+ temp = 0.5 * Log1p(temp+temp*x/(1-x))
+ default:
+ temp = 0.5 * Log1p((x+x)/(1-x))
+ }
+ if sign {
+ temp = -temp
+ }
+ return temp
+}
diff --git a/src/pkg/math/copysign.go b/src/pkg/math/copysign.go
new file mode 100644
index 0000000..6b4cc2a
--- /dev/null
+++ b/src/pkg/math/copysign.go
@@ -0,0 +1,15 @@
+// Copyright 2010 The Go Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+package math
+
+
+// Copysign(x, y) returns a value with the magnitude
+// of x and the sign of y.
+func Copysign(x, y float64) float64 {
+ if x < 0 && y > 0 || x > 0 && y < 0 {
+ return -x
+ }
+ return x
+}
diff --git a/src/pkg/math/erf.go b/src/pkg/math/erf.go
new file mode 100644
index 0000000..b9a945c
--- /dev/null
+++ b/src/pkg/math/erf.go
@@ -0,0 +1,340 @@
+// Copyright 2010 The Go Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+package math
+
+
+/*
+ Floating-point error function and complementary error function.
+*/
+
+// The original C code and the long comment below are
+// from FreeBSD's /usr/src/lib/msun/src/s_erf.c and
+// came with this notice. The go code is a simplified
+// version of the original C.
+//
+// ====================================================
+// Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+//
+// Developed at SunPro, a Sun Microsystems, Inc. business.
+// Permission to use, copy, modify, and distribute this
+// software is freely granted, provided that this notice
+// is preserved.
+// ====================================================
+//
+//
+// double erf(double x)
+// double erfc(double x)
+// x
+// 2 |\
+// erf(x) = --------- | exp(-t*t)dt
+// sqrt(pi) \|
+// 0
+//
+// erfc(x) = 1-erf(x)
+// Note that
+// erf(-x) = -erf(x)
+// erfc(-x) = 2 - erfc(x)
+//
+// Method:
+// 1. For |x| in [0, 0.84375]
+// erf(x) = x + x*R(x^2)
+// erfc(x) = 1 - erf(x) if x in [-.84375,0.25]
+// = 0.5 + ((0.5-x)-x*R) if x in [0.25,0.84375]
+// where R = P/Q where P is an odd poly of degree 8 and
+// Q is an odd poly of degree 10.
+// -57.90
+// | R - (erf(x)-x)/x | <= 2
+//
+//
+// Remark. The formula is derived by noting
+// erf(x) = (2/sqrt(pi))*(x - x^3/3 + x^5/10 - x^7/42 + ....)
+// and that
+// 2/sqrt(pi) = 1.128379167095512573896158903121545171688
+// is close to one. The interval is chosen because the fix
+// point of erf(x) is near 0.6174 (i.e., erf(x)=x when x is
+// near 0.6174), and by some experiment, 0.84375 is chosen to
+// guarantee the error is less than one ulp for erf.
+//
+// 2. For |x| in [0.84375,1.25], let s = |x| - 1, and
+// c = 0.84506291151 rounded to single (24 bits)
+// erf(x) = sign(x) * (c + P1(s)/Q1(s))
+// erfc(x) = (1-c) - P1(s)/Q1(s) if x > 0
+// 1+(c+P1(s)/Q1(s)) if x < 0
+// |P1/Q1 - (erf(|x|)-c)| <= 2**-59.06
+// Remark: here we use the taylor series expansion at x=1.
+// erf(1+s) = erf(1) + s*Poly(s)
+// = 0.845.. + P1(s)/Q1(s)
+// That is, we use rational approximation to approximate
+// erf(1+s) - (c = (single)0.84506291151)
+// Note that |P1/Q1|< 0.078 for x in [0.84375,1.25]
+// where
+// P1(s) = degree 6 poly in s
+// Q1(s) = degree 6 poly in s
+//
+// 3. For x in [1.25,1/0.35(~2.857143)],
+// erfc(x) = (1/x)*exp(-x*x-0.5625+R1/S1)
+// erf(x) = 1 - erfc(x)
+// where
+// R1(z) = degree 7 poly in z, (z=1/x^2)
+// S1(z) = degree 8 poly in z
+//
+// 4. For x in [1/0.35,28]
+// erfc(x) = (1/x)*exp(-x*x-0.5625+R2/S2) if x > 0
+// = 2.0 - (1/x)*exp(-x*x-0.5625+R2/S2) if -6<x<0
+// = 2.0 - tiny (if x <= -6)
+// erf(x) = sign(x)*(1.0 - erfc(x)) if x < 6, else
+// erf(x) = sign(x)*(1.0 - tiny)
+// where
+// R2(z) = degree 6 poly in z, (z=1/x^2)
+// S2(z) = degree 7 poly in z
+//
+// Note1:
+// To compute exp(-x*x-0.5625+R/S), let s be a single
+// precision number and s := x; then
+// -x*x = -s*s + (s-x)*(s+x)
+// exp(-x*x-0.5626+R/S) =
+// exp(-s*s-0.5625)*exp((s-x)*(s+x)+R/S);
+// Note2:
+// Here 4 and 5 make use of the asymptotic series
+// exp(-x*x)
+// erfc(x) ~ ---------- * ( 1 + Poly(1/x^2) )
+// x*sqrt(pi)
+// We use rational approximation to approximate
+// g(s)=f(1/x^2) = log(erfc(x)*x) - x*x + 0.5625
+// Here is the error bound for R1/S1 and R2/S2
+// |R1/S1 - f(x)| < 2**(-62.57)
+// |R2/S2 - f(x)| < 2**(-61.52)
+//
+// 5. For inf > x >= 28
+// erf(x) = sign(x) *(1 - tiny) (raise inexact)
+// erfc(x) = tiny*tiny (raise underflow) if x > 0
+// = 2 - tiny if x<0
+//
+// 7. Special case:
+// erf(0) = 0, erf(inf) = 1, erf(-inf) = -1,
+// erfc(0) = 1, erfc(inf) = 0, erfc(-inf) = 2,
+// erfc/erf(NaN) is NaN
+
+const (
+ erx = 8.45062911510467529297e-01 // 0x3FEB0AC160000000
+ // Coefficients for approximation to erf in [0, 0.84375]
+ efx = 1.28379167095512586316e-01 // 0x3FC06EBA8214DB69
+ efx8 = 1.02703333676410069053e+00 // 0x3FF06EBA8214DB69
+ pp0 = 1.28379167095512558561e-01 // 0x3FC06EBA8214DB68
+ pp1 = -3.25042107247001499370e-01 // 0xBFD4CD7D691CB913
+ pp2 = -2.84817495755985104766e-02 // 0xBF9D2A51DBD7194F
+ pp3 = -5.77027029648944159157e-03 // 0xBF77A291236668E4
+ pp4 = -2.37630166566501626084e-05 // 0xBEF8EAD6120016AC
+ qq1 = 3.97917223959155352819e-01 // 0x3FD97779CDDADC09
+ qq2 = 6.50222499887672944485e-02 // 0x3FB0A54C5536CEBA
+ qq3 = 5.08130628187576562776e-03 // 0x3F74D022C4D36B0F
+ qq4 = 1.32494738004321644526e-04 // 0x3F215DC9221C1A10
+ qq5 = -3.96022827877536812320e-06 // 0xBED09C4342A26120
+ // Coefficients for approximation to erf in [0.84375, 1.25]
+ pa0 = -2.36211856075265944077e-03 // 0xBF6359B8BEF77538
+ pa1 = 4.14856118683748331666e-01 // 0x3FDA8D00AD92B34D
+ pa2 = -3.72207876035701323847e-01 // 0xBFD7D240FBB8C3F1
+ pa3 = 3.18346619901161753674e-01 // 0x3FD45FCA805120E4
+ pa4 = -1.10894694282396677476e-01 // 0xBFBC63983D3E28EC
+ pa5 = 3.54783043256182359371e-02 // 0x3FA22A36599795EB
+ pa6 = -2.16637559486879084300e-03 // 0xBF61BF380A96073F
+ qa1 = 1.06420880400844228286e-01 // 0x3FBB3E6618EEE323
+ qa2 = 5.40397917702171048937e-01 // 0x3FE14AF092EB6F33
+ qa3 = 7.18286544141962662868e-02 // 0x3FB2635CD99FE9A7
+ qa4 = 1.26171219808761642112e-01 // 0x3FC02660E763351F
+ qa5 = 1.36370839120290507362e-02 // 0x3F8BEDC26B51DD1C
+ qa6 = 1.19844998467991074170e-02 // 0x3F888B545735151D
+ // Coefficients for approximation to erfc in [1.25, 1/0.35]
+ ra0 = -9.86494403484714822705e-03 // 0xBF843412600D6435
+ ra1 = -6.93858572707181764372e-01 // 0xBFE63416E4BA7360
+ ra2 = -1.05586262253232909814e+01 // 0xC0251E0441B0E726
+ ra3 = -6.23753324503260060396e+01 // 0xC04F300AE4CBA38D
+ ra4 = -1.62396669462573470355e+02 // 0xC0644CB184282266
+ ra5 = -1.84605092906711035994e+02 // 0xC067135CEBCCABB2
+ ra6 = -8.12874355063065934246e+01 // 0xC054526557E4D2F2
+ ra7 = -9.81432934416914548592e+00 // 0xC023A0EFC69AC25C
+ sa1 = 1.96512716674392571292e+01 // 0x4033A6B9BD707687
+ sa2 = 1.37657754143519042600e+02 // 0x4061350C526AE721
+ sa3 = 4.34565877475229228821e+02 // 0x407B290DD58A1A71
+ sa4 = 6.45387271733267880336e+02 // 0x40842B1921EC2868
+ sa5 = 4.29008140027567833386e+02 // 0x407AD02157700314
+ sa6 = 1.08635005541779435134e+02 // 0x405B28A3EE48AE2C
+ sa7 = 6.57024977031928170135e+00 // 0x401A47EF8E484A93
+ sa8 = -6.04244152148580987438e-02 // 0xBFAEEFF2EE749A62
+ // Coefficients for approximation to erfc in [1/.35, 28]
+ rb0 = -9.86494292470009928597e-03 // 0xBF84341239E86F4A
+ rb1 = -7.99283237680523006574e-01 // 0xBFE993BA70C285DE
+ rb2 = -1.77579549177547519889e+01 // 0xC031C209555F995A
+ rb3 = -1.60636384855821916062e+02 // 0xC064145D43C5ED98
+ rb4 = -6.37566443368389627722e+02 // 0xC083EC881375F228
+ rb5 = -1.02509513161107724954e+03 // 0xC09004616A2E5992
+ rb6 = -4.83519191608651397019e+02 // 0xC07E384E9BDC383F
+ sb1 = 3.03380607434824582924e+01 // 0x403E568B261D5190
+ sb2 = 3.25792512996573918826e+02 // 0x40745CAE221B9F0A
+ sb3 = 1.53672958608443695994e+03 // 0x409802EB189D5118
+ sb4 = 3.19985821950859553908e+03 // 0x40A8FFB7688C246A
+ sb5 = 2.55305040643316442583e+03 // 0x40A3F219CEDF3BE6
+ sb6 = 4.74528541206955367215e+02 // 0x407DA874E79FE763
+ sb7 = -2.24409524465858183362e+01 // 0xC03670E242712D62
+)
+
+// Erf(x) returns the error function of x.
+//
+// Special cases are:
+// Erf(+Inf) = 1
+// Erf(-Inf) = -1
+// Erf(NaN) = NaN
+func Erf(x float64) float64 {
+ const (
+ VeryTiny = 2.848094538889218e-306 // 0x0080000000000000
+ Small = 1.0 / (1 << 28) // 2^-28
+ )
+ // special cases
+ // TODO(rsc): Remove manual inlining of IsNaN, IsInf
+ // when compiler does it for us
+ switch {
+ case x != x: // IsNaN(x):
+ return NaN()
+ case x > MaxFloat64: // IsInf(x, 1):
+ return 1
+ case x < -MaxFloat64: // IsInf(x, -1):
+ return -1
+ }
+ sign := false
+ if x < 0 {
+ x = -x
+ sign = true
+ }
+ if x < 0.84375 { // |x| < 0.84375
+ var temp float64
+ if x < Small { // |x| < 2^-28
+ if x < VeryTiny {
+ temp = 0.125 * (8.0*x + efx8*x) // avoid underflow
+ } else {
+ temp = x + efx*x
+ }
+ } else {
+ z := x * x
+ r := pp0 + z*(pp1+z*(pp2+z*(pp3+z*pp4)))
+ s := 1 + z*(qq1+z*(qq2+z*(qq3+z*(qq4+z*qq5))))
+ y := r / s
+ temp = x + x*y
+ }
+ if sign {
+ return -temp
+ }
+ return temp
+ }
+ if x < 1.25 { // 0.84375 <= |x| < 1.25
+ s := x - 1
+ P := pa0 + s*(pa1+s*(pa2+s*(pa3+s*(pa4+s*(pa5+s*pa6)))))
+ Q := 1 + s*(qa1+s*(qa2+s*(qa3+s*(qa4+s*(qa5+s*qa6)))))
+ if sign {
+ return -erx - P/Q
+ }
+ return erx + P/Q
+ }
+ if x >= 6 { // inf > |x| >= 6
+ if sign {
+ return -1
+ }
+ return 1
+ }
+ s := 1 / (x * x)
+ var R, S float64
+ if x < 1/0.35 { // |x| < 1 / 0.35 ~ 2.857143
+ R = ra0 + s*(ra1+s*(ra2+s*(ra3+s*(ra4+s*(ra5+s*(ra6+s*ra7))))))
+ S = 1 + s*(sa1+s*(sa2+s*(sa3+s*(sa4+s*(sa5+s*(sa6+s*(sa7+s*sa8)))))))
+ } else { // |x| >= 1 / 0.35 ~ 2.857143
+ R = rb0 + s*(rb1+s*(rb2+s*(rb3+s*(rb4+s*(rb5+s*rb6)))))
+ S = 1 + s*(sb1+s*(sb2+s*(sb3+s*(sb4+s*(sb5+s*(sb6+s*sb7))))))
+ }
+ z := Float64frombits(Float64bits(x) & 0xffffffff00000000) // pseudo-single (20-bit) precison x
+ r := Exp(-z*z-0.5625) * Exp((z-x)*(z+x)+R/S)
+ if sign {
+ return r/x - 1
+ }
+ return 1 - r/x
+}
+
+// Erfc(x) returns the complementary error function of x.
+//
+// Special cases are:
+// Erf(+Inf) = 0
+// Erf(-Inf) = 2
+// Erf(NaN) = NaN
+func Erfc(x float64) float64 {
+ const Tiny = 1.0 / (1 << 56) // 2^-56
+ // special cases
+ // TODO(rsc): Remove manual inlining of IsNaN, IsInf
+ // when compiler does it for us
+ switch {
+ case x != x: // IsNaN(x):
+ return NaN()
+ case x > MaxFloat64: // IsInf(x, 1):
+ return 0
+ case x < -MaxFloat64: // IsInf(x, -1):
+ return 2
+ }
+ sign := false
+ if x < 0 {
+ x = -x
+ sign = true
+ }
+ if x < 0.84375 { // |x| < 0.84375
+ var temp float64
+ if x < Tiny { // |x| < 2^-56
+ temp = x
+ } else {
+ z := x * x
+ r := pp0 + z*(pp1+z*(pp2+z*(pp3+z*pp4)))
+ s := 1 + z*(qq1+z*(qq2+z*(qq3+z*(qq4+z*qq5))))
+ y := r / s
+ if x < 0.25 { // |x| < 1/4
+ temp = x + x*y
+ } else {
+ temp = 0.5 + (x*y + (x - 0.5))
+ }
+ }
+ if sign {
+ return 1 + temp
+ }
+ return 1 - temp
+ }
+ if x < 1.25 { // 0.84375 <= |x| < 1.25
+ s := x - 1
+ P := pa0 + s*(pa1+s*(pa2+s*(pa3+s*(pa4+s*(pa5+s*pa6)))))
+ Q := 1 + s*(qa1+s*(qa2+s*(qa3+s*(qa4+s*(qa5+s*qa6)))))
+ if sign {
+ return 1 + erx + P/Q
+ }
+ return 1 - erx - P/Q
+
+ }
+ if x < 28 { // |x| < 28
+ s := 1 / (x * x)
+ var R, S float64
+ if x < 1/0.35 { // |x| < 1 / 0.35 ~ 2.857143
+ R = ra0 + s*(ra1+s*(ra2+s*(ra3+s*(ra4+s*(ra5+s*(ra6+s*ra7))))))
+ S = 1 + s*(sa1+s*(sa2+s*(sa3+s*(sa4+s*(sa5+s*(sa6+s*(sa7+s*sa8)))))))
+ } else { // |x| >= 1 / 0.35 ~ 2.857143
+ if sign && x > 6 {
+ return 2 // x < -6
+ }
+ R = rb0 + s*(rb1+s*(rb2+s*(rb3+s*(rb4+s*(rb5+s*rb6)))))
+ S = 1 + s*(sb1+s*(sb2+s*(sb3+s*(sb4+s*(sb5+s*(sb6+s*sb7))))))
+ }
+ z := Float64frombits(Float64bits(x) & 0xffffffff00000000) // pseudo-single (20-bit) precison x
+ r := Exp(-z*z-0.5625) * Exp((z-x)*(z+x)+R/S)
+ if sign {
+ return 2 - r/x
+ }
+ return r / x
+ }
+ if sign {
+ return 2
+ }
+ return 0
+}
diff --git a/src/pkg/math/exp.go b/src/pkg/math/exp.go
index bc02fda..5ea58c0 100644
--- a/src/pkg/math/exp.go
+++ b/src/pkg/math/exp.go
@@ -86,7 +86,7 @@
// Special cases are:
// Exp(+Inf) = +Inf
// Exp(NaN) = NaN
-// Very large values overflow to -Inf or +Inf.
+// Very large values overflow to 0 or +Inf.
// Very small values underflow to 1.
func Exp(x float64) float64 {
const (
diff --git a/src/pkg/math/expm1.go b/src/pkg/math/expm1.go
new file mode 100644
index 0000000..9c516fb
--- /dev/null
+++ b/src/pkg/math/expm1.go
@@ -0,0 +1,238 @@
+// Copyright 2010 The Go Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+package math
+
+
+// The original C code, the long comment, and the constants
+// below are from FreeBSD's /usr/src/lib/msun/src/s_expm1.c
+// and came with this notice. The go code is a simplified
+// version of the original C.
+//
+// ====================================================
+// Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+//
+// Developed at SunPro, a Sun Microsystems, Inc. business.
+// Permission to use, copy, modify, and distribute this
+// software is freely granted, provided that this notice
+// is preserved.
+// ====================================================
+//
+// expm1(x)
+// Returns exp(x)-1, the exponential of x minus 1.
+//
+// Method
+// 1. Argument reduction:
+// Given x, find r and integer k such that
+//
+// x = k*ln2 + r, |r| <= 0.5*ln2 ~ 0.34658
+//
+// Here a correction term c will be computed to compensate
+// the error in r when rounded to a floating-point number.
+//
+// 2. Approximating expm1(r) by a special rational function on
+// the interval [0,0.34658]:
+// Since
+// r*(exp(r)+1)/(exp(r)-1) = 2+ r^2/6 - r^4/360 + ...
+// we define R1(r*r) by
+// r*(exp(r)+1)/(exp(r)-1) = 2+ r^2/6 * R1(r*r)
+// That is,
+// R1(r**2) = 6/r *((exp(r)+1)/(exp(r)-1) - 2/r)
+// = 6/r * ( 1 + 2.0*(1/(exp(r)-1) - 1/r))
+// = 1 - r^2/60 + r^4/2520 - r^6/100800 + ...
+// We use a special Reme algorithm on [0,0.347] to generate
+// a polynomial of degree 5 in r*r to approximate R1. The
+// maximum error of this polynomial approximation is bounded
+// by 2**-61. In other words,
+// R1(z) ~ 1.0 + Q1*z + Q2*z**2 + Q3*z**3 + Q4*z**4 + Q5*z**5
+// where Q1 = -1.6666666666666567384E-2,
+// Q2 = 3.9682539681370365873E-4,
+// Q3 = -9.9206344733435987357E-6,
+// Q4 = 2.5051361420808517002E-7,
+// Q5 = -6.2843505682382617102E-9;
+// (where z=r*r, and the values of Q1 to Q5 are listed below)
+// with error bounded by
+// | 5 | -61
+// | 1.0+Q1*z+...+Q5*z - R1(z) | <= 2
+// | |
+//
+// expm1(r) = exp(r)-1 is then computed by the following
+// specific way which minimize the accumulation rounding error:
+// 2 3
+// r r [ 3 - (R1 + R1*r/2) ]
+// expm1(r) = r + --- + --- * [--------------------]
+// 2 2 [ 6 - r*(3 - R1*r/2) ]
+//
+// To compensate the error in the argument reduction, we use
+// expm1(r+c) = expm1(r) + c + expm1(r)*c
+// ~ expm1(r) + c + r*c
+// Thus c+r*c will be added in as the correction terms for
+// expm1(r+c). Now rearrange the term to avoid optimization
+// screw up:
+// ( 2 2 )
+// ({ ( r [ R1 - (3 - R1*r/2) ] ) } r )
+// expm1(r+c)~r - ({r*(--- * [--------------------]-c)-c} - --- )
+// ({ ( 2 [ 6 - r*(3 - R1*r/2) ] ) } 2 )
+// ( )
+//
+// = r - E
+// 3. Scale back to obtain expm1(x):
+// From step 1, we have
+// expm1(x) = either 2^k*[expm1(r)+1] - 1
+// = or 2^k*[expm1(r) + (1-2^-k)]
+// 4. Implementation notes:
+// (A). To save one multiplication, we scale the coefficient Qi
+// to Qi*2^i, and replace z by (x^2)/2.
+// (B). To achieve maximum accuracy, we compute expm1(x) by
+// (i) if x < -56*ln2, return -1.0, (raise inexact if x!=inf)
+// (ii) if k=0, return r-E
+// (iii) if k=-1, return 0.5*(r-E)-0.5
+// (iv) if k=1 if r < -0.25, return 2*((r+0.5)- E)
+// else return 1.0+2.0*(r-E);
+// (v) if (k<-2||k>56) return 2^k(1-(E-r)) - 1 (or exp(x)-1)
+// (vi) if k <= 20, return 2^k((1-2^-k)-(E-r)), else
+// (vii) return 2^k(1-((E+2^-k)-r))
+//
+// Special cases:
+// expm1(INF) is INF, expm1(NaN) is NaN;
+// expm1(-INF) is -1, and
+// for finite argument, only expm1(0)=0 is exact.
+//
+// Accuracy:
+// according to an error analysis, the error is always less than
+// 1 ulp (unit in the last place).
+//
+// Misc. info.
+// For IEEE double
+// if x > 7.09782712893383973096e+02 then expm1(x) overflow
+//
+// Constants:
+// The hexadecimal values are the intended ones for the following
+// constants. The decimal values may be used, provided that the
+// compiler will convert from decimal to binary accurately enough
+// to produce the hexadecimal values shown.
+//
+
+// Expm1 returns e^x - 1, the base-e exponential of x minus 1.
+// It is more accurate than Exp(x) - 1 when x is near zero.
+//
+// Special cases are:
+// Expm1(+Inf) = +Inf
+// Expm1(-Inf) = -1
+// Expm1(NaN) = NaN
+// Very large values overflow to -1 or +Inf.
+func Expm1(x float64) float64 {
+ const (
+ Othreshold = 7.09782712893383973096e+02 // 0x40862E42FEFA39EF
+ Ln2X56 = 3.88162421113569373274e+01 // 0x4043687a9f1af2b1
+ Ln2HalfX3 = 1.03972077083991796413e+00 // 0x3ff0a2b23f3bab73
+ Ln2Half = 3.46573590279972654709e-01 // 0x3fd62e42fefa39ef
+ Ln2Hi = 6.93147180369123816490e-01 // 0x3fe62e42fee00000
+ Ln2Lo = 1.90821492927058770002e-10 // 0x3dea39ef35793c76
+ InvLn2 = 1.44269504088896338700e+00 // 0x3ff71547652b82fe
+ Tiny = 1.0 / (1 << 54) // 2^-54 = 0x3c90000000000000
+ // scaled coefficients related to expm1
+ Q1 = -3.33333333333331316428e-02 // 0xBFA11111111110F4
+ Q2 = 1.58730158725481460165e-03 // 0x3F5A01A019FE5585
+ Q3 = -7.93650757867487942473e-05 // 0xBF14CE199EAADBB7
+ Q4 = 4.00821782732936239552e-06 // 0x3ED0CFCA86E65239
+ Q5 = -2.01099218183624371326e-07 // 0xBE8AFDB76E09C32D
+ )
+
+ // special cases
+ // TODO(rsc): Remove manual inlining of IsNaN, IsInf
+ // when compiler does it for us
+ switch {
+ case x > MaxFloat64 || x != x: // IsInf(x, 1) || IsNaN(x):
+ return x
+ case x < -MaxFloat64: // IsInf(x, -1):
+ return -1
+ }
+
+ absx := x
+ sign := false
+ if x < 0 {
+ absx = -absx
+ sign = true
+ }
+
+ // filter out huge argument
+ if absx >= Ln2X56 { // if |x| >= 56 * ln2
+ if absx >= Othreshold { // if |x| >= 709.78...
+ return Inf(1) // overflow
+ }
+ if sign {
+ return -1 // x < -56*ln2, return -1.0
+ }
+ }
+
+ // argument reduction
+ var c float64
+ var k int
+ if absx > Ln2Half { // if |x| > 0.5 * ln2
+ var hi, lo float64
+ if absx < Ln2HalfX3 { // and |x| < 1.5 * ln2
+ if !sign {
+ hi = x - Ln2Hi
+ lo = Ln2Lo
+ k = 1
+ } else {
+ hi = x + Ln2Hi
+ lo = -Ln2Lo
+ k = -1
+ }
+ } else {
+ if !sign {
+ k = int(InvLn2*x + 0.5)
+ } else {
+ k = int(InvLn2*x - 0.5)
+ }
+ t := float64(k)
+ hi = x - t*Ln2Hi // t * Ln2Hi is exact here
+ lo = t * Ln2Lo
+ }
+ x = hi - lo
+ c = (hi - x) - lo
+ } else if absx < Tiny { // when |x| < 2^-54, return x
+ return x
+ } else {
+ k = 0
+ }
+
+ // x is now in primary range
+ hfx := 0.5 * x
+ hxs := x * hfx
+ r1 := 1 + hxs*(Q1+hxs*(Q2+hxs*(Q3+hxs*(Q4+hxs*Q5))))
+ t := 3 - r1*hfx
+ e := hxs * ((r1 - t) / (6.0 - x*t))
+ if k != 0 {
+ e = (x*(e-c) - c)
+ e -= hxs
+ switch {
+ case k == -1:
+ return 0.5*(x-e) - 0.5
+ case k == 1:
+ if x < -0.25 {
+ return -2 * (e - (x + 0.5))
+ }
+ return 1 + 2*(x-e)
+ case k <= -2 || k > 56: // suffice to return exp(x)-1
+ y := 1 - (e - x)
+ y = Float64frombits(Float64bits(y) + uint64(k)<<52) // add k to y's exponent
+ return y - 1
+ }
+ if k < 20 {
+ t := Float64frombits(0x3ff0000000000000 - (0x20000000000000 >> uint(k))) // t=1-2^-k
+ y := t - (e - x)
+ y = Float64frombits(Float64bits(y) + uint64(k)<<52) // add k to y's exponent
+ return y
+ }
+ t := Float64frombits(uint64((0x3ff - k) << 52)) // 2^-k
+ y := x - (e + t)
+ y += 1
+ y = Float64frombits(Float64bits(y) + uint64(k)<<52) // add k to y's exponent
+ return y
+ }
+ return x - (x*e - hxs) // c is 0
+}
diff --git a/src/pkg/math/fabs.go b/src/pkg/math/fabs.go
index ca78f3b..fcddb85 100644
--- a/src/pkg/math/fabs.go
+++ b/src/pkg/math/fabs.go
@@ -5,6 +5,11 @@
package math
// Fabs returns the absolute value of x.
+//
+// Special cases are:
+// Fabs(+Inf) = +Inf
+// Fabs(-Inf) = +Inf
+// Fabs(NaN) = NaN
func Fabs(x float64) float64 {
if x < 0 {
return -x
diff --git a/src/pkg/math/floor.go b/src/pkg/math/floor.go
index c5e496d..9270ba6 100644
--- a/src/pkg/math/floor.go
+++ b/src/pkg/math/floor.go
@@ -1,4 +1,4 @@
-// Copyright 2009 The Go Authors. All rights reserved.
+// Copyright 2009-2010 The Go Authors. All rights reserved.
// Use of this source code is governed by a BSD-style
// license that can be found in the LICENSE file.
@@ -12,6 +12,8 @@
// Floor(-Inf) = -Inf
// Floor(NaN) = NaN
func Floor(x float64) float64 {
+ // TODO(rsc): Remove manual inlining of IsNaN, IsInf
+ // when compiler does it for us
if x != x || x > MaxFloat64 || x < -MaxFloat64 { // IsNaN(x) || IsInf(x, 0)
return x
}
@@ -33,3 +35,19 @@
// Ceil(-Inf) = -Inf
// Ceil(NaN) = NaN
func Ceil(x float64) float64 { return -Floor(-x) }
+
+// Trunc returns the integer value of x.
+//
+// Special cases are:
+// Trunc(+Inf) = +Inf
+// Trunc(-Inf) = -Inf
+// Trunc(NaN) = NaN
+func Trunc(x float64) float64 {
+ // TODO(rsc): Remove manual inlining of IsNaN, IsInf
+ // when compiler does it for us
+ if x != x || x > MaxFloat64 || x < -MaxFloat64 { // IsNaN(x) || IsInf(x, 0)
+ return x
+ }
+ d, _ := Modf(x)
+ return d
+}
diff --git a/src/pkg/math/floor_386.s b/src/pkg/math/floor_386.s
index 2c5ff27..a4ae9d2 100644
--- a/src/pkg/math/floor_386.s
+++ b/src/pkg/math/floor_386.s
@@ -8,7 +8,7 @@
FSTCW -2(SP) // save old Control Word
MOVW -2(SP), AX
ANDW $0xf3ff, AX
- ORW $0x0800, AX // Rounding Control set to +Inf
+ ORW $0x0800, AX // Rounding Control set to +Inf
MOVW AX, -4(SP) // store new Control Word
FLDCW -4(SP) // load new Control Word
FRNDINT // F0=Ceil(x)
@@ -22,10 +22,23 @@
FSTCW -2(SP) // save old Control Word
MOVW -2(SP), AX
ANDW $0xf3ff, AX
- ORW $0x0400, AX // Rounding Control set to -Inf
+ ORW $0x0400, AX // Rounding Control set to -Inf
MOVW AX, -4(SP) // store new Control Word
FLDCW -4(SP) // load new Control Word
- FRNDINT // F0=floor(x)
+ FRNDINT // F0=Floor(x)
+ FLDCW -2(SP) // load old Control Word
+ FMOVDP F0, r+8(FP)
+ RET
+
+// func Trunc(x float64) float64
+TEXT ·Trunc(SB),7,$0
+ FMOVD x+0(FP), F0 // F0=x
+ FSTCW -2(SP) // save old Control Word
+ MOVW -2(SP), AX
+ ORW $0x0c00, AX // Rounding Control set to truncate
+ MOVW AX, -4(SP) // store new Control Word
+ FLDCW -4(SP) // load new Control Word
+ FRNDINT // F0=Trunc(x)
FLDCW -2(SP) // load old Control Word
FMOVDP F0, r+8(FP)
RET
diff --git a/src/pkg/math/floor_decl.go b/src/pkg/math/floor_decl.go
index 09f5646..7da4201 100644
--- a/src/pkg/math/floor_decl.go
+++ b/src/pkg/math/floor_decl.go
@@ -6,3 +6,4 @@
func Ceil(x float64) float64
func Floor(x float64) float64
+func Trunc(x float64) float64
diff --git a/src/pkg/math/fmod_386.s b/src/pkg/math/fmod_386.s
new file mode 100644
index 0000000..eb37bef
--- /dev/null
+++ b/src/pkg/math/fmod_386.s
@@ -0,0 +1,15 @@
+// Copyright 2010 The Go Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+// func Fmod(x, y float64) float64
+TEXT ·Fmod(SB),7,$0
+ FMOVD y+8(FP), F0 // F0=y
+ FMOVD x+0(FP), F0 // F0=x, F1=y
+ FPREM // F0=reduced_x, F1=y
+ FSTSW AX // AX=status word
+ ANDW $0x0400, AX
+ JNE -3(PC) // jump if reduction incomplete
+ FMOVDP F0, F1 // F0=x-q*y
+ FMOVDP F0, r+16(FP)
+ RET
diff --git a/src/pkg/math/fmod_decl.go b/src/pkg/math/fmod_decl.go
new file mode 100644
index 0000000..8d97cdf
--- /dev/null
+++ b/src/pkg/math/fmod_decl.go
@@ -0,0 +1,7 @@
+// Copyright 2010 The Go Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+package math
+
+func Fmod(x, y float64) float64
diff --git a/src/pkg/math/hypot_386.s b/src/pkg/math/hypot_386.s
index 73e4f57..212bb74 100644
--- a/src/pkg/math/hypot_386.s
+++ b/src/pkg/math/hypot_386.s
@@ -23,7 +23,7 @@
FTST // compare F0 to 0
FSTSW AX
ANDW $0x4000, AX
- JNE 10(PC) // jump if F0 = 0
+ JNE 10(PC) // jump if F0 = 0
FXCHD F0, F1 // F0=y (smaller), F1=x (larger)
FDIVD F1, F0 // F0=y(=y/x), F1=x
FMULD F0, F0 // F0=y*y, F1=x
diff --git a/src/pkg/math/log1p.go b/src/pkg/math/log1p.go
new file mode 100644
index 0000000..87a8b0e
--- /dev/null
+++ b/src/pkg/math/log1p.go
@@ -0,0 +1,200 @@
+// Copyright 2010 The Go Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+package math
+
+
+// The original C code, the long comment, and the constants
+// below are from FreeBSD's /usr/src/lib/msun/src/s_log1p.c
+// and came with this notice. The go code is a simplified
+// version of the original C.
+//
+// ====================================================
+// Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+//
+// Developed at SunPro, a Sun Microsystems, Inc. business.
+// Permission to use, copy, modify, and distribute this
+// software is freely granted, provided that this notice
+// is preserved.
+// ====================================================
+//
+//
+// double log1p(double x)
+//
+// Method :
+// 1. Argument Reduction: find k and f such that
+// 1+x = 2^k * (1+f),
+// where sqrt(2)/2 < 1+f < sqrt(2) .
+//
+// Note. If k=0, then f=x is exact. However, if k!=0, then f
+// may not be representable exactly. In that case, a correction
+// term is need. Let u=1+x rounded. Let c = (1+x)-u, then
+// log(1+x) - log(u) ~ c/u. Thus, we proceed to compute log(u),
+// and add back the correction term c/u.
+// (Note: when x > 2**53, one can simply return log(x))
+//
+// 2. Approximation of log1p(f).
+// Let s = f/(2+f) ; based on log(1+f) = log(1+s) - log(1-s)
+// = 2s + 2/3 s**3 + 2/5 s**5 + .....,
+// = 2s + s*R
+// We use a special Reme algorithm on [0,0.1716] to generate
+// a polynomial of degree 14 to approximate R The maximum error
+// of this polynomial approximation is bounded by 2**-58.45. In
+// other words,
+// 2 4 6 8 10 12 14
+// R(z) ~ Lp1*s +Lp2*s +Lp3*s +Lp4*s +Lp5*s +Lp6*s +Lp7*s
+// (the values of Lp1 to Lp7 are listed in the program)
+// a-0.2929nd
+// | 2 14 | -58.45
+// | Lp1*s +...+Lp7*s - R(z) | <= 2
+// | |
+// Note that 2s = f - s*f = f - hfsq + s*hfsq, where hfsq = f*f/2.
+// In order to guarantee error in log below 1ulp, we compute log
+// by
+// log1p(f) = f - (hfsq - s*(hfsq+R)).
+//
+// 3. Finally, log1p(x) = k*ln2 + log1p(f).
+// = k*ln2_hi+(f-(hfsq-(s*(hfsq+R)+k*ln2_lo)))
+// Here ln2 is split into two floating point number:
+// ln2_hi + ln2_lo,
+// where n*ln2_hi is always exact for |n| < 2000.
+//
+// Special cases:
+// log1p(x) is NaN with signal if x < -1 (including -INF) ;
+// log1p(+INF) is +INF; log1p(-1) is -INF with signal;
+// log1p(NaN) is that NaN with no signal.
+//
+// Accuracy:
+// according to an error analysis, the error is always less than
+// 1 ulp (unit in the last place).
+//
+// Constants:
+// The hexadecimal values are the intended ones for the following
+// constants. The decimal values may be used, provided that the
+// compiler will convert from decimal to binary accurately enough
+// to produce the hexadecimal values shown.
+//
+// Note: Assuming log() return accurate answer, the following
+// algorithm can be used to compute log1p(x) to within a few ULP:
+//
+// u = 1+x;
+// if(u==1.0) return x ; else
+// return log(u)*(x/(u-1.0));
+//
+// See HP-15C Advanced Functions Handbook, p.193.
+
+// Log1p returns the natural logarithm of 1 plus its argument x.
+// It is more accurate than Log(1 + x) when x is near zero.
+//
+// Special cases are:
+// Log1p(+Inf) = +Inf
+// Log1p(-1) = -Inf
+// Log1p(x < -1) = NaN
+// Log1p(NaN) = NaN
+func Log1p(x float64) float64 {
+ const (
+ Sqrt2M1 = 4.142135623730950488017e-01 // Sqrt(2)-1 = 0x3fda827999fcef34
+ Sqrt2HalfM1 = -2.928932188134524755992e-01 // Sqrt(2)/2-1 = 0xbfd2bec333018866
+ Small = 1.0 / (1 << 29) // 2^-29 = 0x3e20000000000000
+ Tiny = 1.0 / (1 << 54) // 2^-54
+ Two53 = 1 << 53 // 2^53
+ Ln2Hi = 6.93147180369123816490e-01 // 3fe62e42fee00000
+ Ln2Lo = 1.90821492927058770002e-10 // 3dea39ef35793c76
+ Lp1 = 6.666666666666735130e-01 // 3FE5555555555593
+ Lp2 = 3.999999999940941908e-01 // 3FD999999997FA04
+ Lp3 = 2.857142874366239149e-01 // 3FD2492494229359
+ Lp4 = 2.222219843214978396e-01 // 3FCC71C51D8E78AF
+ Lp5 = 1.818357216161805012e-01 // 3FC7466496CB03DE
+ Lp6 = 1.531383769920937332e-01 // 3FC39A09D078C69F
+ Lp7 = 1.479819860511658591e-01 // 3FC2F112DF3E5244
+ )
+
+ // special cases
+ // TODO(rsc): Remove manual inlining of IsNaN, IsInf
+ // when compiler does it for us
+ switch {
+ case x < -1 || x != x: // x < -1 || IsNaN(x): // includes -Inf
+ return NaN()
+ case x == -1:
+ return Inf(-1)
+ case x > MaxFloat64: // IsInf(x, 1):
+ return Inf(1)
+ }
+
+ absx := x
+ if absx < 0 {
+ absx = -absx
+ }
+
+ var f float64
+ var iu uint64
+ k := 1
+ if absx < Sqrt2M1 { // |x| < Sqrt(2)-1
+ if absx < Small { // |x| < 2^-29
+ if absx < Tiny { // |x| < 2^-54
+ return x
+ }
+ return x - x*x*0.5
+ }
+ if x > Sqrt2HalfM1 { // Sqrt(2)/2-1 < x
+ // (Sqrt(2)/2-1) < x < (Sqrt(2)-1)
+ k = 0
+ f = x
+ iu = 1
+ }
+ }
+ var c float64
+ if k != 0 {
+ var u float64
+ if absx < Two53 { // 1<<53
+ u = 1.0 + x
+ iu = Float64bits(u)
+ k = int((iu >> 52) - 1023)
+ if k > 0 {
+ c = 1.0 - (u - x)
+ } else {
+ c = x - (u - 1.0) // correction term
+ c /= u
+ }
+ } else {
+ u = x
+ iu = Float64bits(u)
+ k = int((iu >> 52) - 1023)
+ c = 0
+ }
+ iu &= 0x000fffffffffffff
+ if iu < 0x0006a09e667f3bcd { // mantissa of Sqrt(2)
+ u = Float64frombits(iu | 0x3ff0000000000000) // normalize u
+ } else {
+ k += 1
+ u = Float64frombits(iu | 0x3fe0000000000000) // normalize u/2
+ iu = (0x0010000000000000 - iu) >> 2
+ }
+ f = u - 1.0 // Sqrt(2)/2 < u < Sqrt(2)
+ }
+ hfsq := 0.5 * f * f
+ var s, R, z float64
+ if iu == 0 { // |f| < 2^-20
+ if f == 0 {
+ if k == 0 {
+ return 0
+ } else {
+ c += float64(k) * Ln2Lo
+ return float64(k)*Ln2Hi + c
+ }
+ }
+ R = hfsq * (1.0 - 0.66666666666666666*f) // avoid division
+ if k == 0 {
+ return f - R
+ }
+ return float64(k)*Ln2Hi - ((R - (float64(k)*Ln2Lo + c)) - f)
+ }
+ s = f / (2.0 + f)
+ z = s * s
+ R = z * (Lp1 + z*(Lp2+z*(Lp3+z*(Lp4+z*(Lp5+z*(Lp6+z*Lp7))))))
+ if k == 0 {
+ return f - (hfsq - s*(hfsq+R))
+ }
+ return float64(k)*Ln2Hi - ((hfsq - (s*(hfsq+R) + (float64(k)*Ln2Lo + c))) - f)
+}
