C/C++: Euler's Math Library

Found use for these functions when writing code for Nvidia CUDA, and not wanting to or not being able to import libraries.

// eulers.h

#ifndef EULERS_H
#define EULERS_H

#ifndef INFINITY
#define INFINITY   ((float)(1e+300 * 1e+300))
#endif

#ifndef NAN
#define NAN        (-(float)(INFINITY * 0.0F))
#endif

// Is INFINITY
int Isinf(double x) {
	union {
		double d;
		unsigned long long i;
	} u = { x };

	// Extract exponent and mantissa
	unsigned long long exp = (u.i >> 52) & 0x7FF;
	unsigned long long mantissa = u.i & ((1ULL << 52) - 1);

	// Check if exponent is all 1's and mantissa is 0
	if (exp == 0x7FF && mantissa == 0) {
		return (u.i >> 63) ? -1 : 1;  // Check sign bit
	}
	return 0;
}

// Is Not a Number
// Returns 1 if x is NaN, 0 otherwise
int Isnan(double x) {
	union {
		double d;
		unsigned long long i;
	} u = { x };

	// Extract exponent and mantissa
	unsigned long long exp = (u.i >> 52) & 0x7FF;
	unsigned long long mantissa = u.i & ((1ULL << 52) - 1);

	// NaN has all 1's in exponent and non-zero mantissa
	return (exp == 0x7FF && mantissa != 0);
}

// Ceiling
double Ceil(double x) {
	// Handle special cases
	if (Isinf(x) || Isnan(x)) return x;

	long long i = (long long)x;
	return (x > 0.0 && x > i) ? i + 1.0 : (double)i;
}

// Returns absolute value of x
double Fabs(double x) {
	// Use bit manipulation to clear sign bit
	union {
		double d;
		unsigned long long i;
	} u = { x };
	u.i &= ~(1ULL << 63);  // Clear sign bit
	return u.d;
}

// Implementation of sqrt(x) using Newton's method
// Valid for x >= 0
double Sqrt(double x) {
	// Handle special cases
	if (x < 0.0) return NAN;
	if (x == 0.0 || x == 1.0 || Isinf(x)) return x;

	// Initial guess: use bit manipulation for a good starting point
	union {
		double d;
		unsigned long long i;
	} u = { x };
	u.i = 0x5fe6eb50c7b537a9 - (u.i >> 1);
	double guess = u.d;

	// Newton's method: repeat until convergence
	for (int i = 0; i < 5; i++) {
		guess = (guess + x / guess) * 0.5;
	}

	return guess;
}

// Implementation of exp(x) using Taylor series
// Valid for any real number x
double Exp(double x) {
	// Handle special cases
	if (x == 0.0) return 1.0;
	if (Isinf(x)) return x > 0 ? INFINITY : 0.0;

	// Use the property exp(x) = (exp(x/n))^n to reduce large values
	int n = 1;
	if (Fabs(x) > 1.0) {
		n = (int)Ceil(Fabs(x));
		x /= n;
	}

	double sum = 1.0;  // First term
	double term = 1.0; // Current term
	double i = 1.0;    // Current iteration

	// Add terms until desired precision reached
	while (Fabs(term) > 1e-15 && i < 100) {
		term *= x / i;
		sum += term;
		i += 1.0;
	}

	// If dividing x by n earlier, need to take result to nth power
	double result = sum;
	for (int j = 1; j < n; j++) {
		result *= sum;
	}

	return result;
}

// Implementation of log(x) using series expansion
// Valid for x > 0
double Log(double x) {
	// Handle special cases
	if (x <= 0.0) return NAN;
	if (x == 1.0) return 0.0;
	if (Isinf(x)) return INFINITY;

	// Use log(x) = 2 * log(sqrt(x))
	// Repeatedly take sqrt until x is close to 1
	int power = 0;
	while (x >= 2.0) {
		x = Sqrt(x);
		power++;
	}
	while (x < 0.5) {
		x = Sqrt(x);
		power--;
	}

	// Now use series log(1+y) = y - y^2/2 + y^3/3 - ...
	// where y = (x-1)/(x+1)
	double y = (x - 1.0) / (x + 1.0);
	double y2 = y * y;
	double sum = y;
	double term = y;
	int i = 1;

	// Add terms until desired precision
	while (Fabs(term) > 1e-15 && i < 100) {
		term *= y2 * (2 * i - 1) / (2 * i + 1);
		sum += term;
		i++;
	}

	// Adjust for the sqrt operations that were done earlier
	return (sum * 2.0) * (1 << power);
}

