/*
 * Copyright (c) 2011-2024 Ross Cunniff
 *
 * 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.
 */

/* OADL libstd - most of the same functionality of C libm/libc, but handles
 * vectors as well as scalars
 */

namespace math {

/******************************************************************************
 * Math constants - single precision, to 10 significant digits
 *****************************************************************************/
const E =               2.718281828;    // e
const LOG2E =           1.442695041;    // log[2](e)
const LOG10E =          .4342944819;    // log[10](e)
const LN2 =             .6931471806;    // log[e](2)
const LN10 =            2.302585093;    // log[e](10)
const PI =              3.141592654;    // pi
const PI_2 =            1.570796327;    // pi/2
const PI_4 =            .7853981634;    // pi/4
const ONE_DIV_PI =      .3183098862;    // 1/pi
const TWO_DIV_PI =      .6366197724;    // 2/pi
const TWO_DIV_SQRTPI =  1.128379167;    // 2/sqrt(pi)
const SQRT2 =           1.414213562;    // sqrt(2)
const SQRT1_2 =         .7071067812;    // 1/sqrt(2)
const EPSILON =         2.384185791e-07; // Minimum increment of 1.0

/******************************************************************************
 * Math constants - double precision, to 18 significant digits
 *****************************************************************************/
const D_E =             2.71828182845904524D;   // e
const D_LOG2E =         1.44269504088896341D;   // log[2](e)
const D_LOG10E =        .434294481903251828D;   // log[10](e)
const D_LN2 =           .693147180559945309D;   // log[e](2)
const D_LN10 =          2.30258509299404568D;   // log[e](10)
const D_PI =            3.14159265358979324D;   // pi
const D_PI_2 =          1.57079632679489662D;   // pi/2
const D_PI_4 =          .785398163397448310D;   // pi/4
const D_1_PI =          .318309886183790672D;   // 1/pi
const D_2_PI =          .636619772367581343D;   // 2/pi
const D_2_SQRTPI =      1.12837916709551257D;   // 2/sqrt(pi)
const D_SQRT2 =         1.41421356237309505D;   // sqrt(2)
const D_SQRT1_2 =       .707106781186547524D;   // 1/sqrt(2)
const D_EPSILON =       4.4408920985006261617e-16D; // Minimum increment of 1.0

/******************************************************************************
 * Math functions
 *
 * Note that all arguments may be array-valued, but all must contain only
 * floating point data (Float, Half, or Double), except for ldexp which
 * requires integral data for the exponent
 *****************************************************************************/
extern acos;      // r = acos(x)      inverse cosine
extern asin;      // r = asin(x)      inverse sine
extern atan;      // r = atan(x)      one-parameter inverse tangent
extern atan2;     // r = atan2(y,x)   two-parameter inverse tangent
extern ceil;      // r = ceil(x)      smallest integer not less than x
extern cos;       // r = cos(x)       cosine
extern cosh;      // r = cosh(x)      hyperbolic cosine
extern exp;       // r = exp(x)       exponential function
extern fabs;      // r = fabs(x)      absolute value
extern floor;     // r = floor(x)     largest integer not greater than x
extern fmod;      // r = fmod(x,y)    floating-point remainder: x - y*(int)(x/y)
extern frexp;     // [m,e] = frexp(x) break number into mantissa and exponent
extern ldexp;     // r = ldexp([m,e]) compose number from mantissa and exponent
extern log;       // r = log(x)       natural logarithm
extern log10;     // r = log10(x)     base-10 logarithm
extern modf;      // [f,i] = modf(x)  integral and fractional parts of x
extern pow;       // r = pow(x,y)     raise x to the power of y, x**y
extern sin;       // r = sin(x)       sine
extern sinh;      // r = sinh(x)      hyperbolic sine
extern sqrt;      // r = sqrt(x)      square root
extern tan;       // r = tan(x)       tangent
extern tanh;      // r = tanh(x)      hyperbolic tangent
extern isinf;     // r = isinf(x)     true if IEEE Inf, false otherwise
extern isnan;     // r = isnan(x)     true if IEE NaN, false otherwise
extern isint;     // r = isint(x)     true if no fractional bits, else false
extern gamma;     // r = gamma(x)     calculate the real gamma of x
extern factorial; // r = factorail(x) calculate the real factorial of x
extern binomial;  // r = binomial(x,y) returns x!/(y!*(x-y)!)
extern hypot;     // r = hypot(v1,v2) returns ((v1-v2)*2).reduce(`+)**0.5
extern vlen;      // r = vlen(vec)    returns (vec*vec).reduce(`+)**0.5
extern polar;     // r = polar(vec)   returns [rad,ang], etc.
extern cartesian; // r = cartesian(vec) returns rad*[ang.sin(),ang.cos()] etc.
}

namespace linalg {
/******************************************************************************
 * linear algebra functions
 *****************************************************************************/
extern matinv;  // r = matinv(m)        mat must be a 2D packed array
extern svd;     // {UDV} = svd(m)       Singular Value Decomposition of 2D m
extern solve;   // x = solve(udv,b)     Solve (udv)x=b for x
extern dot;       // r = dot(v1,v2)     returns (v1*v2).reduce(`+)
extern cross;     // r = cross(v1,v2)   returns 3-dimensional cross product
extern unitize;   // r = unitize(vec)   returns (v1**2).reduce(`+)**0.5
extern ucross;    // r = ucross(v1,v2)  returns unit 3-dimensional cross product
}

