dec64_math.c

/*
dec64_math.c
Elementary functions for DEC64.

dec64.com
2022-09-07
Public Domain

No warranty.

This file is a placeholder. It should be replaced with functions that are
faster and more accurate.
*/

#include <stdlib.h>
#include <stdlib.h>
#include "dec64.h"
#include "dec64_math.h"

#define D_E              0x6092A113D8D574F0LL
#define D_HALF                        0x5FFLL
#define D_HALF_PI        0x37CE4F32BB21A6F0LL
#define D_NHALF_PI       0xC831B0CD44DE59F0LL
#define D_NPI            0x9063619A89BCB4F0LL
#define D_PI             0x6F9C9E6576434CF0LL
#define D_2PI            0x165286144ADA42F1LL

#define FAC              93

dec64 factorials[FAC] = {
    (1LL << 8) + 0,
    (1LL << 8) + 0,
    (2LL << 8) + 0,
    (6LL << 8) + 0,
    (24LL << 8) + 0,
    (120LL << 8) + 0,
    (720LL << 8) + 0,
    (5040LL << 8) + 0,
    (40320LL << 8) + 0,
    (362880LL << 8) + 0,
    (3628800LL << 8) + 0,
    (39916800LL << 8) + 0,
    (479001600LL << 8) + 0,
    (6227020800LL << 8) + 0,
    (87178291200LL << 8) + 0,
    (1307674368000LL << 8) + 0,
    (20922789888000LL << 8) + 0,
    (355687428096000LL << 8) + 0,
    (6402373705728000LL << 8) + 0,
    (12164510040883200LL << 8) + 1,
    (24329020081766400LL << 8) + 2,
    (5109094217170944LL << 8) + 4,        /* last exact */
    (11240007277776077LL << 8) + 5,
    (25852016738884977LL << 8) + 6,
    (6204484017332394LL << 8) + 8,
    (15511210043330986LL << 8) + 9,
    (4032914611266056LL << 8) + 11,
    (10888869450418352LL << 8) + 12,
    (30488834461171386LL << 8) + 13,
    (8841761993739702LL << 8) + 15,
    (26525285981219106LL << 8) + 16,
    (8222838654177923LL << 8) + 18,
    (26313083693369353LL << 8) + 19,
    (8683317618811886LL << 8) + 21,
    (29523279903960414LL << 8) + 22,
    (10333147966386145LL << 8) + 24,
    (3719933267899012LL << 8) + 26,
    (13763753091226345LL << 8) + 27,
    (5230226174666011LL << 8) + 29,
    (20397882081197443LL << 8) + 30,
    (8159152832478977LL << 8) + 32,
    (33452526613163807LL << 8) + 33,
    (14050061177528799LL << 8) + 35,
    (6041526306337384LL << 8) + 37,
    (26582715747884488LL << 8) + 38,
    (11962222086548019LL << 8) + 40,
    (5502622159812089LL << 8) + 42,
    (25862324151116818LL << 8) + 43,
    (12413915592536073LL << 8) + 45,
    (6082818640342676LL << 8) + 47,
    (30414093201713378LL << 8) + 48,
    (15511187532873823LL << 8) + 50,
    (8065817517094388LL << 8) + 52,
    (4274883284060026LL << 8) + 54,
    (23084369733924138LL << 8) + 55,
    (12696403353658276LL << 8) + 57,
    (7109985878048635LL << 8) + 59,
    (4052691950487722LL << 8) + 61,
    (23505613312828786LL << 8) + 62,
    (13868311854568984LL << 8) + 64,
    (8320987112741390LL << 8) + 66,
    (5075802138772248LL << 8) + 68,
    (31469973260387938LL << 8) + 69,
    (19826083154044401LL << 8) + 71,
    (12688693218588416LL << 8) + 73,
    (8247650592082471LL << 8) + 75,
    (5443449390774431LL << 8) + 77,
    (3647111091818869LL << 8) + 79,
    (24800355424368306LL << 8) + 80,
    (17112245242814131LL << 8) + 82,
    (11978571669969892LL << 8) + 84,
    (8504785885678623LL << 8) + 86,
    (6123445837688609LL << 8) + 88,
    (4470115461512684LL << 8) + 90,
    (33078854415193864LL << 8) + 91,
    (24809140811395398LL << 8) + 93,
    (18854947016660503LL << 8) + 95,
    (14518309202828587LL << 8) + 97,
    (11324281178206298LL << 8) + 99,
    (8946182130782975LL << 8) + 101,
    (7156945704626380LL << 8) + 103,
    (5797126020747368LL << 8) + 105,
    (4753643337012842LL << 8) + 107,
    (3945523969720659LL << 8) + 109,
    (33142401345653533LL << 8) + 110,
    (28171041143805503LL << 8) + 112,
    (24227095383672732LL << 8) + 114,
    (21077572983795277LL << 8) + 116,
    (18548264225739844LL << 8) + 118,
    (16507955160908461LL << 8) + 120,
    (14857159644817615LL << 8) + 122,
    (1352001527678403LL << 8) + 124,
    (12438414054641307LL << 8) + 126
};