#endif

Revised to be more IEEE compliant. This may work better on certain systems.

// eulers2.h

#ifndef EULERS_H
#define EULERS_H

// Define IEEE 754 double-precision bits for Infinity and NaN as constants:
static const union { unsigned long long __u; double __d; } __infinity_union = { 0x7FF0000000000000ULL };
static const union { unsigned long long __u; double __d; } __nan_union = { 0x7FF8000000000001ULL };

// Then define macros if they are not already defined (shouldn't be):
#ifndef INFINITY
#define INFINITY (__infinity_union.__d)
#endif

#ifndef NAN
#define NAN (__nan_union.__d)
#endif

// Check for infinity: returns +1 for +inf, -1 for -inf, 0 otherwise
int Isinf(double x) {
	union { double d; unsigned long long u; } u = { x };
	unsigned long long exp = (u.u >> 52) & 0x7FFULL;
	unsigned long long frac = u.u & 0x000FFFFFFFFFFFFFULL;

	if (exp == 0x7FF && frac == 0ULL) {
		return (u.u >> 63) ? -1 : 1;
	}
	return 0;
}

// Check for NaN: returns 1 if NaN, 0 otherwise
int Isnan(double x) {
	union { double d; unsigned long long u; } u = { x };
	unsigned long long exp = (u.u >> 52) & 0x7FFULL;
	unsigned long long frac = u.u & 0x000FFFFFFFFFFFFFULL;

	// NaN: exponent all 1's and non-zero fraction
	return (exp == 0x7FF && frac != 0ULL);
}

// Ceiling function: returns the smallest integer >= x
double Ceil(double x) {
	if (Isinf(x) || Isnan(x)) return x;

	long long i = (long long)x;
	// If x > 0 and x not an integer, ceil is i+1; else just i.
	return (x > 0.0 && x > (double)i) ? (double)(i + 1) : (double)i;
}

// Absolute value
double Fabs(double x) {
	union { double d; unsigned long long u; } u = { x };
	u.u &= ~(1ULL << 63);
	return u.d;
}

// Square root using Newton's method
double Sqrt(double x) {
	if (x < 0.0) return NAN;
	if (x == 0.0 || x == 1.0 || Isinf(x)) return x;

	// Initial guess: a clever bit trick (approximate)
	union { double d; unsigned long long i; } u = { x };
	// Initial guess: from Quake's fast inverse sqrt approximation
	u.i = 0x5fe6eb50c7b537a9ULL - (u.i >> 1);
	double guess = u.d;

	double old_guess;
	// Iterate until convergence
	do {
		old_guess = guess;
		guess = 0.5 * (guess + x / guess);
	} while (Fabs(guess - old_guess) > 1e-15 * guess);

	return guess;
}

// Exp(x): has range reduction
double Exp(double x) {
	// Handle special cases
	if (Isnan(x)) return NAN;
	if (Isinf(x)) return (x > 0.0) ? INFINITY : 0.0;
	if (x == 0.0) return 1.0;

	// Constants
	const double LN2 = 0.69314718055994530942; // Approximation of ln(2)

	// For large |x|, break down x = n*ln(2) + r, where r in [-ln(2)/2, ln(2)/2]
	// This gives us exp(x) = 2^n * exp(r), with |r| <= 0.3466 approx.
	int n = (int)(x / LN2);
	double r = x - n * LN2;

	// If x is very small (like negative large), then might have a large negative n.
	// For example, x = -100 might give a large negative n. This is still fine.

	// Now use a Taylor series expansion for exp(r) around 0:
	// exp(r) = 1 + r + r^2/2! + r^3/3! + ...
	double term = 1.0;
	double sum = 1.0;
	double i = 1.0;
	while (Fabs(term) > 1e-15 && i < 100) {
		term *= r / i;
		sum += term;
		i += 1.0;
	}

	// Now multiply by 2^n
	double result = sum;
	if (n > 0) {
		// multiply result by 2^n
		for (int j = 0; j < n; j++) {
			result *= 2.0;
			if (Isinf(result)) break; // avoid overflow
		}
	}
	else if (n < 0) {
		// divide result by 2^-n
		for (int j = 0; j < -n; j++) {
			result *= 0.5;
			if (result == 0.0) break; // underflow
		}
	}

	return result;
}

