Add `math.signbit`; Add `math.gamma` based on http://netlib.sandia.gov/cephes/cprob/gamma.c

core/math/math.odin

@@ -1642,6 +1642,28 @@ nextafter :: proc{
 	nextafter_f64, nextafter_f64le, nextafter_f64be,
 }
 
+signbit_f16 :: proc "contextless" (x: f16) -> bool {
+	return (transmute(u16)x)&(1<<15) != 0
+}
+signbit_f32 :: proc "contextless" (x: f32) -> bool {
+	return (transmute(u32)x)&(1<<31) != 0
+}
+signbit_f64 :: proc "contextless" (x: f64) -> bool {
+	return (transmute(u64)x)&(1<<63) != 0
+}
+signbit_f16le :: proc "contextless" (x: f16le) -> bool { return signbit_f16(f16(x)) }
+signbit_f32le :: proc "contextless" (x: f32le) -> bool { return signbit_f32(f32(x)) }
+signbit_f64le :: proc "contextless" (x: f64le) -> bool { return signbit_f64(f64(x)) }
+signbit_f16be :: proc "contextless" (x: f16be) -> bool { return signbit_f16(f16(x)) }
+signbit_f32be :: proc "contextless" (x: f32be) -> bool { return signbit_f32(f32(x)) }
+signbit_f64be :: proc "contextless" (x: f64be) -> bool { return signbit_f64(f64(x)) }
+
+signbit :: proc{
+	signbit_f16, signbit_f16le, signbit_f16be,
+	signbit_f32, signbit_f32le, signbit_f32be,
+	signbit_f64, signbit_f64le, signbit_f64be,
+}
+
 
 F16_DIG        :: 3
 F16_EPSILON    :: 0.00097656

core/math/math_gamma.odin

