108 lines
2.4 KiB
C
108 lines
2.4 KiB
C
// Part of JoeyLib, Copyright (c) 2026 Scott Duensing, under the MIT License.
|
|
// See LICENSE.md for license information.
|
|
// SPDX-License-Identifier: MIT
|
|
//
|
|
// Tiny libm replacement on the m68k-atari-mint toolchain, built into the
|
|
// ST's libjoey.a. Only the functions libxmp-lite once called are implemented.
|
|
// Accuracy targets: a few hundred parts per million, plenty for
|
|
// note-frequency math.
|
|
|
|
#include "math.h"
|
|
|
|
double ceil(double x) {
|
|
long i = (long)x;
|
|
if ((double)i < x) {
|
|
i++;
|
|
}
|
|
return (double)i;
|
|
}
|
|
|
|
|
|
double floor(double x) {
|
|
long i = (long)x;
|
|
if ((double)i > x) {
|
|
i--;
|
|
}
|
|
return (double)i;
|
|
}
|
|
|
|
|
|
// Natural log via 2*atanh((x-1)/(x+1)) after a power-of-two range
|
|
// reduction into [1, 2). The atanh series converges fast on that
|
|
// interval.
|
|
static double logSeries(double m) {
|
|
double y;
|
|
double y2;
|
|
|
|
y = (m - 1.0) / (m + 1.0);
|
|
y2 = y * y;
|
|
return 2.0 * y * (1.0 + y2 * (1.0 / 3.0 + y2 * (1.0 / 5.0 + y2 * (1.0 / 7.0 + y2 * (1.0 / 9.0)))));
|
|
}
|
|
|
|
|
|
double log(double x) {
|
|
int k;
|
|
|
|
if (x <= 0.0) {
|
|
return 0.0;
|
|
}
|
|
k = 0;
|
|
while (x >= 2.0) {
|
|
x *= 0.5;
|
|
k++;
|
|
}
|
|
while (x < 1.0) {
|
|
x *= 2.0;
|
|
k--;
|
|
}
|
|
return (double)k * M_LN2 + logSeries(x);
|
|
}
|
|
|
|
|
|
// 2^f for f in [0, 1) -- 5-term polynomial good to ~1e-6.
|
|
static double exp2Frac(double f) {
|
|
return 1.0 + f * (0.69314718 + f * (0.24022650 + f * (0.05550410 + f * (0.00961813 + f * 0.00133335))));
|
|
}
|
|
|
|
|
|
double pow(double base, double exp) {
|
|
double k;
|
|
long ki;
|
|
double frac;
|
|
double result;
|
|
|
|
if (base <= 0.0) {
|
|
return 0.0;
|
|
}
|
|
if (exp == 0.0) {
|
|
return 1.0;
|
|
}
|
|
k = exp * (log(base) / M_LN2);
|
|
ki = (long)k;
|
|
if (k < 0.0 && (double)ki != k) {
|
|
ki--;
|
|
}
|
|
frac = k - (double)ki;
|
|
result = exp2Frac(frac);
|
|
while (ki > 0) {
|
|
result *= 2.0;
|
|
ki--;
|
|
}
|
|
while (ki < 0) {
|
|
result *= 0.5;
|
|
ki++;
|
|
}
|
|
return result;
|
|
}
|
|
|
|
|
|
double trunc(double x) {
|
|
return (double)(long)x;
|
|
}
|
|
|
|
|
|
float ceilf (float x) { return (float)ceil ((double)x); }
|
|
float floorf(float x) { return (float)floor((double)x); }
|
|
float logf (float x) { return (float)log ((double)x); }
|
|
float powf (float b, float e) { return (float)pow ((double)b, (double)e); }
|
|
float truncf(float x) { return (float)(long)x; }
|