From 554fd8c5195424bdbcabf5de30fdc183aba391bd Mon Sep 17 00:00:00 2001 From: upstream source tree Date: Sun, 15 Mar 2015 20:14:05 -0400 Subject: obtained gcc-4.6.4.tar.bz2 from upstream website; verified gcc-4.6.4.tar.bz2.sig; imported gcc-4.6.4 source tree from verified upstream tarball. downloading a git-generated archive based on the 'upstream' tag should provide you with a source tree that is binary identical to the one extracted from the above tarball. if you have obtained the source via the command 'git clone', however, do note that line-endings of files in your working directory might differ from line-endings of the respective files in the upstream repository. --- libgo/go/math/acosh.go | 62 + libgo/go/math/all_test.go | 2737 +++++++++++++++++++++++++++++++++++++++++++ libgo/go/math/asin.go | 50 + libgo/go/math/asinh.go | 72 ++ libgo/go/math/atan.go | 62 + libgo/go/math/atan2.go | 71 ++ libgo/go/math/atanh.go | 79 ++ libgo/go/math/bits.go | 59 + libgo/go/math/cbrt.go | 79 ++ libgo/go/math/const.go | 53 + libgo/go/math/copysign.go | 12 + libgo/go/math/erf.go | 340 ++++++ libgo/go/math/exp.go | 14 + libgo/go/math/exp2.go | 10 + libgo/go/math/exp_port.go | 192 +++ libgo/go/math/exp_test.go | 10 + libgo/go/math/expm1.go | 238 ++++ libgo/go/math/fabs.go | 21 + libgo/go/math/fdim.go | 29 + libgo/go/math/floor.go | 53 + libgo/go/math/fmod.go | 48 + libgo/go/math/frexp.go | 33 + libgo/go/math/gamma.go | 188 +++ libgo/go/math/hypot.go | 41 + libgo/go/math/hypot_port.go | 63 + libgo/go/math/hypot_test.go | 9 + libgo/go/math/j0.go | 433 +++++++ libgo/go/math/j1.go | 426 +++++++ libgo/go/math/jn.go | 310 +++++ libgo/go/math/ldexp.go | 45 + libgo/go/math/lgamma.go | 350 ++++++ libgo/go/math/log.go | 123 ++ libgo/go/math/log10.go | 13 + libgo/go/math/log10_decl.go | 8 + libgo/go/math/log1p.go | 200 ++++ libgo/go/math/logb.go | 54 + libgo/go/math/modf.go | 33 + libgo/go/math/nextafter.go | 29 + libgo/go/math/pow.go | 139 +++ libgo/go/math/pow10.go | 30 + libgo/go/math/remainder.go | 85 ++ libgo/go/math/signbit.go | 10 + libgo/go/math/sin.go | 66 ++ libgo/go/math/sincos.go | 13 + libgo/go/math/sinh.go | 68 ++ libgo/go/math/sqrt.go | 28 + libgo/go/math/sqrt_decl.go | 7 + libgo/go/math/sqrt_port.go | 143 +++ libgo/go/math/sqrt_test.go | 9 + libgo/go/math/tan.go | 65 + libgo/go/math/tanh.go | 28 + libgo/go/math/unsafe.go | 21 + 52 files changed, 7361 insertions(+) create mode 100644 libgo/go/math/acosh.go create mode 100644 libgo/go/math/all_test.go create mode 100644 libgo/go/math/asin.go create mode 100644 libgo/go/math/asinh.go create mode 100644 libgo/go/math/atan.go create mode 100644 libgo/go/math/atan2.go create mode 100644 libgo/go/math/atanh.go create mode 100644 libgo/go/math/bits.go create mode 100644 libgo/go/math/cbrt.go create mode 100644 libgo/go/math/const.go create mode 100644 libgo/go/math/copysign.go create mode 100644 libgo/go/math/erf.go create mode 100644 libgo/go/math/exp.go create mode 100644 libgo/go/math/exp2.go create mode 100644 libgo/go/math/exp_port.go create mode 100644 libgo/go/math/exp_test.go create mode 100644 libgo/go/math/expm1.go create mode 100644 libgo/go/math/fabs.go create mode 100644 libgo/go/math/fdim.go create mode 100644 libgo/go/math/floor.go create mode 100644 libgo/go/math/fmod.go create mode 100644 libgo/go/math/frexp.go create mode 100644 libgo/go/math/gamma.go create mode 100644 libgo/go/math/hypot.go create mode 100644 libgo/go/math/hypot_port.go create mode 100644 libgo/go/math/hypot_test.go create mode 100644 libgo/go/math/j0.go create mode 100644 libgo/go/math/j1.go create mode 100644 libgo/go/math/jn.go create mode 100644 libgo/go/math/ldexp.go create mode 100644 libgo/go/math/lgamma.go create mode 100644 libgo/go/math/log.go create mode 100644 libgo/go/math/log10.go create mode 100644 libgo/go/math/log10_decl.go create mode 100644 libgo/go/math/log1p.go create mode 100644 libgo/go/math/logb.go create mode 100644 libgo/go/math/modf.go create mode 100644 libgo/go/math/nextafter.go create mode 100644 libgo/go/math/pow.go create mode 100644 libgo/go/math/pow10.go create mode 100644 libgo/go/math/remainder.go create mode 100644 libgo/go/math/signbit.go create mode 100644 libgo/go/math/sin.go create mode 100644 libgo/go/math/sincos.go create mode 100644 libgo/go/math/sinh.go create mode 100644 libgo/go/math/sqrt.go create mode 100644 libgo/go/math/sqrt_decl.go create mode 100644 libgo/go/math/sqrt_port.go create mode 100644 libgo/go/math/sqrt_test.go create mode 100644 libgo/go/math/tan.go create mode 100644 libgo/go/math/tanh.go create mode 100644 libgo/go/math/unsafe.go (limited to 'libgo/go/math') diff --git a/libgo/go/math/acosh.go b/libgo/go/math/acosh.go new file mode 100644 index 000000000..d8067c065 --- /dev/null +++ b/libgo/go/math/acosh.go @@ -0,0 +1,62 @@ +// 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 + ) + // TODO(rsc): Remove manual inlining of IsNaN + // when compiler does it for us + // first case is special case + switch { + case x < 1 || x != x: // 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/libgo/go/math/all_test.go b/libgo/go/math/all_test.go new file mode 100644 index 000000000..d2a7d411e --- /dev/null +++ b/libgo/go/math/all_test.go @@ -0,0 +1,2737 @@ +// Copyright 2009 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_test + +import ( + "fmt" + . "math" + "runtime" + "testing" +) + +var vf = []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, +} +// 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). +var acos = []float64{ + 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.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.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 atan2 = []float64{ + 1.1088291730037004444527075e+00, + 9.1218183188715804018797795e-01, + 1.5984772603216203736068915e+00, + 2.0352918654092086637227327e+00, + 8.0391819139044720267356014e-01, + 1.2861075249894661588866752e+00, + 1.0889904479131695712182587e+00, + 1.3044821793397925293797357e+00, + 1.3902530903455392306872261e+00, + 2.2859857424479142655411058e+00, +} +var cbrt = []float64{ + 1.7075799841925094446722675e+00, + 1.9779982212970353936691498e+00, + -6.5177429017779910853339447e-01, + -1.7111838886544019873338113e+00, + 2.1279920909827937423960472e+00, + 1.4303536770460741452312367e+00, + 1.7357021059106154902341052e+00, + 1.3972633462554328350552916e+00, + 1.2221149580905388454977636e+00, + -2.0556003730500069110343596e+00, +} +var ceil = []float64{ + 5.0000000000000000e+00, + 8.0000000000000000e+00, + 0.0000000000000000e+00, + -5.0000000000000000e+00, + 1.0000000000000000e+01, + 3.0000000000000000e+00, + 6.0000000000000000e+00, + 3.0000000000000000e+00, + 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.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 exp2 = []float64{ + 3.1537839463286288034313104e+01, + 2.1361549283756232296144849e+02, + 8.2537402562185562902577219e-01, + 3.1021158628740294833424229e-02, + 7.9581744110252191462569661e+02, + 7.6019905892596359262696423e+00, + 3.7506882048388096973183084e+01, + 6.6250893439173561733216375e+00, + 3.5438267900243941544605339e+00, + 2.4281533133513300984289196e-03, +} +var fabs = []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 fdim = []float64{ + 4.9790119248836735e+00, + 7.7388724745781045e+00, + 0.0000000000000000e+00, + 0.0000000000000000e+00, + 9.6362937071984173e+00, + 2.9263772392439646e+00, + 5.2290834314593066e+00, + 2.7279399104360102e+00, + 1.8253080916808550e+00, + 0.0000000000000000e+00, +} +var floor = []float64{ + 4.0000000000000000e+00, + 7.0000000000000000e+00, + -1.0000000000000000e+00, + -6.0000000000000000e+00, + 9.0000000000000000e+00, + 2.0000000000000000e+00, + 5.0000000000000000e+00, + 2.0000000000000000e+00, + 1.0000000000000000e+00, + -9.0000000000000000e+00, +} +var fmod = []float64{ + 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, +} + +type fi struct { + f float64 + i int +} + +var frexp = []fi{ + {6.2237649061045918750e-01, 3}, + {9.6735905932226306250e-01, 3}, + {-5.5376011438400318000e-01, -1}, + {-6.2632545228388436250e-01, 3}, + {6.02268356699901081250e-01, 4}, + {7.3159430981099115000e-01, 2}, + {6.5363542893241332500e-01, 3}, + {6.8198497760900255000e-01, 2}, + {9.1265404584042750000e-01, 1}, + {-5.4287029803597508250e-01, 4}, +} +var gamma = []float64{ + 2.3254348370739963835386613898e+01, + 2.991153837155317076427529816e+03, + -4.561154336726758060575129109e+00, + 7.719403468842639065959210984e-01, + 1.6111876618855418534325755566e+05, + 1.8706575145216421164173224946e+00, + 3.4082787447257502836734201635e+01, + 1.579733951448952054898583387e+00, + 9.3834586598354592860187267089e-01, + -2.093995902923148389186189429e-05, +} +var j0 = []float64{ + -1.8444682230601672018219338e-01, + 2.27353668906331975435892e-01, + 9.809259936157051116270273e-01, + -1.741170131426226587841181e-01, + -2.1389448451144143352039069e-01, + -2.340905848928038763337414e-01, + -1.0029099691890912094586326e-01, + -1.5466726714884328135358907e-01, + 3.252650187653420388714693e-01, + -8.72218484409407250005360235e-03, +} +var j1 = []float64{ + -3.251526395295203422162967e-01, + 1.893581711430515718062564e-01, + -1.3711761352467242914491514e-01, + 3.287486536269617297529617e-01, + 1.3133899188830978473849215e-01, + 3.660243417832986825301766e-01, + -3.4436769271848174665420672e-01, + 4.329481396640773768835036e-01, + 5.8181350531954794639333955e-01, + -2.7030574577733036112996607e-01, +} +var j2 = []float64{ + 5.3837518920137802565192769e-02, + -1.7841678003393207281244667e-01, + 9.521746934916464142495821e-03, + 4.28958355470987397983072e-02, + 2.4115371837854494725492872e-01, + 4.842458532394520316844449e-01, + -3.142145220618633390125946e-02, + 4.720849184745124761189957e-01, + 3.122312022520957042957497e-01, + 7.096213118930231185707277e-02, +} +var jM3 = []float64{ + -3.684042080996403091021151e-01, + 2.8157665936340887268092661e-01, + 4.401005480841948348343589e-04, + 3.629926999056814081597135e-01, + 3.123672198825455192489266e-02, + -2.958805510589623607540455e-01, + -3.2033177696533233403289416e-01, + -2.592737332129663376736604e-01, + -1.0241334641061485092351251e-01, + -2.3762660886100206491674503e-01, +} +var lgamma = []fi{ + {3.146492141244545774319734e+00, 1}, + {8.003414490659126375852113e+00, 1}, + {1.517575735509779707488106e+00, -1}, + {-2.588480028182145853558748e-01, 1}, + {1.1989897050205555002007985e+01, 1}, + {6.262899811091257519386906e-01, 1}, + {3.5287924899091566764846037e+00, 1}, + {4.5725644770161182299423372e-01, 1}, + {-6.363667087767961257654854e-02, 1}, + {-1.077385130910300066425564e+01, -1}, +} +var log = []float64{ + 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 logb = []float64{ + 2.0000000000000000e+00, + 2.0000000000000000e+00, + -2.0000000000000000e+00, + 2.0000000000000000e+00, + 3.0000000000000000e+00, + 1.0000000000000000e+00, + 2.0000000000000000e+00, + 1.0000000000000000e+00, + 0.0000000000000000e+00, + 3.0000000000000000e+00, +} +var log10 = []float64{ + 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 log2 = []float64{ + 2.3158594707062190618898251e+00, + 2.9521233862883917703341018e+00, + -1.8526669502700329984917062e+00, + 2.3249844127278861543568029e+00, + 3.268478366538305087466309e+00, + 1.5491157592596970278166492e+00, + 2.3865580889631732407886495e+00, + 1.447811865817085365540347e+00, + 8.6813999540425116282815557e-01, + 3.118679457227342224364709e+00, +} +var modf = [][2]float64{ + {4.0000000000000000e+00, 9.7901192488367350108546816e-01}, + {7.0000000000000000e+00, 7.3887247457810456552351752e-01}, + {0.0000000000000000e+00, -2.7688005719200159404635997e-01}, + {-5.0000000000000000e+00, -1.060361827107492160848778e-02}, + {9.0000000000000000e+00, 6.3629370719841737980004837e-01}, + {2.0000000000000000e+00, 9.2637723924396464525443662e-01}, + {5.0000000000000000e+00, 2.2908343145930665230025625e-01}, + {2.0000000000000000e+00, 7.2793991043601025126008608e-01}, + {1.0000000000000000e+00, 8.2530809168085506044576505e-01}, + {-8.0000000000000000e+00, -6.8592476857560136238589621e-01}, +} +var nextafter = []float64{ + 4.97901192488367438926388786e+00, + 7.73887247457810545370193722e+00, + -2.7688005719200153853520874e-01, + -5.01060361827107403343006808e+00, + 9.63629370719841915615688777e+00, + 2.92637723924396508934364647e+00, + 5.22908343145930754047867595e+00, + 2.72793991043601069534929593e+00, + 1.82530809168085528249036997e+00, + -8.68592476857559958602905681e+00, +} +var pow = []float64{ + 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 remainder = []float64{ + 4.197615023265299782906368e-02, + 2.261127525421895434476482e+00, + 3.231794108794261433104108e-02, + -2.120723654214984321697556e-02, + 3.637062928015826201999516e-01, + 1.220868282268106064236690e+00, + -4.581668629186133046005125e-01, + -9.117596417440410050403443e-01, + 8.734595415957246977711748e-01, + 1.314075231424398637614104e+00, +} +var signbit = []bool{ + false, + false, + true, + true, + false, + false, + false, + false, + false, + true, +} +var sin = []float64{ + -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.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.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.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.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, +} +var y0 = []float64{ + -3.053399153780788357534855e-01, + 1.7437227649515231515503649e-01, + -8.6221781263678836910392572e-01, + -3.100664880987498407872839e-01, + 1.422200649300982280645377e-01, + 4.000004067997901144239363e-01, + -3.3340749753099352392332536e-01, + 4.5399790746668954555205502e-01, + 4.8290004112497761007536522e-01, + 2.7036697826604756229601611e-01, +} +var y1 = []float64{ + 0.15494213737457922210218611, + -0.2165955142081145245075746, + -2.4644949631241895201032829, + 0.1442740489541836405154505, + 0.2215379960518984777080163, + 0.3038800915160754150565448, + 0.0691107642452362383808547, + 0.2380116417809914424860165, + -0.20849492979459761009678934, + 0.0242503179793232308250804, +} +var y2 = []float64{ + 0.3675780219390303613394936, + -0.23034826393250119879267257, + -16.939677983817727205631397, + 0.367653980523052152867791, + -0.0962401471767804440353136, + -0.1923169356184851105200523, + 0.35984072054267882391843766, + -0.2794987252299739821654982, + -0.7113490692587462579757954, + -0.2647831587821263302087457, +} +var yM3 = []float64{ + -0.14035984421094849100895341, + -0.097535139617792072703973, + 242.25775994555580176377379, + -0.1492267014802818619511046, + 0.26148702629155918694500469, + 0.56675383593895176530394248, + -0.206150264009006981070575, + 0.64784284687568332737963658, + 1.3503631555901938037008443, + 0.1461869756579956803341844, +} + +// arguments and expected results for special cases +var vfacosSC = []float64{ + -Pi, + 1, + Pi, + NaN(), +} +var acosSC = []float64{ + NaN(), + 0, + NaN(), + NaN(), +} + +var vfacoshSC = []float64{ + Inf(-1), + 0.5, + 1, + Inf(1), + NaN(), +} +var acoshSC = []float64{ + NaN(), + NaN(), + 0, + Inf(1), + NaN(), +} + +var vfasinSC = []float64{ + -Pi, + Copysign(0, -1), + 0, + Pi, + NaN(), +} +var asinSC = []float64{ + NaN(), + Copysign(0, -1), + 0, + NaN(), + NaN(), +} + +var vfasinhSC = []float64{ + Inf(-1), + Copysign(0, -1), + 0, + Inf(1), + NaN(), +} +var asinhSC = []float64{ + Inf(-1), + Copysign(0, -1), + 0, + Inf(1), + NaN(), +} + +var vfatanSC = []float64{ + Inf(-1), + Copysign(0, -1), + 0, + Inf(1), + NaN(), +} +var atanSC = []float64{ + -Pi / 2, + Copysign(0, -1), + 0, + Pi / 2, + NaN(), +} + +var vfatanhSC = []float64{ + Inf(-1), + -Pi, + -1, + Copysign(0, -1), + 0, + 1, + Pi, + Inf(1), + NaN(), +} +var atanhSC = []float64{ + NaN(), + NaN(), + Inf(-1), + Copysign(0, -1), + 0, + Inf(1), + NaN(), + NaN(), + NaN(), +} +var vfatan2SC = [][2]float64{ + {Inf(-1), Inf(-1)}, + {Inf(-1), -Pi}, + {Inf(-1), 0}, + {Inf(-1), +Pi}, + {Inf(-1), Inf(1)}, + {Inf(-1), NaN()}, + {-Pi, Inf(-1)}, + {-Pi, 0}, + {-Pi, Inf(1)}, + {-Pi, NaN()}, + {Copysign(0, -1), Inf(-1)}, + {Copysign(0, -1), -Pi}, + {Copysign(0, -1), Copysign(0, -1)}, + {Copysign(0, -1), 0}, + {Copysign(0, -1), +Pi}, + {Copysign(0, -1), Inf(1)}, + {Copysign(0, -1), NaN()}, + {0, Inf(-1)}, + {0, -Pi}, + {0, Copysign(0, -1)}, + {0, 0}, + {0, +Pi}, + {0, Inf(1)}, + {0, NaN()}, + {+Pi, Inf(-1)}, + {+Pi, 0}, + {+Pi, Inf(1)}, + {+Pi, NaN()}, + {Inf(1), Inf(-1)}, + {Inf(1), -Pi}, + {Inf(1), 0}, + {Inf(1), +Pi}, + {Inf(1), Inf(1)}, + {Inf(1), NaN()}, + {NaN(), NaN()}, +} +var atan2SC = []float64{ + -3 * Pi / 4, // atan2(-Inf, -Inf) + -Pi / 2, // atan2(-Inf, -Pi) + -Pi / 2, // atan2(-Inf, +0) + -Pi / 2, // atan2(-Inf, +Pi) + -Pi / 4, // atan2(-Inf, +Inf) + NaN(), // atan2(-Inf, NaN) + -Pi, // atan2(-Pi, -Inf) + -Pi / 2, // atan2(-Pi, +0) + Copysign(0, -1), // atan2(-Pi, Inf) + NaN(), // atan2(-Pi, NaN) + -Pi, // atan2(-0, -Inf) + -Pi, // atan2(-0, -Pi) + -Pi, // atan2(-0, -0) + Copysign(0, -1), // atan2(-0, +0) + Copysign(0, -1), // atan2(-0, +Pi) + Copysign(0, -1), // atan2(-0, +Inf) + NaN(), // atan2(-0, NaN) + Pi, // atan2(+0, -Inf) + Pi, // atan2(+0, -Pi) + Pi, // atan2(+0, -0) + 0, // atan2(+0, +0) + 0, // atan2(+0, +Pi) + 0, // atan2(+0, +Inf) + NaN(), // atan2(+0, NaN) + Pi, // atan2(+Pi, -Inf) + Pi / 2, // atan2(+Pi, +0) + 0, // atan2(+Pi, +Inf) + NaN(), // atan2(+Pi, NaN) + 3 * Pi / 4, // atan2(+Inf, -Inf) + Pi / 2, // atan2(+Inf, -Pi) + Pi / 2, // atan2(+Inf, +0) + Pi / 2, // atan2(+Inf, +Pi) + Pi / 4, // atan2(+Inf, +Inf) + NaN(), // atan2(+Inf, NaN) + NaN(), // atan2(NaN, NaN) +} + +var vfcbrtSC = []float64{ + Inf(-1), + Copysign(0, -1), + 0, + Inf(1), + NaN(), +} +var cbrtSC = []float64{ + Inf(-1), + Copysign(0, -1), + 0, + Inf(1), + NaN(), +} + +var vfceilSC = []float64{ + Inf(-1), + Copysign(0, -1), + 0, + Inf(1), + NaN(), +} +var ceilSC = []float64{ + Inf(-1), + Copysign(0, -1), + 0, + Inf(1), + NaN(), +} + +var vfcopysignSC = []float64{ + Inf(-1), + Inf(1), + NaN(), +} +var copysignSC = []float64{ + Inf(-1), + Inf(-1), + NaN(), +} + +var vfcosSC = []float64{ + Inf(-1), + Inf(1), + NaN(), +} +var cosSC = []float64{ + NaN(), + NaN(), + NaN(), +} + +var vfcoshSC = []float64{ + Inf(-1), + Copysign(0, -1), + 0, + Inf(1), + NaN(), +} +var coshSC = []float64{ + Inf(1), + 1, + 1, + Inf(1), + NaN(), +} + +var vferfSC = []float64{ + Inf(-1), + Copysign(0, -1), + 0, + Inf(1), + NaN(), +} +var erfSC = []float64{ + -1, + Copysign(0, -1), + 0, + 1, + NaN(), +} + +var vferfcSC = []float64{ + Inf(-1), + Inf(1), + NaN(), +} +var erfcSC = []float64{ + 2, + 0, + NaN(), +} + +var vfexpSC = []float64{ + Inf(-1), + -2000, + 2000, + Inf(1), + NaN(), +} +var expSC = []float64{ + 0, + 0, + Inf(1), + Inf(1), + NaN(), +} + +var vfexpm1SC = []float64{ + Inf(-1), + Copysign(0, -1), + 0, + Inf(1), + NaN(), +} +var expm1SC = []float64{ + -1, + Copysign(0, -1), + 0, + Inf(1), + NaN(), +} + +var vffabsSC = []float64{ + Inf(-1), + Copysign(0, -1), + 0, + Inf(1), + NaN(), +} +var fabsSC = []float64{ + Inf(1), + 0, + 0, + Inf(1), + NaN(), +} + +var vffmodSC = [][2]float64{ + {Inf(-1), Inf(-1)}, + {Inf(-1), -Pi}, + {Inf(-1), 0}, + {Inf(-1), Pi}, + {Inf(-1), Inf(1)}, + {Inf(-1), NaN()}, + {-Pi, Inf(-1)}, + {-Pi, 0}, + {-Pi, Inf(1)}, + {-Pi, NaN()}, + {Copysign(0, -1), Inf(-1)}, + {Copysign(0, -1), 0}, + {Copysign(0, -1), Inf(1)}, + {Copysign(0, -1), NaN()}, + {0, Inf(-1)}, + {0, 0}, + {0, Inf(1)}, + {0, NaN()}, + {Pi, Inf(-1)}, + {Pi, 0}, + {Pi, Inf(1)}, + {Pi, NaN()}, + {Inf(1), Inf(-1)}, + {Inf(1), -Pi}, + {Inf(1), 0}, + {Inf(1), Pi}, + {Inf(1), Inf(1)}, + {Inf(1), NaN()}, + {NaN(), Inf(-1)}, + {NaN(), -Pi}, + {NaN(), 0}, + {NaN(), Pi}, + {NaN(), Inf(1)}, + {NaN(), NaN()}, +} +var fmodSC = []float64{ + NaN(), // fmod(-Inf, -Inf) + NaN(), // fmod(-Inf, -Pi) + NaN(), // fmod(-Inf, 0) + NaN(), // fmod(-Inf, Pi) + NaN(), // fmod(-Inf, +Inf) + NaN(), // fmod(-Inf, NaN) + -Pi, // fmod(-Pi, -Inf) + NaN(), // fmod(-Pi, 0) + -Pi, // fmod(-Pi, +Inf) + NaN(), // fmod(-Pi, NaN) + Copysign(0, -1), // fmod(-0, -Inf) + NaN(), // fmod(-0, 0) + Copysign(0, -1), // fmod(-0, Inf) + NaN(), // fmod(-0, NaN) + 0, // fmod(0, -Inf) + NaN(), // fmod(0, 0) + 0, // fmod(0, +Inf) + NaN(), // fmod(0, NaN) + Pi, // fmod(Pi, -Inf) + NaN(), // fmod(Pi, 0) + Pi, // fmod(Pi, +Inf) + NaN(), // fmod(Pi, NaN) + NaN(), // fmod(+Inf, -Inf) + NaN(), // fmod(+Inf, -Pi) + NaN(), // fmod(+Inf, 0) + NaN(), // fmod(+Inf, Pi) + NaN(), // fmod(+Inf, +Inf) + NaN(), // fmod(+Inf, NaN) + NaN(), // fmod(NaN, -Inf) + NaN(), // fmod(NaN, -Pi) + NaN(), // fmod(NaN, 0) + NaN(), // fmod(NaN, Pi) + NaN(), // fmod(NaN, +Inf) + NaN(), // fmod(NaN, NaN) +} + +var vffrexpSC = []float64{ + Inf(-1), + Copysign(0, -1), + 0, + Inf(1), + NaN(), +} +var frexpSC = []fi{ + {Inf(-1), 0}, + {Copysign(0, -1), 0}, + {0, 0}, + {Inf(1), 0}, + {NaN(), 0}, +} + +var vfgammaSC = []float64{ + Inf(-1), + -3, + Copysign(0, -1), + 0, + Inf(1), + NaN(), +} +var gammaSC = []float64{ + Inf(-1), + Inf(1), + Inf(1), + Inf(1), + Inf(1), + NaN(), +} + +var vfhypotSC = [][2]float64{ + {Inf(-1), Inf(-1)}, + {Inf(-1), 0}, + {Inf(-1), Inf(1)}, + {Inf(-1), NaN()}, + {Copysign(0, -1), Copysign(0, -1)}, + {Copysign(0, -1), 0}, + {0, Copysign(0, -1)}, + {0, 0}, // +0, +0 + {0, Inf(-1)}, + {0, Inf(1)}, + {0, NaN()}, + {Inf(1), Inf(-1)}, + {Inf(1), 0}, + {Inf(1), Inf(1)}, + {Inf(1), NaN()}, + {NaN(), Inf(-1)}, + {NaN(), 0}, + {NaN(), Inf(1)}, + {NaN(), NaN()}, +} +var hypotSC = []float64{ + Inf(1), + Inf(1), + Inf(1), + Inf(1), + 0, + 0, + 0, + 0, + Inf(1), + Inf(1), + NaN(), + Inf(1), + Inf(1), + Inf(1), + Inf(1), + Inf(1), + NaN(), + Inf(1), + NaN(), +} + +var vfilogbSC = []float64{ + Inf(-1), + 0, + Inf(1), + NaN(), +} +var ilogbSC = []int{ + MaxInt32, + MinInt32, + MaxInt32, + MaxInt32, +} + +var vfj0SC = []float64{ + Inf(-1), + 0, + Inf(1), + NaN(), +} +var j0SC = []float64{ + 0, + 1, + 0, + NaN(), +} +var j1SC = []float64{ + 0, + 0, + 0, + NaN(), +} +var j2SC = []float64{ + 0, + 0, + 0, + NaN(), +} +var jM3SC = []float64{ + 0, + 0, + 0, + NaN(), +} + +var vfldexpSC = []fi{ + {0, 0}, + {0, -1075}, + {0, 1024}, + {Copysign(0, -1), 0}, + {Copysign(0, -1), -1075}, + {Copysign(0, -1), 1024}, + {Inf(1), 0}, + {Inf(1), -1024}, + {Inf(-1), 0}, + {Inf(-1), -1024}, + {NaN(), -1024}, +} +var ldexpSC = []float64{ + 0, + 0, + 0, + Copysign(0, -1), + Copysign(0, -1), + Copysign(0, -1), + Inf(1), + Inf(1), + Inf(-1), + Inf(-1), + NaN(), +} + +var vflgammaSC = []float64{ + Inf(-1), + -3, + 0, + 1, + 2, + Inf(1), + NaN(), +} +var lgammaSC = []fi{ + {Inf(-1), 1}, + {Inf(1), 1}, + {Inf(1), 1}, + {0, 1}, + {0, 1}, + {Inf(1), 1}, + {NaN(), 1}, +} + +var vflogSC = []float64{ + Inf(-1), + -Pi, + Copysign(0, -1), + 0, + 1, + Inf(1), + NaN(), +} +var logSC = []float64{ + NaN(), + NaN(), + Inf(-1), + Inf(-1), + 0, + Inf(1), + NaN(), +} + +var vflogbSC = []float64{ + Inf(-1), + 0, + Inf(1), + NaN(), +} +var logbSC = []float64{ + Inf(1), + Inf(-1), + Inf(1), + NaN(), +} + +var vflog1pSC = []float64{ + Inf(-1), + -Pi, + -1, + Copysign(0, -1), + 0, + Inf(1), + NaN(), +} +var log1pSC = []float64{ + NaN(), + NaN(), + Inf(-1), + Copysign(0, -1), + 0, + Inf(1), + NaN(), +} + +var vfmodfSC = []float64{ + Inf(-1), + Inf(1), + NaN(), +} +var modfSC = [][2]float64{ + {Inf(-1), NaN()}, // [2]float64{Copysign(0, -1), Inf(-1)}, + {Inf(1), NaN()}, // [2]float64{0, Inf(1)}, + {NaN(), NaN()}, +} + +var vfnextafterSC = [][2]float64{ + {0, NaN()}, + {NaN(), 0}, + {NaN(), NaN()}, +} +var nextafterSC = []float64{ + NaN(), + NaN(), + NaN(), +} + +var vfpowSC = [][2]float64{ + {Inf(-1), -Pi}, + {Inf(-1), -3}, + {Inf(-1), Copysign(0, -1)}, + {Inf(-1), 0}, + {Inf(-1), 1}, + {Inf(-1), 3}, + {Inf(-1), Pi}, + {Inf(-1), NaN()}, + + {-Pi, Inf(-1)}, + {-Pi, -Pi}, + {-Pi, Copysign(0, -1)}, + {-Pi, 0}, + {-Pi, 1}, + {-Pi, Pi}, + {-Pi, Inf(1)}, + {-Pi, NaN()}, + + {-1, Inf(-1)}, + {-1, Inf(1)}, + {-1, NaN()}, + {-1 / 2, Inf(-1)}, + {-1 / 2, Inf(1)}, + {Copysign(0, -1), Inf(-1)}, + {Copysign(0, -1), -Pi}, + {Copysign(0, -1), -3}, + {Copysign(0, -1), 3}, + {Copysign(0, -1), Pi}, + {Copysign(0, -1), Inf(1)}, + + {0, Inf(-1)}, + {0, -Pi}, + {0, -3}, + {0, Copysign(0, -1)}, + {0, 0}, + {0, 3}, + {0, Pi}, + {0, Inf(1)}, + {0, NaN()}, + + {1 / 2, Inf(-1)}, + {1 / 2, Inf(1)}, + {1, Inf(-1)}, + {1, Inf(1)}, + {1, NaN()}, + + {Pi, Inf(-1)}, + {Pi, Copysign(0, -1)}, + {Pi, 0}, + {Pi, 1}, + {Pi, Inf(1)}, + {Pi, NaN()}, + {Inf(1), -Pi}, + {Inf(1), Copysign(0, -1)}, + {Inf(1), 0}, + {Inf(1), 1}, + {Inf(1), Pi}, + {Inf(1), NaN()}, + {NaN(), -Pi}, + {NaN(), Copysign(0, -1)}, + {NaN(), 0}, + {NaN(), 1}, + {NaN(), Pi}, + {NaN(), NaN()}, +} +var powSC = []float64{ + 0, // pow(-Inf, -Pi) + Copysign(0, -1), // pow(-Inf, -3) + 1, // pow(-Inf, -0) + 1, // pow(-Inf, +0) + Inf(-1), // pow(-Inf, 1) + Inf(-1), // pow(-Inf, 3) + Inf(1), // pow(-Inf, Pi) + NaN(), // pow(-Inf, NaN) + 0, // pow(-Pi, -Inf) + NaN(), // pow(-Pi, -Pi) + 1, // pow(-Pi, -0) + 1, // pow(-Pi, +0) + -Pi, // pow(-Pi, 1) + NaN(), // pow(-Pi, Pi) + Inf(1), // pow(-Pi, +Inf) + NaN(), // pow(-Pi, NaN) + 1, // pow(-1, -Inf) IEEE 754-2008 + 1, // pow(-1, +Inf) IEEE 754-2008 + NaN(), // pow(-1, NaN) + Inf(1), // pow(-1/2, -Inf) + 0, // pow(-1/2, +Inf) + Inf(1), // pow(-0, -Inf) + Inf(1), // pow(-0, -Pi) + Inf(-1), // pow(-0, -3) IEEE 754-2008 + Copysign(0, -1), // pow(-0, 3) IEEE 754-2008 + 0, // pow(-0, +Pi) + 0, // pow(-0, +Inf) + Inf(1), // pow(+0, -Inf) + Inf(1), // pow(+0, -Pi) + Inf(1), // pow(+0, -3) + 1, // pow(+0, -0) + 1, // pow(+0, +0) + 0, // pow(+0, 3) + 0, // pow(+0, +Pi) + 0, // pow(+0, +Inf) + NaN(), // pow(+0, NaN) + Inf(1), // pow(1/2, -Inf) + 0, // pow(1/2, +Inf) + 1, // pow(1, -Inf) IEEE 754-2008 + 1, // pow(1, +Inf) IEEE 754-2008 + 1, // pow(1, NaN) IEEE 754-2008 + 0, // pow(+Pi, -Inf) + 1, // pow(+Pi, -0) + 1, // pow(+Pi, +0) + Pi, // pow(+Pi, 1) + Inf(1), // pow(+Pi, +Inf) + NaN(), // pow(+Pi, NaN) + 0, // pow(+Inf, -Pi) + 1, // pow(+Inf, -0) + 1, // pow(+Inf, +0) + Inf(1), // pow(+Inf, 1) + Inf(1), // pow(+Inf, Pi) + NaN(), // pow(+Inf, NaN) + NaN(), // pow(NaN, -Pi) + 1, // pow(NaN, -0) + 1, // pow(NaN, +0) + NaN(), // pow(NaN, 1) + NaN(), // pow(NaN, +Pi) + NaN(), // pow(NaN, NaN) +} + +var vfsignbitSC = []float64{ + Inf(-1), + Copysign(0, -1), + 0, + Inf(1), + NaN(), +} +var signbitSC = []bool{ + true, + true, + false, + false, + false, +} + +var vfsinSC = []float64{ + Inf(-1), + Copysign(0, -1), + 0, + Inf(1), + NaN(), +} +var sinSC = []float64{ + NaN(), + Copysign(0, -1), + 0, + NaN(), + NaN(), +} + +var vfsinhSC = []float64{ + Inf(-1), + Copysign(0, -1), + 0, + Inf(1), + NaN(), +} +var sinhSC = []float64{ + Inf(-1), + Copysign(0, -1), + 0, + Inf(1), + NaN(), +} + +var vfsqrtSC = []float64{ + Inf(-1), + -Pi, + Copysign(0, -1), + 0, + Inf(1), + NaN(), +} +var sqrtSC = []float64{ + NaN(), + NaN(), + Copysign(0, -1), + 0, + Inf(1), + NaN(), +} + +var vftanhSC = []float64{ + Inf(-1), + Copysign(0, -1), + 0, + Inf(1), + NaN(), +} +var tanhSC = []float64{ + -1, + Copysign(0, -1), + 0, + 1, + NaN(), +} + +var vfy0SC = []float64{ + Inf(-1), + 0, + Inf(1), + NaN(), +} +var y0SC = []float64{ + NaN(), + Inf(-1), + 0, + NaN(), +} +var y1SC = []float64{ + NaN(), + Inf(-1), + 0, + NaN(), +} +var y2SC = []float64{ + NaN(), + Inf(-1), + 0, + NaN(), +} +var yM3SC = []float64{ + NaN(), + Inf(1), + 0, + NaN(), +} + +// arguments and expected results for boundary cases +const ( + SmallestNormalFloat64 = 2.2250738585072014e-308 // 2**-1022 + LargestSubnormalFloat64 = SmallestNormalFloat64 - SmallestNonzeroFloat64 +) + +var vffrexpBC = []float64{ + SmallestNormalFloat64, + LargestSubnormalFloat64, + SmallestNonzeroFloat64, + MaxFloat64, + -SmallestNormalFloat64, + -LargestSubnormalFloat64, + -SmallestNonzeroFloat64, + -MaxFloat64, +} +var frexpBC = []fi{ + {0.5, -1021}, + {0.99999999999999978, -1022}, + {0.5, -1073}, + {0.99999999999999989, 1024}, + {-0.5, -1021}, + {-0.99999999999999978, -1022}, + {-0.5, -1073}, + {-0.99999999999999989, 1024}, +} + +var vfldexpBC = []fi{ + {SmallestNormalFloat64, -52}, + {LargestSubnormalFloat64, -51}, + {SmallestNonzeroFloat64, 1074}, + {MaxFloat64, -(1023 + 1074)}, + {1, -1075}, + {-1, -1075}, + {1, 1024}, + {-1, 1024}, +} +var ldexpBC = []float64{ + SmallestNonzeroFloat64, + 1e-323, // 2**-1073 + 1, + 1e-323, // 2**-1073 + 0, + Copysign(0, -1), + Inf(1), + Inf(-1), +} + +var logbBC = []float64{ + -1022, + -1023, + -1074, + 1023, + -1022, + -1023, + -1074, + 1023, +} + +func tolerance(a, b, e float64) bool { + d := a - b + if d < 0 { + d = -d + } + + if a != 0 { + e = e * a + if e < 0 { + e = -e + } + } + return d < e +} +func kindaclose(a, b float64) bool { return tolerance(a, b, 1e-8) } +func close(a, b float64) bool { return tolerance(a, b, 1e-14) } +func veryclose(a, b float64) bool { return tolerance(a, b, 4e-16) } +func soclose(a, b, e float64) bool { return tolerance(a, b, e) } +func alike(a, b float64) bool { + switch { + case IsNaN(a) && IsNaN(b): + return true + case a == b: + return Signbit(a) == Signbit(b) + } + return false +} + +func TestAcos(t *testing.T) { + for i := 0; i < len(vf); i++ { + a := vf[i] / 10 + if f := Acos(a); !close(acos[i], f) { + t.Errorf("Acos(%g) = %g, want %g", a, f, acos[i]) + } + } + for i := 0; i < len(vfacosSC); i++ { + if f := Acos(vfacosSC[i]); !alike(acosSC[i], f) { + t.Errorf("Acos(%g) = %g, want %g", vfacosSC[i], f, acosSC[i]) + } + } +} + +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", 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", vfacoshSC[i], f, acoshSC[i]) + } + } +} + +func TestAsin(t *testing.T) { + for i := 0; i < len(vf); i++ { + a := vf[i] / 10 + if f := Asin(a); !veryclose(asin[i], f) { + t.Errorf("Asin(%g) = %g, want %g", a, f, asin[i]) + } + } + for i := 0; i < len(vfasinSC); i++ { + if f := Asin(vfasinSC[i]); !alike(asinSC[i], f) { + t.Errorf("Asin(%g) = %g, want %g", vfasinSC[i], f, asinSC[i]) + } + } +} + +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", 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", 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) { + t.Errorf("Atan(%g) = %g, want %g", vf[i], f, atan[i]) + } + } + for i := 0; i < len(vfatanSC); i++ { + if f := Atan(vfatanSC[i]); !alike(atanSC[i], f) { + t.Errorf("Atan(%g) = %g, want %g", vfatanSC[i], f, atanSC[i]) + } + } +} + +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", 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", vfatanhSC[i], f, atanhSC[i]) + } + } +} + +func TestAtan2(t *testing.T) { + for i := 0; i < len(vf); i++ { + if f := Atan2(10, vf[i]); !veryclose(atan2[i], f) { + t.Errorf("Atan2(10, %g) = %g, want %g", vf[i], f, atan2[i]) + } + } + for i := 0; i < len(vfatan2SC); i++ { + if f := Atan2(vfatan2SC[i][0], vfatan2SC[i][1]); !alike(atan2SC[i], f) { + t.Errorf("Atan2(%g, %g) = %g, want %g", vfatan2SC[i][0], vfatan2SC[i][1], f, atan2SC[i]) + } + } +} + +func TestCbrt(t *testing.T) { + for i := 0; i < len(vf); i++ { + if f := Cbrt(vf[i]); !veryclose(cbrt[i], f) { + t.Errorf("Cbrt(%g) = %g, want %g", vf[i], f, cbrt[i]) + } + } + for i := 0; i < len(vfcbrtSC); i++ { + if f := Cbrt(vfcbrtSC[i]); !alike(cbrtSC[i], f) { + t.Errorf("Cbrt(%g) = %g, want %g", vfcbrtSC[i], f, cbrtSC[i]) + } + } +} + +func TestCeil(t *testing.T) { + for i := 0; i < len(vf); i++ { + if f := Ceil(vf[i]); ceil[i] != f { + t.Errorf("Ceil(%g) = %g, want %g", vf[i], f, ceil[i]) + } + } + for i := 0; i < len(vfceilSC); i++ { + if f := Ceil(vfceilSC[i]); !alike(ceilSC[i], f) { + t.Errorf("Ceil(%g) = %g, want %g", vfceilSC[i], f, ceilSC[i]) + } + } +} + +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", vf[i], f, copysign[i]) + } + } + for i := 0; i < len(vf); i++ { + if f := Copysign(vf[i], 1); -copysign[i] != f { + t.Errorf("Copysign(%g, 1) = %g, want %g", 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", 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", vf[i], f, cos[i]) + } + } + for i := 0; i < len(vfcosSC); i++ { + if f := Cos(vfcosSC[i]); !alike(cosSC[i], f) { + t.Errorf("Cos(%g) = %g, want %g", vfcosSC[i], f, cosSC[i]) + } + } +} + +func TestCosh(t *testing.T) { + for i := 0; i < len(vf); i++ { + if f := Cosh(vf[i]); !close(cosh[i], f) { + t.Errorf("Cosh(%g) = %g, want %g", vf[i], f, cosh[i]) + } + } + for i := 0; i < len(vfcoshSC); i++ { + if f := Cosh(vfcoshSC[i]); !alike(coshSC[i], f) { + t.Errorf("Cosh(%g) = %g, want %g", vfcoshSC[i], f, coshSC[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", 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", 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", a, f, erfc[i]) + } + } + for i := 0; i < len(vferfcSC); i++ { + if f := Erfc(vferfcSC[i]); !alike(erfcSC[i], f) { + t.Errorf("Erfc(%g) = %g, want %g", vferfcSC[i], f, erfcSC[i]) + } + } +} + +func TestExp(t *testing.T) { + testExp(t, Exp, "Exp") + testExp(t, ExpGo, "ExpGo") +} + +func testExp(t *testing.T, Exp func(float64) float64, name string) { + for i := 0; i < len(vf); i++ { + if f := Exp(vf[i]); !close(exp[i], f) { + t.Errorf("%s(%g) = %g, want %g", name, vf[i], f, exp[i]) + } + } + for i := 0; i < len(vfexpSC); i++ { + if f := Exp(vfexpSC[i]); !alike(expSC[i], f) { + t.Errorf("%s(%g) = %g, want %g", name, vfexpSC[i], f, expSC[i]) + } + } +} + +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(%g) = %g, want %g", a, f, expm1[i]) + } + } + for i := 0; i < len(vfexpm1SC); i++ { + if f := Expm1(vfexpm1SC[i]); !alike(expm1SC[i], f) { + t.Errorf("Expm1(%g) = %g, want %g", vfexpm1SC[i], f, expm1SC[i]) + } + } +} + +func TestExp2(t *testing.T) { + testExp2(t, Exp2, "Exp2") + testExp2(t, Exp2Go, "Exp2Go") +} + +func testExp2(t *testing.T, Exp2 func(float64) float64, name string) { + for i := 0; i < len(vf); i++ { + if f := Exp2(vf[i]); !close(exp2[i], f) { + t.Errorf("%s(%g) = %g, want %g", name, vf[i], f, exp2[i]) + } + } + for i := 0; i < len(vfexpSC); i++ { + if f := Exp2(vfexpSC[i]); !alike(expSC[i], f) { + t.Errorf("%s(%g) = %g, want %g", name, vfexpSC[i], f, expSC[i]) + } + } + for n := -1074; n < 1024; n++ { + f := Exp2(float64(n)) + vf := Ldexp(1, n) + if f != vf { + t.Errorf("%s(%d) = %g, want %g", name, n, f, vf) + } + } +} + +func TestFabs(t *testing.T) { + for i := 0; i < len(vf); i++ { + if f := Fabs(vf[i]); fabs[i] != f { + t.Errorf("Fabs(%g) = %g, want %g", vf[i], f, fabs[i]) + } + } + for i := 0; i < len(vffabsSC); i++ { + if f := Fabs(vffabsSC[i]); !alike(fabsSC[i], f) { + t.Errorf("Fabs(%g) = %g, want %g", vffabsSC[i], f, fabsSC[i]) + } + } +} + +func TestFdim(t *testing.T) { + for i := 0; i < len(vf); i++ { + if f := Fdim(vf[i], 0); fdim[i] != f { + t.Errorf("Fdim(%g, %g) = %g, want %g", vf[i], 0.0, f, fdim[i]) + } + } +} + +func TestFloor(t *testing.T) { + for i := 0; i < len(vf); i++ { + if f := Floor(vf[i]); floor[i] != f { + t.Errorf("Floor(%g) = %g, want %g", vf[i], f, floor[i]) + } + } + for i := 0; i < len(vfceilSC); i++ { + if f := Floor(vfceilSC[i]); !alike(ceilSC[i], f) { + t.Errorf("Floor(%g) = %g, want %g", vfceilSC[i], f, ceilSC[i]) + } + } +} + +func TestFmax(t *testing.T) { + for i := 0; i < len(vf); i++ { + if f := Fmax(vf[i], ceil[i]); ceil[i] != f { + t.Errorf("Fmax(%g, %g) = %g, want %g", vf[i], ceil[i], f, ceil[i]) + } + } +} + +func TestFmin(t *testing.T) { + for i := 0; i < len(vf); i++ { + if f := Fmin(vf[i], floor[i]); floor[i] != f { + t.Errorf("Fmin(%g, %g) = %g, want %g", vf[i], floor[i], f, floor[i]) + } + } +} + +func TestFmod(t *testing.T) { + for i := 0; i < len(vf); i++ { + if f := Fmod(10, vf[i]); fmod[i] != f { + t.Errorf("Fmod(10, %g) = %g, want %g", vf[i], f, fmod[i]) + } + } + for i := 0; i < len(vffmodSC); i++ { + if f := Fmod(vffmodSC[i][0], vffmodSC[i][1]); !alike(fmodSC[i], f) { + t.Errorf("Fmod(%g, %g) = %g, want %g", vffmodSC[i][0], vffmodSC[i][1], f, fmodSC[i]) + } + } +} + +func TestFrexp(t *testing.T) { + for i := 0; i < len(vf); i++ { + if f, j := Frexp(vf[i]); !veryclose(frexp[i].f, f) || frexp[i].i != j { + t.Errorf("Frexp(%g) = %g, %d, want %g, %d", vf[i], f, j, frexp[i].f, frexp[i].i) + } + } + for i := 0; i < len(vffrexpSC); i++ { + if f, j := Frexp(vffrexpSC[i]); !alike(frexpSC[i].f, f) || frexpSC[i].i != j { + t.Errorf("Frexp(%g) = %g, %d, want %g, %d", vffrexpSC[i], f, j, frexpSC[i].f, frexpSC[i].i) + } + } + for i := 0; i < len(vffrexpBC); i++ { + if f, j := Frexp(vffrexpBC[i]); !alike(frexpBC[i].f, f) || frexpBC[i].i != j { + t.Errorf("Frexp(%g) = %g, %d, want %g, %d", vffrexpBC[i], f, j, frexpBC[i].f, frexpBC[i].i) + } + } +} + +func TestGamma(t *testing.T) { + for i := 0; i < len(vf); i++ { + if f := Gamma(vf[i]); !close(gamma[i], f) { + t.Errorf("Gamma(%g) = %g, want %g", vf[i], f, gamma[i]) + } + } + for i := 0; i < len(vfgammaSC); i++ { + if f := Gamma(vfgammaSC[i]); !alike(gammaSC[i], f) { + t.Errorf("Gamma(%g) = %g, want %g", vfgammaSC[i], f, gammaSC[i]) + } + } +} + +func TestHypot(t *testing.T) { + for i := 0; i < len(vf); i++ { + a := Fabs(1e200 * tanh[i] * Sqrt(2)) + if f := Hypot(1e200*tanh[i], 1e200*tanh[i]); !veryclose(a, f) { + t.Errorf("Hypot(%g, %g) = %g, want %g", 1e200*tanh[i], 1e200*tanh[i], f, a) + } + } + for i := 0; i < len(vfhypotSC); i++ { + if f := Hypot(vfhypotSC[i][0], vfhypotSC[i][1]); !alike(hypotSC[i], f) { + t.Errorf("Hypot(%g, %g) = %g, want %g", vfhypotSC[i][0], vfhypotSC[i][1], f, hypotSC[i]) + } + } +} + +func TestIlogb(t *testing.T) { + for i := 0; i < len(vf); i++ { + a := frexp[i].i - 1 // adjust because fr in the interval [½, 1) + if e := Ilogb(vf[i]); a != e { + t.Errorf("Ilogb(%g) = %d, want %d", vf[i], e, a) + } + } + for i := 0; i < len(vflogbSC); i++ { + if e := Ilogb(vflogbSC[i]); ilogbSC[i] != e { + t.Errorf("Ilogb(%g) = %d, want %d", vflogbSC[i], e, ilogbSC[i]) + } + } + for i := 0; i < len(vffrexpBC); i++ { + if e := Ilogb(vffrexpBC[i]); int(logbBC[i]) != e { + t.Errorf("Ilogb(%g) = %d, want %d", vffrexpBC[i], e, int(logbBC[i])) + } + } +} + +func TestJ0(t *testing.T) { + for i := 0; i < len(vf); i++ { + if f := J0(vf[i]); !soclose(j0[i], f, 4e-14) { + t.Errorf("J0(%g) = %g, want %g", vf[i], f, j0[i]) + } + } + for i := 0; i < len(vfj0SC); i++ { + if f := J0(vfj0SC[i]); !alike(j0SC[i], f) { + t.Errorf("J0(%g) = %g, want %g", vfj0SC[i], f, j0SC[i]) + } + } +} + +func TestJ1(t *testing.T) { + for i := 0; i < len(vf); i++ { + if f := J1(vf[i]); !close(j1[i], f) { + t.Errorf("J1(%g) = %g, want %g", vf[i], f, j1[i]) + } + } + for i := 0; i < len(vfj0SC); i++ { + if f := J1(vfj0SC[i]); !alike(j1SC[i], f) { + t.Errorf("J1(%g) = %g, want %g", vfj0SC[i], f, j1SC[i]) + } + } +} + +func TestJn(t *testing.T) { + for i := 0; i < len(vf); i++ { + if f := Jn(2, vf[i]); !close(j2[i], f) { + t.Errorf("Jn(2, %g) = %g, want %g", vf[i], f, j2[i]) + } + if f := Jn(-3, vf[i]); !close(jM3[i], f) { + t.Errorf("Jn(-3, %g) = %g, want %g", vf[i], f, jM3[i]) + } + } + for i := 0; i < len(vfj0SC); i++ { + if f := Jn(2, vfj0SC[i]); !alike(j2SC[i], f) { + t.Errorf("Jn(2, %g) = %g, want %g", vfj0SC[i], f, j2SC[i]) + } + if f := Jn(-3, vfj0SC[i]); !alike(jM3SC[i], f) { + t.Errorf("Jn(-3, %g) = %g, want %g", vfj0SC[i], f, jM3SC[i]) + } + } +} + +func TestLdexp(t *testing.T) { + for i := 0; i < len(vf); i++ { + if f := Ldexp(frexp[i].f, frexp[i].i); !veryclose(vf[i], f) { + t.Errorf("Ldexp(%g, %d) = %g, want %g", frexp[i].f, frexp[i].i, f, vf[i]) + } + } + for i := 0; i < len(vffrexpSC); i++ { + if f := Ldexp(frexpSC[i].f, frexpSC[i].i); !alike(vffrexpSC[i], f) { + t.Errorf("Ldexp(%g, %d) = %g, want %g", frexpSC[i].f, frexpSC[i].i, f, vffrexpSC[i]) + } + } + for i := 0; i < len(vfldexpSC); i++ { + if f := Ldexp(vfldexpSC[i].f, vfldexpSC[i].i); !alike(ldexpSC[i], f) { + t.Errorf("Ldexp(%g, %d) = %g, want %g", vfldexpSC[i].f, vfldexpSC[i].i, f, ldexpSC[i]) + } + } + for i := 0; i < len(vffrexpBC); i++ { + if f := Ldexp(frexpBC[i].f, frexpBC[i].i); !alike(vffrexpBC[i], f) { + t.Errorf("Ldexp(%g, %d) = %g, want %g", frexpBC[i].f, frexpBC[i].i, f, vffrexpBC[i]) + } + } + for i := 0; i < len(vfldexpBC); i++ { + if f := Ldexp(vfldexpBC[i].f, vfldexpBC[i].i); !alike(ldexpBC[i], f) { + t.Errorf("Ldexp(%g, %d) = %g, want %g", vfldexpBC[i].f, vfldexpBC[i].i, f, ldexpBC[i]) + } + } +} + +func TestLgamma(t *testing.T) { + for i := 0; i < len(vf); i++ { + if f, s := Lgamma(vf[i]); !close(lgamma[i].f, f) || lgamma[i].i != s { + t.Errorf("Lgamma(%g) = %g, %d, want %g, %d", vf[i], f, s, lgamma[i].f, lgamma[i].i) + } + } + for i := 0; i < len(vflgammaSC); i++ { + if f, s := Lgamma(vflgammaSC[i]); !alike(lgammaSC[i].f, f) || lgammaSC[i].i != s { + t.Errorf("Lgamma(%g) = %g, %d, want %g, %d", vflgammaSC[i], f, s, lgammaSC[i].f, lgammaSC[i].i) + } + } +} + +func TestLog(t *testing.T) { + for i := 0; i < len(vf); i++ { + a := Fabs(vf[i]) + if f := Log(a); log[i] != f { + t.Errorf("Log(%g) = %g, want %g", a, f, log[i]) + } + } + if f := Log(10); f != Ln10 { + t.Errorf("Log(%g) = %g, want %g", 10.0, f, Ln10) + } + for i := 0; i < len(vflogSC); i++ { + if f := Log(vflogSC[i]); !alike(logSC[i], f) { + t.Errorf("Log(%g) = %g, want %g", vflogSC[i], f, logSC[i]) + } + } +} + +func TestLogb(t *testing.T) { + for i := 0; i < len(vf); i++ { + if f := Logb(vf[i]); logb[i] != f { + t.Errorf("Logb(%g) = %g, want %g", vf[i], f, logb[i]) + } + } + for i := 0; i < len(vflogbSC); i++ { + if f := Logb(vflogbSC[i]); !alike(logbSC[i], f) { + t.Errorf("Logb(%g) = %g, want %g", vflogbSC[i], f, logbSC[i]) + } + } + for i := 0; i < len(vffrexpBC); i++ { + if e := Logb(vffrexpBC[i]); !alike(logbBC[i], e) { + t.Errorf("Ilogb(%g) = %g, want %g", vffrexpBC[i], e, logbBC[i]) + } + } +} + +func TestLog10(t *testing.T) { + for i := 0; i < len(vf); i++ { + a := Fabs(vf[i]) + if f := Log10(a); !veryclose(log10[i], f) { + t.Errorf("Log10(%g) = %g, want %g", a, f, log10[i]) + } + } + if f := Log10(E); f != Log10E { + t.Errorf("Log10(%g) = %g, want %g", E, f, Log10E) + } + for i := 0; i < len(vflogSC); i++ { + if f := Log10(vflogSC[i]); !alike(logSC[i], f) { + t.Errorf("Log10(%g) = %g, want %g", vflogSC[i], f, logSC[i]) + } + } +} + +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", a, f, log1p[i]) + } + } + a := 9.0 + if f := Log1p(a); f != Ln10 { + t.Errorf("Log1p(%g) = %g, want %g", 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", vflog1pSC[i], f, log1pSC[i]) + } + } +} + +func TestLog2(t *testing.T) { + for i := 0; i < len(vf); i++ { + a := Fabs(vf[i]) + if f := Log2(a); !veryclose(log2[i], f) { + t.Errorf("Log2(%g) = %g, want %g", a, f, log2[i]) + } + } + if f := Log2(E); f != Log2E { + t.Errorf("Log2(%g) = %g, want %g", E, f, Log2E) + } + for i := 0; i < len(vflogSC); i++ { + if f := Log2(vflogSC[i]); !alike(logSC[i], f) { + t.Errorf("Log2(%g) = %g, want %g", vflogSC[i], f, logSC[i]) + } + } +} + +func TestModf(t *testing.T) { + for i := 0; i < len(vf); i++ { + if f, g := Modf(vf[i]); !veryclose(modf[i][0], f) || !veryclose(modf[i][1], g) { + t.Errorf("Modf(%g) = %g, %g, want %g, %g", vf[i], f, g, modf[i][0], modf[i][1]) + } + } + for i := 0; i < len(vfmodfSC); i++ { + if f, g := Modf(vfmodfSC[i]); !alike(modfSC[i][0], f) || !alike(modfSC[i][1], g) { + t.Errorf("Modf(%g) = %g, %g, want %g, %g", vfmodfSC[i], f, g, modfSC[i][0], modfSC[i][1]) + } + } +} + +func TestNextafter(t *testing.T) { + for i := 0; i < len(vf); i++ { + if f := Nextafter(vf[i], 10); nextafter[i] != f { + t.Errorf("Nextafter(%g, %g) = %g want %g", vf[i], 10.0, f, nextafter[i]) + } + } + for i := 0; i < len(vfmodfSC); i++ { + if f := Nextafter(vfnextafterSC[i][0], vfnextafterSC[i][1]); !alike(nextafterSC[i], f) { + t.Errorf("Nextafter(%g, %g) = %g want %g", vfnextafterSC[i][0], vfnextafterSC[i][1], f, nextafterSC[i]) + } + } +} + +func TestPow(t *testing.T) { + for i := 0; i < len(vf); i++ { + if f := Pow(10, vf[i]); !close(pow[i], f) { + t.Errorf("Pow(10, %g) = %g, want %g", vf[i], f, pow[i]) + } + } + for i := 0; i < len(vfpowSC); i++ { + if f := Pow(vfpowSC[i][0], vfpowSC[i][1]); !alike(powSC[i], f) { + t.Errorf("Pow(%g, %g) = %g, want %g", vfpowSC[i][0], vfpowSC[i][1], f, powSC[i]) + } + } +} + +func TestRemainder(t *testing.T) { + for i := 0; i < len(vf); i++ { + if f := Remainder(10, vf[i]); remainder[i] != f { + t.Errorf("Remainder(10, %g) = %g, want %g", vf[i], f, remainder[i]) + } + } + for i := 0; i < len(vffmodSC); i++ { + if f := Remainder(vffmodSC[i][0], vffmodSC[i][1]); !alike(fmodSC[i], f) { + t.Errorf("Remainder(%g, %g) = %g, want %g", vffmodSC[i][0], vffmodSC[i][1], f, fmodSC[i]) + } + } +} + +func TestSignbit(t *testing.T) { + for i := 0; i < len(vf); i++ { + if f := Signbit(vf[i]); signbit[i] != f { + t.Errorf("Signbit(%g) = %t, want %t", vf[i], f, signbit[i]) + } + } + for i := 0; i < len(vfsignbitSC); i++ { + if f := Signbit(vfsignbitSC[i]); signbitSC[i] != f { + t.Errorf("Signbit(%g) = %t, want %t", vfsignbitSC[i], f, signbitSC[i]) + } + } +} +func TestSin(t *testing.T) { + for i := 0; i < len(vf); i++ { + if f := Sin(vf[i]); !close(sin[i], f) { + t.Errorf("Sin(%g) = %g, want %g", vf[i], f, sin[i]) + } + } + for i := 0; i < len(vfsinSC); i++ { + if f := Sin(vfsinSC[i]); !alike(sinSC[i], f) { + t.Errorf("Sin(%g) = %g, want %g", vfsinSC[i], f, sinSC[i]) + } + } +} + +func TestSincos(t *testing.T) { + for i := 0; i < len(vf); i++ { + if s, c := Sincos(vf[i]); !close(sin[i], s) || !close(cos[i], c) { + t.Errorf("Sincos(%g) = %g, %g want %g, %g", vf[i], s, c, sin[i], cos[i]) + } + } +} + +func TestSinh(t *testing.T) { + for i := 0; i < len(vf); i++ { + if f := Sinh(vf[i]); !close(sinh[i], f) { + t.Errorf("Sinh(%g) = %g, want %g", vf[i], f, sinh[i]) + } + } + for i := 0; i < len(vfsinhSC); i++ { + if f := Sinh(vfsinhSC[i]); !alike(sinhSC[i], f) { + t.Errorf("Sinh(%g) = %g, want %g", vfsinhSC[i], f, sinhSC[i]) + } + } +} + +func TestSqrt(t *testing.T) { + for i := 0; i < len(vf); i++ { + a := Fabs(vf[i]) + if f := SqrtGo(a); sqrt[i] != f { + t.Errorf("SqrtGo(%g) = %g, want %g", a, f, sqrt[i]) + } + a = Fabs(vf[i]) + if f := Sqrt(a); sqrt[i] != f { + t.Errorf("Sqrt(%g) = %g, want %g", a, f, sqrt[i]) + } + } + for i := 0; i < len(vfsqrtSC); i++ { + if f := SqrtGo(vfsqrtSC[i]); !alike(sqrtSC[i], f) { + t.Errorf("SqrtGo(%g) = %g, want %g", vfsqrtSC[i], f, sqrtSC[i]) + } + if f := Sqrt(vfsqrtSC[i]); !alike(sqrtSC[i], f) { + t.Errorf("Sqrt(%g) = %g, want %g", vfsqrtSC[i], f, sqrtSC[i]) + } + } +} + +func TestTan(t *testing.T) { + for i := 0; i < len(vf); i++ { + if f := Tan(vf[i]); !close(tan[i], f) { + t.Errorf("Tan(%g) = %g, want %g", vf[i], f, tan[i]) + } + } + // same special cases as Sin + for i := 0; i < len(vfsinSC); i++ { + if f := Tan(vfsinSC[i]); !alike(sinSC[i], f) { + t.Errorf("Tan(%g) = %g, want %g", vfsinSC[i], f, sinSC[i]) + } + } + + // Make sure portable Tan(Pi/2) doesn't panic (it used to). + // The portable implementation returns NaN. + // Assembly implementations might not, + // because Pi/2 is not exactly representable. + if runtime.GOARCH != "386" { + if f := Tan(Pi / 2); !alike(f, NaN()) { + t.Errorf("Tan(%g) = %g, want %g", Pi/2, f, NaN()) + } + } +} + +func TestTanh(t *testing.T) { + for i := 0; i < len(vf); i++ { + if f := Tanh(vf[i]); !veryclose(tanh[i], f) { + t.Errorf("Tanh(%g) = %g, want %g", vf[i], f, tanh[i]) + } + } + for i := 0; i < len(vftanhSC); i++ { + if f := Tanh(vftanhSC[i]); !alike(tanhSC[i], f) { + t.Errorf("Tanh(%g) = %g, want %g", vftanhSC[i], f, tanhSC[i]) + } + } +} + +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", 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", vfceilSC[i], f, ceilSC[i]) + } + } +} + +func TestY0(t *testing.T) { + for i := 0; i < len(vf); i++ { + a := Fabs(vf[i]) + if f := Y0(a); !close(y0[i], f) { + t.Errorf("Y0(%g) = %g, want %g", a, f, y0[i]) + } + } + for i := 0; i < len(vfy0SC); i++ { + if f := Y0(vfy0SC[i]); !alike(y0SC[i], f) { + t.Errorf("Y0(%g) = %g, want %g", vfy0SC[i], f, y0SC[i]) + } + } +} + +func TestY1(t *testing.T) { + for i := 0; i < len(vf); i++ { + a := Fabs(vf[i]) + if f := Y1(a); !soclose(y1[i], f, 2e-14) { + t.Errorf("Y1(%g) = %g, want %g", a, f, y1[i]) + } + } + for i := 0; i < len(vfy0SC); i++ { + if f := Y1(vfy0SC[i]); !alike(y1SC[i], f) { + t.Errorf("Y1(%g) = %g, want %g", vfy0SC[i], f, y1SC[i]) + } + } +} + +func TestYn(t *testing.T) { + for i := 0; i < len(vf); i++ { + a := Fabs(vf[i]) + if f := Yn(2, a); !close(y2[i], f) { + t.Errorf("Yn(2, %g) = %g, want %g", a, f, y2[i]) + } + if f := Yn(-3, a); !close(yM3[i], f) { + t.Errorf("Yn(-3, %g) = %g, want %g", a, f, yM3[i]) + } + } + for i := 0; i < len(vfy0SC); i++ { + if f := Yn(2, vfy0SC[i]); !alike(y2SC[i], f) { + t.Errorf("Yn(2, %g) = %g, want %g", vfy0SC[i], f, y2SC[i]) + } + if f := Yn(-3, vfy0SC[i]); !alike(yM3SC[i], f) { + t.Errorf("Yn(-3, %g) = %g, want %g", vfy0SC[i], f, yM3SC[i]) + } + } +} + +// Check that math functions of high angle values +// return similar results to low angle values +func TestLargeCos(t *testing.T) { + large := float64(100000 * Pi) + for i := 0; i < len(vf); i++ { + f1 := Cos(vf[i]) + f2 := Cos(vf[i] + large) + if !kindaclose(f1, f2) { + t.Errorf("Cos(%g) = %g, want %g", vf[i]+large, f2, f1) + } + } +} + +func TestLargeSin(t *testing.T) { + large := float64(100000 * Pi) + for i := 0; i < len(vf); i++ { + f1 := Sin(vf[i]) + f2 := Sin(vf[i] + large) + if !kindaclose(f1, f2) { + t.Errorf("Sin(%g) = %g, want %g", vf[i]+large, f2, f1) + } + } +} + +func TestLargeSincos(t *testing.T) { + large := float64(100000 * Pi) + for i := 0; i < len(vf); i++ { + f1, g1 := Sincos(vf[i]) + f2, g2 := Sincos(vf[i] + large) + if !kindaclose(f1, f2) || !kindaclose(g1, g2) { + t.Errorf("Sincos(%g) = %g, %g, want %g, %g", vf[i]+large, f2, g2, f1, g1) + } + } +} + +func TestLargeTan(t *testing.T) { + large := float64(100000 * Pi) + for i := 0; i < len(vf); i++ { + f1 := Tan(vf[i]) + f2 := Tan(vf[i] + large) + if !kindaclose(f1, f2) { + t.Errorf("Tan(%g) = %g, want %g", vf[i]+large, f2, f1) + } + } +} + +// Check that math constants are accepted by compiler +// and have right value (assumes strconv.Atof works). +// http://code.google.com/p/go/issues/detail?id=201 + +type floatTest struct { + val interface{} + name string + str string +} + +var floatTests = []floatTest{ + {float64(MaxFloat64), "MaxFloat64", "1.7976931348623157e+308"}, + {float64(SmallestNonzeroFloat64), "SmallestNonzeroFloat64", "5e-324"}, + {float32(MaxFloat32), "MaxFloat32", "3.4028235e+38"}, + {float32(SmallestNonzeroFloat32), "SmallestNonzeroFloat32", "1e-45"}, +} + +func TestFloatMinMax(t *testing.T) { + for _, tt := range floatTests { + s := fmt.Sprint(tt.val) + if s != tt.str { + t.Errorf("Sprint(%v) = %s, want %s", tt.name, s, tt.str) + } + } +} + +// Benchmarks + +func BenchmarkAcos(b *testing.B) { + for i := 0; i < b.N; i++ { + Acos(.5) + } +} + +func BenchmarkAcosh(b *testing.B) { + for i := 0; i < b.N; i++ { + Acosh(1.5) + } +} + +func BenchmarkAsin(b *testing.B) { + for i := 0; i < b.N; i++ { + Asin(.5) + } +} + +func BenchmarkAsinh(b *testing.B) { + for i := 0; i < b.N; i++ { + 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 BenchmarkAtan2(b *testing.B) { + for i := 0; i < b.N; i++ { + Atan2(.5, 1) + } +} + +func BenchmarkCbrt(b *testing.B) { + for i := 0; i < b.N; i++ { + Cbrt(10) + } +} + +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 BenchmarkExpGo(b *testing.B) { + for i := 0; i < b.N; i++ { + ExpGo(.5) + } +} + +func BenchmarkExpm1(b *testing.B) { + for i := 0; i < b.N; i++ { + Expm1(.5) + } +} + +func BenchmarkExp2(b *testing.B) { + for i := 0; i < b.N; i++ { + Exp2(.5) + } +} + +func BenchmarkExp2Go(b *testing.B) { + for i := 0; i < b.N; i++ { + Exp2Go(.5) + } +} + +func BenchmarkFabs(b *testing.B) { + for i := 0; i < b.N; i++ { + Fabs(.5) + } +} + +func BenchmarkFdim(b *testing.B) { + for i := 0; i < b.N; i++ { + Fdim(10, 3) + } +} + +func BenchmarkFloor(b *testing.B) { + for i := 0; i < b.N; i++ { + Floor(.5) + } +} + +func BenchmarkFmax(b *testing.B) { + for i := 0; i < b.N; i++ { + Fmax(10, 3) + } +} + +func BenchmarkFmin(b *testing.B) { + for i := 0; i < b.N; i++ { + Fmin(10, 3) + } +} + +func BenchmarkFmod(b *testing.B) { + for i := 0; i < b.N; i++ { + Fmod(10, 3) + } +} + +func BenchmarkFrexp(b *testing.B) { + for i := 0; i < b.N; i++ { + Frexp(8) + } +} + +func BenchmarkGamma(b *testing.B) { + for i := 0; i < b.N; i++ { + Gamma(2.5) + } +} + +func BenchmarkHypot(b *testing.B) { + for i := 0; i < b.N; i++ { + Hypot(3, 4) + } +} + +func BenchmarkHypotGo(b *testing.B) { + for i := 0; i < b.N; i++ { + HypotGo(3, 4) + } +} + +func BenchmarkIlogb(b *testing.B) { + for i := 0; i < b.N; i++ { + Ilogb(.5) + } +} + +func BenchmarkJ0(b *testing.B) { + for i := 0; i < b.N; i++ { + J0(2.5) + } +} + +func BenchmarkJ1(b *testing.B) { + for i := 0; i < b.N; i++ { + J1(2.5) + } +} + +func BenchmarkJn(b *testing.B) { + for i := 0; i < b.N; i++ { + Jn(2, 2.5) + } +} + +func BenchmarkLdexp(b *testing.B) { + for i := 0; i < b.N; i++ { + Ldexp(.5, 2) + } +} + +func BenchmarkLgamma(b *testing.B) { + for i := 0; i < b.N; i++ { + Lgamma(2.5) + } +} + +func BenchmarkLog(b *testing.B) { + for i := 0; i < b.N; i++ { + Log(.5) + } +} + +func BenchmarkLogb(b *testing.B) { + for i := 0; i < b.N; i++ { + Logb(.5) + } +} + +func BenchmarkLog1p(b *testing.B) { + for i := 0; i < b.N; i++ { + Log1p(.5) + } +} + +func BenchmarkLog10(b *testing.B) { + for i := 0; i < b.N; i++ { + Log10(.5) + } +} + +func BenchmarkLog2(b *testing.B) { + for i := 0; i < b.N; i++ { + Log2(.5) + } +} + +func BenchmarkModf(b *testing.B) { + for i := 0; i < b.N; i++ { + Modf(1.5) + } +} + +func BenchmarkNextafter(b *testing.B) { + for i := 0; i < b.N; i++ { + Nextafter(.5, 1) + } +} + +func BenchmarkPowInt(b *testing.B) { + for i := 0; i < b.N; i++ { + Pow(2, 2) + } +} + +func BenchmarkPowFrac(b *testing.B) { + for i := 0; i < b.N; i++ { + Pow(2.5, 1.5) + } +} + +func BenchmarkRemainder(b *testing.B) { + for i := 0; i < b.N; i++ { + Remainder(10, 3) + } +} + +func BenchmarkSignbit(b *testing.B) { + for i := 0; i < b.N; i++ { + Signbit(2.5) + } +} + +func BenchmarkSin(b *testing.B) { + for i := 0; i < b.N; i++ { + Sin(.5) + } +} + +func BenchmarkSincos(b *testing.B) { + for i := 0; i < b.N; i++ { + Sincos(.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) + } +} + +func BenchmarkY0(b *testing.B) { + for i := 0; i < b.N; i++ { + Y0(2.5) + } +} + +func BenchmarkY1(b *testing.B) { + for i := 0; i < b.N; i++ { + Y1(2.5) + } +} + +func BenchmarkYn(b *testing.B) { + for i := 0; i < b.N; i++ { + Yn(2, 2.5) + } +} diff --git a/libgo/go/math/asin.go b/libgo/go/math/asin.go new file mode 100644 index 000000000..3bace8ff1 --- /dev/null +++ b/libgo/go/math/asin.go @@ -0,0 +1,50 @@ +// Copyright 2009 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 arcsine and arccosine. + + They are implemented by computing the arctangent + after appropriate range reduction. +*/ + +// Asin returns the arcsine of x. +// +// Special cases are: +// Asin(±0) = ±0 +// Asin(x) = NaN if x < -1 or x > 1 +func Asin(x float64) float64 { + if x == 0 { + return x // special case + } + sign := false + if x < 0 { + x = -x + sign = true + } + if x > 1 { + return NaN() // special case + } + + temp := Sqrt(1 - x*x) + if x > 0.7 { + temp = Pi/2 - satan(temp/x) + } else { + temp = satan(x / temp) + } + + if sign { + temp = -temp + } + return temp +} + +// 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/libgo/go/math/asinh.go b/libgo/go/math/asinh.go new file mode 100644 index 000000000..90dcd27ab --- /dev/null +++ b/libgo/go/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/libgo/go/math/atan.go b/libgo/go/math/atan.go new file mode 100644 index 000000000..9d4ec2f72 --- /dev/null +++ b/libgo/go/math/atan.go @@ -0,0 +1,62 @@ +// Copyright 2009 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 arctangent. + + Atan returns the value of the arctangent of its + argument in the range [-pi/2,pi/2]. + There are no error returns. + Coefficients are #5077 from Hart & Cheney. (19.56D) +*/ + +// xatan evaluates a series valid in the +// range [-0.414...,+0.414...]. (tan(pi/8)) +func xatan(arg float64) float64 { + const ( + P4 = .161536412982230228262e2 + P3 = .26842548195503973794141e3 + P2 = .11530293515404850115428136e4 + P1 = .178040631643319697105464587e4 + P0 = .89678597403663861959987488e3 + Q4 = .5895697050844462222791e2 + Q3 = .536265374031215315104235e3 + Q2 = .16667838148816337184521798e4 + Q1 = .207933497444540981287275926e4 + Q0 = .89678597403663861962481162e3 + ) + sq := arg * arg + value := ((((P4*sq+P3)*sq+P2)*sq+P1)*sq + P0) + value = value / (((((sq+Q4)*sq+Q3)*sq+Q2)*sq+Q1)*sq + Q0) + return value * arg +} + +// satan reduces its argument (known to be positive) +// to the range [0,0.414...] and calls xatan. +func satan(arg float64) float64 { + if arg < Sqrt2-1 { + return xatan(arg) + } + if arg > Sqrt2+1 { + return Pi/2 - xatan(1/arg) + } + return Pi/4 + xatan((arg-1)/(arg+1)) +} + +// Atan returns the arctangent of x. +// +// Special cases are: +// Atan(±0) = ±0 +// Atan(±Inf) = ±Pi/2 +func Atan(x float64) float64 { + if x == 0 { + return x + } + if x > 0 { + return satan(x) + } + return -satan(-x) +} diff --git a/libgo/go/math/atan2.go b/libgo/go/math/atan2.go new file mode 100644 index 000000000..49d4bdd71 --- /dev/null +++ b/libgo/go/math/atan2.go @@ -0,0 +1,71 @@ +// Copyright 2009 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 + +// Atan2 returns the arc tangent of y/x, using +// the signs of the two to determine the quadrant +// of the return value. +// +// Special cases are (in order): +// Atan2(y, NaN) = NaN +// Atan2(NaN, x) = NaN +// Atan2(+0, x>=0) = +0 +// Atan2(-0, x>=0) = -0 +// Atan2(+0, x<=-0) = +Pi +// Atan2(-0, x<=-0) = -Pi +// Atan2(y>0, 0) = +Pi/2 +// Atan2(y<0, 0) = -Pi/2 +// Atan2(+Inf, +Inf) = +Pi/4 +// Atan2(-Inf, +Inf) = -Pi/4 +// Atan2(+Inf, -Inf) = 3Pi/4 +// Atan2(-Inf, -Inf) = -3Pi/4 +// Atan2(y, +Inf) = 0 +// Atan2(y>0, -Inf) = +Pi +// Atan2(y<0, -Inf) = -Pi +// Atan2(+Inf, x) = +Pi/2 +// Atan2(-Inf, x) = -Pi/2 +func Atan2(y, x float64) float64 { + // TODO(rsc): Remove manual inlining of IsNaN, IsInf + // when compiler does it for us + // special cases + switch { + case y != y || x != x: // IsNaN(y) || IsNaN(x): + return NaN() + case y == 0: + if x >= 0 && !Signbit(x) { + return Copysign(0, y) + } + return Copysign(Pi, y) + case x == 0: + return Copysign(Pi/2, y) + case x < -MaxFloat64 || x > MaxFloat64: // IsInf(x, 0): + if x > MaxFloat64 { // IsInf(x, 1) { + switch { + case y < -MaxFloat64 || y > MaxFloat64: // IsInf(y, -1) || IsInf(y, 1): + return Copysign(Pi/4, y) + default: + return Copysign(0, y) + } + } + switch { + case y < -MaxFloat64 || y > MaxFloat64: // IsInf(y, -1) || IsInf(y, 1): + return Copysign(3*Pi/4, y) + default: + return Copysign(Pi, y) + } + case y < -MaxFloat64 || y > MaxFloat64: //IsInf(y, 0): + return Copysign(Pi/2, y) + } + + // Call atan and determine the quadrant. + q := Atan(y / x) + if x < 0 { + if q <= 0 { + return q + Pi + } + return q - Pi + } + return q +} diff --git a/libgo/go/math/atanh.go b/libgo/go/math/atanh.go new file mode 100644 index 000000000..6aecb7b3b --- /dev/null +++ b/libgo/go/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/libgo/go/math/bits.go b/libgo/go/math/bits.go new file mode 100644 index 000000000..a1dca3ed6 --- /dev/null +++ b/libgo/go/math/bits.go @@ -0,0 +1,59 @@ +// Copyright 2009 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 + +const ( + uvnan = 0x7FF0000000000001 + uvinf = 0x7FF0000000000000 + uvneginf = 0xFFF0000000000000 + mask = 0x7FF + shift = 64 - 11 - 1 + bias = 1023 +) + +// Inf returns positive infinity if sign >= 0, negative infinity if sign < 0. +func Inf(sign int) float64 { + var v uint64 + if sign >= 0 { + v = uvinf + } else { + v = uvneginf + } + return Float64frombits(v) +} + +// NaN returns an IEEE 754 ``not-a-number'' value. +func NaN() float64 { return Float64frombits(uvnan) } + +// IsNaN returns whether f is an IEEE 754 ``not-a-number'' value. +func IsNaN(f float64) (is bool) { + // IEEE 754 says that only NaNs satisfy f != f. + // To avoid the floating-point hardware, could use: + // x := Float64bits(f); + // return uint32(x>>shift)&mask == mask && x != uvinf && x != uvneginf + return f != f +} + +// IsInf returns whether f is an infinity, according to sign. +// If sign > 0, IsInf returns whether f is positive infinity. +// If sign < 0, IsInf returns whether f is negative infinity. +// If sign == 0, IsInf returns whether f is either infinity. +func IsInf(f float64, sign int) bool { + // Test for infinity by comparing against maximum float. + // To avoid the floating-point hardware, could use: + // x := Float64bits(f); + // return sign >= 0 && x == uvinf || sign <= 0 && x == uvneginf; + return sign >= 0 && f > MaxFloat64 || sign <= 0 && f < -MaxFloat64 +} + +// normalize returns a normal number y and exponent exp +// satisfying x == y × 2**exp. It assumes x is finite and non-zero. +func normalize(x float64) (y float64, exp int) { + const SmallestNormal = 2.2250738585072014e-308 // 2**-1022 + if Fabs(x) < SmallestNormal { + return x * (1 << 52), -52 + } + return x, 0 +} diff --git a/libgo/go/math/cbrt.go b/libgo/go/math/cbrt.go new file mode 100644 index 000000000..d2b7e910b --- /dev/null +++ b/libgo/go/math/cbrt.go @@ -0,0 +1,79 @@ +// Copyright 2009 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 algorithm is based in part on "Optimal Partitioning of + Newton's Method for Calculating Roots", by Gunter Meinardus + and G. D. Taylor, Mathematics of Computation © 1980 American + Mathematical Society. + (http://www.jstor.org/stable/2006387?seq=9, accessed 11-Feb-2010) +*/ + +// Cbrt returns the cube root of its argument. +// +// Special cases are: +// Cbrt(±0) = ±0 +// Cbrt(±Inf) = ±Inf +// Cbrt(NaN) = NaN +func Cbrt(x float64) float64 { + const ( + A1 = 1.662848358e-01 + A2 = 1.096040958e+00 + A3 = 4.105032829e-01 + A4 = 5.649335816e-01 + B1 = 2.639607233e-01 + B2 = 8.699282849e-01 + B3 = 1.629083358e-01 + B4 = 2.824667908e-01 + C1 = 4.190115298e-01 + C2 = 6.904625373e-01 + C3 = 6.46502159e-02 + C4 = 1.412333954e-01 + ) + // TODO(rsc): Remove manual inlining of IsNaN, IsInf + // when compiler does it for us + // special cases + switch { + case x == 0 || x != x || x < -MaxFloat64 || x > MaxFloat64: // x == 0 || IsNaN(x) || IsInf(x, 0): + return x + } + sign := false + if x < 0 { + x = -x + sign = true + } + // Reduce argument + f, e := Frexp(x) + m := e % 3 + if m > 0 { + m -= 3 + e -= m // e is multiple of 3 + } + f = Ldexp(f, m) // 0.125 <= f < 1.0 + + // Estimate cube root + switch m { + case 0: // 0.5 <= f < 1.0 + f = A1*f + A2 - A3/(A4+f) + case -1: // 0.25 <= f < 0.5 + f = B1*f + B2 - B3/(B4+f) + default: // 0.125 <= f < 0.25 + f = C1*f + C2 - C3/(C4+f) + } + y := Ldexp(f, e/3) // e/3 = exponent of cube root + + // Iterate + s := y * y * y + t := s + x + y *= (t + x) / (s + t) + // Reiterate + s = (y*y*y - x) / x + y -= y * (((14.0/81.0)*s-(2.0/9.0))*s + (1.0 / 3.0)) * s + if sign { + y = -y + } + return y +} diff --git a/libgo/go/math/const.go b/libgo/go/math/const.go new file mode 100644 index 000000000..b53527a4f --- /dev/null +++ b/libgo/go/math/const.go @@ -0,0 +1,53 @@ +// Copyright 2009 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. + +// The math package provides basic constants and mathematical functions. +package math + +// Mathematical constants. +// Reference: http://www.research.att.com/~njas/sequences/Axxxxxx +const ( + E = 2.71828182845904523536028747135266249775724709369995957496696763 // A001113 + Pi = 3.14159265358979323846264338327950288419716939937510582097494459 // A000796 + Phi = 1.61803398874989484820458683436563811772030917980576286213544862 // A001622 + + Sqrt2 = 1.41421356237309504880168872420969807856967187537694807317667974 // A002193 + SqrtE = 1.64872127070012814684865078781416357165377610071014801157507931 // A019774 + SqrtPi = 1.77245385090551602729816748334114518279754945612238712821380779 // A002161 + SqrtPhi = 1.27201964951406896425242246173749149171560804184009624861664038 // A139339 + + Ln2 = 0.693147180559945309417232121458176568075500134360255254120680009 // A002162 + Log2E = 1 / Ln2 + Ln10 = 2.30258509299404568401799145468436420760110148862877297603332790 // A002392 + Log10E = 1 / Ln10 +) + +// Floating-point limit values. +// Max is the largest finite value representable by the type. +// SmallestNonzero is the smallest positive, non-zero value representable by the type. +const ( + MaxFloat32 = 3.40282346638528859811704183484516925440e+38 /* 2**127 * (2**24 - 1) / 2**23 */ + SmallestNonzeroFloat32 = 1.401298464324817070923729583289916131280e-45 /* 1 / 2**(127 - 1 + 23) */ + + MaxFloat64 = 1.797693134862315708145274237317043567981e+308 /* 2**1023 * (2**53 - 1) / 2**52 */ + SmallestNonzeroFloat64 = 4.940656458412465441765687928682213723651e-324 /* 1 / 2**(1023 - 1 + 52) */ +) + +// Integer limit values. +const ( + MaxInt8 = 1<<7 - 1 + MinInt8 = -1 << 7 + MaxInt16 = 1<<15 - 1 + MinInt16 = -1 << 15 + MaxInt32 = 1<<31 - 1 + MinInt32 = -1 << 31 + MaxInt64 = 1<<63 - 1 + MinInt64 = -1 << 63 + MaxUint8 = 1<<8 - 1 + MaxUint16 = 1<<16 - 1 + MaxUint32 = 1<<32 - 1 + MaxUint64 = 1<<64 - 1 +) + +// BUG(rsc): The manual should define the special cases for all of these functions. diff --git a/libgo/go/math/copysign.go b/libgo/go/math/copysign.go new file mode 100644 index 000000000..ee65456a1 --- /dev/null +++ b/libgo/go/math/copysign.go @@ -0,0 +1,12 @@ +// 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 { + const sign = 1 << 63 + return Float64frombits(Float64bits(x)&^sign | Float64bits(y)&sign) +} diff --git a/libgo/go/math/erf.go b/libgo/go/math/erf.go new file mode 100644 index 000000000..b60899933 --- /dev/null +++ b/libgo/go/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: +// Erfc(+Inf) = 0 +// Erfc(-Inf) = 2 +// Erfc(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/libgo/go/math/exp.go b/libgo/go/math/exp.go new file mode 100644 index 000000000..c519c2cb6 --- /dev/null +++ b/libgo/go/math/exp.go @@ -0,0 +1,14 @@ +// Copyright 2009 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 + +// Exp returns e**x, the base-e exponential of x. +// +// Special cases are: +// Exp(+Inf) = +Inf +// Exp(NaN) = NaN +// Very large values overflow to 0 or +Inf. +// Very small values underflow to 1. +func Exp(x float64) float64 { return expGo(x) } diff --git a/libgo/go/math/exp2.go b/libgo/go/math/exp2.go new file mode 100644 index 000000000..1cface9d3 --- /dev/null +++ b/libgo/go/math/exp2.go @@ -0,0 +1,10 @@ +// 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 + +// Exp2 returns 2**x, the base-2 exponential of x. +// +// Special cases are the same as Exp. +func Exp2(x float64) float64 { return exp2Go(x) } diff --git a/libgo/go/math/exp_port.go b/libgo/go/math/exp_port.go new file mode 100644 index 000000000..071420c24 --- /dev/null +++ b/libgo/go/math/exp_port.go @@ -0,0 +1,192 @@ +// Copyright 2009 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_exp.c +// and came with this notice. The go code is a simplified +// version of the original C. +// +// ==================================================== +// Copyright (C) 2004 by Sun Microsystems, Inc. All rights reserved. +// +// Permission to use, copy, modify, and distribute this +// software is freely granted, provided that this notice +// is preserved. +// ==================================================== +// +// +// exp(x) +// Returns the exponential of x. +// +// Method +// 1. Argument reduction: +// Reduce x to an r so that |r| <= 0.5*ln2 ~ 0.34658. +// Given x, find r and integer k such that +// +// x = k*ln2 + r, |r| <= 0.5*ln2. +// +// Here r will be represented as r = hi-lo for better +// accuracy. +// +// 2. Approximation of exp(r) by a special rational function on +// the interval [0,0.34658]: +// Write +// R(r**2) = r*(exp(r)+1)/(exp(r)-1) = 2 + r*r/6 - r**4/360 + ... +// We use a special Remes algorithm on [0,0.34658] to generate +// a polynomial of degree 5 to approximate R. The maximum error +// of this polynomial approximation is bounded by 2**-59. In +// other words, +// R(z) ~ 2.0 + P1*z + P2*z**2 + P3*z**3 + P4*z**4 + P5*z**5 +// (where z=r*r, and the values of P1 to P5 are listed below) +// and +// | 5 | -59 +// | 2.0+P1*z+...+P5*z - R(z) | <= 2 +// | | +// The computation of exp(r) thus becomes +// 2*r +// exp(r) = 1 + ------- +// R - r +// r*R1(r) +// = 1 + r + ----------- (for better accuracy) +// 2 - R1(r) +// where +// 2 4 10 +// R1(r) = r - (P1*r + P2*r + ... + P5*r ). +// +// 3. Scale back to obtain exp(x): +// From step 1, we have +// exp(x) = 2**k * exp(r) +// +// Special cases: +// exp(INF) is INF, exp(NaN) is NaN; +// exp(-INF) is 0, and +// for finite argument, only exp(0)=1 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 exp(x) overflow +// if x < -7.45133219101941108420e+02 then exp(x) underflow +// +// 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. + +// Exp returns e**x, the base-e exponential of x. +// +// Special cases are: +// Exp(+Inf) = +Inf +// Exp(NaN) = NaN +// Very large values overflow to 0 or +Inf. +// Very small values underflow to 1. +func expGo(x float64) float64 { + const ( + Ln2Hi = 6.93147180369123816490e-01 + Ln2Lo = 1.90821492927058770002e-10 + Log2e = 1.44269504088896338700e+00 + + Overflow = 7.09782712893383973096e+02 + Underflow = -7.45133219101941108420e+02 + NearZero = 1.0 / (1 << 28) // 2**-28 + ) + + // TODO(rsc): Remove manual inlining of IsNaN, IsInf + // when compiler does it for us + // special cases + switch { + case x != x || x > MaxFloat64: // IsNaN(x) || IsInf(x, 1): + return x + case x < -MaxFloat64: // IsInf(x, -1): + return 0 + case x > Overflow: + return Inf(1) + case x < Underflow: + return 0 + case -NearZero < x && x < NearZero: + return 1 + x + } + + // reduce; computed as r = hi - lo for extra precision. + var k int + switch { + case x < 0: + k = int(Log2e*x - 0.5) + case x > 0: + k = int(Log2e*x + 0.5) + } + hi := x - float64(k)*Ln2Hi + lo := float64(k) * Ln2Lo + + // compute + return exp(hi, lo, k) +} + +// Exp2 returns 2**x, the base-2 exponential of x. +// +// Special cases are the same as Exp. +func exp2Go(x float64) float64 { + const ( + Ln2Hi = 6.93147180369123816490e-01 + Ln2Lo = 1.90821492927058770002e-10 + + Overflow = 1.0239999999999999e+03 + Underflow = -1.0740e+03 + ) + + // TODO: remove manual inlining of IsNaN and IsInf + // when compiler does it for us + // special cases + switch { + case x != x || x > MaxFloat64: // IsNaN(x) || IsInf(x, 1): + return x + case x < -MaxFloat64: // IsInf(x, -1): + return 0 + case x > Overflow: + return Inf(1) + case x < Underflow: + return 0 + } + + // argument reduction; x = r×lg(e) + k with |r| ≤ ln(2)/2. + // computed as r = hi - lo for extra precision. + var k int + switch { + case x > 0: + k = int(x + 0.5) + case x < 0: + k = int(x - 0.5) + } + t := x - float64(k) + hi := t * Ln2Hi + lo := -t * Ln2Lo + + // compute + return exp(hi, lo, k) +} + +// exp returns e**r × 2**k where r = hi - lo and |r| ≤ ln(2)/2. +func exp(hi, lo float64, k int) float64 { + const ( + P1 = 1.66666666666666019037e-01 /* 0x3FC55555; 0x5555553E */ + P2 = -2.77777777770155933842e-03 /* 0xBF66C16C; 0x16BEBD93 */ + P3 = 6.61375632143793436117e-05 /* 0x3F11566A; 0xAF25DE2C */ + P4 = -1.65339022054652515390e-06 /* 0xBEBBBD41; 0xC5D26BF1 */ + P5 = 4.13813679705723846039e-08 /* 0x3E663769; 0x72BEA4D0 */ + ) + + r := hi - lo + t := r * r + c := r - t*(P1+t*(P2+t*(P3+t*(P4+t*P5)))) + y := 1 - ((lo - (r*c)/(2-c)) - hi) + // TODO(rsc): make sure Ldexp can handle boundary k + return Ldexp(y, k) +} diff --git a/libgo/go/math/exp_test.go b/libgo/go/math/exp_test.go new file mode 100644 index 000000000..7381fd5ad --- /dev/null +++ b/libgo/go/math/exp_test.go @@ -0,0 +1,10 @@ +// 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 + +// Make expGo and exp2Go available for testing. + +func ExpGo(x float64) float64 { return expGo(x) } +func Exp2Go(x float64) float64 { return exp2Go(x) } diff --git a/libgo/go/math/expm1.go b/libgo/go/math/expm1.go new file mode 100644 index 000000000..35100caa4 --- /dev/null +++ b/libgo/go/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/libgo/go/math/fabs.go b/libgo/go/math/fabs.go new file mode 100644 index 000000000..343123126 --- /dev/null +++ b/libgo/go/math/fabs.go @@ -0,0 +1,21 @@ +// Copyright 2009 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 + +// Fabs returns the absolute value of x. +// +// Special cases are: +// Fabs(+Inf) = +Inf +// Fabs(-Inf) = +Inf +// Fabs(NaN) = NaN +func Fabs(x float64) float64 { + switch { + case x < 0: + return -x + case x == 0: + return 0 // return correctly fabs(-0) + } + return x +} diff --git a/libgo/go/math/fdim.go b/libgo/go/math/fdim.go new file mode 100644 index 000000000..18993137a --- /dev/null +++ b/libgo/go/math/fdim.go @@ -0,0 +1,29 @@ +// 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 + +// Fdim returns the maximum of x-y or 0. +func Fdim(x, y float64) float64 { + if x > y { + return x - y + } + return 0 +} + +// Fmax returns the larger of x or y. +func Fmax(x, y float64) float64 { + if x > y { + return x + } + return y +} + +// Fmin returns the smaller of x or y. +func Fmin(x, y float64) float64 { + if x < y { + return x + } + return y +} diff --git a/libgo/go/math/floor.go b/libgo/go/math/floor.go new file mode 100644 index 000000000..b22b94ad6 --- /dev/null +++ b/libgo/go/math/floor.go @@ -0,0 +1,53 @@ +// 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. + +package math + + +// Floor returns the greatest integer value less than or equal to x. +// +// Special cases are: +// Floor(+Inf) = +Inf +// 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 == 0 || x != x || x > MaxFloat64 || x < -MaxFloat64 { // x == 0 || IsNaN(x) || IsInf(x, 0) + return x + } + if x < 0 { + d, fract := Modf(-x) + if fract != 0.0 { + d = d + 1 + } + return -d + } + d, _ := Modf(x) + return d +} + +// Ceil returns the least integer value greater than or equal to x. +// +// Special cases are: +// Ceil(+Inf) = +Inf +// 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 == 0 || x != x || x > MaxFloat64 || x < -MaxFloat64 { // x == 0 || IsNaN(x) || IsInf(x, 0) + return x + } + d, _ := Modf(x) + return d +} diff --git a/libgo/go/math/fmod.go b/libgo/go/math/fmod.go new file mode 100644 index 000000000..fc57f7483 --- /dev/null +++ b/libgo/go/math/fmod.go @@ -0,0 +1,48 @@ +// 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. + +package math + + +/* + Floating-point mod function. +*/ + +// Fmod returns the floating-point remainder of x/y. +// The magnitude of the result is less than y and its +// sign agrees with that of x. +// +// Special cases are: +// if x is not finite, Fmod returns NaN +// if y is 0 or NaN, Fmod returns NaN +func Fmod(x, y float64) float64 { + // TODO(rsc): Remove manual inlining of IsNaN, IsInf + // when compiler does it for us. + if y == 0 || x > MaxFloat64 || x < -MaxFloat64 || x != x || y != y { // y == 0 || IsInf(x, 0) || IsNaN(x) || IsNan(y) + return NaN() + } + if y < 0 { + y = -y + } + + yfr, yexp := Frexp(y) + sign := false + r := x + if x < 0 { + r = -x + sign = true + } + + for r >= y { + rfr, rexp := Frexp(r) + if rfr < yfr { + rexp = rexp - 1 + } + r = r - Ldexp(y, rexp-yexp) + } + if sign { + r = -r + } + return r +} diff --git a/libgo/go/math/frexp.go b/libgo/go/math/frexp.go new file mode 100644 index 000000000..867b78f36 --- /dev/null +++ b/libgo/go/math/frexp.go @@ -0,0 +1,33 @@ +// Copyright 2009 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 + +// Frexp breaks f into a normalized fraction +// and an integral power of two. +// It returns frac and exp satisfying f == frac × 2**exp, +// with the absolute value of frac in the interval [½, 1). +// +// Special cases are: +// Frexp(±0) = ±0, 0 +// Frexp(±Inf) = ±Inf, 0 +// Frexp(NaN) = NaN, 0 +func Frexp(f float64) (frac float64, exp int) { + // TODO(rsc): Remove manual inlining of IsNaN, IsInf + // when compiler does it for us + // special cases + switch { + case f == 0: + return f, 0 // correctly return -0 + case f < -MaxFloat64 || f > MaxFloat64 || f != f: // IsInf(f, 0) || IsNaN(f): + return f, 0 + } + f, exp = normalize(f) + x := Float64bits(f) + exp += int((x>>shift)&mask) - bias + 1 + x &^= mask << shift + x |= (-1 + bias) << shift + frac = Float64frombits(x) + return +} diff --git a/libgo/go/math/gamma.go b/libgo/go/math/gamma.go new file mode 100644 index 000000000..73ca0e53a --- /dev/null +++ b/libgo/go/math/gamma.go @@ -0,0 +1,188 @@ +// 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 http://netlib.sandia.gov/cephes/cprob/gamma.c. +// The go code is a simplified version of the original C. +// +// tgamma.c +// +// Gamma function +// +// SYNOPSIS: +// +// double x, y, tgamma(); +// extern int signgam; +// +// y = tgamma( x ); +// +// DESCRIPTION: +// +// Returns gamma function of the argument. The result is +// correctly signed, and the sign (+1 or -1) is also +// returned in a global (extern) variable named signgam. +// This variable is also filled in by the logarithmic gamma +// function lgamma(). +// +// Arguments |x| <= 34 are reduced by recurrence and the function +// approximated by a rational function of degree 6/7 in the +// interval (2,3). Large arguments are handled by Stirling's +// formula. Large negative arguments are made positive using +// a reflection formula. +// +// ACCURACY: +// +// Relative error: +// arithmetic domain # trials peak rms +// DEC -34, 34 10000 1.3e-16 2.5e-17 +// IEEE -170,-33 20000 2.3e-15 3.3e-16 +// IEEE -33, 33 20000 9.4e-16 2.2e-16 +// IEEE 33, 171.6 20000 2.3e-15 3.2e-16 +// +// Error for arguments outside the test range will be larger +// owing to error amplification by the exponential function. +// +// Cephes Math Library Release 2.8: June, 2000 +// Copyright 1984, 1987, 1989, 1992, 2000 by Stephen L. Moshier +// +// The readme file at http://netlib.sandia.gov/cephes/ says: +// Some software in this archive may be from the book _Methods and +// Programs for Mathematical Functions_ (Prentice-Hall or Simon & Schuster +// International, 1989) or from the Cephes Mathematical Library, a +// commercial product. In either event, it is copyrighted by the author. +// What you see here may be used freely but it comes with no support or +// guarantee. +// +// The two known misprints in the book are repaired here in the +// source listings for the gamma function and the incomplete beta +// integral. +// +// Stephen L. Moshier +// moshier@na-net.ornl.gov + +var _P = []float64{ + 1.60119522476751861407e-04, + 1.19135147006586384913e-03, + 1.04213797561761569935e-02, + 4.76367800457137231464e-02, + 2.07448227648435975150e-01, + 4.94214826801497100753e-01, + 9.99999999999999996796e-01, +} +var _Q = []float64{ + -2.31581873324120129819e-05, + 5.39605580493303397842e-04, + -4.45641913851797240494e-03, + 1.18139785222060435552e-02, + 3.58236398605498653373e-02, + -2.34591795718243348568e-01, + 7.14304917030273074085e-02, + 1.00000000000000000320e+00, +} +var _S = []float64{ + 7.87311395793093628397e-04, + -2.29549961613378126380e-04, + -2.68132617805781232825e-03, + 3.47222221605458667310e-03, + 8.33333333333482257126e-02, +} + +// Gamma function computed by Stirling's formula. +// The polynomial is valid for 33 <= x <= 172. +func stirling(x float64) float64 { + const ( + SqrtTwoPi = 2.506628274631000502417 + MaxStirling = 143.01608 + ) + w := 1 / x + w = 1 + w*((((_S[0]*w+_S[1])*w+_S[2])*w+_S[3])*w+_S[4]) + y := Exp(x) + if x > MaxStirling { // avoid Pow() overflow + v := Pow(x, 0.5*x-0.25) + y = v * (v / y) + } else { + y = Pow(x, x-0.5) / y + } + y = SqrtTwoPi * y * w + return y +} + +// Gamma(x) returns the Gamma function of x. +// +// Special cases are: +// Gamma(Inf) = Inf +// Gamma(-Inf) = -Inf +// Gamma(NaN) = NaN +// Large values overflow to +Inf. +// Negative integer values equal ±Inf. +func Gamma(x float64) float64 { + const Euler = 0.57721566490153286060651209008240243104215933593992 // A001620 + // special cases + switch { + case x < -MaxFloat64 || x != x: // IsInf(x, -1) || IsNaN(x): + return x + case x < -170.5674972726612 || x > 171.61447887182298: + return Inf(1) + } + q := Fabs(x) + p := Floor(q) + if q > 33 { + if x >= 0 { + return stirling(x) + } + signgam := 1 + if ip := int(p); ip&1 == 0 { + signgam = -1 + } + z := q - p + if z > 0.5 { + p = p + 1 + z = q - p + } + z = q * Sin(Pi*z) + if z == 0 { + return Inf(signgam) + } + z = Pi / (Fabs(z) * stirling(q)) + return float64(signgam) * z + } + + // Reduce argument + z := 1.0 + for x >= 3 { + x = x - 1 + z = z * x + } + for x < 0 { + if x > -1e-09 { + goto small + } + z = z / x + x = x + 1 + } + for x < 2 { + if x < 1e-09 { + goto small + } + z = z / x + x = x + 1 + } + + if x == 2 { + return z + } + + x = x - 2 + p = (((((x*_P[0]+_P[1])*x+_P[2])*x+_P[3])*x+_P[4])*x+_P[5])*x + _P[6] + q = ((((((x*_Q[0]+_Q[1])*x+_Q[2])*x+_Q[3])*x+_Q[4])*x+_Q[5])*x+_Q[6])*x + _Q[7] + return z * p / q + +small: + if x == 0 { + return Inf(1) + } + return z / ((1 + Euler*x) * x) +} diff --git a/libgo/go/math/hypot.go b/libgo/go/math/hypot.go new file mode 100644 index 000000000..ecd115d9e --- /dev/null +++ b/libgo/go/math/hypot.go @@ -0,0 +1,41 @@ +// 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 + +/* + Hypot -- sqrt(p*p + q*q), but overflows only if the result does. +*/ + +// Hypot computes Sqrt(p*p + q*q), taking care to avoid +// unnecessary overflow and underflow. +// +// Special cases are: +// Hypot(p, q) = +Inf if p or q is infinite +// Hypot(p, q) = NaN if p or q is NaN +func Hypot(p, q float64) float64 { + // TODO(rsc): Remove manual inlining of IsNaN, IsInf + // when compiler does it for us + // special cases + switch { + case p < -MaxFloat64 || p > MaxFloat64 || q < -MaxFloat64 || q > MaxFloat64: // IsInf(p, 0) || IsInf(q, 0): + return Inf(1) + case p != p || q != q: // IsNaN(p) || IsNaN(q): + return NaN() + } + if p < 0 { + p = -p + } + if q < 0 { + q = -q + } + if p < q { + p, q = q, p + } + if p == 0 { + return 0 + } + q = q / p + return p * Sqrt(1+q*q) +} diff --git a/libgo/go/math/hypot_port.go b/libgo/go/math/hypot_port.go new file mode 100644 index 000000000..27f335ba2 --- /dev/null +++ b/libgo/go/math/hypot_port.go @@ -0,0 +1,63 @@ +// 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. + +package math + +/* + Hypot -- sqrt(p*p + q*q), but overflows only if the result does. + See: + Cleve Moler and Donald Morrison, + Replacing Square Roots by Pythagorean Sums + IBM Journal of Research and Development, + Vol. 27, Number 6, pp. 577-581, Nov. 1983 +*/ + +// Hypot computes Sqrt(p*p + q*q), taking care to avoid +// unnecessary overflow and underflow. +// +// Special cases are: +// Hypot(p, q) = +Inf if p or q is infinite +// Hypot(p, q) = NaN if p or q is NaN +func hypotGo(p, q float64) float64 { + // TODO(rsc): Remove manual inlining of IsNaN, IsInf + // when compiler does it for us + // special cases + switch { + case p < -MaxFloat64 || p > MaxFloat64 || q < -MaxFloat64 || q > MaxFloat64: // IsInf(p, 0) || IsInf(q, 0): + return Inf(1) + case p != p || q != q: // IsNaN(p) || IsNaN(q): + return NaN() + } + if p < 0 { + p = -p + } + if q < 0 { + q = -q + } + + if p < q { + p, q = q, p + } + + if p == 0 { + return 0 + } + + pfac := p + q = q / p + r := q + p = 1 + for { + r = r * r + s := r + 4 + if s == 4 { + return p * pfac + } + r = r / s + p = p + 2*r*p + q = q * r + r = q / p + } + panic("unreachable") +} diff --git a/libgo/go/math/hypot_test.go b/libgo/go/math/hypot_test.go new file mode 100644 index 000000000..85ce1d404 --- /dev/null +++ b/libgo/go/math/hypot_test.go @@ -0,0 +1,9 @@ +// 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 + +// Make hypotGo available for testing. + +func HypotGo(x, y float64) float64 { return hypotGo(x, y) } diff --git a/libgo/go/math/j0.go b/libgo/go/math/j0.go new file mode 100644 index 000000000..5aaf4ab9c --- /dev/null +++ b/libgo/go/math/j0.go @@ -0,0 +1,433 @@ +// 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 + +/* + Bessel function of the first and second kinds of order zero. +*/ + +// The original C code and the long comment below are +// from FreeBSD's /usr/src/lib/msun/src/e_j0.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_j0(x), __ieee754_y0(x) +// Bessel function of the first and second kinds of order zero. +// Method -- j0(x): +// 1. For tiny x, we use j0(x) = 1 - x**2/4 + x**4/64 - ... +// 2. Reduce x to |x| since j0(x)=j0(-x), and +// for x in (0,2) +// j0(x) = 1-z/4+ z**2*R0/S0, where z = x*x; +// (precision: |j0-1+z/4-z**2R0/S0 |<2**-63.67 ) +// for x in (2,inf) +// j0(x) = sqrt(2/(pi*x))*(p0(x)*cos(x0)-q0(x)*sin(x0)) +// where x0 = x-pi/4. It is better to compute sin(x0),cos(x0) +// as follow: +// cos(x0) = cos(x)cos(pi/4)+sin(x)sin(pi/4) +// = 1/sqrt(2) * (cos(x) + sin(x)) +// sin(x0) = sin(x)cos(pi/4)-cos(x)sin(pi/4) +// = 1/sqrt(2) * (sin(x) - cos(x)) +// (To avoid cancellation, use +// sin(x) +- cos(x) = -cos(2x)/(sin(x) -+ cos(x)) +// to compute the worse one.) +// +// 3 Special cases +// j0(nan)= nan +// j0(0) = 1 +// j0(inf) = 0 +// +// Method -- y0(x): +// 1. For x<2. +// Since +// y0(x) = 2/pi*(j0(x)*(ln(x/2)+Euler) + x**2/4 - ...) +// therefore y0(x)-2/pi*j0(x)*ln(x) is an even function. +// We use the following function to approximate y0, +// y0(x) = U(z)/V(z) + (2/pi)*(j0(x)*ln(x)), z= x**2 +// where +// U(z) = u00 + u01*z + ... + u06*z**6 +// V(z) = 1 + v01*z + ... + v04*z**4 +// with absolute approximation error bounded by 2**-72. +// Note: For tiny x, U/V = u0 and j0(x)~1, hence +// y0(tiny) = u0 + (2/pi)*ln(tiny), (choose tiny<2**-27) +// 2. For x>=2. +// y0(x) = sqrt(2/(pi*x))*(p0(x)*cos(x0)+q0(x)*sin(x0)) +// where x0 = x-pi/4. It is better to compute sin(x0),cos(x0) +// by the method mentioned above. +// 3. Special cases: y0(0)=-inf, y0(x<0)=NaN, y0(inf)=0. +// + +// J0 returns the order-zero Bessel function of the first kind. +// +// Special cases are: +// J0(±Inf) = 0 +// J0(0) = 1 +// J0(NaN) = NaN +func J0(x float64) float64 { + const ( + Huge = 1e300 + TwoM27 = 1.0 / (1 << 27) // 2**-27 0x3e40000000000000 + TwoM13 = 1.0 / (1 << 13) // 2**-13 0x3f20000000000000 + Two129 = 1 << 129 // 2**129 0x4800000000000000 + // R0/S0 on [0, 2] + R02 = 1.56249999999999947958e-02 // 0x3F8FFFFFFFFFFFFD + R03 = -1.89979294238854721751e-04 // 0xBF28E6A5B61AC6E9 + R04 = 1.82954049532700665670e-06 // 0x3EBEB1D10C503919 + R05 = -4.61832688532103189199e-09 // 0xBE33D5E773D63FCE + S01 = 1.56191029464890010492e-02 // 0x3F8FFCE882C8C2A4 + S02 = 1.16926784663337450260e-04 // 0x3F1EA6D2DD57DBF4 + S03 = 5.13546550207318111446e-07 // 0x3EA13B54CE84D5A9 + S04 = 1.16614003333790000205e-09 // 0x3E1408BCF4745D8F + ) + // TODO(rsc): Remove manual inlining of IsNaN, IsInf + // when compiler does it for us + // special cases + switch { + case x != x: // IsNaN(x) + return x + case x < -MaxFloat64 || x > MaxFloat64: // IsInf(x, 0): + return 0 + case x == 0: + return 1 + } + + if x < 0 { + x = -x + } + if x >= 2 { + s, c := Sincos(x) + ss := s - c + cc := s + c + + // make sure x+x does not overflow + if x < MaxFloat64/2 { + z := -Cos(x + x) + if s*c < 0 { + cc = z / ss + } else { + ss = z / cc + } + } + + // j0(x) = 1/sqrt(pi) * (P(0,x)*cc - Q(0,x)*ss) / sqrt(x) + // y0(x) = 1/sqrt(pi) * (P(0,x)*ss + Q(0,x)*cc) / sqrt(x) + + var z float64 + if x > Two129 { // |x| > ~6.8056e+38 + z = (1 / SqrtPi) * cc / Sqrt(x) + } else { + u := pzero(x) + v := qzero(x) + z = (1 / SqrtPi) * (u*cc - v*ss) / Sqrt(x) + } + return z // |x| >= 2.0 + } + if x < TwoM13 { // |x| < ~1.2207e-4 + if x < TwoM27 { + return 1 // |x| < ~7.4506e-9 + } + return 1 - 0.25*x*x // ~7.4506e-9 < |x| < ~1.2207e-4 + } + z := x * x + r := z * (R02 + z*(R03+z*(R04+z*R05))) + s := 1 + z*(S01+z*(S02+z*(S03+z*S04))) + if x < 1 { + return 1 + z*(-0.25+(r/s)) // |x| < 1.00 + } + u := 0.5 * x + return (1+u)*(1-u) + z*(r/s) // 1.0 < |x| < 2.0 +} + +// Y0 returns the order-zero Bessel function of the second kind. +// +// Special cases are: +// Y0(+Inf) = 0 +// Y0(0) = -Inf +// Y0(x < 0) = NaN +// Y0(NaN) = NaN +func Y0(x float64) float64 { + const ( + TwoM27 = 1.0 / (1 << 27) // 2**-27 0x3e40000000000000 + Two129 = 1 << 129 // 2**129 0x4800000000000000 + U00 = -7.38042951086872317523e-02 // 0xBFB2E4D699CBD01F + U01 = 1.76666452509181115538e-01 // 0x3FC69D019DE9E3FC + U02 = -1.38185671945596898896e-02 // 0xBF8C4CE8B16CFA97 + U03 = 3.47453432093683650238e-04 // 0x3F36C54D20B29B6B + U04 = -3.81407053724364161125e-06 // 0xBECFFEA773D25CAD + U05 = 1.95590137035022920206e-08 // 0x3E5500573B4EABD4 + U06 = -3.98205194132103398453e-11 // 0xBDC5E43D693FB3C8 + V01 = 1.27304834834123699328e-02 // 0x3F8A127091C9C71A + V02 = 7.60068627350353253702e-05 // 0x3F13ECBBF578C6C1 + V03 = 2.59150851840457805467e-07 // 0x3E91642D7FF202FD + V04 = 4.41110311332675467403e-10 // 0x3DFE50183BD6D9EF + ) + // TODO(rsc): Remove manual inlining of IsNaN, IsInf + // when compiler does it for us + // special cases + switch { + case x < 0 || x != x: // x < 0 || IsNaN(x): + return NaN() + case x > MaxFloat64: // IsInf(x, 1): + return 0 + case x == 0: + return Inf(-1) + } + + if x >= 2 { // |x| >= 2.0 + + // y0(x) = sqrt(2/(pi*x))*(p0(x)*sin(x0)+q0(x)*cos(x0)) + // where x0 = x-pi/4 + // Better formula: + // cos(x0) = cos(x)cos(pi/4)+sin(x)sin(pi/4) + // = 1/sqrt(2) * (sin(x) + cos(x)) + // sin(x0) = sin(x)cos(3pi/4)-cos(x)sin(3pi/4) + // = 1/sqrt(2) * (sin(x) - cos(x)) + // To avoid cancellation, use + // sin(x) +- cos(x) = -cos(2x)/(sin(x) -+ cos(x)) + // to compute the worse one. + + s, c := Sincos(x) + ss := s - c + cc := s + c + + // j0(x) = 1/sqrt(pi) * (P(0,x)*cc - Q(0,x)*ss) / sqrt(x) + // y0(x) = 1/sqrt(pi) * (P(0,x)*ss + Q(0,x)*cc) / sqrt(x) + + // make sure x+x does not overflow + if x < MaxFloat64/2 { + z := -Cos(x + x) + if s*c < 0 { + cc = z / ss + } else { + ss = z / cc + } + } + var z float64 + if x > Two129 { // |x| > ~6.8056e+38 + z = (1 / SqrtPi) * ss / Sqrt(x) + } else { + u := pzero(x) + v := qzero(x) + z = (1 / SqrtPi) * (u*ss + v*cc) / Sqrt(x) + } + return z // |x| >= 2.0 + } + if x <= TwoM27 { + return U00 + (2/Pi)*Log(x) // |x| < ~7.4506e-9 + } + z := x * x + u := U00 + z*(U01+z*(U02+z*(U03+z*(U04+z*(U05+z*U06))))) + v := 1 + z*(V01+z*(V02+z*(V03+z*V04))) + return u/v + (2/Pi)*J0(x)*Log(x) // ~7.4506e-9 < |x| < 2.0 +} + +// The asymptotic expansions of pzero is +// 1 - 9/128 s**2 + 11025/98304 s**4 - ..., where s = 1/x. +// For x >= 2, We approximate pzero by +// pzero(x) = 1 + (R/S) +// where R = pR0 + pR1*s**2 + pR2*s**4 + ... + pR5*s**10 +// S = 1 + pS0*s**2 + ... + pS4*s**10 +// and +// | pzero(x)-1-R/S | <= 2 ** ( -60.26) + +// for x in [inf, 8]=1/[0,0.125] +var p0R8 = [6]float64{ + 0.00000000000000000000e+00, // 0x0000000000000000 + -7.03124999999900357484e-02, // 0xBFB1FFFFFFFFFD32 + -8.08167041275349795626e+00, // 0xC02029D0B44FA779 + -2.57063105679704847262e+02, // 0xC07011027B19E863 + -2.48521641009428822144e+03, // 0xC0A36A6ECD4DCAFC + -5.25304380490729545272e+03, // 0xC0B4850B36CC643D +} +var p0S8 = [5]float64{ + 1.16534364619668181717e+02, // 0x405D223307A96751 + 3.83374475364121826715e+03, // 0x40ADF37D50596938 + 4.05978572648472545552e+04, // 0x40E3D2BB6EB6B05F + 1.16752972564375915681e+05, // 0x40FC810F8F9FA9BD + 4.76277284146730962675e+04, // 0x40E741774F2C49DC +} + +// for x in [8,4.5454]=1/[0.125,0.22001] +var p0R5 = [6]float64{ + -1.14125464691894502584e-11, // 0xBDA918B147E495CC + -7.03124940873599280078e-02, // 0xBFB1FFFFE69AFBC6 + -4.15961064470587782438e+00, // 0xC010A370F90C6BBF + -6.76747652265167261021e+01, // 0xC050EB2F5A7D1783 + -3.31231299649172967747e+02, // 0xC074B3B36742CC63 + -3.46433388365604912451e+02, // 0xC075A6EF28A38BD7 +} +var p0S5 = [5]float64{ + 6.07539382692300335975e+01, // 0x404E60810C98C5DE + 1.05125230595704579173e+03, // 0x40906D025C7E2864 + 5.97897094333855784498e+03, // 0x40B75AF88FBE1D60 + 9.62544514357774460223e+03, // 0x40C2CCB8FA76FA38 + 2.40605815922939109441e+03, // 0x40A2CC1DC70BE864 +} + +// for x in [4.547,2.8571]=1/[0.2199,0.35001] +var p0R3 = [6]float64{ + -2.54704601771951915620e-09, // 0xBE25E1036FE1AA86 + -7.03119616381481654654e-02, // 0xBFB1FFF6F7C0E24B + -2.40903221549529611423e+00, // 0xC00345B2AEA48074 + -2.19659774734883086467e+01, // 0xC035F74A4CB94E14 + -5.80791704701737572236e+01, // 0xC04D0A22420A1A45 + -3.14479470594888503854e+01, // 0xC03F72ACA892D80F +} +var p0S3 = [5]float64{ + 3.58560338055209726349e+01, // 0x4041ED9284077DD3 + 3.61513983050303863820e+02, // 0x40769839464A7C0E + 1.19360783792111533330e+03, // 0x4092A66E6D1061D6 + 1.12799679856907414432e+03, // 0x40919FFCB8C39B7E + 1.73580930813335754692e+02, // 0x4065B296FC379081 +} + +// for x in [2.8570,2]=1/[0.3499,0.5] +var p0R2 = [6]float64{ + -8.87534333032526411254e-08, // 0xBE77D316E927026D + -7.03030995483624743247e-02, // 0xBFB1FF62495E1E42 + -1.45073846780952986357e+00, // 0xBFF736398A24A843 + -7.63569613823527770791e+00, // 0xC01E8AF3EDAFA7F3 + -1.11931668860356747786e+01, // 0xC02662E6C5246303 + -3.23364579351335335033e+00, // 0xC009DE81AF8FE70F +} +var p0S2 = [5]float64{ + 2.22202997532088808441e+01, // 0x40363865908B5959 + 1.36206794218215208048e+02, // 0x4061069E0EE8878F + 2.70470278658083486789e+02, // 0x4070E78642EA079B + 1.53875394208320329881e+02, // 0x40633C033AB6FAFF + 1.46576176948256193810e+01, // 0x402D50B344391809 +} + +func pzero(x float64) float64 { + var p [6]float64 + var q [5]float64 + if x >= 8 { + p = p0R8 + q = p0S8 + } else if x >= 4.5454 { + p = p0R5 + q = p0S5 + } else if x >= 2.8571 { + p = p0R3 + q = p0S3 + } else if x >= 2 { + p = p0R2 + q = p0S2 + } + z := 1 / (x * x) + r := p[0] + z*(p[1]+z*(p[2]+z*(p[3]+z*(p[4]+z*p[5])))) + s := 1 + z*(q[0]+z*(q[1]+z*(q[2]+z*(q[3]+z*q[4])))) + return 1 + r/s +} + +// For x >= 8, the asymptotic expansions of qzero is +// -1/8 s + 75/1024 s**3 - ..., where s = 1/x. +// We approximate pzero by +// qzero(x) = s*(-1.25 + (R/S)) +// where R = qR0 + qR1*s**2 + qR2*s**4 + ... + qR5*s**10 +// S = 1 + qS0*s**2 + ... + qS5*s**12 +// and +// | qzero(x)/s +1.25-R/S | <= 2**(-61.22) + +// for x in [inf, 8]=1/[0,0.125] +var q0R8 = [6]float64{ + 0.00000000000000000000e+00, // 0x0000000000000000 + 7.32421874999935051953e-02, // 0x3FB2BFFFFFFFFE2C + 1.17682064682252693899e+01, // 0x402789525BB334D6 + 5.57673380256401856059e+02, // 0x40816D6315301825 + 8.85919720756468632317e+03, // 0x40C14D993E18F46D + 3.70146267776887834771e+04, // 0x40E212D40E901566 +} +var q0S8 = [6]float64{ + 1.63776026895689824414e+02, // 0x406478D5365B39BC + 8.09834494656449805916e+03, // 0x40BFA2584E6B0563 + 1.42538291419120476348e+05, // 0x4101665254D38C3F + 8.03309257119514397345e+05, // 0x412883DA83A52B43 + 8.40501579819060512818e+05, // 0x4129A66B28DE0B3D + -3.43899293537866615225e+05, // 0xC114FD6D2C9530C5 +} + +// for x in [8,4.5454]=1/[0.125,0.22001] +var q0R5 = [6]float64{ + 1.84085963594515531381e-11, // 0x3DB43D8F29CC8CD9 + 7.32421766612684765896e-02, // 0x3FB2BFFFD172B04C + 5.83563508962056953777e+00, // 0x401757B0B9953DD3 + 1.35111577286449829671e+02, // 0x4060E3920A8788E9 + 1.02724376596164097464e+03, // 0x40900CF99DC8C481 + 1.98997785864605384631e+03, // 0x409F17E953C6E3A6 +} +var q0S5 = [6]float64{ + 8.27766102236537761883e+01, // 0x4054B1B3FB5E1543 + 2.07781416421392987104e+03, // 0x40A03BA0DA21C0CE + 1.88472887785718085070e+04, // 0x40D267D27B591E6D + 5.67511122894947329769e+04, // 0x40EBB5E397E02372 + 3.59767538425114471465e+04, // 0x40E191181F7A54A0 + -5.35434275601944773371e+03, // 0xC0B4EA57BEDBC609 +} + +// for x in [4.547,2.8571]=1/[0.2199,0.35001] +var q0R3 = [6]float64{ + 4.37741014089738620906e-09, // 0x3E32CD036ADECB82 + 7.32411180042911447163e-02, // 0x3FB2BFEE0E8D0842 + 3.34423137516170720929e+00, // 0x400AC0FC61149CF5 + 4.26218440745412650017e+01, // 0x40454F98962DAEDD + 1.70808091340565596283e+02, // 0x406559DBE25EFD1F + 1.66733948696651168575e+02, // 0x4064D77C81FA21E0 +} +var q0S3 = [6]float64{ + 4.87588729724587182091e+01, // 0x40486122BFE343A6 + 7.09689221056606015736e+02, // 0x40862D8386544EB3 + 3.70414822620111362994e+03, // 0x40ACF04BE44DFC63 + 6.46042516752568917582e+03, // 0x40B93C6CD7C76A28 + 2.51633368920368957333e+03, // 0x40A3A8AAD94FB1C0 + -1.49247451836156386662e+02, // 0xC062A7EB201CF40F +} + +// for x in [2.8570,2]=1/[0.3499,0.5] +var q0R2 = [6]float64{ + 1.50444444886983272379e-07, // 0x3E84313B54F76BDB + 7.32234265963079278272e-02, // 0x3FB2BEC53E883E34 + 1.99819174093815998816e+00, // 0x3FFFF897E727779C + 1.44956029347885735348e+01, // 0x402CFDBFAAF96FE5 + 3.16662317504781540833e+01, // 0x403FAA8E29FBDC4A + 1.62527075710929267416e+01, // 0x403040B171814BB4 +} +var q0S2 = [6]float64{ + 3.03655848355219184498e+01, // 0x403E5D96F7C07AED + 2.69348118608049844624e+02, // 0x4070D591E4D14B40 + 8.44783757595320139444e+02, // 0x408A664522B3BF22 + 8.82935845112488550512e+02, // 0x408B977C9C5CC214 + 2.12666388511798828631e+02, // 0x406A95530E001365 + -5.31095493882666946917e+00, // 0xC0153E6AF8B32931 +} + +func qzero(x float64) float64 { + var p, q [6]float64 + if x >= 8 { + p = q0R8 + q = q0S8 + } else if x >= 4.5454 { + p = q0R5 + q = q0S5 + } else if x >= 2.8571 { + p = q0R3 + q = q0S3 + } else if x >= 2 { + p = q0R2 + q = q0S2 + } + z := 1 / (x * x) + r := p[0] + z*(p[1]+z*(p[2]+z*(p[3]+z*(p[4]+z*p[5])))) + s := 1 + z*(q[0]+z*(q[1]+z*(q[2]+z*(q[3]+z*(q[4]+z*q[5]))))) + return (-0.125 + r/s) / x +} diff --git a/libgo/go/math/j1.go b/libgo/go/math/j1.go new file mode 100644 index 000000000..278162e9d --- /dev/null +++ b/libgo/go/math/j1.go @@ -0,0 +1,426 @@ +// 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 + +/* + Bessel function of the first and second kinds of order one. +*/ + +// The original C code and the long comment below are +// from FreeBSD's /usr/src/lib/msun/src/e_j1.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_j1(x), __ieee754_y1(x) +// Bessel function of the first and second kinds of order one. +// Method -- j1(x): +// 1. For tiny x, we use j1(x) = x/2 - x**3/16 + x**5/384 - ... +// 2. Reduce x to |x| since j1(x)=-j1(-x), and +// for x in (0,2) +// j1(x) = x/2 + x*z*R0/S0, where z = x*x; +// (precision: |j1/x - 1/2 - R0/S0 |<2**-61.51 ) +// for x in (2,inf) +// j1(x) = sqrt(2/(pi*x))*(p1(x)*cos(x1)-q1(x)*sin(x1)) +// y1(x) = sqrt(2/(pi*x))*(p1(x)*sin(x1)+q1(x)*cos(x1)) +// where x1 = x-3*pi/4. It is better to compute sin(x1),cos(x1) +// as follow: +// cos(x1) = cos(x)cos(3pi/4)+sin(x)sin(3pi/4) +// = 1/sqrt(2) * (sin(x) - cos(x)) +// sin(x1) = sin(x)cos(3pi/4)-cos(x)sin(3pi/4) +// = -1/sqrt(2) * (sin(x) + cos(x)) +// (To avoid cancellation, use +// sin(x) +- cos(x) = -cos(2x)/(sin(x) -+ cos(x)) +// to compute the worse one.) +// +// 3 Special cases +// j1(nan)= nan +// j1(0) = 0 +// j1(inf) = 0 +// +// Method -- y1(x): +// 1. screen out x<=0 cases: y1(0)=-inf, y1(x<0)=NaN +// 2. For x<2. +// Since +// y1(x) = 2/pi*(j1(x)*(ln(x/2)+Euler)-1/x-x/2+5/64*x**3-...) +// therefore y1(x)-2/pi*j1(x)*ln(x)-1/x is an odd function. +// We use the following function to approximate y1, +// y1(x) = x*U(z)/V(z) + (2/pi)*(j1(x)*ln(x)-1/x), z= x**2 +// where for x in [0,2] (abs err less than 2**-65.89) +// U(z) = U0[0] + U0[1]*z + ... + U0[4]*z**4 +// V(z) = 1 + v0[0]*z + ... + v0[4]*z**5 +// Note: For tiny x, 1/x dominate y1 and hence +// y1(tiny) = -2/pi/tiny, (choose tiny<2**-54) +// 3. For x>=2. +// y1(x) = sqrt(2/(pi*x))*(p1(x)*sin(x1)+q1(x)*cos(x1)) +// where x1 = x-3*pi/4. It is better to compute sin(x1),cos(x1) +// by method mentioned above. + +// J1 returns the order-one Bessel function of the first kind. +// +// Special cases are: +// J1(±Inf) = 0 +// J1(NaN) = NaN +func J1(x float64) float64 { + const ( + TwoM27 = 1.0 / (1 << 27) // 2**-27 0x3e40000000000000 + Two129 = 1 << 129 // 2**129 0x4800000000000000 + // R0/S0 on [0, 2] + R00 = -6.25000000000000000000e-02 // 0xBFB0000000000000 + R01 = 1.40705666955189706048e-03 // 0x3F570D9F98472C61 + R02 = -1.59955631084035597520e-05 // 0xBEF0C5C6BA169668 + R03 = 4.96727999609584448412e-08 // 0x3E6AAAFA46CA0BD9 + S01 = 1.91537599538363460805e-02 // 0x3F939D0B12637E53 + S02 = 1.85946785588630915560e-04 // 0x3F285F56B9CDF664 + S03 = 1.17718464042623683263e-06 // 0x3EB3BFF8333F8498 + S04 = 5.04636257076217042715e-09 // 0x3E35AC88C97DFF2C + S05 = 1.23542274426137913908e-11 // 0x3DAB2ACFCFB97ED8 + ) + // TODO(rsc): Remove manual inlining of IsNaN, IsInf + // when compiler does it for us + // special cases + switch { + case x != x: // IsNaN(x) + return x + case x < -MaxFloat64 || x > MaxFloat64 || x == 0: // IsInf(x, 0) || x == 0: + return 0 + } + + sign := false + if x < 0 { + x = -x + sign = true + } + if x >= 2 { + s, c := Sincos(x) + ss := -s - c + cc := s - c + + // make sure x+x does not overflow + if x < MaxFloat64/2 { + z := Cos(x + x) + if s*c > 0 { + cc = z / ss + } else { + ss = z / cc + } + } + + // j1(x) = 1/sqrt(pi) * (P(1,x)*cc - Q(1,x)*ss) / sqrt(x) + // y1(x) = 1/sqrt(pi) * (P(1,x)*ss + Q(1,x)*cc) / sqrt(x) + + var z float64 + if x > Two129 { + z = (1 / SqrtPi) * cc / Sqrt(x) + } else { + u := pone(x) + v := qone(x) + z = (1 / SqrtPi) * (u*cc - v*ss) / Sqrt(x) + } + if sign { + return -z + } + return z + } + if x < TwoM27 { // |x|<2**-27 + return 0.5 * x // inexact if x!=0 necessary + } + z := x * x + r := z * (R00 + z*(R01+z*(R02+z*R03))) + s := 1.0 + z*(S01+z*(S02+z*(S03+z*(S04+z*S05)))) + r *= x + z = 0.5*x + r/s + if sign { + return -z + } + return z +} + +// Y1 returns the order-one Bessel function of the second kind. +// +// Special cases are: +// Y1(+Inf) = 0 +// Y1(0) = -Inf +// Y1(x < 0) = NaN +// Y1(NaN) = NaN +func Y1(x float64) float64 { + const ( + TwoM54 = 1.0 / (1 << 54) // 2**-54 0x3c90000000000000 + Two129 = 1 << 129 // 2**129 0x4800000000000000 + U00 = -1.96057090646238940668e-01 // 0xBFC91866143CBC8A + U01 = 5.04438716639811282616e-02 // 0x3FA9D3C776292CD1 + U02 = -1.91256895875763547298e-03 // 0xBF5F55E54844F50F + U03 = 2.35252600561610495928e-05 // 0x3EF8AB038FA6B88E + U04 = -9.19099158039878874504e-08 // 0xBE78AC00569105B8 + V00 = 1.99167318236649903973e-02 // 0x3F94650D3F4DA9F0 + V01 = 2.02552581025135171496e-04 // 0x3F2A8C896C257764 + V02 = 1.35608801097516229404e-06 // 0x3EB6C05A894E8CA6 + V03 = 6.22741452364621501295e-09 // 0x3E3ABF1D5BA69A86 + V04 = 1.66559246207992079114e-11 // 0x3DB25039DACA772A + ) + // TODO(rsc): Remove manual inlining of IsNaN, IsInf + // when compiler does it for us + // special cases + switch { + case x < 0 || x != x: // x < 0 || IsNaN(x): + return NaN() + case x > MaxFloat64: // IsInf(x, 1): + return 0 + case x == 0: + return Inf(-1) + } + + if x >= 2 { + s, c := Sincos(x) + ss := -s - c + cc := s - c + + // make sure x+x does not overflow + if x < MaxFloat64/2 { + z := Cos(x + x) + if s*c > 0 { + cc = z / ss + } else { + ss = z / cc + } + } + // y1(x) = sqrt(2/(pi*x))*(p1(x)*sin(x0)+q1(x)*cos(x0)) + // where x0 = x-3pi/4 + // Better formula: + // cos(x0) = cos(x)cos(3pi/4)+sin(x)sin(3pi/4) + // = 1/sqrt(2) * (sin(x) - cos(x)) + // sin(x0) = sin(x)cos(3pi/4)-cos(x)sin(3pi/4) + // = -1/sqrt(2) * (cos(x) + sin(x)) + // To avoid cancellation, use + // sin(x) +- cos(x) = -cos(2x)/(sin(x) -+ cos(x)) + // to compute the worse one. + + var z float64 + if x > Two129 { + z = (1 / SqrtPi) * ss / Sqrt(x) + } else { + u := pone(x) + v := qone(x) + z = (1 / SqrtPi) * (u*ss + v*cc) / Sqrt(x) + } + return z + } + if x <= TwoM54 { // x < 2**-54 + return -(2 / Pi) / x + } + z := x * x + u := U00 + z*(U01+z*(U02+z*(U03+z*U04))) + v := 1 + z*(V00+z*(V01+z*(V02+z*(V03+z*V04)))) + return x*(u/v) + (2/Pi)*(J1(x)*Log(x)-1/x) +} + +// For x >= 8, the asymptotic expansions of pone is +// 1 + 15/128 s**2 - 4725/2**15 s**4 - ..., where s = 1/x. +// We approximate pone by +// pone(x) = 1 + (R/S) +// where R = pr0 + pr1*s**2 + pr2*s**4 + ... + pr5*s**10 +// S = 1 + ps0*s**2 + ... + ps4*s**10 +// and +// | pone(x)-1-R/S | <= 2**(-60.06) + +// for x in [inf, 8]=1/[0,0.125] +var p1R8 = [6]float64{ + 0.00000000000000000000e+00, // 0x0000000000000000 + 1.17187499999988647970e-01, // 0x3FBDFFFFFFFFFCCE + 1.32394806593073575129e+01, // 0x402A7A9D357F7FCE + 4.12051854307378562225e+02, // 0x4079C0D4652EA590 + 3.87474538913960532227e+03, // 0x40AE457DA3A532CC + 7.91447954031891731574e+03, // 0x40BEEA7AC32782DD +} +var p1S8 = [5]float64{ + 1.14207370375678408436e+02, // 0x405C8D458E656CAC + 3.65093083420853463394e+03, // 0x40AC85DC964D274F + 3.69562060269033463555e+04, // 0x40E20B8697C5BB7F + 9.76027935934950801311e+04, // 0x40F7D42CB28F17BB + 3.08042720627888811578e+04, // 0x40DE1511697A0B2D +} + +// for x in [8,4.5454] = 1/[0.125,0.22001] +var p1R5 = [6]float64{ + 1.31990519556243522749e-11, // 0x3DAD0667DAE1CA7D + 1.17187493190614097638e-01, // 0x3FBDFFFFE2C10043 + 6.80275127868432871736e+00, // 0x401B36046E6315E3 + 1.08308182990189109773e+02, // 0x405B13B9452602ED + 5.17636139533199752805e+02, // 0x40802D16D052D649 + 5.28715201363337541807e+02, // 0x408085B8BB7E0CB7 +} +var p1S5 = [5]float64{ + 5.92805987221131331921e+01, // 0x404DA3EAA8AF633D + 9.91401418733614377743e+02, // 0x408EFB361B066701 + 5.35326695291487976647e+03, // 0x40B4E9445706B6FB + 7.84469031749551231769e+03, // 0x40BEA4B0B8A5BB15 + 1.50404688810361062679e+03, // 0x40978030036F5E51 +} + +// for x in[4.5453,2.8571] = 1/[0.2199,0.35001] +var p1R3 = [6]float64{ + 3.02503916137373618024e-09, // 0x3E29FC21A7AD9EDD + 1.17186865567253592491e-01, // 0x3FBDFFF55B21D17B + 3.93297750033315640650e+00, // 0x400F76BCE85EAD8A + 3.51194035591636932736e+01, // 0x40418F489DA6D129 + 9.10550110750781271918e+01, // 0x4056C3854D2C1837 + 4.85590685197364919645e+01, // 0x4048478F8EA83EE5 +} +var p1S3 = [5]float64{ + 3.47913095001251519989e+01, // 0x40416549A134069C + 3.36762458747825746741e+02, // 0x40750C3307F1A75F + 1.04687139975775130551e+03, // 0x40905B7C5037D523 + 8.90811346398256432622e+02, // 0x408BD67DA32E31E9 + 1.03787932439639277504e+02, // 0x4059F26D7C2EED53 +} + +// for x in [2.8570,2] = 1/[0.3499,0.5] +var p1R2 = [6]float64{ + 1.07710830106873743082e-07, // 0x3E7CE9D4F65544F4 + 1.17176219462683348094e-01, // 0x3FBDFF42BE760D83 + 2.36851496667608785174e+00, // 0x4002F2B7F98FAEC0 + 1.22426109148261232917e+01, // 0x40287C377F71A964 + 1.76939711271687727390e+01, // 0x4031B1A8177F8EE2 + 5.07352312588818499250e+00, // 0x40144B49A574C1FE +} +var p1S2 = [5]float64{ + 2.14364859363821409488e+01, // 0x40356FBD8AD5ECDC + 1.25290227168402751090e+02, // 0x405F529314F92CD5 + 2.32276469057162813669e+02, // 0x406D08D8D5A2DBD9 + 1.17679373287147100768e+02, // 0x405D6B7ADA1884A9 + 8.36463893371618283368e+00, // 0x4020BAB1F44E5192 +} + +func pone(x float64) float64 { + var p [6]float64 + var q [5]float64 + if x >= 8 { + p = p1R8 + q = p1S8 + } else if x >= 4.5454 { + p = p1R5 + q = p1S5 + } else if x >= 2.8571 { + p = p1R3 + q = p1S3 + } else if x >= 2 { + p = p1R2 + q = p1S2 + } + z := 1 / (x * x) + r := p[0] + z*(p[1]+z*(p[2]+z*(p[3]+z*(p[4]+z*p[5])))) + s := 1.0 + z*(q[0]+z*(q[1]+z*(q[2]+z*(q[3]+z*q[4])))) + return 1 + r/s +} + +// For x >= 8, the asymptotic expansions of qone is +// 3/8 s - 105/1024 s**3 - ..., where s = 1/x. +// We approximate qone by +// qone(x) = s*(0.375 + (R/S)) +// where R = qr1*s**2 + qr2*s**4 + ... + qr5*s**10 +// S = 1 + qs1*s**2 + ... + qs6*s**12 +// and +// | qone(x)/s -0.375-R/S | <= 2**(-61.13) + +// for x in [inf, 8] = 1/[0,0.125] +var q1R8 = [6]float64{ + 0.00000000000000000000e+00, // 0x0000000000000000 + -1.02539062499992714161e-01, // 0xBFBA3FFFFFFFFDF3 + -1.62717534544589987888e+01, // 0xC0304591A26779F7 + -7.59601722513950107896e+02, // 0xC087BCD053E4B576 + -1.18498066702429587167e+04, // 0xC0C724E740F87415 + -4.84385124285750353010e+04, // 0xC0E7A6D065D09C6A +} +var q1S8 = [6]float64{ + 1.61395369700722909556e+02, // 0x40642CA6DE5BCDE5 + 7.82538599923348465381e+03, // 0x40BE9162D0D88419 + 1.33875336287249578163e+05, // 0x4100579AB0B75E98 + 7.19657723683240939863e+05, // 0x4125F65372869C19 + 6.66601232617776375264e+05, // 0x412457D27719AD5C + -2.94490264303834643215e+05, // 0xC111F9690EA5AA18 +} + +// for x in [8,4.5454] = 1/[0.125,0.22001] +var q1R5 = [6]float64{ + -2.08979931141764104297e-11, // 0xBDB6FA431AA1A098 + -1.02539050241375426231e-01, // 0xBFBA3FFFCB597FEF + -8.05644828123936029840e+00, // 0xC0201CE6CA03AD4B + -1.83669607474888380239e+02, // 0xC066F56D6CA7B9B0 + -1.37319376065508163265e+03, // 0xC09574C66931734F + -2.61244440453215656817e+03, // 0xC0A468E388FDA79D +} +var q1S5 = [6]float64{ + 8.12765501384335777857e+01, // 0x405451B2FF5A11B2 + 1.99179873460485964642e+03, // 0x409F1F31E77BF839 + 1.74684851924908907677e+04, // 0x40D10F1F0D64CE29 + 4.98514270910352279316e+04, // 0x40E8576DAABAD197 + 2.79480751638918118260e+04, // 0x40DB4B04CF7C364B + -4.71918354795128470869e+03, // 0xC0B26F2EFCFFA004 +} + +// for x in [4.5454,2.8571] = 1/[0.2199,0.35001] ??? +var q1R3 = [6]float64{ + -5.07831226461766561369e-09, // 0xBE35CFA9D38FC84F + -1.02537829820837089745e-01, // 0xBFBA3FEB51AEED54 + -4.61011581139473403113e+00, // 0xC01270C23302D9FF + -5.78472216562783643212e+01, // 0xC04CEC71C25D16DA + -2.28244540737631695038e+02, // 0xC06C87D34718D55F + -2.19210128478909325622e+02, // 0xC06B66B95F5C1BF6 +} +var q1S3 = [6]float64{ + 4.76651550323729509273e+01, // 0x4047D523CCD367E4 + 6.73865112676699709482e+02, // 0x40850EEBC031EE3E + 3.38015286679526343505e+03, // 0x40AA684E448E7C9A + 5.54772909720722782367e+03, // 0x40B5ABBAA61D54A6 + 1.90311919338810798763e+03, // 0x409DBC7A0DD4DF4B + -1.35201191444307340817e+02, // 0xC060E670290A311F +} + +// for x in [2.8570,2] = 1/[0.3499,0.5] +var q1R2 = [6]float64{ + -1.78381727510958865572e-07, // 0xBE87F12644C626D2 + -1.02517042607985553460e-01, // 0xBFBA3E8E9148B010 + -2.75220568278187460720e+00, // 0xC006048469BB4EDA + -1.96636162643703720221e+01, // 0xC033A9E2C168907F + -4.23253133372830490089e+01, // 0xC04529A3DE104AAA + -2.13719211703704061733e+01, // 0xC0355F3639CF6E52 +} +var q1S2 = [6]float64{ + 2.95333629060523854548e+01, // 0x403D888A78AE64FF + 2.52981549982190529136e+02, // 0x406F9F68DB821CBA + 7.57502834868645436472e+02, // 0x4087AC05CE49A0F7 + 7.39393205320467245656e+02, // 0x40871B2548D4C029 + 1.55949003336666123687e+02, // 0x40637E5E3C3ED8D4 + -4.95949898822628210127e+00, // 0xC013D686E71BE86B +} + +func qone(x float64) float64 { + var p, q [6]float64 + if x >= 8 { + p = q1R8 + q = q1S8 + } else if x >= 4.5454 { + p = q1R5 + q = q1S5 + } else if x >= 2.8571 { + p = q1R3 + q = q1S3 + } else if x >= 2 { + p = q1R2 + q = q1S2 + } + z := 1 / (x * x) + r := p[0] + z*(p[1]+z*(p[2]+z*(p[3]+z*(p[4]+z*p[5])))) + s := 1 + z*(q[0]+z*(q[1]+z*(q[2]+z*(q[3]+z*(q[4]+z*q[5]))))) + return (0.375 + r/s) / x +} diff --git a/libgo/go/math/jn.go b/libgo/go/math/jn.go new file mode 100644 index 000000000..9024af3c2 --- /dev/null +++ b/libgo/go/math/jn.go @@ -0,0 +1,310 @@ +// 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 + +/* + Bessel function of the first and second kinds of order n. +*/ + +// The original C code and the long comment below are +// from FreeBSD's /usr/src/lib/msun/src/e_jn.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_jn(n, x), __ieee754_yn(n, x) +// floating point Bessel's function of the 1st and 2nd kind +// of order n +// +// Special cases: +// y0(0)=y1(0)=yn(n,0) = -inf with division by zero signal; +// y0(-ve)=y1(-ve)=yn(n,-ve) are NaN with invalid signal. +// Note 2. About jn(n,x), yn(n,x) +// For n=0, j0(x) is called, +// for n=1, j1(x) is called, +// for nx, a continued fraction approximation to +// j(n,x)/j(n-1,x) is evaluated and then backward +// recursion is used starting from a supposed value +// for j(n,x). The resulting value of j(0,x) is +// compared with the actual value to correct the +// supposed value of j(n,x). +// +// yn(n,x) is similar in all respects, except +// that forward recursion is used for all +// values of n>1. + +// Jn returns the order-n Bessel function of the first kind. +// +// Special cases are: +// Jn(n, ±Inf) = 0 +// Jn(n, NaN) = NaN +func Jn(n int, x float64) float64 { + const ( + TwoM29 = 1.0 / (1 << 29) // 2**-29 0x3e10000000000000 + Two302 = 1 << 302 // 2**302 0x52D0000000000000 + ) + // TODO(rsc): Remove manual inlining of IsNaN, IsInf + // when compiler does it for us + // special cases + switch { + case x != x: // IsNaN(x) + return x + case x < -MaxFloat64 || x > MaxFloat64: // IsInf(x, 0): + return 0 + } + // J(-n, x) = (-1)**n * J(n, x), J(n, -x) = (-1)**n * J(n, x) + // Thus, J(-n, x) = J(n, -x) + + if n == 0 { + return J0(x) + } + if x == 0 { + return 0 + } + if n < 0 { + n, x = -n, -x + } + if n == 1 { + return J1(x) + } + sign := false + if x < 0 { + x = -x + if n&1 == 1 { + sign = true // odd n and negative x + } + } + var b float64 + if float64(n) <= x { + // Safe to use J(n+1,x)=2n/x *J(n,x)-J(n-1,x) + if x >= Two302 { // x > 2**302 + + // (x >> n**2) + // Jn(x) = cos(x-(2n+1)*pi/4)*sqrt(2/x*pi) + // Yn(x) = sin(x-(2n+1)*pi/4)*sqrt(2/x*pi) + // Let s=sin(x), c=cos(x), + // xn=x-(2n+1)*pi/4, sqt2 = sqrt(2),then + // + // n sin(xn)*sqt2 cos(xn)*sqt2 + // ---------------------------------- + // 0 s-c c+s + // 1 -s-c -c+s + // 2 -s+c -c-s + // 3 s+c c-s + + var temp float64 + switch n & 3 { + case 0: + temp = Cos(x) + Sin(x) + case 1: + temp = -Cos(x) + Sin(x) + case 2: + temp = -Cos(x) - Sin(x) + case 3: + temp = Cos(x) - Sin(x) + } + b = (1 / SqrtPi) * temp / Sqrt(x) + } else { + b = J1(x) + for i, a := 1, J0(x); i < n; i++ { + a, b = b, b*(float64(i+i)/x)-a // avoid underflow + } + } + } else { + if x < TwoM29 { // x < 2**-29 + // x is tiny, return the first Taylor expansion of J(n,x) + // J(n,x) = 1/n!*(x/2)**n - ... + + if n > 33 { // underflow + b = 0 + } else { + temp := x * 0.5 + b = temp + a := 1.0 + for i := 2; i <= n; i++ { + a *= float64(i) // a = n! + b *= temp // b = (x/2)**n + } + b /= a + } + } else { + // use backward recurrence + // x x**2 x**2 + // J(n,x)/J(n-1,x) = ---- ------ ------ ..... + // 2n - 2(n+1) - 2(n+2) + // + // 1 1 1 + // (for large x) = ---- ------ ------ ..... + // 2n 2(n+1) 2(n+2) + // -- - ------ - ------ - + // x x x + // + // Let w = 2n/x and h=2/x, then the above quotient + // is equal to the continued fraction: + // 1 + // = ----------------------- + // 1 + // w - ----------------- + // 1 + // w+h - --------- + // w+2h - ... + // + // To determine how many terms needed, let + // Q(0) = w, Q(1) = w(w+h) - 1, + // Q(k) = (w+k*h)*Q(k-1) - Q(k-2), + // When Q(k) > 1e4 good for single + // When Q(k) > 1e9 good for double + // When Q(k) > 1e17 good for quadruple + + // determine k + w := float64(n+n) / x + h := 2 / x + q0 := w + z := w + h + q1 := w*z - 1 + k := 1 + for q1 < 1e9 { + k += 1 + z += h + q0, q1 = q1, z*q1-q0 + } + m := n + n + t := 0.0 + for i := 2 * (n + k); i >= m; i -= 2 { + t = 1 / (float64(i)/x - t) + } + a := t + b = 1 + // estimate log((2/x)**n*n!) = n*log(2/x)+n*ln(n) + // Hence, if n*(log(2n/x)) > ... + // single 8.8722839355e+01 + // double 7.09782712893383973096e+02 + // long double 1.1356523406294143949491931077970765006170e+04 + // then recurrent value may overflow and the result is + // likely underflow to zero + + tmp := float64(n) + v := 2 / x + tmp = tmp * Log(Fabs(v*tmp)) + if tmp < 7.09782712893383973096e+02 { + for i := n - 1; i > 0; i-- { + di := float64(i + i) + a, b = b, b*di/x-a + di -= 2 + } + } else { + for i := n - 1; i > 0; i-- { + di := float64(i + i) + a, b = b, b*di/x-a + di -= 2 + // scale b to avoid spurious overflow + if b > 1e100 { + a /= b + t /= b + b = 1 + } + } + } + b = t * J0(x) / b + } + } + if sign { + return -b + } + return b +} + +// Yn returns the order-n Bessel function of the second kind. +// +// Special cases are: +// Yn(n, +Inf) = 0 +// Yn(n > 0, 0) = -Inf +// Yn(n < 0, 0) = +Inf if n is odd, -Inf if n is even +// Y1(n, x < 0) = NaN +// Y1(n, NaN) = NaN +func Yn(n int, x float64) float64 { + const Two302 = 1 << 302 // 2**302 0x52D0000000000000 + // TODO(rsc): Remove manual inlining of IsNaN, IsInf + // when compiler does it for us + // special cases + switch { + case x < 0 || x != x: // x < 0 || IsNaN(x): + return NaN() + case x > MaxFloat64: // IsInf(x, 1) + return 0 + } + + if n == 0 { + return Y0(x) + } + if x == 0 { + if n < 0 && n&1 == 1 { + return Inf(1) + } + return Inf(-1) + } + sign := false + if n < 0 { + n = -n + if n&1 == 1 { + sign = true // sign true if n < 0 && |n| odd + } + } + if n == 1 { + if sign { + return -Y1(x) + } + return Y1(x) + } + var b float64 + if x >= Two302 { // x > 2**302 + // (x >> n**2) + // Jn(x) = cos(x-(2n+1)*pi/4)*sqrt(2/x*pi) + // Yn(x) = sin(x-(2n+1)*pi/4)*sqrt(2/x*pi) + // Let s=sin(x), c=cos(x), + // xn=x-(2n+1)*pi/4, sqt2 = sqrt(2),then + // + // n sin(xn)*sqt2 cos(xn)*sqt2 + // ---------------------------------- + // 0 s-c c+s + // 1 -s-c -c+s + // 2 -s+c -c-s + // 3 s+c c-s + + var temp float64 + switch n & 3 { + case 0: + temp = Sin(x) - Cos(x) + case 1: + temp = -Sin(x) - Cos(x) + case 2: + temp = -Sin(x) + Cos(x) + case 3: + temp = Sin(x) + Cos(x) + } + b = (1 / SqrtPi) * temp / Sqrt(x) + } else { + a := Y0(x) + b = Y1(x) + // quit if b is -inf + for i := 1; i < n && b >= -MaxFloat64; i++ { // for i := 1; i < n && !IsInf(b, -1); i++ { + a, b = b, (float64(i+i)/x)*b-a + } + } + if sign { + return -b + } + return b +} diff --git a/libgo/go/math/ldexp.go b/libgo/go/math/ldexp.go new file mode 100644 index 000000000..96c95cad4 --- /dev/null +++ b/libgo/go/math/ldexp.go @@ -0,0 +1,45 @@ +// Copyright 2009 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 + +// Ldexp is the inverse of Frexp. +// It returns frac × 2**exp. +// +// Special cases are: +// Ldexp(±0, exp) = ±0 +// Ldexp(±Inf, exp) = ±Inf +// Ldexp(NaN, exp) = NaN +func Ldexp(frac float64, exp int) float64 { + // TODO(rsc): Remove manual inlining of IsNaN, IsInf + // when compiler does it for us + // special cases + switch { + case frac == 0: + return frac // correctly return -0 + case frac < -MaxFloat64 || frac > MaxFloat64 || frac != frac: // IsInf(frac, 0) || IsNaN(frac): + return frac + } + frac, e := normalize(frac) + exp += e + x := Float64bits(frac) + exp += int(x>>shift)&mask - bias + if exp < -1074 { + return Copysign(0, frac) // underflow + } + if exp > 1023 { // overflow + if frac < 0 { + return Inf(-1) + } + return Inf(1) + } + var m float64 = 1 + if exp < -1022 { // denormal + exp += 52 + m = 1.0 / (1 << 52) // 2**-52 + } + x &^= mask << shift + x |= uint64(exp+bias) << shift + return m * Float64frombits(x) +} diff --git a/libgo/go/math/lgamma.go b/libgo/go/math/lgamma.go new file mode 100644 index 000000000..dc30f468f --- /dev/null +++ b/libgo/go/math/lgamma.go @@ -0,0 +1,350 @@ +// 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 logarithm of the Gamma function. +*/ + +// The original C code and the long comment below are +// from FreeBSD's /usr/src/lib/msun/src/e_lgamma_r.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_lgamma_r(x, signgamp) +// Reentrant version of the logarithm of the Gamma function +// with user provided pointer for the sign of Gamma(x). +// +// Method: +// 1. Argument Reduction for 0 < x <= 8 +// Since gamma(1+s)=s*gamma(s), for x in [0,8], we may +// reduce x to a number in [1.5,2.5] by +// lgamma(1+s) = log(s) + lgamma(s) +// for example, +// lgamma(7.3) = log(6.3) + lgamma(6.3) +// = log(6.3*5.3) + lgamma(5.3) +// = log(6.3*5.3*4.3*3.3*2.3) + lgamma(2.3) +// 2. Polynomial approximation of lgamma around its +// minimum (ymin=1.461632144968362245) to maintain monotonicity. +// On [ymin-0.23, ymin+0.27] (i.e., [1.23164,1.73163]), use +// Let z = x-ymin; +// lgamma(x) = -1.214862905358496078218 + z**2*poly(z) +// poly(z) is a 14 degree polynomial. +// 2. Rational approximation in the primary interval [2,3] +// We use the following approximation: +// s = x-2.0; +// lgamma(x) = 0.5*s + s*P(s)/Q(s) +// with accuracy +// |P/Q - (lgamma(x)-0.5s)| < 2**-61.71 +// Our algorithms are based on the following observation +// +// zeta(2)-1 2 zeta(3)-1 3 +// lgamma(2+s) = s*(1-Euler) + --------- * s - --------- * s + ... +// 2 3 +// +// where Euler = 0.5772156649... is the Euler constant, which +// is very close to 0.5. +// +// 3. For x>=8, we have +// lgamma(x)~(x-0.5)log(x)-x+0.5*log(2pi)+1/(12x)-1/(360x**3)+.... +// (better formula: +// lgamma(x)~(x-0.5)*(log(x)-1)-.5*(log(2pi)-1) + ...) +// Let z = 1/x, then we approximation +// f(z) = lgamma(x) - (x-0.5)(log(x)-1) +// by +// 3 5 11 +// w = w0 + w1*z + w2*z + w3*z + ... + w6*z +// where +// |w - f(z)| < 2**-58.74 +// +// 4. For negative x, since (G is gamma function) +// -x*G(-x)*G(x) = pi/sin(pi*x), +// we have +// G(x) = pi/(sin(pi*x)*(-x)*G(-x)) +// since G(-x) is positive, sign(G(x)) = sign(sin(pi*x)) for x<0 +// Hence, for x<0, signgam = sign(sin(pi*x)) and +// lgamma(x) = log(|Gamma(x)|) +// = log(pi/(|x*sin(pi*x)|)) - lgamma(-x); +// Note: one should avoid computing pi*(-x) directly in the +// computation of sin(pi*(-x)). +// +// 5. Special Cases +// lgamma(2+s) ~ s*(1-Euler) for tiny s +// lgamma(1)=lgamma(2)=0 +// lgamma(x) ~ -log(x) for tiny x +// lgamma(0) = lgamma(inf) = inf +// lgamma(-integer) = +-inf +// +// + +// Lgamma returns the natural logarithm and sign (-1 or +1) of Gamma(x). +// +// Special cases are: +// Lgamma(+Inf) = +Inf +// Lgamma(0) = +Inf +// Lgamma(-integer) = +Inf +// Lgamma(-Inf) = -Inf +// Lgamma(NaN) = NaN +func Lgamma(x float64) (lgamma float64, sign int) { + const ( + Ymin = 1.461632144968362245 + Two52 = 1 << 52 // 0x4330000000000000 ~4.5036e+15 + Two53 = 1 << 53 // 0x4340000000000000 ~9.0072e+15 + Two58 = 1 << 58 // 0x4390000000000000 ~2.8823e+17 + Tiny = 1.0 / (1 << 70) // 0x3b90000000000000 ~8.47033e-22 + A0 = 7.72156649015328655494e-02 // 0x3FB3C467E37DB0C8 + A1 = 3.22467033424113591611e-01 // 0x3FD4A34CC4A60FAD + A2 = 6.73523010531292681824e-02 // 0x3FB13E001A5562A7 + A3 = 2.05808084325167332806e-02 // 0x3F951322AC92547B + A4 = 7.38555086081402883957e-03 // 0x3F7E404FB68FEFE8 + A5 = 2.89051383673415629091e-03 // 0x3F67ADD8CCB7926B + A6 = 1.19270763183362067845e-03 // 0x3F538A94116F3F5D + A7 = 5.10069792153511336608e-04 // 0x3F40B6C689B99C00 + A8 = 2.20862790713908385557e-04 // 0x3F2CF2ECED10E54D + A9 = 1.08011567247583939954e-04 // 0x3F1C5088987DFB07 + A10 = 2.52144565451257326939e-05 // 0x3EFA7074428CFA52 + A11 = 4.48640949618915160150e-05 // 0x3F07858E90A45837 + Tc = 1.46163214496836224576e+00 // 0x3FF762D86356BE3F + Tf = -1.21486290535849611461e-01 // 0xBFBF19B9BCC38A42 + // Tt = -(tail of Tf) + Tt = -3.63867699703950536541e-18 // 0xBC50C7CAA48A971F + T0 = 4.83836122723810047042e-01 // 0x3FDEF72BC8EE38A2 + T1 = -1.47587722994593911752e-01 // 0xBFC2E4278DC6C509 + T2 = 6.46249402391333854778e-02 // 0x3FB08B4294D5419B + T3 = -3.27885410759859649565e-02 // 0xBFA0C9A8DF35B713 + T4 = 1.79706750811820387126e-02 // 0x3F9266E7970AF9EC + T5 = -1.03142241298341437450e-02 // 0xBF851F9FBA91EC6A + T6 = 6.10053870246291332635e-03 // 0x3F78FCE0E370E344 + T7 = -3.68452016781138256760e-03 // 0xBF6E2EFFB3E914D7 + T8 = 2.25964780900612472250e-03 // 0x3F6282D32E15C915 + T9 = -1.40346469989232843813e-03 // 0xBF56FE8EBF2D1AF1 + T10 = 8.81081882437654011382e-04 // 0x3F4CDF0CEF61A8E9 + T11 = -5.38595305356740546715e-04 // 0xBF41A6109C73E0EC + T12 = 3.15632070903625950361e-04 // 0x3F34AF6D6C0EBBF7 + T13 = -3.12754168375120860518e-04 // 0xBF347F24ECC38C38 + T14 = 3.35529192635519073543e-04 // 0x3F35FD3EE8C2D3F4 + U0 = -7.72156649015328655494e-02 // 0xBFB3C467E37DB0C8 + U1 = 6.32827064025093366517e-01 // 0x3FE4401E8B005DFF + U2 = 1.45492250137234768737e+00 // 0x3FF7475CD119BD6F + U3 = 9.77717527963372745603e-01 // 0x3FEF497644EA8450 + U4 = 2.28963728064692451092e-01 // 0x3FCD4EAEF6010924 + U5 = 1.33810918536787660377e-02 // 0x3F8B678BBF2BAB09 + V1 = 2.45597793713041134822e+00 // 0x4003A5D7C2BD619C + V2 = 2.12848976379893395361e+00 // 0x40010725A42B18F5 + V3 = 7.69285150456672783825e-01 // 0x3FE89DFBE45050AF + V4 = 1.04222645593369134254e-01 // 0x3FBAAE55D6537C88 + V5 = 3.21709242282423911810e-03 // 0x3F6A5ABB57D0CF61 + S0 = -7.72156649015328655494e-02 // 0xBFB3C467E37DB0C8 + S1 = 2.14982415960608852501e-01 // 0x3FCB848B36E20878 + S2 = 3.25778796408930981787e-01 // 0x3FD4D98F4F139F59 + S3 = 1.46350472652464452805e-01 // 0x3FC2BB9CBEE5F2F7 + S4 = 2.66422703033638609560e-02 // 0x3F9B481C7E939961 + S5 = 1.84028451407337715652e-03 // 0x3F5E26B67368F239 + S6 = 3.19475326584100867617e-05 // 0x3F00BFECDD17E945 + R1 = 1.39200533467621045958e+00 // 0x3FF645A762C4AB74 + R2 = 7.21935547567138069525e-01 // 0x3FE71A1893D3DCDC + R3 = 1.71933865632803078993e-01 // 0x3FC601EDCCFBDF27 + R4 = 1.86459191715652901344e-02 // 0x3F9317EA742ED475 + R5 = 7.77942496381893596434e-04 // 0x3F497DDACA41A95B + R6 = 7.32668430744625636189e-06 // 0x3EDEBAF7A5B38140 + W0 = 4.18938533204672725052e-01 // 0x3FDACFE390C97D69 + W1 = 8.33333333333329678849e-02 // 0x3FB555555555553B + W2 = -2.77777777728775536470e-03 // 0xBF66C16C16B02E5C + W3 = 7.93650558643019558500e-04 // 0x3F4A019F98CF38B6 + W4 = -5.95187557450339963135e-04 // 0xBF4380CB8C0FE741 + W5 = 8.36339918996282139126e-04 // 0x3F4B67BA4CDAD5D1 + W6 = -1.63092934096575273989e-03 // 0xBF5AB89D0B9E43E4 + ) + // TODO(rsc): Remove manual inlining of IsNaN, IsInf + // when compiler does it for us + // special cases + sign = 1 + switch { + case x != x: // IsNaN(x): + lgamma = x + return + case x < -MaxFloat64 || x > MaxFloat64: // IsInf(x, 0): + lgamma = x + return + case x == 0: + lgamma = Inf(1) + return + } + + neg := false + if x < 0 { + x = -x + neg = true + } + + if x < Tiny { // if |x| < 2**-70, return -log(|x|) + if neg { + sign = -1 + } + lgamma = -Log(x) + return + } + var nadj float64 + if neg { + if x >= Two52 { // |x| >= 2**52, must be -integer + lgamma = Inf(1) + return + } + t := sinPi(x) + if t == 0 { + lgamma = Inf(1) // -integer + return + } + nadj = Log(Pi / Fabs(t*x)) + if t < 0 { + sign = -1 + } + } + + switch { + case x == 1 || x == 2: // purge off 1 and 2 + lgamma = 0 + return + case x < 2: // use lgamma(x) = lgamma(x+1) - log(x) + var y float64 + var i int + if x <= 0.9 { + lgamma = -Log(x) + switch { + case x >= (Ymin - 1 + 0.27): // 0.7316 <= x <= 0.9 + y = 1 - x + i = 0 + case x >= (Ymin - 1 - 0.27): // 0.2316 <= x < 0.7316 + y = x - (Tc - 1) + i = 1 + default: // 0 < x < 0.2316 + y = x + i = 2 + } + } else { + lgamma = 0 + switch { + case x >= (Ymin + 0.27): // 1.7316 <= x < 2 + y = 2 - x + i = 0 + case x >= (Ymin - 0.27): // 1.2316 <= x < 1.7316 + y = x - Tc + i = 1 + default: // 0.9 < x < 1.2316 + y = x - 1 + i = 2 + } + } + switch i { + case 0: + z := y * y + p1 := A0 + z*(A2+z*(A4+z*(A6+z*(A8+z*A10)))) + p2 := z * (A1 + z*(A3+z*(A5+z*(A7+z*(A9+z*A11))))) + p := y*p1 + p2 + lgamma += (p - 0.5*y) + case 1: + z := y * y + w := z * y + p1 := T0 + w*(T3+w*(T6+w*(T9+w*T12))) // parallel comp + p2 := T1 + w*(T4+w*(T7+w*(T10+w*T13))) + p3 := T2 + w*(T5+w*(T8+w*(T11+w*T14))) + p := z*p1 - (Tt - w*(p2+y*p3)) + lgamma += (Tf + p) + case 2: + p1 := y * (U0 + y*(U1+y*(U2+y*(U3+y*(U4+y*U5))))) + p2 := 1 + y*(V1+y*(V2+y*(V3+y*(V4+y*V5)))) + lgamma += (-0.5*y + p1/p2) + } + case x < 8: // 2 <= x < 8 + i := int(x) + y := x - float64(i) + p := y * (S0 + y*(S1+y*(S2+y*(S3+y*(S4+y*(S5+y*S6)))))) + q := 1 + y*(R1+y*(R2+y*(R3+y*(R4+y*(R5+y*R6))))) + lgamma = 0.5*y + p/q + z := 1.0 // Lgamma(1+s) = Log(s) + Lgamma(s) + switch i { + case 7: + z *= (y + 6) + fallthrough + case 6: + z *= (y + 5) + fallthrough + case 5: + z *= (y + 4) + fallthrough + case 4: + z *= (y + 3) + fallthrough + case 3: + z *= (y + 2) + lgamma += Log(z) + } + case x < Two58: // 8 <= x < 2**58 + t := Log(x) + z := 1 / x + y := z * z + w := W0 + z*(W1+y*(W2+y*(W3+y*(W4+y*(W5+y*W6))))) + lgamma = (x-0.5)*(t-1) + w + default: // 2**58 <= x <= Inf + lgamma = x * (Log(x) - 1) + } + if neg { + lgamma = nadj - lgamma + } + return +} + +// sinPi(x) is a helper function for negative x +func sinPi(x float64) float64 { + const ( + Two52 = 1 << 52 // 0x4330000000000000 ~4.5036e+15 + Two53 = 1 << 53 // 0x4340000000000000 ~9.0072e+15 + ) + if x < 0.25 { + return -Sin(Pi * x) + } + + // argument reduction + z := Floor(x) + var n int + if z != x { // inexact + x = Fmod(x, 2) + n = int(x * 4) + } else { + if x >= Two53 { // x must be even + x = 0 + n = 0 + } else { + if x < Two52 { + z = x + Two52 // exact + } + n = int(1 & Float64bits(z)) + x = float64(n) + n <<= 2 + } + } + switch n { + case 0: + x = Sin(Pi * x) + case 1, 2: + x = Cos(Pi * (0.5 - x)) + case 3, 4: + x = Sin(Pi * (1 - x)) + case 5, 6: + x = -Cos(Pi * (x - 1.5)) + default: + x = Sin(Pi * (x - 2)) + } + return -x +} diff --git a/libgo/go/math/log.go b/libgo/go/math/log.go new file mode 100644 index 000000000..39d94512d --- /dev/null +++ b/libgo/go/math/log.go @@ -0,0 +1,123 @@ +// Copyright 2009 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 logarithm. +*/ + +// The original C code, the long comment, and the constants +// below are from FreeBSD's /usr/src/lib/msun/src/e_log.c +// and came with this notice. The go code is a simpler +// 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_log(x) +// Return the logrithm of x +// +// Method : +// 1. Argument Reduction: find k and f such that +// x = 2**k * (1+f), +// where sqrt(2)/2 < 1+f < sqrt(2) . +// +// 2. Approximation of log(1+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) ~ L1*s +L2*s +L3*s +L4*s +L5*s +L6*s +L7*s +// (the values of L1 to L7 are listed in the program) and +// | 2 14 | -58.45 +// | L1*s +...+L7*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 +// log(1+f) = f - s*(f - R) (if f is not too large) +// log(1+f) = f - (hfsq - s*(hfsq+R)). (better accuracy) +// +// 3. Finally, log(x) = k*Ln2 + log(1+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: +// log(x) is NaN with signal if x < 0 (including -INF) ; +// log(+INF) is +INF; log(0) is -INF with signal; +// log(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. + +// Log returns the natural logarithm of x. +// +// Special cases are: +// Log(+Inf) = +Inf +// Log(0) = -Inf +// Log(x < 0) = NaN +// Log(NaN) = NaN +func Log(x float64) float64 { + const ( + Ln2Hi = 6.93147180369123816490e-01 /* 3fe62e42 fee00000 */ + Ln2Lo = 1.90821492927058770002e-10 /* 3dea39ef 35793c76 */ + L1 = 6.666666666666735130e-01 /* 3FE55555 55555593 */ + L2 = 3.999999999940941908e-01 /* 3FD99999 9997FA04 */ + L3 = 2.857142874366239149e-01 /* 3FD24924 94229359 */ + L4 = 2.222219843214978396e-01 /* 3FCC71C5 1D8E78AF */ + L5 = 1.818357216161805012e-01 /* 3FC74664 96CB03DE */ + L6 = 1.531383769920937332e-01 /* 3FC39A09 D078C69F */ + L7 = 1.479819860511658591e-01 /* 3FC2F112 DF3E5244 */ + ) + + // TODO(rsc): Remove manual inlining of IsNaN, IsInf + // when compiler does it for us + // special cases + switch { + case x != x || x > MaxFloat64: // IsNaN(x) || IsInf(x, 1): + return x + case x < 0: + return NaN() + case x == 0: + return Inf(-1) + } + + // reduce + f1, ki := Frexp(x) + if f1 < Sqrt2/2 { + f1 *= 2 + ki-- + } + f := f1 - 1 + k := float64(ki) + + // compute + s := f / (2 + f) + s2 := s * s + s4 := s2 * s2 + t1 := s2 * (L1 + s4*(L3+s4*(L5+s4*L7))) + t2 := s4 * (L2 + s4*(L4+s4*L6)) + R := t1 + t2 + hfsq := 0.5 * f * f + return k*Ln2Hi - ((hfsq - (s*(hfsq+R) + k*Ln2Lo)) - f) +} diff --git a/libgo/go/math/log10.go b/libgo/go/math/log10.go new file mode 100644 index 000000000..6d18baae2 --- /dev/null +++ b/libgo/go/math/log10.go @@ -0,0 +1,13 @@ +// Copyright 2009 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 + +// Log10 returns the decimal logarithm of x. +// The special cases are the same as for Log. +func Log10(x float64) float64 { return Log(x) * (1 / Ln10) } + +// Log2 returns the binary logarithm of x. +// The special cases are the same as for Log. +func Log2(x float64) float64 { return Log(x) * (1 / Ln2) } diff --git a/libgo/go/math/log10_decl.go b/libgo/go/math/log10_decl.go new file mode 100644 index 000000000..5aec94e1c --- /dev/null +++ b/libgo/go/math/log10_decl.go @@ -0,0 +1,8 @@ +// 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 Log10(x float64) float64 +func Log2(x float64) float64 diff --git a/libgo/go/math/log1p.go b/libgo/go/math/log1p.go new file mode 100644 index 000000000..e1fc275d0 --- /dev/null +++ b/libgo/go/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) +} diff --git a/libgo/go/math/logb.go b/libgo/go/math/logb.go new file mode 100644 index 000000000..072281ddf --- /dev/null +++ b/libgo/go/math/logb.go @@ -0,0 +1,54 @@ +// 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 + +// Logb(x) returns the binary exponent of x. +// +// Special cases are: +// Logb(±Inf) = +Inf +// Logb(0) = -Inf +// Logb(NaN) = NaN +func Logb(x float64) float64 { + // TODO(rsc): Remove manual inlining of IsNaN, IsInf + // when compiler does it for us + // special cases + switch { + case x == 0: + return Inf(-1) + case x < -MaxFloat64 || x > MaxFloat64: // IsInf(x, 0): + return Inf(1) + case x != x: // IsNaN(x): + return x + } + return float64(ilogb(x)) +} + +// Ilogb(x) returns the binary exponent of x as an integer. +// +// Special cases are: +// Ilogb(±Inf) = MaxInt32 +// Ilogb(0) = MinInt32 +// Ilogb(NaN) = MaxInt32 +func Ilogb(x float64) int { + // TODO(rsc): Remove manual inlining of IsNaN, IsInf + // when compiler does it for us + // special cases + switch { + case x == 0: + return MinInt32 + case x != x: // IsNaN(x): + return MaxInt32 + case x < -MaxFloat64 || x > MaxFloat64: // IsInf(x, 0): + return MaxInt32 + } + return ilogb(x) +} + +// logb returns the binary exponent of x. It assumes x is finite and +// non-zero. +func ilogb(x float64) int { + x, exp := normalize(x) + return int((Float64bits(x)>>shift)&mask) - bias + exp +} diff --git a/libgo/go/math/modf.go b/libgo/go/math/modf.go new file mode 100644 index 000000000..315174b70 --- /dev/null +++ b/libgo/go/math/modf.go @@ -0,0 +1,33 @@ +// Copyright 2009 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 + +// Modf returns integer and fractional floating-point numbers +// that sum to f. Both values have the same sign as f. +// +// Special cases are: +// Modf(+Inf) = +Inf, NaN +// Modf(-Inf) = -Inf, NaN +// Modf(NaN) = NaN, NaN +func Modf(f float64) (int float64, frac float64) { + if f < 1 { + if f < 0 { + int, frac = Modf(-f) + return -int, -frac + } + return 0, f + } + + x := Float64bits(f) + e := uint(x>>shift)&mask - bias + + // Keep the top 12+e bits, the integer part; clear the rest. + if e < 64-12 { + x &^= 1<<(64-12-e) - 1 + } + int = Float64frombits(x) + frac = f - int + return +} diff --git a/libgo/go/math/nextafter.go b/libgo/go/math/nextafter.go new file mode 100644 index 000000000..86114340c --- /dev/null +++ b/libgo/go/math/nextafter.go @@ -0,0 +1,29 @@ +// 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 + +// Nextafter returns the next representable value after x towards y. +// If x == y, then x is returned. +// +// Special cases are: +// Nextafter(NaN, y) = NaN +// Nextafter(x, NaN) = NaN +func Nextafter(x, y float64) (r float64) { + // TODO(rsc): Remove manual inlining of IsNaN + // when compiler does it for us + switch { + case x != x || y != y: // IsNaN(x) || IsNaN(y): // special case + r = NaN() + case x == y: + r = x + case x == 0: + r = Copysign(Float64frombits(1), y) + case (y > x) == (x > 0): + r = Float64frombits(Float64bits(x) + 1) + default: + r = Float64frombits(Float64bits(x) - 1) + } + return r +} diff --git a/libgo/go/math/pow.go b/libgo/go/math/pow.go new file mode 100644 index 000000000..06b107401 --- /dev/null +++ b/libgo/go/math/pow.go @@ -0,0 +1,139 @@ +// Copyright 2009 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 isOddInt(x float64) bool { + xi, xf := Modf(x) + return xf == 0 && int64(xi)&1 == 1 +} + +// Special cases taken from FreeBSD's /usr/src/lib/msun/src/e_pow.c +// updated by IEEE Std. 754-2008 "Section 9.2.1 Special values". + +// Pow returns x**y, the base-x exponential of y. +// +// Special cases are (in order): +// Pow(x, ±0) = 1 for any x +// Pow(1, y) = 1 for any y +// Pow(x, 1) = x for any x +// Pow(NaN, y) = NaN +// Pow(x, NaN) = NaN +// Pow(±0, y) = ±Inf for y an odd integer < 0 +// Pow(±0, -Inf) = +Inf +// Pow(±0, +Inf) = +0 +// Pow(±0, y) = +Inf for finite y < 0 and not an odd integer +// Pow(±0, y) = ±0 for y an odd integer > 0 +// Pow(±0, y) = +0 for finite y > 0 and not an odd integer +// Pow(-1, ±Inf) = 1 +// Pow(x, +Inf) = +Inf for |x| > 1 +// Pow(x, -Inf) = +0 for |x| > 1 +// Pow(x, +Inf) = +0 for |x| < 1 +// Pow(x, -Inf) = +Inf for |x| < 1 +// Pow(+Inf, y) = +Inf for y > 0 +// Pow(+Inf, y) = +0 for y < 0 +// Pow(-Inf, y) = Pow(-0, -y) +// Pow(x, y) = NaN for finite x < 0 and finite non-integer y +func Pow(x, y float64) float64 { + // TODO(rsc): Remove manual inlining of IsNaN, IsInf + // when compiler does it for us + switch { + case y == 0 || x == 1: + return 1 + case y == 1: + return x + case y == 0.5: + return Sqrt(x) + case y == -0.5: + return 1 / Sqrt(x) + case x != x || y != y: // IsNaN(x) || IsNaN(y): + return NaN() + case x == 0: + switch { + case y < 0: + if isOddInt(y) { + return Copysign(Inf(1), x) + } + return Inf(1) + case y > 0: + if isOddInt(y) { + return x + } + return 0 + } + case y > MaxFloat64 || y < -MaxFloat64: // IsInf(y, 0): + switch { + case x == -1: + return 1 + case (Fabs(x) < 1) == IsInf(y, 1): + return 0 + default: + return Inf(1) + } + case x > MaxFloat64 || x < -MaxFloat64: // IsInf(x, 0): + if IsInf(x, -1) { + return Pow(1/x, -y) // Pow(-0, -y) + } + switch { + case y < 0: + return 0 + case y > 0: + return Inf(1) + } + } + + absy := y + flip := false + if absy < 0 { + absy = -absy + flip = true + } + yi, yf := Modf(absy) + if yf != 0 && x < 0 { + return NaN() + } + if yi >= 1<<63 { + return Exp(y * Log(x)) + } + + // ans = a1 * 2**ae (= 1 for now). + a1 := 1.0 + ae := 0 + + // ans *= x**yf + if yf != 0 { + if yf > 0.5 { + yf-- + yi++ + } + a1 = Exp(yf * Log(x)) + } + + // ans *= x**yi + // by multiplying in successive squarings + // of x according to bits of yi. + // accumulate powers of two into exp. + x1, xe := Frexp(x) + for i := int64(yi); i != 0; i >>= 1 { + if i&1 == 1 { + a1 *= x1 + ae += xe + } + x1 *= x1 + xe <<= 1 + if x1 < .5 { + x1 += x1 + xe-- + } + } + + // ans = a1*2**ae + // if flip { ans = 1 / ans } + // but in the opposite order + if flip { + a1 = 1 / a1 + ae = -ae + } + return Ldexp(a1, ae) +} diff --git a/libgo/go/math/pow10.go b/libgo/go/math/pow10.go new file mode 100644 index 000000000..bda2e824e --- /dev/null +++ b/libgo/go/math/pow10.go @@ -0,0 +1,30 @@ +// Copyright 2009 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 + +// This table might overflow 127-bit exponent representations. +// In that case, truncate it after 1.0e38. +var pow10tab [70]float64 + +// Pow10 returns 10**e, the base-10 exponential of e. +func Pow10(e int) float64 { + if e < 0 { + return 1 / Pow10(-e) + } + if e < len(pow10tab) { + return pow10tab[e] + } + m := e / 2 + return Pow10(m) * Pow10(e-m) +} + +func init() { + pow10tab[0] = 1.0e0 + pow10tab[1] = 1.0e1 + for i := 2; i < len(pow10tab); i++ { + m := i / 2 + pow10tab[i] = pow10tab[m] * pow10tab[i-m] + } +} diff --git a/libgo/go/math/remainder.go b/libgo/go/math/remainder.go new file mode 100644 index 000000000..be8724c7f --- /dev/null +++ b/libgo/go/math/remainder.go @@ -0,0 +1,85 @@ +// 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 and the the comment below are from +// FreeBSD's /usr/src/lib/msun/src/e_remainder.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_remainder(x,y) +// Return : +// returns x REM y = x - [x/y]*y as if in infinite +// precision arithmetic, where [x/y] is the (infinite bit) +// integer nearest x/y (in half way cases, choose the even one). +// Method : +// Based on fmod() returning x - [x/y]chopped * y exactly. + +// Remainder returns the IEEE 754 floating-point remainder of x/y. +// +// Special cases are: +// Remainder(x, NaN) = NaN +// Remainder(NaN, y) = NaN +// Remainder(Inf, y) = NaN +// Remainder(x, 0) = NaN +// Remainder(x, Inf) = x +func Remainder(x, y float64) float64 { + const ( + Tiny = 4.45014771701440276618e-308 // 0x0020000000000000 + HalfMax = MaxFloat64 / 2 + ) + // TODO(rsc): Remove manual inlining of IsNaN, IsInf + // when compiler does it for us + // special cases + switch { + case x != x || y != y || x < -MaxFloat64 || x > MaxFloat64 || y == 0: // IsNaN(x) || IsNaN(y) || IsInf(x, 0) || y == 0: + return NaN() + case y < -MaxFloat64 || y > MaxFloat64: // IsInf(y): + return x + } + sign := false + if x < 0 { + x = -x + sign = true + } + if y < 0 { + y = -y + } + if x == y { + return 0 + } + if y <= HalfMax { + x = Fmod(x, y+y) // now x < 2y + } + if y < Tiny { + if x+x > y { + x -= y + if x+x >= y { + x -= y + } + } + } else { + yHalf := 0.5 * y + if x > yHalf { + x -= y + if x >= yHalf { + x -= y + } + } + } + if sign { + x = -x + } + return x +} diff --git a/libgo/go/math/signbit.go b/libgo/go/math/signbit.go new file mode 100644 index 000000000..670cc1a66 --- /dev/null +++ b/libgo/go/math/signbit.go @@ -0,0 +1,10 @@ +// 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 + +// Signbit returns true if x is negative or negative zero. +func Signbit(x float64) bool { + return Float64bits(x)&(1<<63) != 0 +} diff --git a/libgo/go/math/sin.go b/libgo/go/math/sin.go new file mode 100644 index 000000000..35220cb3e --- /dev/null +++ b/libgo/go/math/sin.go @@ -0,0 +1,66 @@ +// Copyright 2009 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 sine and cosine. + + Coefficients are #5077 from Hart & Cheney. (18.80D) +*/ + +func sinus(x float64, quad int) float64 { + const ( + P0 = .1357884097877375669092680e8 + P1 = -.4942908100902844161158627e7 + P2 = .4401030535375266501944918e6 + P3 = -.1384727249982452873054457e5 + P4 = .1459688406665768722226959e3 + Q0 = .8644558652922534429915149e7 + Q1 = .4081792252343299749395779e6 + Q2 = .9463096101538208180571257e4 + Q3 = .1326534908786136358911494e3 + ) + if x < 0 { + x = -x + quad = quad + 2 + } + x = x * (2 / Pi) /* underflow? */ + var y float64 + if x > 32764 { + var e float64 + e, y = Modf(x) + e = e + float64(quad) + f, _ := Modf(0.25 * e) + quad = int(e - 4*f) + } else { + k := int32(x) + y = x - float64(k) + quad = (quad + int(k)) & 3 + } + + if quad&1 != 0 { + y = 1 - y + } + if quad > 1 { + y = -y + } + + yy := y * y + temp1 := ((((P4*yy+P3)*yy+P2)*yy+P1)*yy + P0) * y + temp2 := ((((yy+Q3)*yy+Q2)*yy+Q1)*yy + Q0) + return temp1 / temp2 +} + +// Cos returns the cosine of x. +func Cos(x float64) float64 { + if x < 0 { + x = -x + } + return sinus(x, 1) +} + +// Sin returns the sine of x. +func Sin(x float64) float64 { return sinus(x, 0) } diff --git a/libgo/go/math/sincos.go b/libgo/go/math/sincos.go new file mode 100644 index 000000000..4c1576bea --- /dev/null +++ b/libgo/go/math/sincos.go @@ -0,0 +1,13 @@ +// 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 + +// Sincos(x) returns Sin(x), Cos(x). +// +// Special conditions are: +// Sincos(+Inf) = NaN, NaN +// Sincos(-Inf) = NaN, NaN +// Sincos(NaN) = NaN, NaN +func Sincos(x float64) (sin, cos float64) { return Sin(x), Cos(x) } diff --git a/libgo/go/math/sinh.go b/libgo/go/math/sinh.go new file mode 100644 index 000000000..23a8719f2 --- /dev/null +++ b/libgo/go/math/sinh.go @@ -0,0 +1,68 @@ +// Copyright 2009 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 hyperbolic sine and cosine. + + The exponential func is called for arguments + greater in magnitude than 0.5. + + A series is used for arguments smaller in magnitude than 0.5. + + Cosh(x) is computed from the exponential func for + all arguments. +*/ + +// Sinh returns the hyperbolic sine of x. +func Sinh(x float64) float64 { + // The coefficients are #2029 from Hart & Cheney. (20.36D) + const ( + P0 = -0.6307673640497716991184787251e+6 + P1 = -0.8991272022039509355398013511e+5 + P2 = -0.2894211355989563807284660366e+4 + P3 = -0.2630563213397497062819489e+2 + Q0 = -0.6307673640497716991212077277e+6 + Q1 = 0.1521517378790019070696485176e+5 + Q2 = -0.173678953558233699533450911e+3 + ) + + sign := false + if x < 0 { + x = -x + sign = true + } + + var temp float64 + switch true { + case x > 21: + temp = Exp(x) / 2 + + case x > 0.5: + temp = (Exp(x) - Exp(-x)) / 2 + + default: + sq := x * x + temp = (((P3*sq+P2)*sq+P1)*sq + P0) * x + temp = temp / (((sq+Q2)*sq+Q1)*sq + Q0) + } + + if sign { + temp = -temp + } + return temp +} + +// Cosh returns the hyperbolic cosine of x. +func Cosh(x float64) float64 { + if x < 0 { + x = -x + } + if x > 21 { + return Exp(x) / 2 + } + return (Exp(x) + Exp(-x)) / 2 +} diff --git a/libgo/go/math/sqrt.go b/libgo/go/math/sqrt.go new file mode 100644 index 000000000..bf6ef64a0 --- /dev/null +++ b/libgo/go/math/sqrt.go @@ -0,0 +1,28 @@ +// Copyright 2009 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 libc_sqrt(float64) float64 __asm__("sqrt") + +// Sqrt returns the square root of x. +// +// Special cases are: +// Sqrt(+Inf) = +Inf +// Sqrt(±0) = ±0 +// Sqrt(x < 0) = NaN +// Sqrt(NaN) = NaN +func Sqrt(x float64) float64 { + // special cases + // TODO(rsc): Remove manual inlining of IsNaN, IsInf + // when compiler does it for us + switch { + case x == 0 || x != x || x > MaxFloat64: // x == 0 || IsNaN(x) || IsInf(x, 1): + return x + case x < 0: + return NaN() + } + + return libc_sqrt(x) +} diff --git a/libgo/go/math/sqrt_decl.go b/libgo/go/math/sqrt_decl.go new file mode 100644 index 000000000..e50774645 --- /dev/null +++ b/libgo/go/math/sqrt_decl.go @@ -0,0 +1,7 @@ +// Copyright 2009 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 Sqrt(x float64) float64 diff --git a/libgo/go/math/sqrt_port.go b/libgo/go/math/sqrt_port.go new file mode 100644 index 000000000..6f35a383d --- /dev/null +++ b/libgo/go/math/sqrt_port.go @@ -0,0 +1,143 @@ +// Copyright 2009 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 square root. +*/ + +// The original C code and the long comment below are +// from FreeBSD's /usr/src/lib/msun/src/e_sqrt.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_sqrt(x) +// Return correctly rounded sqrt. +// ----------------------------------------- +// | Use the hardware sqrt if you have one | +// ----------------------------------------- +// Method: +// Bit by bit method using integer arithmetic. (Slow, but portable) +// 1. Normalization +// Scale x to y in [1,4) with even powers of 2: +// find an integer k such that 1 <= (y=x*2**(2k)) < 4, then +// sqrt(x) = 2**k * sqrt(y) +// 2. Bit by bit computation +// Let q = sqrt(y) truncated to i bit after binary point (q = 1), +// i 0 +// i+1 2 +// s = 2*q , and y = 2 * ( y - q ). (1) +// i i i i +// +// To compute q from q , one checks whether +// i+1 i +// +// -(i+1) 2 +// (q + 2 ) <= y. (2) +// i +// -(i+1) +// If (2) is false, then q = q ; otherwise q = q + 2 . +// i+1 i i+1 i +// +// With some algebric manipulation, it is not difficult to see +// that (2) is equivalent to +// -(i+1) +// s + 2 <= y (3) +// i i +// +// The advantage of (3) is that s and y can be computed by +// i i +// the following recurrence formula: +// if (3) is false +// +// s = s , y = y ; (4) +// i+1 i i+1 i +// +// otherwise, +// -i -(i+1) +// s = s + 2 , y = y - s - 2 (5) +// i+1 i i+1 i i +// +// One may easily use induction to prove (4) and (5). +// Note. Since the left hand side of (3) contain only i+2 bits, +// it does not necessary to do a full (53-bit) comparison +// in (3). +// 3. Final rounding +// After generating the 53 bits result, we compute one more bit. +// Together with the remainder, we can decide whether the +// result is exact, bigger than 1/2ulp, or less than 1/2ulp +// (it will never equal to 1/2ulp). +// The rounding mode can be detected by checking whether +// huge + tiny is equal to huge, and whether huge - tiny is +// equal to huge for some floating point number "huge" and "tiny". +// +// +// Notes: Rounding mode detection omitted. The constants "mask", "shift", +// and "bias" are found in src/pkg/math/bits.go + +// Sqrt returns the square root of x. +// +// Special cases are: +// Sqrt(+Inf) = +Inf +// Sqrt(±0) = ±0 +// Sqrt(x < 0) = NaN +// Sqrt(NaN) = NaN +func sqrtGo(x float64) float64 { + // special cases + // TODO(rsc): Remove manual inlining of IsNaN, IsInf + // when compiler does it for us + switch { + case x == 0 || x != x || x > MaxFloat64: // x == 0 || IsNaN(x) || IsInf(x, 1): + return x + case x < 0: + return NaN() + } + ix := Float64bits(x) + // normalize x + exp := int((ix >> shift) & mask) + if exp == 0 { // subnormal x + for ix&1<>= 1 // exp = exp/2, exponent of square root + // generate sqrt(x) bit by bit + ix <<= 1 + var q, s uint64 // q = sqrt(x) + r := uint64(1 << (shift + 1)) // r = moving bit from MSB to LSB + for r != 0 { + t := s + r + if t <= ix { + s = t + r + ix -= t + q += r + } + ix <<= 1 + r >>= 1 + } + // final rounding + if ix != 0 { // remainder, result not exact + q += q & 1 // round according to extra bit + } + ix = q>>1 + uint64(exp-1+bias)< 21 { + return -1 + } + return -Sinh(x) / Cosh(x) + } + if x > 21 { + return 1 + } + return Sinh(x) / Cosh(x) +} diff --git a/libgo/go/math/unsafe.go b/libgo/go/math/unsafe.go new file mode 100644 index 000000000..5ae67420f --- /dev/null +++ b/libgo/go/math/unsafe.go @@ -0,0 +1,21 @@ +// Copyright 2009 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 + +import "unsafe" + +// Float32bits returns the IEEE 754 binary representation of f. +func Float32bits(f float32) uint32 { return *(*uint32)(unsafe.Pointer(&f)) } + +// Float32frombits returns the floating point number corresponding +// to the IEEE 754 binary representation b. +func Float32frombits(b uint32) float32 { return *(*float32)(unsafe.Pointer(&b)) } + +// Float64bits returns the IEEE 754 binary representation of f. +func Float64bits(f float64) uint64 { return *(*uint64)(unsafe.Pointer(&f)) } + +// Float64frombits returns the floating point number corresponding +// the IEEE 754 binary representation b. +func Float64frombits(b uint64) float64 { return *(*float64)(unsafe.Pointer(&b)) } -- cgit v1.2.3