@@ -0,0 +1,226 @@
+package math
+
+// The original C code, the long comment, and the constants
+// below are from http://netlib.sandia.gov/cephes/cprob/gamma.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
+//   [email protected]
+
+// Gamma function computed by Stirling's formula.
+// The pair of results must be multiplied together to get the actual answer.
+// The multiplication is left to the caller so that, if careful, the caller can avoid
+// infinity for 172 <= x <= 180.
+// The polynomial is valid for 33 <= x <= 172; larger values are only used
+// in reciprocal and produce denormalized floats. The lower precision there
+// masks any imprecision in the polynomial.
+@(private="file")
+stirling :: proc "contextless" (x: f64) -> (f64, f64) {
+	@(static) gamS := [?]f64{
+		7.87311395793093628397e-04,
+		-2.29549961613378126380e-04,
+		-2.68132617805781232825e-03,
+		3.47222221605458667310e-03,
+		8.33333333333482257126e-02,
+	}
+	
+	if x > 200 {
+		return inf_f64(1), 1
+	}
+	SQRT_TWO_PI :: 2.506628274631000502417
+	MAX_STIRLING :: 143.01608
+	w := 1 / x
+	w = 1 + w*((((gamS[0]*w+gamS[1])*w+gamS[2])*w+gamS[3])*w+gamS[4])
+	y1 := exp(x)
+	y2 := 1.0
+	if x > MAX_STIRLING { // avoid pow() overflow
+		v := pow(x, 0.5*x-0.25)
+		y1, y2 = v, v/y1
+	} else {
+		y1 = pow(x, x-0.5) / y1
+	}
+	return y1, SQRT_TWO_PI * w * y2
+}
+
+gamma_f64 :: proc "contextless" (x: f64) -> f64 {
+	is_neg_int :: proc "contextless" (x: f64) -> bool {
+		if x < 0 {
+			_, xf := modf(x)
+			return xf == 0
+		}
+		return false
+	}
+	
+	@(static) gamP := [?]f64{
+		1.60119522476751861407e-04,
+		1.19135147006586384913e-03,
+		1.04213797561761569935e-02,
+		4.76367800457137231464e-02,
+		2.07448227648435975150e-01,
+		4.94214826801497100753e-01,
+		9.99999999999999996796e-01,
+	}
+	@(static) gamQ := [?]f64{
+		-2.31581873324120129819e-05,
+		5.39605580493303397842e-04,
+		-4.45641913851797240494e-03,
+		1.18139785222060435552e-02,
+		3.58236398605498653373e-02,
+		-2.34591795718243348568e-01,
+		7.14304917030273074085e-02,
+		1.00000000000000000320e+00,
+	}
+
+	
+	EULER :: 0.57721566490153286060651209008240243104215933593992 // A001620
+	
+	switch {
+	case is_neg_int(x) || is_inf(x, -1) || is_nan(x):
+		return nan_f64()
+	case is_inf(x, 1):
+		return inf_f64(1)
+	case x == 0:
+		if signbit(x) {
+			return inf_f64(-1)
+		}
+		return inf_f64(1)
+	}
+	
+	x := x
+	q := abs(x)
+	p := floor(q)
+	if q > 33 {
+		if x >= 0 {
+			y1, y2 := stirling(x)
+			return y1 * y2
+		}
+		// Note: x is negative but (checked above) not a negative integer,
+		// so x must be small enough to be in range for conversion to i64.
+		// If |x| were >= 2⁶³ it would have to be an integer.
+		signgam := 1
+		if ip := i64(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_f64(signgam)
+		}
+		sq1, sq2 := stirling(q)
+		absz := abs(z)
+		d := absz * sq1 * sq2
+		if is_inf(d, 0) {
+			z = PI / absz / sq1 / sq2
+		} else {
+			z = PI / d
+		}
+		return f64(signgam) * z
+	}
+
+	// Reduce argument
+	z := 1.0
+	for x >= 3 {
+		x = x - 1
+		z = z * x
+	}
+	for x < 0 {
+		if x > -1e-09 {
+			if x == 0 {
+				return inf_f64(1)
+			}
+			return z / ((1 + EULER*x) * x)
+		}
+		z = z / x
+		x = x + 1
+	}
+	for x < 2 {
+		if x < 1e-09 {
+			if x == 0 {
+				return inf_f64(1)
+			}
+			return z / ((1 + EULER*x) * x)
+		}
+		z = z / x
+		x = x + 1
+	}
+
+	if x == 2 {
+		return z
+	}
+
+	x = x - 2
+	p = (((((x*gamP[0]+gamP[1])*x+gamP[2])*x+gamP[3])*x+gamP[4])*x+gamP[5])*x + gamP[6]
+	q = ((((((x*gamQ[0]+gamQ[1])*x+gamQ[2])*x+gamQ[3])*x+gamQ[4])*x+gamQ[5])*x+gamQ[6])*x + gamQ[7]
+	return z * p / q
+}
+
+
+gamma_f16   :: proc "contextless" (x: f16)   -> f16   { return f16(gamma_f64(f64(x))) }
+gamma_f16le :: proc "contextless" (x: f16le) -> f16le { return f16le(gamma_f64(f64(x))) }
+gamma_f16be :: proc "contextless" (x: f16be) -> f16be { return f16be(gamma_f64(f64(x))) }
+gamma_f32   :: proc "contextless" (x: f32)   -> f32   { return f32(gamma_f64(f64(x))) }
+gamma_f32le :: proc "contextless" (x: f32le) -> f32le { return f32le(gamma_f64(f64(x))) }
+gamma_f32be :: proc "contextless" (x: f32be) -> f32be { return f32be(gamma_f64(f64(x))) }
+gamma_f64le :: proc "contextless" (x: f64le) -> f64le { return f64le(gamma_f64(f64(x))) }
+gamma_f64be :: proc "contextless" (x: f64be) -> f64be { return f64be(gamma_f64(f64(x))) }
+
+gamma :: proc{
+	gamma_f16, gamma_f16le, gamma_f16be,
+	gamma_f32, gamma_f32le, gamma_f32be,
+	gamma_f64, gamma_f64le, gamma_f64be,
+}