/* fp_api.c
** Convenience functions for Owl extensions that use floating point
*/

#include <stdint.h>
#include <math.h>
#include <limits.h>
#include "ovm.h"

// Convert a double to an owl rational
//
word mkrat_approx(double xin)
{
    const double eps = 1.0e-014; // "near zero" termination limit for remainder

    double x = fabs(xin);
    if (x > (double )LONG_MAX) return IFALSE; // not representable with int64_t
    if ((uint64_t )x > (uint64_t )LONG_MAX) return IFALSE; // not representable with int64_t
    if (x != x) return IFALSE; // NaN

    uint64_t a2 = (uint64_t )x; // floor, since x is positive; fits in int64_t

    uint64_t pnm1 = 1u; // p[n-1] initialized to p[1]
    uint64_t qnm1 = 0u; // q[n-1] initialized to q[1]
    uint64_t pn = a2;   // p[n  ] initialized to p[2]
    uint64_t qn = 1u;   // q[n  ] initialized to q[2]

    double r = x - (double )a2; // the remainder of the continuing fraction

    while (r > eps)
    {
        // next a[n+1] and p[n+1] / q[n+1]
        uint64_t anp1 = (uint64_t )(1.0 / r);
        uint64_t pnp1 = anp1 * pn + pnm1;
        uint64_t qnp1 = anp1 * qn + qnm1;

        // test for overflow and exceeding int64_t limits
        if (   anp1 == 0   // "can't" happen, but protects against divide by zero
            || ((anp1 * pn) / anp1) != pn  // multiply overflowed u64
            || ((anp1 * qn) / anp1) != qn  // multiply overflowed u64
            || (anp1 * pn) > LONG_MAX      // multiply overflowed u63
            || (anp1 * qn) > LONG_MAX      // multiply overflowed u63
            || pnp1 > LONG_MAX  // sum too large
            || qnp1 > LONG_MAX  // sum too large
        ) {
            // XXX - semi-convergents could be added here:
            // for (K from an/2 to an-1) K * pn + pnm1, K * qn + qnm1
            // but this requires weeding out semi-convergents that are not best
            // approximations, and perhaps not worth the trouble
            break;
        }
        pnm1 = pn;
        pn = pnp1;
        qnm1 = qn;
        qn = qnp1;
        r = (1.0 / r) - (double )anp1;
    }
    pn = (xin < 0) ? -pn : pn; // apply sign of input
    return mkrat(pn, qn);
}

// Convert an owl rational to a double
//
double cdouble(word x)
{
    if (allocp(x))
    {
        if (is_type(*(word *)x,TRAT))
        {
            int64_t p = cnum(car(x));
            int64_t q = cnum(cdr(x));
            return (double )((long double )p / (long double )q);
        }
        if (header(x) == NUMHDR || header(x) == NUMNHDR)
        {
            return (double )cnum(x);
        }
    }
    else if (is_type(x,TNUM)||is_type(x,TNUMN))
    {
        return (double )cnum(x);
    }
    return 0.0/0.0; // NaN
}