/* ext.c - a test for mkrat() and mkrat_approx()
**
Copyright (c) 2024 Doug Currie

Permission is hereby granted, free of charge, to any person obtaining a copy
of this software and associated documentation files (the "Software"), to deal
in the Software without restriction, including without limitation the rights
to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
copies of the Software, and to permit persons to whom the Software is
furnished to do so, subject to the following conditions:

The above copyright notice and this permission notice shall be included in
all copies or substantial portions of the Software.

THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
THE SOFTWARE.
**
*/

#include "ovm.h"

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
//#include <float.h>
#include <math.h>
#include <limits.h>

static int verbose = 0;

static int64_t rat_p(word r)
{
    //fprintf(stderr,"rat_p \n");
    if (allocp(r))
    {
        if (is_type(*(word *)r,TRAT))
        {
            return cnum(car(r));
        }
    }
    if (is_type(r,TNUM)||is_type(r,TNUMN))
    {
        return cnum(r);
    }
    return 0;
}

static int64_t rat_q(word r)
{
    //fprintf(stderr,"rat_q \n");
    if (allocp(r))
    {
        if (is_type(*(word *)r,TRAT))
        {
            return cnum(cdr(r));
        }
    }
    if (is_type(r,TNUM)||is_type(r,TNUMN))
    {
        return 1;
    }
    return 0;
}

static double randdbl(void)
{
    uint64_t r53 = (((uint64_t)(rand()) << 26) ^ (uint64_t)rand()) & ((1uLL << 53) - 1uLL);
    return (double )r53 / 9007199254740991.0; // 2^53 - 1
}

static inline uint64_t extractbits(double x)
{
    union { double x; uint64_t u;} m;
    m.x = x;
    return m.u;
}

// This tests both the accuracy of mkrat_approx
// and that mkrat_approx <- roundtrips -> cdouble

static int64_t qerrors[4] = {0};
static int64_t max_err = 0;

static void test_init(unsigned int seed)
{
    memcpy(qerrors, (int64_t[]) {0, 0, 0, 0}, sizeof qerrors);
    max_err = 0;
    srand(seed); // make the test repeatable
}
static void test_step(void)
{
    double u = randdbl();
    int exp = 64; // -32
    while (1)
    {
        double x = ldexp(u, exp); // u * 2^exp
        while (x <= 1.0e-014) x *= 2.0; // get above the eps for mkrat
        while (x >= (double )LONG_MAX) x /= 2.0; // get below the limit for mkrat
        word r = mkrat_approx(x);
        uint64_t xb = extractbits(x);
        double c = cdouble(r);
        uint64_t cb = extractbits(c);
        if ((xb & 0x8000000000000000uL) != (cb & 0x8000000000000000uL))
        {
            // mismatch signs
            printf("Oops sign %d %a ne %a = %lld / %lld\n", exp, x, c, rat_p(r), rat_q(r));
        }
        if ((xb & 0x7ff0000000000000uL) != (cb & 0x7ff0000000000000uL))
        {
            // mismatch exponents
            printf("Oops expt %d %a ne %a = %lld / %lld\n", exp, x, c, rat_p(r), rat_q(r));
        }
        int64_t err = (xb & 0x000fffffffffffffuL) - (cb & 0x000fffffffffffffuL);
        if (err != 0)
        {
            // mismatch mantissas
            int64_t ulp = (err < 0) ? -err : err;
            if      (ulp <= 1) qerrors[1] += 1;
            else if (ulp <= 3) qerrors[2] += 1;
            else               qerrors[3] += 1;
            if (ulp > max_err) max_err = err;
        }
        else
        {
            qerrors[0] += 1;
        }
        if (exp < 0)
        {
            exp = -exp;
        }
        else if (exp > 0)
        {
            exp = -(exp / 2);
        }
        else
        {
            break; // exp was 0
        }
    }
}

static word test_report(void)
{
    printf("Correct: %lld Errors by ULP: %lld,%lld,%lld Max err: %lld\n"
            , qerrors[0], qerrors[1], qerrors[2], qerrors[3], max_err);
    return onum(max_err, 1);
}

word
prim_custom(int op, word a, word b, word c)
{
    switch (op)
    {
        case 100:
        {
            printf("Hello from C!\n");
            return ITRUE;
        }
        case 101:
        {
            double d = (double )((long double )cnum(a) / (long double )cnum(b));
            word r = mkrat_approx(d);
            //if (r == IFALSE) printf("Made rat: %lld  / %lld -> #f\n", cnum(a), cnum(b));
            //else printf("Made rat: %lld / %lld\n", rat_p(r), rat_q(r));
            printf("Made rat: %lld  / %lld -> %a ->\n", cnum(a), cnum(b), d);
            return r;
        }
        case 102:
        {
            test_init((unsigned int )cnum(a)); // arg is seed
            return ITRUE;
        }
        case 103:
        {
            test_step();
            return ITRUE;
        }
        case 104:
        {
            return test_report();
        }
        case 105:
        {
            //printf("C onum(%lld, %lld) -> %lld\n", cnum(a), cnum(b), cnum(onum(cnum(a), cnum(b))));
            return onum(cnum(a), cnum(b));
        }
        case 106:
        {
            return mkrat(cnum(a), cnum(b));
        }
        case 107:
        {
            double x = cdouble(IFALSE); // pass in #f to get a NaN
            if (x != x) printf("ok cdouble(#f) -> NaN\n");
            else printf("ng cdouble(#f) -> %g\n", x);
            break;
        }
    }
    return IFALSE;
}