-
-
Save LingDong-/7e4c4cae5cbbc44400a05fba65f06f23 to your computer and use it in GitHub Desktop.
simple, fast, accurate natural log approximation without <math.h>
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
// ln.c | |
// | |
// simple, fast, accurate natural log approximation | |
// when without <math.h> | |
// featuring * floating point bit level hacking, | |
// * x=m*2^p => ln(x)=ln(m)+ln(2)p, | |
// * Remez algorithm | |
// by Lingdong Huang, 2020. Public domain. | |
// ============================================ | |
float ln(float x) { | |
unsigned int bx = * (unsigned int *) (&x); | |
unsigned int ex = bx >> 23; | |
signed int t = (signed int)ex-(signed int)127; | |
unsigned int s = (t < 0) ? (-t) : t; | |
bx = 1065353216 | (bx & 8388607); | |
x = * (float *) (&bx); | |
return -1.49278+(2.11263+(-0.729104+0.10969*x)*x)*x+0.6931471806*t; | |
} | |
// done. | |
// ============================================ | |
// Exact same function, with added comments: | |
float natural_log(float x) { | |
// ASSUMING: | |
// - non-denormalized numbers i.e. x > 2^−126 | |
// - integer is 32 bit. float is IEEE 32 bit. | |
// INSPIRED BY: | |
// - https://stackoverflow.com/a/44232045 | |
// - http://mathonweb.com/help_ebook/html/algorithms.htm#ln | |
// - https://en.wikipedia.org/wiki/Fast_inverse_square_root | |
// FORMULA: | |
// x = m * 2^p => | |
// ln(x) = ln(m) + ln(2)p, | |
// first normalize the value to between 1.0 and 2.0 | |
// assuming normalized IEEE float | |
// sign exp frac | |
// 0b 0 [00000000] 00000000000000000000000 | |
// value = (-1)^s * M * 2 ^ (exp-127) | |
// | |
// exp = 127 for x = 1, | |
// so 2^(exp-127) is the multiplier | |
// evil floating point bit level hacking | |
unsigned int bx = * (unsigned int *) (&x); | |
// extract exp, since x>0, sign bit must be 0 | |
unsigned int ex = bx >> 23; | |
signed int t = (signed int)ex-(signed int)127; | |
unsigned int s = (t < 0) ? (-t) : t; | |
// reinterpret back to float | |
// 127 << 23 = 1065353216 | |
// 0b11111111111111111111111 = 8388607 | |
bx = 1065353216 | (bx & 8388607); | |
x = * (float *) (&bx); | |
// use remez algorithm to find approximation between [1,2] | |
// - see this answer https://stackoverflow.com/a/44232045 | |
// - or this usage of C++/boost's remez implementation | |
// https://computingandrecording.wordpress.com/2017/04/24/ | |
// e.g. | |
// boost::math::tools::remez_minimax<double> approx( | |
// [](const double& x) { return log(x); }, | |
// 4, 0, 1, 2, false, 0, 0, 64); | |
// | |
// 4th order is: | |
// { -1.74178, 2.82117, -1.46994, 0.447178, -0.0565717 } | |
// | |
// 3rd order is: | |
// { -1.49278, 2.11263, -0.729104, 0.10969 } | |
return | |
/* less accurate */ | |
-1.49278+(2.11263+(-0.729104+0.10969*x)*x)*x | |
/* OR more accurate */ | |
// -1.7417939+(2.8212026+(-1.4699568+(0.44717955-0.056570851*x)*x)*x)*x | |
/* compensate for the ln(2)s. ln(2)=0.6931471806 */ | |
+ 0.6931471806*t; | |
} |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
Dude youre a genuius. I need an exp function as well