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 b10df65300db92b0fbaf0cc80138281b52cb44fd..f22808429164187e445185306dae5a2fab903180 100644 (file)
--- a/src/pkg/math/Makefile
+++ b/src/pkg/math/Makefile
@@ -15,6 +15,7 @@ OFILES_386=\
        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=\
        $(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 (file)
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 9490e06664d003c0242f4492aa9b4bd0b555df5e..59fdd9e1a7a06c65d0fd373e453a2f86d61ffdd8 100644 (file)
--- a/src/pkg/math/all_test.go
+++ b/src/pkg/math/all_test.go
@@ -22,41 +22,83 @@ var vf = []float64{
        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 @@ var ceil = []float64{
        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 @@ var floor = []float64{
        -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 @@ var asinSC = []float64{
        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 @@ var atanSC = []float64{
        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 @@ var ceilSC = []float64{
        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 @@ var expSC = []float64{
        Inf(1),
        NaN(),
 }
+var expm1SC = []float64{
+       -1,
+       Inf(1),
+       NaN(),
+}
 
 var vffmodSC = [][2]float64{
        [2]float64{Inf(-1), Inf(-1)},
@@ -359,6 +566,21 @@ var logSC = []float64{
        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 alike(a, b float64) bool {
 
 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 TestAcos(t *testing.T) {
        }
 }
 
+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 TestAsin(t *testing.T) {
        }
 }
 
+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 TestAtan(t *testing.T) {
        }
 }
 
+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 TestCeil(t *testing.T) {
        }
 }
 
+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 TestExp(t *testing.T) {
        }
 }
 
+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 TestFloor(t *testing.T) {
 
 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 TestLog10(t *testing.T) {
        }
 }
 
+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 TestTanh(t *testing.T) {
        }
 }
 
+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 TestLargeCos(t *testing.T) {
        }
 }
 
-
 func TestLargeTan(t *testing.T) {
        large := float64(100000 * Pi)
        for i := 0; i < len(vf); i++ {
@@ -758,9 +1124,15 @@ func TestFloatMinMax(t *testing.T) {
 
 // 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 BenchmarkAsin(b *testing.B) {
        }
 }
 
-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 BenchmarkPowFrac(b *testing.B) {
        }
 }
 
+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 a9df6631133cfc5c28fa95e5af081636e9132711..63cac1c4bc4f7c7d9dbf7102d12b4042985469bd 100644 (file)
--- a/src/pkg/math/asin.go
+++ b/src/pkg/math/asin.go
@@ -6,13 +6,16 @@ package math
 
 
 /*
-       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 @@ func Asin(x float64) float64 {
 }
 
 // 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 (file)
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 (file)
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 (file)
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 (file)
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 bc02fda10ccc1e6e1de3e6d2f4fac1e1acf1e736..5ea58c0fb391d302d285745e3c42f6015e0c069f 100644 (file)
--- a/src/pkg/math/exp.go
+++ b/src/pkg/math/exp.go
@@ -86,7 +86,7 @@ package math
 // 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 (file)
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 ca78f3b59f259c232c2d6aad81aa6f8921a84e43..fcddb8510019a60bcd3a8a642eb44ff55a0b8935 100644 (file)
--- 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 c5e496d8fab7b3c8c9c982801573b20f7a7ee336..9270ba6aa5145c07551e3991e404b0257f4a0a0d 100644 (file)
--- 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 @@ package math
 //     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 @@ func Floor(x float64) float64 {
 //     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 2c5ff270a1445cdce23855698f23b3d878b2daf4..a4ae9d2ebad48af3baad6ead41a80bd827be3a24 100644 (file)
--- a/src/pkg/math/floor_386.s
+++ b/src/pkg/math/floor_386.s
@@ -8,7 +8,7 @@ TEXT ·Ceil(SB),7,$0
        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 @@ TEXT ·Floor(SB),7,$0
        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 09f5646e3e2a0b738743ae20cec96382873a11bf..7da42017944d688dde053d98990401f85fdf96ca 100644 (file)
--- a/src/pkg/math/floor_decl.go
+++ b/src/pkg/math/floor_decl.go
@@ -6,3 +6,4 @@ package math
 
 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 (file)
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 (file)
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 73e4f577c32499e7cab864f7a3653dcd962a20d9..212bb74753f00650aa0010d79cdd776f40a0e17b 100644 (file)
--- a/src/pkg/math/hypot_386.s
+++ b/src/pkg/math/hypot_386.s
@@ -23,7 +23,7 @@ TEXT ·Hypot(SB),7,$0
        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 (file)
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)
+}