Odin Documentation

core/math

math

Types

1

Float_Class

Float_Class :: enum int { Normal = 0, // an ordinary nonzero floating point value Subnormal = 1, // a subnormal floating point value Zero = 2, // zero Neg_Zero = 3, // the negative zero NaN = 4, // Not-A-Number (NaN) Inf = 5, // positive infinity Neg_Inf = 6, // negative infinity }Source

Constants

79

PI

PI :: 3.14159265358979323846264338327950288Source

Procedures

468

_trig_reduce_f64

_trig_reduce_f64 :: proc(x: f64) -> (j: u64, z: f64)Source

_trig_reduce_f64 implements Payne-Hanek range reduction by Pi/4 for x > 0. It returns the integer part mod 8 (j) and the fractional part (z) of x / (Pi/4). The implementation is based on: "ARGUMENT REDUCTION FOR HUGE ARGUMENTS: Good to the Last Bit" K. C. Ng et al, March 24, 1992 The simulated multi-precision calculation of x*B uses 64-bit integer arithmetic.

atan

atan :: proc(x: T) -> (T)Source

Return the arc tangent of x, in radians. Defined on the domain of [-∞, ∞] with a range of [-π/2, π/2]

is_inf_f16

is_inf_f16 :: proc(x: f16, sign: int) -> (bool)Source

is_inf reports whether f is an infinity, according to sign. If sign > 0, is_inf reports whether f is positive infinity. If sign < 0, is_inf reports whether f is negative infinity. If sign == 0, is_inf reports whether f is either infinity.

remap

remap :: proc(old_value: T, old_min: T, old_max: T, new_min: T, new_max: T) -> (x: T)Source

Procedure Groups

44

atan2

atan2 :: proc{atan2_f64, atan2_f32, atan2_f16, atan2_f64le, atan2_f64be, atan2_f32le, atan2_f32be, atan2_f16le, atan2_f16be}Source

Return the arc tangent of y/x in radians. Defined on the domain [-∞, ∞] for x and y with a range of [-π, π]

Special cases:

atan2(y, NaN)     = NaN
	atan2(NaN, x)     = NaN
	atan2(+0, x>=0)   = + 0
	atan2(-0, x>=0)   = - 0
	atan2(+0, x<=-0)  = + π
	atan2(-0, x<=-0)  = - π
	atan2(y>0, 0)     = + π/2
	atan2(y<0, 0)     = - π/2
	atan2(+∞, +∞)     = + π/4
	atan2(-∞, +∞)     = - π/4
	atan2(+∞, -∞)     =   3π/4
	atan2(-∞, -∞)     = - 3π/4
	atan2(y, +∞)      =   0
	atan2(y>0, -∞)    = + π
	atan2(y<0, -∞)    = - π
	atan2(+∞, x)      = + π/2
	atan2(-∞, x)      = - π/2

erf

erf :: proc{erf_f16, erf_f16le, erf_f16be, erf_f32, erf_f32le, erf_f32be, erf_f64}Source

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.

==================================================== 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 + xR(x2) erfc(x) = 1 - erf(x) if x in [-.84375,0.25] = 0.5 + ((0.5-x)-xR) 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 - x3/3 + x5/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) + sPoly(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(-xx-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(-xx-0.5625+R2/S2) if x > 0 = 2.0 - (1/x)exp(-xx-0.5625+R2/S2) if -6<x<0 = 2.0 - tiny (if x <= -6) erf(x) = sign(x)(1.0 - erfc(x)) if x < 6, else erf(x) = sign(x)(1.0 - tiny) where R2(z) = degree 6 poly in z, (z=1/x**2) S2(z) = degree 7 poly in z

Note1: To compute exp(-xx-0.5625+R/S), let s be a single precision number and s := x; then -xx = -ss + (s-x)(s+x) exp(-xx-0.5626+R/S) = exp(-ss-0.5625)exp((s-x)(s+x)+R/S); Note2: Here 4 and 5 make use of the asymptotic series exp(-xx) erfc(x) ~ ---------- ( 1 + Poly(1/x2) ) x*sqrt(pi) We use rational approximation to approximate g(s)=f(1/x2) = log(erfc(x)x) - xx + 0.5625 Here is the error bound for R1/S1 and R2/S2 |R1/S1 - f(x)| < 2(-62.57) |R2/S2 - f(x)| < 2(-61.52)

5. For inf > x >= 28 erf(x) = sign(x) (1 - tiny) (raise inexact) erfc(x) = tinytiny (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

log1p

log1p :: proc{log1p_f16, log1p_f32, log1p_f64, log1p_f16le, log1p_f16be, log1p_f32le, log1p_f32be, log1p_f64le, log1p_f64be}Source

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 s3 + 2/5 s5 + ....., = 2s + sR 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) ~ Lp1s +Lp2s +Lp3s +Lp4s +Lp5s +Lp6s +Lp7s (the values of Lp1 to Lp7 are listed in the program) and | 2 14 | -58.45 | Lp1s +...+Lp7s - R(z) | <= 2

Note that 2s = f - sf = f - hfsq + shfsq, where hfsq = ff/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) = kln2 + log1p(f). = kln2_hi+(f-(hfsq-(s(hfsq+R)+kln2_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.

sincos

sincos :: proc{sincos_f16, sincos_f16le, sincos_f16be, sincos_f32, sincos_f32le, sincos_f32be, sincos_f64, sincos_f64le, sincos_f64be}Source

The original C code, the long comment, and the constants below were from http://netlib.sandia.gov/cephes/cmath/sin.c, available from http://www.netlib.org/cephes/cmath.tgz. The go code is a simplified version of the original C.

sin.c

Circular sine

SYNOPSIS:

double x, y, sin(); y = sin( x );

DESCRIPTION:

Range reduction is into intervals of pi/4. The reduction error is nearly eliminated by contriving an extended precision modular arithmetic.

Two polynomial approximating functions are employed. Between 0 and pi/4 the sine is approximated by x + x3 P(x2). Between pi/4 and pi/2 the cosine is represented as 1 - x2 Q(x2).

ACCURACY:

Relative error: arithmetic domain # trials peak rms DEC 0, 10 150000 3.0e-17 7.8e-18 IEEE -1.07e9,+1.07e9 130000 2.1e-16 5.4e-17

Partial loss of accuracy begins to occur at x = 230 = 1.074e9. The loss is not gradual, but jumps suddenly to about 1 part in 10e7. Results may be meaningless for x > 249 = 5.6e14.

cos.c

Circular cosine

SYNOPSIS:

double x, y, cos(); y = cos( x );

DESCRIPTION:

Range reduction is into intervals of pi/4. The reduction error is nearly eliminated by contriving an extended precision modular arithmetic.

Two polynomial approximating functions are employed. Between 0 and pi/4 the cosine is approximated by 1 - x2 Q(x2). Between pi/4 and pi/2 the sine is represented as x + x3 P(x2).

ACCURACY:

Relative error: arithmetic domain # trials peak rms IEEE -1.07e9,+1.07e9 130000 2.1e-16 5.4e-17 DEC 0,+1.07e9 17000 3.0e-17 7.2e-18

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

Reference search

Find anything

Documentation preferences

Settings

System theme variants

Used only while Theme is set to System.