namespace std {

/******************************************************************************
 * ctype functions
 *****************************************************************************/
extern isalnum;   // r = isalnum(c)
extern isalpha;   // r = isalpha(c)
extern iscntrl;   // r = iscntrl(c)
extern isdigit;   // r = isdigit(c)
extern isgraph;   // r = isgraph(c)
extern islower;   // r = islower(c)
extern isprint;   // r = isprint(c)
extern ispunct;   // r = ispunct(c)
extern isspace;   // r = isspace(c)
extern isupper;   // r = isupper(c)
extern isxdigit;  // r = isxdigit(c)
extern tolower;   // r = tolower(c)
extern toupper;   // r = toupper(c)

/******************************************************************************
 * string functions
 *****************************************************************************/
extern strstr;    // r = strstr(big,little) find little in big
extern strtok;    // r = strtok(str, sepchars) tokenize a string

}

// Public versions of the calls
default public
    acos = math::acos,  // r = x.acos()      inverse cosine
    asin = math::asin,  // r = x.asin()      inverse sine
    atan = math::atan,  // r = x.atan()      one-parameter inverse tangent
    atan2 = math::atan2, // r = y.atan2(x)   two-parameter inverse tangent
    ceil = math::ceil,  // r = x.ceil()      smallest integer not less than x
    cos = math::cos,    // r = x.cos()       cosine
    cosh = math::cosh,  // r = x.cosh()      hyperbolic cosine
    exp = math::exp,    // r = x.exp()       exponential function
    fabs = math::fabs,  // r = x.fabs()      absolute value
    floor = math::floor,// r = x.floor()     largest integer not greater than x
    fmod = math::fmod,  // r = x.fmod(y)     float remainder: x-y*(int)(x/y)
    frexp = math::frexp, // {m,e} = x.frexp() break number into mant. and expon.
    ldexp = math::ldexp, // r = m.ldexp(e)   make number from mant. and expon.
    log = math::log,    // r = x.log()       natural logarithm
    log10 = math::log10, // r = x.log10()    base-10 logarithm
    modf = math::modf,  // {f,i} = x.modf()  integral and fractional parts of x
    pow = math::pow,    // r = x.pow(y)      raise x to the power of y, x**y
    sin = math::sin,    // r = x.sin()       sine
    sinh = math::sinh,  // r = x.sinh()      hyperbolic sine
    sqrt = math::sqrt,  // r = x.sqrt()      square root
    tan = math::tan,    // r = x.tan()       tangent
    tanh = math::tanh,  // r = x.tanh()      hyperbolic tangent
    isinf = math::isinf,// r = x.isinf()     true if IEEE Inf, false otherwise
    isnan = math::isnan,// r = x.isnan()     true if IEEE NaN, false otherwise
    isint = math::isint,// r = x.isint()     true if no frac. bits, else false
    gamma = math::gamma,// r = x.gamma()
    factorial = math::factorial, // r = x.factorial()
    binomial = math::binomial,  // r = x.binomial(y) returns x!/(y!*(x-y)!)
    hypot = math::hypot, // r = v1.hypot(v2) returns ((v1-v2)*2).reduce(`+)**0.5
    vlen = math::vlen,   // r = vlen(vec)     returns (vec*vec).reduce(`+)**0.5
    polar = math::polar, // r = vec.polar()    returns [rad,ang], etc.
    cartesian = math::cartesian, // r = vec.cartesian() returns
                                 //              rad*[ang.sin(),ang.cos()] etc.

    matinv = linalg::matinv, // r = m.matinv()    mat must be a 2D packed array
    svd = linalg::svd,       // {UDV} = m.svd()   Singular Value Decomposition
    solve = linalg::solve,   // x = udv.solve(b)  Solve (udv)x=b for x
    dot = linalg::dot,       // r = dot(v1,v2)    returns (v1*v2).reduce(`+)
    cross = linalg::cross,   // r = cross(v1,v2)  returns n-dim cross product
    unitize = linalg::unitize, // r = unitize(vec)  returns vec/vec.vlen()
    ucross = linalg::ucross, // r = ucross(v1,v2) returns unit n-dim cross prod

    isalnum = std::isalnum,   // r = c.isalnum()
    isalpha = std::isalpha,   // r = c.isalpha()
    iscntrl = std::iscntrl,   // r = c.iscntrl()
    isdigit = std::isdigit,   // r = c.isdigit()
    isgraph = std::isgraph,   // r = c.isgraph()
    islower = std::islower,   // r = c.islower()
    isprint = std::isprint,   // r = c.isprint()
    ispunct = std::ispunct,   // r = c.ispunct()
    isspace = std::isspace,   // r = c.isspace()
    isupper = std::isupper,   // r = c.isupper()
    isxdigit = std::isxdigit, // r = c.isxdigit()
    tolower = std::tolower,   // r = c.tolower()
    toupper = std::toupper,   // r = c.toupper()

    strstr = std::strstr,     // r = big.strstr(little) find little in big
    strtok = std::strtok;     // r = str.strtok(seps) split str based on seps

using extern "libstd";