From: Charles L. Dorian Date: Tue, 2 Feb 2010 06:21:40 +0000 (-0800) Subject: math: add functions; update tests and special cases X-Git-Tag: weekly.2010-02-04~29 X-Git-Url: http://www.git.cypherpunks.su/?a=commitdiff_plain;h=a0690b69dad53ebeebc181749ede6d298272d120;p=gostls13.git 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 b10df65300..f228084291 100644 --- 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 index 0000000000..1f0d3f3380 --- /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 9490e06664..59fdd9e1a7 100644 --- 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 a9df663113..63cac1c4bc 100644 --- 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 index 0000000000..b38bbd78d0 --- /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 0000000000..72ae2a60f6 --- /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 0000000000..6b4cc2a4cf --- /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 0000000000..b9a945ce4b --- /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 >= 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 bc02fda10c..5ea58c0fb3 100644 --- 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 index 0000000000..9c516fb4ff --- /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 ca78f3b59f..fcddb85100 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 c5e496d8fa..9270ba6aa5 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 @@ 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 2c5ff270a1..a4ae9d2eba 100644 --- 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 09f5646e3e..7da4201794 100644 --- 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 index 0000000000..eb37bef406 --- /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 0000000000..8d97cdf4a0 --- /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 73e4f577c3..212bb74753 100644 --- 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 index 0000000000..87a8b0e221 --- /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) +}