/*
 * fputest.c - exercise and sanity-check a 387 FPU (old-compiler friendly, C89)
 *
 * Build (hardware FPU):
 *   gcc -O2 -o fputest fputest.c -lm
 *
 * Run:
 *   ./fputest        normal tests + timing
 *   ./fputest -x     also test FP exception delivery (expects SIGFPE)
 *
 * Exit code 0 = all tests passed, 1 = at least one failure.
 */
#include <stdio.h>
#include <math.h>
#include <time.h>
#include <string.h>
#include <signal.h>
#include <stdlib.h>
#include <unistd.h>

#ifdef USE_FPU_CONTROL
#include <fpu_control.h>
#endif

static int failures = 0;

/* volatile so the compiler can't fold everything at compile time
   and skip the FPU entirely */
static volatile double vtwo = 2.0, vone = 1.0, vzero = 0.0, v3 = 3.0;

static void check(const char *name, double got, double expect, double tol)
{
    double diff = fabs(got - expect);
    if (diff <= tol) {
        printf("  PASS  %-14s got %.15g\n", name, got);
    } else {
        printf("  FAIL  %-14s got %.15g expected %.15g (diff %g)\n",
               name, got, expect, diff);
        failures++;
    }
}

static void test_basic(void)
{
    printf("Basic arithmetic:\n");
    check("add", vone + vtwo, 3.0, 0.0);
    check("mul", vtwo * v3, 6.0, 0.0);
    check("div 1/3*3", (vone / v3) * v3, 1.0, 1e-15);
    check("sqrt(2)", sqrt(vtwo), 1.41421356237309505, 1e-15);
}

static void test_transcendental(void)
{
    double pi;
    printf("Transcendental functions:\n");
    pi = atan(vone) * 4.0;
    check("pi=4*atan(1)", pi, 3.14159265358979324, 1e-15);
    check("sin(pi/6)", sin(pi / 6.0), 0.5, 1e-15);
    check("cos(pi/3)", cos(pi / 3.0), 0.5, 1e-15);
    check("tan(pi/4)", tan(pi / 4.0), 1.0, 1e-15);
    check("exp(1)", exp(vone), 2.71828182845904524, 1e-15);
    check("log(e)", log(exp(vone)), 1.0, 1e-15);
    check("log10(1000)", log10(1000.0 * vone), 3.0, 1e-14);
    check("pow(2,0.5)", pow(vtwo, 0.5), 1.41421356237309505, 1e-15);
}

static void test_precision(void)
{
    /* 80-bit extended precision check: only true x87 hardware
       gives the eps of roughly 1.08e-19 for long double */
    long double e = 1.0L, prev = 1.0L;
    printf("Extended precision:\n");
    while (1.0L + e != 1.0L) {
        prev = e;
        e = e / 2.0L;
    }
    printf("  long double epsilon = %Lg (expect ~1.08e-19 on x87)\n", prev);
    if (prev > 1.2e-19L || prev < 1.0e-19L) {
        printf("  FAIL  not 64-bit mantissa\n");
        failures++;
    } else {
        printf("  PASS\n");
    }
}

static void test_series(void)
{
    /* Leibniz series for pi/4 and a harmonic sum: long chains of
       dependent FP ops, good at catching flaky hardware */
    long i;
    double s = 0.0, sign = 1.0, h = 0.0;
    printf("Series (100000 iterations):\n");
    for (i = 0; i < 100000L; i++) {
        s += sign / (2.0 * i + 1.0);
        sign = -sign;
        h += 1.0 / (i + 1.0);
    }
    check("leibniz*4", s * 4.0, 3.14159265358979324, 2e-5);
    check("harmonic", h, 12.0901461298634, 1e-9);
}

static void test_special(void)
{
    double inf, nan_;
    printf("Special values:\n");
    inf = vone / vzero;   /* masked by default: gives +inf */
    nan_ = vzero / vzero; /* gives NaN */
    if (inf > 1e308) printf("  PASS  1/0 = inf\n");
    else { printf("  FAIL  1/0 = %g\n", inf); failures++; }
    if (nan_ != nan_) printf("  PASS  0/0 = NaN\n");
    else { printf("  FAIL  0/0 = %g\n", nan_); failures++; }
}

static void test_speed(void)
{
    clock_t t0, t1;
    long i, n = 2000000L;
    double x = 1.0, secs;
    printf("Timing (%ld mul+add+div+sqrt iterations):\n", n);
    t0 = clock();
    for (i = 0; i < n; i++) {
        x = sqrt(x * 1.0000001 + 0.5) / 1.0000001;
    }
    t1 = clock();
    secs = (double)(t1 - t0) / CLOCKS_PER_SEC;
    printf("  result %.12g, %.2f s (hardware should be several times\n"
           "  faster than emulation: compare with 'no387' boot option)\n",
           x, secs);
}

static void on_sigfpe(int sig)
{
    static const char msg[] = "  PASS  got SIGFPE: exception delivery works\n";
    (void)sig;
    write(1, msg, sizeof(msg) - 1);
    _exit(failures ? 1 : 0);
}

static void test_exception(void)
{
#ifdef USE_FPU_CONTROL
    fpu_control_t cw;
    printf("Exception delivery (unmasking divide-by-zero):\n");
    signal(SIGFPE, on_sigfpe);
    _FPU_GETCW(cw);
    cw &= ~_FPU_MASK_ZM;
    _FPU_SETCW(cw);
    printf("  dividing by zero now...\n");
    fflush(stdout);
    printf("  value: %f\n", vone / vzero);
    /* the x87 reports errors on the *next* FP instruction */
    printf("  value: %f\n", vone + vone);
    printf("  FAIL  no SIGFPE received (ERROR# wiring / IRQ13 problem?)\n");
    failures++;
#else
    printf("Exception test not compiled in. Rebuild with -DUSE_FPU_CONTROL\n");
#endif
}

int main(int argc, char **argv)
{
    int do_exc = (argc > 1 && strcmp(argv[1], "-x") == 0);

    test_basic();
    test_transcendental();
    test_precision();
    test_series();
    test_special();
    test_speed();
    if (do_exc)
        test_exception();

    printf("\n%s (%d failure%s)\n", failures ? "FAILED" : "ALL PASSED",
           failures, failures == 1 ? "" : "s");
    return failures ? 1 : 0;
}