dec64 dec64_acos(dec64 slope) {
    dec64 result = dec64_subtract(
        D_HALF_PI,
        dec64_asin(slope)
    );
    return result;
}

dec64 dec64_asin(dec64 slope) {
    if (dec64_is_equal(slope, DEC64_ONE) == DEC64_TRUE) {
        return D_HALF_PI;
    }
    if (dec64_is_equal(slope, DEC64_NEGATIVE_ONE) == DEC64_TRUE) {
        return D_NHALF_PI;
    }
    if (
        dec64_is_nan(slope) == DEC64_TRUE ||
        dec64_is_less(DEC64_ONE, dec64_abs(slope)) == DEC64_TRUE
    ) {
        return DEC64_NULL;
    }
    dec64 bottom = DEC64_TWO;
    dec64 factor = slope;
    dec64 x2 = dec64_multiply(slope, slope);
    dec64 result = factor;
    while (1) {
        factor = dec64_divide(
            dec64_multiply(
                dec64_multiply(dec64_add(DEC64_NEGATIVE_ONE, bottom), x2),
                factor
            ),
            bottom
        );
        dec64 progress = dec64_add(
            result,
            dec64_divide(factor, dec64_add(DEC64_ONE, bottom))
        );
        if (result == progress) {
            break;
        }
        result = progress;
        bottom = dec64_add(bottom, DEC64_TWO);
    }
    return result;
}

dec64 dec64_atan(dec64 slope) {
    return dec64_asin(
        dec64_divide(
            slope,
            dec64_sqrt(dec64_add(DEC64_ONE, dec64_multiply(slope, slope)))
        )
    );
}

dec64 dec64_atan2(dec64 y, dec64 x) {
    if (dec64_is_zero(x) == DEC64_TRUE) {
        if (dec64_is_zero(y) == DEC64_TRUE) {
            return DEC64_NULL;
        } else if (y < 0) {
            return D_NHALF_PI;
        } else {
            return D_HALF_PI;
        }
    } else {
        dec64 atan = dec64_atan(dec64_divide(y, x));
        if (x < 0) {
            if (y < 0) {
                return dec64_subtract(atan, D_HALF_PI);
            } else {
                return dec64_add(atan, D_HALF_PI);
            }
        } else {
            return atan;
        }
    }
}

dec64 dec64_cos(dec64 radians) {
    return dec64_sin(dec64_add(radians, D_HALF_PI));
}

dec64 dec64_exp(dec64 exponent) {
        dec64 result = dec64_add(DEC64_ONE, exponent);
        dec64 divisor = DEC64_TWO;
        dec64 term = exponent;
        while (1) {
            term = dec64_divide(
                dec64_multiply(term, exponent),
                divisor
            );
            dec64 progress = dec64_add(result, term);
            if (result == progress) {
                break;
            }
            result = progress;
            divisor = dec64_add(DEC64_ONE, divisor);
        }
        return result;
}