// Log(x): scaling into [0.5, 2) then using series
double Log(double x) {
	if (x < 0.0) return NAN;
	if (x == 0.0) return -INFINITY;
	if (Isinf(x)) return INFINITY;
	if (Isnan(x)) return NAN;
	if (x == 1.0) return 0.0;

	const double LN2 = 0.69314718055994530942;

	// Scale x into [0.5, 2) by adjusting an integer k such that:
	// x = 2^k * y, with y in [0.5, 2)
	int k = 0;
	while (x >= 2.0) {
		x *= 0.5;
		k++;
	}
	while (x < 0.5) {
		x *= 2.0;
		k--;
	}

	// Now use y = (x-1)/(x+1)
	// log(x) = 2 * [ y + y^3/3 + y^5/5 + ... ] for x in (0.5,2)
	double y = (x - 1.0) / (x + 1.0);
	double y2 = y * y;
	double term = y;
	double sum = term;
	int i = 1;
	// Add higher-order terms
	// The standard series: log(x) = 2*(y + y^3/3 + y^5/5 + ...)
	// Stop when terms get very small.
	do {
		i += 2;
		term *= y2;
		double add = term / i;
		if (Fabs(add) < 1e-15) break;
		sum += add;
	} while (i < 200);

	double result = 2.0 * sum + k * LN2;
	return result;
}

#endif

Now, use the above functions!

#include "eulers2.h"

double Tanh(double x) // Hyperbolic Tangent function
{
	double n = Exp(-x);
	double p = Exp(x);
	return (p - n) / (p + n);
}

double TanhPrime(double x) // derivative of Tanh
{
	double n = Exp(-x);
	double p = Exp(x);
	return 1.0 - ((p - n) / (p + n)) * ((p - n) / (p + n)); // this is simply: 1.0 - (tanh(x) * tanh(x))
}

double Sigmoid(double x) // and oldie but goodie, ReLU/ELU has replaced it for the most part
{
	return 1.0 / (1.0 + Exp(-x));
}

double SigmoidPrime(double x) // derivative of Sigmoid
{
	double n = Exp(-x);
	return (1.0 / (1.0 + n)) * (1.0 - (1.0 / (1.0 + n))); // this is simply: Sigmoid(x) * (1.0 - Sigmoid(x))
}
double ELU(double x, double alpha = 1.0) // Exponential Linear Unit function
{
	return x >= 0 ? x : (alpha * (Exp(x) - 1.0));
}

double ELUPrime(double x, double alpha = 1.0) // derivative of ELU; make sure to use same alpha value as passed to ELU()
{
	return x >= 0 ? 1.0 : (alpha * Exp(x));
}

float ConvertSoftmaxLossToAccuracy(float loss, int numClasses)
{
	// Calculate the maximum possible loss (random prediction)
	// For random prediction, probability is 1/numClasses for each class
	float maxLoss = -Log(1.0F / numClasses);

	// Calculate accuracy: 1 when loss is 0, 0 when loss equals maxLoss
	// Using a linear interpolation between these points
	float accuracy = 1.0F - (loss / maxLoss);

	// Clamp the result between 0 and 1 to handle edge cases
	accuracy = accuracy > 1.0F ? 1.0F : accuracy;
	accuracy = accuracy < 0.0F ? 0.0F : accuracy;
	return accuracy;
}

float ConvertBinaryLossToAccuracy(float loss)
{
	return Exp(-loss);
}

Comment