// Tiny libm replacement used by libxmp-lite on the m68k-atari-mint
// toolchain. Only the functions libxmp-lite calls 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; }