dec64 dec64_raise(dec64 coefficient, dec64 exponent) {
    if (dec64_is_zero(exponent) == DEC64_TRUE) {
        return DEC64_ONE;
    }

// Adjust for a negative exponent.

    if (exponent < 0) {
        coefficient = dec64_divide(DEC64_ONE, coefficient);
        exponent = dec64_neg(exponent);
    }
    if (dec64_is_nan(coefficient) == DEC64_TRUE) {
        return DEC64_NULL;
    }
    if (dec64_is_zero(coefficient) == DEC64_TRUE) {
        return 0;
    }

// If the exponent is an integer, then use the squaring method.

    if (exponent > 0 && dec64_exponent(exponent) == 0) {
        dec64 aux = DEC64_ONE;
        int64 n = dec64_coefficient(exponent);
        if (n <= 1) {
            return coefficient;
        }
        while (n > 1) {
            if ((n & 1) != 0) {
                aux = dec64_multiply(aux, coefficient);
            }
            coefficient = dec64_multiply(coefficient, coefficient);
            n /= 2;
        }
        return (n == 1)
            ? dec64_multiply(aux, coefficient)
            : aux;
    }

// Otherwise do it the hard way.

    return dec64_exp(dec64_multiply(
        dec64_log(coefficient),
        exponent
    ));
}

dec64 dec64_factorial(dec64 x) {
    dec64 n = dec64_normal(x);
    int64 c = dec64_coefficient(n);
    if (c >= 0 && c < FAC && dec64_exponent(n) == 0) {
        return factorials[c];
    }
    return DEC64_NULL;
}

dec64 dec64_log(dec64 x) {
    if (x <= 0 || dec64_is_nan(x) == DEC64_TRUE) {
        return DEC64_NULL;
    }
    if (dec64_is_equal(x, DEC64_ONE) == DEC64_TRUE) {
        return DEC64_ZERO;
    }
    if (dec64_is_equal(x, D_HALF) == DEC64_TRUE) {
        return dec64_new(-6931471805599453, -16);
    }
    if (x == D_E) {
        return DEC64_ONE;
    }
    dec64 y = dec64_divide(dec64_add(DEC64_NEGATIVE_ONE, x), x);
    dec64 factor = y;
    dec64 result = factor;
    dec64 divisor = DEC64_TWO;

    while (1) {
        factor = dec64_multiply(factor, y);
        dec64 progress = dec64_add(
            result,
            dec64_divide(factor, divisor)
        );
        if (result == progress || progress == DEC64_NULL) {
            break;
        }
        result = progress;
        divisor = dec64_add(DEC64_ONE, divisor);
    }
    return result;
}

/*
    The seed variables contain the random number generator's state.
    They can be set by dec64_seed.
*/

static uint64 seed_0 = D_E;
static uint64 seed_1 = D_2PI;

dec64 dec64_random() {
/*
    Return a number between 0 and 1 containing 16 randomy digits.
    It uses xorshift128+.
*/
    while (1) {
        uint64 s1 = seed_0;
        uint64 s0 = seed_1;
        s1 ^= s1 << 23;
        s1 ^= s0 ^ (s0 >> 5) ^ (s1 >> 18);
        seed_0 = s0;
        seed_1 = s1;
        uint64 mantissa = (s1 + s0) >> 10;
/*
    mantissa contains an integer between 0 and 18014398509481983.
    If it is less than or equal to 9999999999999999 then we are done.
*/
        if (mantissa <= 9999999999999999LL) {
            return dec64_new(mantissa, -16);
        }
    }
}

