diff options
Diffstat (limited to '05')
| -rw-r--r-- | 05/main.c | 10 | ||||
| -rw-r--r-- | 05/math.h | 409 | ||||
| -rw-r--r-- | 05/stdc_common.h | 1 | ||||
| -rw-r--r-- | 05/tokenize.b | 7 |
4 files changed, 415 insertions, 12 deletions
@@ -1,15 +1,9 @@ #define _STDLIB_DEBUG +#include <math.h> #include <stdio.h> -#include <string.h> -#include <stdlib.h> -#include <stddef.h> -#include <ctype.h> -#include <locale.h> int main(int argc, char **argv) { - setlocale(LC_ALL, "C"); - struct lconv *c = localeconv(); - printf("%s\n",c->negative_sign); + printf("%.16g\n", cosh(10)); return 0; } diff --git a/05/math.h b/05/math.h new file mode 100644 index 0000000..7b029ac --- /dev/null +++ b/05/math.h @@ -0,0 +1,409 @@ +#ifndef _MATH_H +#define _MATH_H + +#include <stdc_common.h> +#define HUGE_VAL _INFINITY // glibc defines HUGE_VAL as infinity (the C standard only requires it to be positive, funnily enough) +#define _NAN (-(_INFINITY-_INFINITY)) +#define _PI 3.141592653589793 +#define _2PI 6.283185307179586 +#define _HALF_PI 1.5707963267948966 +#define _THREE_HALVES_PI 4.71238898038469 + +// NOTE: these functions are not IEEE 754-compliant (the C standard doesn't require them to be), but they're pretty good + +double frexp(double value, int *exp) { + if (value == 0) { + *exp = 0; + return 0; + } + unsigned long u = *(unsigned long *)&value, significand; + *exp = ((u >> 52) & 0x7ff) - 1022; + // replace exponent with 1022 + u &= 0x800fffffffffffff; + u |= 0x3fe0000000000000; + return *(double *)&u; +} + +double ldexp(double x, int exp) { + int e; + double y = frexp(x, &e); + // since x = y * 2^e, x * 2^exp = y * 2^(e+exp) + exp += e; + if (exp < -1022) return 0; + if (exp > 1023) return _INFINITY; + unsigned long pow2 = (unsigned long)(exp + 1023) << 52; + return y * *(double *)&pow2; +} + +double floor(double x) { + if (x >= 0.0) { + if (x > 1073741824.0 * 1073741824.0) + return x; // floats this big must be integers + return (unsigned long)x; + } else { + if (x < -1073741824.0 * 1073741824.0) + return x; // floats this big must be integers + double i = (long)x; + if (x == i) return x; + return i - 1.0; + } +} + +double ceil(double x) { + double f = floor(x); + if (x == f) return f; + return f + 1.; +} + +double fabs(double x) { + // this is better than x >= 0 ? x : -x because it deals with -0 properly + unsigned long u = *(unsigned long *)&x; + u &= 0x7fffffffffffffff; + return *(double *)&u; +} + +double fmod(double x, double y) { + if (y == 0.0) { + errno = EDOM; + return 0.0; + } + return x - (floor(x / y) * y); +} + +double _sin_taylor(double x) { + double i; + double term = x; + // taylor expansion for sin: x - x³/3! + x⁵/5! - ... + + // https://en.wikipedia.org/wiki/Kahan_summation_algorithm + double prev = -1.0; + double sum = 0.0; + double c = 0.0; + for (i = 0.0; i < 100.0 && sum != prev; ++i) { + prev = sum; + double y = term - c; + double t = sum + y; + c = (t - sum) - y; + sum = t; + term *= -(x * x) / ((2.0*i+2.0)*(2.0*i+3.0)); + } + return sum; +} + +double _cos_taylor(double x) { + double i; + double term = 1.0; + // taylor expansion for cos: 1 - x²/2! + x⁴/4! - ... + + // https://en.wikipedia.org/wiki/Kahan_summation_algorithm + double prev = -1.0; + double sum = 0.0; + double c = 0.0; + for (i = 0.0; i < 100.0 && sum != prev; ++i) { + prev = sum; + double y = term - c; + double t = sum + y; + c = (t - sum) - y; + sum = t; + term *= -(x * x) / ((2.0*i+1.0)*(2.0*i+2.0)); + } + return sum; +} + +double sin(double x) { + x = fmod(x, 2.0*_PI); + // the Taylor series works best for small inputs. so, provide _sin_taylor with a value in the range [0,π/2] + if (x < _HALF_PI) + return _sin_taylor(x); + if (x < _PI) + return _sin_taylor(_PI - x); + if (x < _THREE_HALVES_PI) + return -_sin_taylor(x - _PI); + return -_sin_taylor(_2PI - x); +} + +double cos(double x) { + x = fmod(x, 2.0*_PI); + // the Taylor series works best for small inputs. so, provide _cos_taylor with a value in the range [0,π/2] + if (x < _HALF_PI) + return _cos_taylor(x); + if (x < _PI) + return -_cos_taylor(_PI - x); + if (x < _THREE_HALVES_PI) + return -_cos_taylor(x - _PI); + return _cos_taylor(_2PI - x); +} + +double tan(double x) { + return sin(x)/cos(x); +} + +// for sqrt and the inverse trigonometric functions, we use Newton's method +// https://en.wikipedia.org/wiki/Newton%27s_method + +double sqrt(double x) { + if (x < 0.0) { + errno = EDOM; + return _NAN; + } + if (x == 0.0) return 0.0; + if (x == _INFINITY) return _INFINITY; + // we want to find the root of: f(t) = t² - x + // f'(t) = 2t + int exp; + double y = frexp(x, &exp); + if (exp & 1) { + y *= 2; + --exp; + } + // newton's method will be slow for very small or very large numbers. + // so we have ensured that + // 0.5 ≤ y < 2 + // and also x = y * 2^exp; sqrt(x) = sqrt(y) * 2^(exp/2) NB: we've ensured that exp is even + + // 7 iterations seems to be more than enough for any number + double t = y; + t = (y / t + t) * 0.5; + t = (y / t + t) * 0.5; + t = (y / t + t) * 0.5; + t = (y / t + t) * 0.5; + t = (y / t + t) * 0.5; + t = (y / t + t) * 0.5; + t = (y / t + t) * 0.5; + + return ldexp(t, exp>>1); + +} + +double _acos_newton(double x) { + // we want to find the root of: f(t) = cos(t) - x + // f'(t) = -sin(t) + double t = _HALF_PI - x; // reasonably good first approximation + double prev_t = -100.0; + int i; + + for (i = 0; i < 100 && prev_t != t; ++i) { + prev_t = t; + t += (cos(t) - x) / sin(t); + } + return t; +} + +double _asin_newton(double x) { + // we want to find the root of: f(t) = sin(t) - x + // f'(t) = cos(t) + double t = x; // reasonably good first approximation + double prev_t = -100.0; + int i; + + for (i = 0; i < 100 && prev_t != t; ++i) { + prev_t = t; + t += (x - sin(t)) / cos(t); + } + return t; +} + +double acos(double x) { + if (x > 1.0 || x < -1.0) { + errno = EDOM; + return _NAN; + } + // Newton's method doesn't work well near -1 and 1, because f(x) / f'(x) is very large. + if (x > 0.8) + return _asin_newton(sqrt(1-x*x)); + if (x < -0.8) + return _PI-_asin_newton(sqrt(1-x*x)); + + return _acos_newton(x); +} + +double asin(double x) { + if (x > 1.0 || x < -1.0) { + errno = EDOM; + return _NAN; + } + // Newton's method doesn't work well near -1 and 1, because f(x) / f'(x) is very large. + if (x > 0.8) + return _acos_newton(sqrt(1.0-x*x)); + if (x < -0.8) + return -_acos_newton(sqrt(1.0-x*x)); + + return _asin_newton(x); +} + +double atan(double x) { + // the formula below breaks for really large inputs; tan(10^20) as a double is indistinguishable from pi/2 anyways + if (x > 1e20) return _HALF_PI; + if (x < -1e20) return -_HALF_PI; + + // we can use a nice trigonometric identity here + return asin(x / sqrt(1+x*x)); +} + +double atan2(double y, double x) { + if (x == 0.0) { + if (y > 0.0) return _HALF_PI; + if (y < 0.0) return -_HALF_PI; + return 0.0; // this is what IEEE 754 does + } + + double a = atan(y/x); + if (x > 0.0) { + return a; + } else if (y > 0.0) { + return a + _PI; + } else { + return a - _PI; + } +} + +double _exp_taylor(double x) { + double i; + double term = 1.0; + // taylor expansion for exp: 1 + x/1! + x²/2! + ... + + // https://en.wikipedia.org/wiki/Kahan_summation_algorithm + double prev = -1.0; + double sum = 0.0; + double c = 0.0; + for (i = 1.0; i < 100.0 && sum != prev; ++i) { + prev = sum; + double y = term - c; + double t = sum + y; + c = (t - sum) - y; + sum = t; + term *= x / i; + } + return sum; +} + +double exp(double x) { + if (x > 709.782712893384) { + errno = ERANGE; + return _INFINITY; + } + if (x == 0.0) return 1; + if (x < -744.4400719213812) + return 0; + int i, e; + double y = frexp(x, &e); + if (e < 1.0) return _exp_taylor(x); + // the taylor series doesn't work well for large x (positive or negative), + // so we use the fact that exp(y*2^e) = exp(y)^(2^e) + double value = _exp_taylor(y); + for (i = 0; i < e; ++i) + value *= value; + return value; +} + +#define _LOG2 0.6931471805599453 + +double log(double x) { + if (x < 0.0) { + errno = EDOM; + return _NAN; + } + if (x == 0.0) return -_INFINITY; + if (x == 1.0) return 0.0; + int e; + double sum; + double a = frexp(x, &e); + // since x = a * 2^e, log(x) = log(a) + log(2^e) = log(a) + e log(2) + sum = e * _LOG2; + // now that a is in [1/2,1), the series log(a) = (a-1) - (a-1)²/2 + (a-1)³/3 - ... converges quickly + + a -= 1; + // https://en.wikipedia.org/wiki/Kahan_summation_algorithm + double prev = HUGE_VAL; + double c = 0; + double term = a; + double i; + for (i = 1.0; i < 100.0 && sum != prev; ++i) { + prev = sum; + double y = term / i - c; + double t = sum + y; + c = (t - sum) - y; + sum = t; + term *= -a; + } + return sum; +} + +#define _INVLOG10 0.43429448190325176 // = 1/log(10) +double log10(double x) { + return log(x) * _INVLOG10; +} + +double modf(double value, double *iptr) { + double m = fmod(value, 1.0); + if (value >= 0.0) { + *iptr = value - m; + return m; + } else if (m == 0.0) { + *iptr = value; + return 0.0; + } else { + *iptr = value - m + 1.0; + return m - 1.0; + } +} + +// double raised to the power of an integer +double _dpowi(double x, unsigned long y) { + double result = 1.0; + if (y & 1) { + --y; + result *= x; + } + if (y > 0) { + double p = _dpowi(x, y >> 1); + result *= p * p; + } + return result; +} + +double pow(double x, double y) { + if (x > 0.0) { + return exp(y * log(x)); + } else if (x < 0.0) { + if (fmod(y, 1.0) != 0) { + errno = EDOM; + return _NAN; + } + if (y > 1.6602069666338597e+19) + return x < -1. ? -_INFINITY : 0.; + if (y < -1.6602069666338597e+19) + return x < -1. ? 0. : -_INFINITY; + return _dpowi(x, (unsigned long)y); + } else { + if (y < 0) { + errno = EDOM; + return _NAN; + } + if (y > 0) { + // 0^x = 0 for x>0 + return 0.; + } + // 0^0 = 1 + return 1.; + } +} + +double cosh(double x) { + double e = exp(x); + return (e + 1./e) * 0.5; +} + +double sinh(double x) { + double e = exp(x); + return (e - 1./e) * 0.5; +} + +double tanh(double x) { + if (x > 20.0) return 1.; + if (x < -20.0) return -1.; + double e = exp(x); + double f = 1./e; + return (e - f) / (e + f); +} +#endif // _MATH_H diff --git a/05/stdc_common.h b/05/stdc_common.h index 9e5d7a6..af30b06 100644 --- a/05/stdc_common.h +++ b/05/stdc_common.h @@ -424,6 +424,7 @@ double strtod(const char *nptr, char **endptr) { sum = t; } + if (sum == _INFINITY) errno = ERANGE; if (endptr) *endptr = nptr; return sum * sign; } diff --git a/05/tokenize.b b/05/tokenize.b index 809c9ff..56838c5 100644 --- a/05/tokenize.b +++ b/05/tokenize.b @@ -287,6 +287,7 @@ function tokenize fraction = strtoi(&in, 10) ; e.g. to turn 35 into .35, multiply by 10^-2 pow10 = p - in + ;putnln_signed(pow10) if pow10 < -400 goto bad_float :float_no_fraction ; construct the number integer + fraction*10^pow10 @@ -303,16 +304,13 @@ function tokenize exponent = new_exponent significand += right_shift(integer, exponent) :float_no_integer + if *1in == 'e goto float_exponent if *1in == 'E goto float_exponent :float_have_significand_and_exponent if significand == 0 goto float_zero normalize_float(&significand, &exponent) - ; putn(significand) - ; putc(32) - ; putn_signed(exponent) - ; putc(10) ; make number round to the nearest representable float roughly (this is what gcc does) ; this fails for 5e-100 probably because of imprecision, but mostly works significand += 15 @@ -321,6 +319,7 @@ function tokenize exponent += 5 exponent += 52 ; 1001010111... => 1.001010111... n = leftmost_1bit(significand) + exponent += n - 52 ; in most cases, this is 0, but sometimes it turns out to be 1. b = 1 < n significand &= ~b data = significand |