dec64 dec64_root(dec64 index, dec64 radicand) {
    dec64 result;
    index = dec64_normal(index);
    if (
        dec64_is_nan(radicand) == DEC64_TRUE
        || dec64_is_zero(index) == DEC64_TRUE
        || index < 0
        || dec64_exponent(index) != 0
        || (
            radicand < 0
            && (dec64_coefficient(index) & 1) == 0
        )
    ) {
        return DEC64_NULL;
    }
    if (dec64_is_zero(radicand) == DEC64_TRUE) {
        return DEC64_ZERO;
    }
    if (index == DEC64_ONE) {
        return radicand;
    }
    if (index == DEC64_TWO) {
        return dec64_sqrt(radicand);
    }
    dec64 index_minus_one = dec64_add(DEC64_NEGATIVE_ONE, index);
    result = DEC64_ONE;
    dec64 prosult = DEC64_NULL;
    while (1) {
        dec64 progress = dec64_divide(
            dec64_add(
                dec64_multiply(result, index_minus_one),
                dec64_divide(
                    radicand,
                    dec64_raise(result, index_minus_one)
                )
            ),
            index
        );
        if (progress == result) {
            return result;
        }
        if (progress == prosult) {
            return dec64_divide(dec64_add(progress, result), DEC64_TWO);
        }
        prosult = result;
        result = progress;
    }
}

void dec64_seed(uint64 part_0, uint64 part_1) {
/*
    Seed the dec64_random function. It takes any 128 bits as the seed value.
    The seed must contain at least one 1 bit.
*/
    seed_0 = part_0;
    seed_1 = part_1;
    if ((seed_0 | seed_1) == 0) {
        seed_1 = 1;
    }
}

dec64 dec64_sin(dec64 radians) {
    while (dec64_is_less(D_PI, radians) == DEC64_TRUE) {
        radians = dec64_subtract(radians, D_PI);
        radians = dec64_subtract(radians, D_PI);
    }
    while (dec64_is_less(radians, D_NPI) == DEC64_TRUE) {
        radians = dec64_add(radians, D_PI);
        radians = dec64_add(radians, D_PI);
    }
    int neg = 0;
    if (radians < 0) {
        radians = dec64_neg(radians);
        neg = -1;
    }
    if (dec64_is_less(D_HALF_PI, radians) == DEC64_TRUE) {
        radians = dec64_subtract(D_PI, radians);
    }
    dec64 result;
    if (radians == D_HALF_PI) {
        result = DEC64_ONE;
    } else {
        dec64 x2 = dec64_multiply(radians, radians);
        dec64 order = DEC64_ONE;
        dec64 term = radians;
        result = term;
        while (1) {
            term = dec64_multiply(term, x2);
            order = dec64_add(DEC64_ONE, order);
            term = dec64_divide(term, order);
            order = dec64_add(DEC64_ONE, order);
            term = dec64_divide(term, order);
            dec64 progress = dec64_subtract(result, term);

            term = dec64_multiply(term, x2);
            order = dec64_add(DEC64_ONE, order);
            term = dec64_divide(term, order);
            order = dec64_add(DEC64_ONE, order);
            term = dec64_divide(term, order);
            progress = dec64_add(progress, term);

            if (progress == result) {
                break;
            }
            result = progress;
        }
    }
    if (neg) {
        result = dec64_neg(result);
    }
    return result;
}

dec64 dec64_sqrt(dec64 radicand) {
    if (dec64_is_nan(radicand) == DEC64_FALSE && radicand >= 0) {
        if (dec64_coefficient(radicand) == 0) {
            return DEC64_ZERO;
        }
        dec64 result = radicand;
        if (dec64_is_less(DEC64_ONE, result)) {
            result = dec64_divide(result, DEC64_TWO);
        }
        while (1) {
            dec64 progress = dec64_divide(
                dec64_add(
                    result,
                    dec64_divide(radicand, result)
                ),
                DEC64_TWO
            );
            if (progress == result) {
                return result;
            }
            result = progress;
        }
        return result;
    } else {
        return DEC64_NULL;
    }
}

dec64 dec64_tan(dec64 radians) {
    return dec64_divide(
        dec64_sin(radians),
        dec64_cos(radians)
    );
}