mirror of
https://github.com/ashinn/chibi-scheme.git
synced 2025-05-18 21:29:19 +02:00
449 lines
13 KiB
C
449 lines
13 KiB
C
/* bignum.h -- header for bignum utilities */
|
|
/* Copyright (c) 2009-2020 Alex Shinn. All rights reserved. */
|
|
/* BSD-style license: http://synthcode.com/license.txt */
|
|
|
|
#ifndef SEXP_BIGNUM_H
|
|
#define SEXP_BIGNUM_H
|
|
|
|
#include "chibi/eval.h"
|
|
|
|
#if SEXP_USE_CUSTOM_LONG_LONGS
|
|
#ifdef PLAN9
|
|
#include <ape/stdint.h>
|
|
#else
|
|
#include <stdint.h>
|
|
#endif
|
|
typedef struct
|
|
{
|
|
uint64_t hi;
|
|
uint64_t lo;
|
|
} sexp_luint_t;
|
|
typedef struct
|
|
{
|
|
int64_t hi;
|
|
uint64_t lo;
|
|
} sexp_lsint_t;
|
|
#elif SEXP_64_BIT
|
|
typedef unsigned int uint128_t __attribute__((mode(TI)));
|
|
typedef int sint128_t __attribute__((mode(TI)));
|
|
typedef uint128_t sexp_luint_t;
|
|
typedef sint128_t sexp_lsint_t;
|
|
#else
|
|
typedef unsigned long long sexp_luint_t;
|
|
typedef long long sexp_lsint_t;
|
|
#endif
|
|
|
|
#if !SEXP_USE_CUSTOM_LONG_LONGS
|
|
|
|
#define sexp_lsint_fits_sint(x) ((sexp_sint_t)x == x)
|
|
#define sexp_luint_fits_uint(x) ((sexp_uint_t)x == x)
|
|
#define lsint_from_sint(v) ((sexp_lsint_t)v)
|
|
#define luint_from_uint(v) ((sexp_luint_t)v)
|
|
#define lsint_to_sint(v) ((sexp_sint_t)v)
|
|
#define luint_to_uint(v) ((sexp_uint_t)v)
|
|
#define lsint_to_sint_hi(v) ((sexp_sint_t) ((v) >> (8*sizeof(sexp_sint_t))))
|
|
#define luint_to_uint_hi(v) ((sexp_uint_t) ((v) >> (8*sizeof(sexp_uint_t))))
|
|
#define lsint_negate(v) (-((sexp_lsint_t)v))
|
|
#define luint_eq(a, b) (((sexp_luint_t)a)==((sexp_luint_t)b))
|
|
#define luint_lt(a, b) (((sexp_luint_t)a)<((sexp_luint_t)b))
|
|
#define lsint_lt_0(a) (((sexp_lsint_t)a)<0)
|
|
#define luint_shl(a, shift) (((sexp_luint_t)a)<<(shift))
|
|
#define luint_shr(a, shift) (((sexp_luint_t)a)>>(shift))
|
|
#define luint_add(a, b) (((sexp_luint_t)a)+((sexp_luint_t)b))
|
|
#define luint_add_uint(a, b) (((sexp_luint_t)a)+((sexp_uint_t)b))
|
|
#define luint_sub(a, b) (((sexp_luint_t)a)-((sexp_luint_t)b))
|
|
#define luint_mul_uint(a, b) (((sexp_luint_t)a)*((sexp_uint_t)b))
|
|
#define lsint_mul_sint(a, b) (((sexp_lsint_t)a)*((sexp_sint_t)b))
|
|
#define luint_div(a, b) (((sexp_luint_t)a)/((sexp_luint_t)b))
|
|
#define luint_div_uint(a, b) (((sexp_luint_t)a)/((sexp_luint_t)b))
|
|
#define luint_and(a, b) (((sexp_luint_t)a)&((sexp_luint_t)b))
|
|
#define luint_is_fixnum(x) (((sexp_luint_t)x)<=SEXP_MAX_FIXNUM)
|
|
#define lsint_is_fixnum(x) ((SEXP_MIN_FIXNUM <= ((sexp_lsint_t)x)) && (((sexp_lsint_t)x) <= SEXP_MAX_FIXNUM))
|
|
|
|
#else
|
|
|
|
static inline int lsint_lt_0(sexp_lsint_t a) {
|
|
return a.hi < 0;
|
|
}
|
|
|
|
static inline int sexp_lsint_fits_sint(sexp_lsint_t x) {
|
|
return x.hi == (((int64_t)x.lo)>>63) && ((sexp_sint_t)x.lo == x.lo);
|
|
}
|
|
static inline int sexp_luint_fits_uint(sexp_luint_t x) {
|
|
return x.hi == 0 && ((sexp_uint_t)x.lo == x.lo);
|
|
}
|
|
|
|
static inline sexp_luint_t luint_from_lsint(sexp_lsint_t v) {
|
|
sexp_luint_t result;
|
|
result.hi = v.hi;
|
|
result.lo = v.lo;
|
|
return result;
|
|
}
|
|
|
|
static inline sexp_lsint_t lsint_from_luint(sexp_luint_t v) {
|
|
sexp_lsint_t result;
|
|
result.hi = v.hi;
|
|
result.lo = v.lo;
|
|
return result;
|
|
}
|
|
|
|
static inline sexp_lsint_t lsint_from_sint(sexp_sint_t v) {
|
|
sexp_lsint_t result;
|
|
result.hi = v >> 63;
|
|
result.lo = v;
|
|
return result;
|
|
}
|
|
|
|
static inline sexp_luint_t luint_from_uint(sexp_uint_t v) {
|
|
sexp_luint_t result;
|
|
result.hi = 0;
|
|
result.lo = v;
|
|
return result;
|
|
}
|
|
|
|
static inline sexp_sint_t lsint_to_sint(sexp_lsint_t v) {
|
|
return v.lo;
|
|
}
|
|
|
|
static inline sexp_uint_t luint_to_uint(sexp_luint_t v) {
|
|
return v.lo;
|
|
}
|
|
|
|
static inline sexp_sint_t lsint_to_sint_hi(sexp_lsint_t v) {
|
|
#if SEXP_64_BIT
|
|
return v.hi;
|
|
#else
|
|
return v.lo >> 32;
|
|
#endif
|
|
}
|
|
|
|
static inline sexp_uint_t luint_to_uint_hi(sexp_luint_t v) {
|
|
#if SEXP_64_BIT
|
|
return v.hi;
|
|
#else
|
|
return v.lo >> 32;
|
|
#endif
|
|
}
|
|
|
|
static inline sexp_lsint_t lsint_negate(sexp_lsint_t v) {
|
|
sexp_luint_t a;
|
|
a.hi = ~v.hi;
|
|
a.lo = ~v.lo;
|
|
|
|
uint64_t aLoLo = a.lo & 0xFFFFFFFF;
|
|
uint64_t aLoHi = a.lo >> 32;
|
|
uint64_t aHiLo = a.hi & 0xFFFFFFFF;
|
|
uint64_t aHiHi = a.hi >> 32;
|
|
|
|
uint64_t carry;
|
|
uint64_t sumLoLo = aLoLo + 1;
|
|
carry = sumLoLo >> 32;
|
|
uint64_t resultLoLo = sumLoLo & 0xFFFFFFFF;
|
|
|
|
uint64_t sumLoHi = aLoHi + carry;
|
|
uint64_t resultLoHi = sumLoHi & 0xFFFFFFFF;
|
|
carry = sumLoHi >> 32;
|
|
|
|
uint64_t sumHiLo = aHiLo + carry;
|
|
uint64_t resultHiLo = sumHiLo & 0xFFFFFFFF;
|
|
carry = sumHiLo >> 32;
|
|
|
|
uint64_t sumHiHi = aHiHi + carry;
|
|
uint64_t resultHiHi = sumHiHi & 0xFFFFFFFF;
|
|
/* carry = sumHiHi >> 32; */
|
|
|
|
sexp_lsint_t result;
|
|
result.hi = (resultHiHi << 32) | resultHiLo;
|
|
result.lo = (resultLoHi << 32) | resultLoLo;
|
|
return result;
|
|
}
|
|
|
|
static inline int luint_eq(sexp_luint_t a, sexp_luint_t b) {
|
|
return (a.hi == b.hi) && (a.lo == b.lo);
|
|
}
|
|
|
|
static inline int luint_lt(sexp_luint_t a, sexp_luint_t b) {
|
|
if (a.hi < b.hi)
|
|
return 1;
|
|
else if (a.hi > b.hi)
|
|
return 0;
|
|
else
|
|
return a.lo < b.lo;
|
|
}
|
|
|
|
static inline sexp_luint_t luint_shl(sexp_luint_t v, size_t shift) {
|
|
if (shift == 0)
|
|
return v;
|
|
sexp_luint_t result;
|
|
if (shift >= 64) {
|
|
result.hi = v.lo << (shift - 64);
|
|
result.lo = 0;
|
|
} else {
|
|
result.hi = (v.hi << shift) | (v.lo >> (64-shift));
|
|
result.lo = v.lo << shift;
|
|
}
|
|
return result;
|
|
}
|
|
|
|
static inline sexp_luint_t luint_shr(sexp_luint_t v, size_t shift) {
|
|
if (shift == 0)
|
|
return v;
|
|
sexp_luint_t result;
|
|
if (shift >= 64) {
|
|
result.hi = 0;
|
|
result.lo = v.hi >> (shift - 64);
|
|
} else {
|
|
result.hi = v.hi >> shift;
|
|
result.lo = (v.lo >> shift) | (v.hi << (64-shift));
|
|
}
|
|
return result;
|
|
}
|
|
|
|
static inline sexp_luint_t luint_add(sexp_luint_t a, sexp_luint_t b) {
|
|
uint64_t aLoLo = a.lo & 0xFFFFFFFF;
|
|
uint64_t aLoHi = a.lo >> 32;
|
|
uint64_t aHiLo = a.hi & 0xFFFFFFFF;
|
|
uint64_t aHiHi = a.hi >> 32;
|
|
uint64_t bLoLo = b.lo & 0xFFFFFFFF;
|
|
uint64_t bLoHi = b.lo >> 32;
|
|
uint64_t bHiLo = b.hi & 0xFFFFFFFF;
|
|
uint64_t bHiHi = b.hi >> 32;
|
|
|
|
uint64_t carry;
|
|
uint64_t sumLoLo = (aLoLo + bLoLo);
|
|
carry = sumLoLo >> 32;
|
|
uint64_t resultLoLo = sumLoLo & 0xFFFFFFFF;
|
|
|
|
uint64_t sumLoHi = (aLoHi + bLoHi) + carry;
|
|
uint64_t resultLoHi = sumLoHi & 0xFFFFFFFF;
|
|
carry = sumLoHi >> 32;
|
|
|
|
uint64_t sumHiLo = (aHiLo + bHiLo) + carry;
|
|
uint64_t resultHiLo = sumHiLo & 0xFFFFFFFF;
|
|
carry = sumHiLo >> 32;
|
|
|
|
uint64_t sumHiHi = (aHiHi + bHiHi) + carry;
|
|
uint64_t resultHiHi = sumHiHi & 0xFFFFFFFF;
|
|
/* carry = sumHiHi >> 32; */
|
|
|
|
sexp_luint_t result;
|
|
result.hi = (resultHiHi << 32) | resultHiLo;
|
|
result.lo = (resultLoHi << 32) | resultLoLo;
|
|
return result;
|
|
}
|
|
|
|
static inline sexp_luint_t luint_add_uint(sexp_luint_t a, sexp_uint_t b) {
|
|
uint64_t aLoLo = a.lo & 0xFFFFFFFF;
|
|
uint64_t aLoHi = a.lo >> 32;
|
|
uint64_t aHiLo = a.hi & 0xFFFFFFFF;
|
|
uint64_t aHiHi = a.hi >> 32;
|
|
uint64_t bLoLo = b & 0xFFFFFFFF;
|
|
uint64_t bLoHi = b >> 32;
|
|
|
|
uint64_t carry;
|
|
uint64_t sumLoLo = (aLoLo + bLoLo);
|
|
carry = sumLoLo >> 32;
|
|
uint64_t resultLoLo = sumLoLo & 0xFFFFFFFF;
|
|
|
|
uint64_t sumLoHi = (aLoHi + bLoHi) + carry;
|
|
uint64_t resultLoHi = sumLoHi & 0xFFFFFFFF;
|
|
carry = sumLoHi >> 32;
|
|
|
|
uint64_t sumHiLo = aHiLo + carry;
|
|
uint64_t resultHiLo = sumHiLo & 0xFFFFFFFF;
|
|
carry = sumHiLo >> 32;
|
|
|
|
uint64_t sumHiHi = aHiHi + carry;
|
|
uint64_t resultHiHi = sumHiHi & 0xFFFFFFFF;
|
|
/* carry = sumHiHi >> 32; */
|
|
|
|
sexp_luint_t result;
|
|
result.hi = (resultHiHi << 32) | resultHiLo;
|
|
result.lo = (resultLoHi << 32) | resultLoLo;
|
|
return result;
|
|
}
|
|
|
|
static inline sexp_luint_t luint_sub(sexp_luint_t a, sexp_luint_t b) {
|
|
sexp_luint_t negB;
|
|
negB.hi = ~b.hi;
|
|
negB.lo = ~b.lo;
|
|
return luint_add(a, luint_add_uint(negB, 1));
|
|
}
|
|
|
|
static inline sexp_luint_t luint_mul_uint(sexp_luint_t a, sexp_uint_t b) {
|
|
uint64_t aLoLo = a.lo & 0xFFFFFFFF;
|
|
uint64_t aLoHi = a.lo >> 32;
|
|
uint64_t aHiLo = a.hi & 0xFFFFFFFF;
|
|
uint64_t aHiHi = a.hi >> 32;
|
|
|
|
uint64_t bLo = b & 0xFFFFFFFF;
|
|
uint64_t bHi = b >> 32;
|
|
|
|
sexp_luint_t resultBLo, resultBHi;
|
|
{
|
|
sexp_luint_t prodLoLo;
|
|
prodLoLo.hi = 0;
|
|
prodLoLo.lo = aLoLo * bLo;
|
|
|
|
sexp_luint_t prodLoHi;
|
|
prodLoHi.hi = (aLoHi * bLo) >> 32;
|
|
prodLoHi.lo = (aLoHi * bLo) << 32;
|
|
|
|
sexp_luint_t prodHiLo;
|
|
prodHiLo.hi = aHiLo * bLo;
|
|
prodHiLo.lo = 0;
|
|
|
|
sexp_luint_t prodHiHi;
|
|
prodHiHi.hi = (aHiHi * bLo) << 32;
|
|
prodHiHi.lo = 0;
|
|
|
|
resultBLo = luint_add(luint_add(luint_add(prodLoLo, prodLoHi), prodHiLo), prodHiHi);
|
|
}
|
|
{
|
|
sexp_luint_t prodLoLo;
|
|
prodLoLo.hi = 0;
|
|
prodLoLo.lo = aLoLo * bHi;
|
|
|
|
sexp_luint_t prodLoHi;
|
|
prodLoHi.hi = (aLoHi * bHi) >> 32;
|
|
prodLoHi.lo = (aLoHi * bHi) << 32;
|
|
|
|
sexp_luint_t prodHiLo;
|
|
prodHiLo.hi = aHiLo * bHi;
|
|
prodHiLo.lo = 0;
|
|
|
|
sexp_luint_t prodHiHi;
|
|
prodHiHi.hi = (aHiHi * bHi) << 32;
|
|
prodHiHi.lo = 0;
|
|
|
|
resultBHi = luint_add(luint_add(luint_add(prodLoLo, prodLoHi), prodHiLo), prodHiHi);
|
|
}
|
|
|
|
sexp_luint_t result = luint_add(resultBLo, luint_shl(resultBHi, 32));
|
|
|
|
return result;
|
|
}
|
|
|
|
static inline sexp_lsint_t lsint_mul_sint(sexp_lsint_t a, sexp_sint_t b) {
|
|
if (lsint_lt_0(a)) {
|
|
sexp_luint_t minusA = luint_from_lsint(lsint_negate(a));
|
|
if (b < 0)
|
|
return lsint_from_luint(luint_mul_uint(minusA, (sexp_uint_t)-b));
|
|
else
|
|
return lsint_negate(lsint_from_luint(luint_mul_uint(minusA, (sexp_uint_t)b)));
|
|
} else {
|
|
if (b < 0)
|
|
return lsint_negate(lsint_from_luint(luint_mul_uint(luint_from_lsint(a), (sexp_uint_t)-b)));
|
|
else
|
|
return lsint_from_luint(luint_mul_uint(luint_from_lsint(a), (sexp_uint_t)b));
|
|
}
|
|
}
|
|
|
|
static inline sexp_luint_t luint_div(sexp_luint_t a, sexp_luint_t b) {
|
|
if (luint_lt(a, b))
|
|
return luint_from_uint(0);
|
|
else if (luint_eq(a, b))
|
|
return luint_from_uint(1);
|
|
|
|
sexp_luint_t quotient = luint_from_uint(0);
|
|
sexp_luint_t remainder = luint_from_uint(0);
|
|
|
|
for (int i = 0; i < 128; i++) {
|
|
quotient = luint_shl(quotient, 1);
|
|
|
|
remainder = luint_shl(remainder, 1);
|
|
remainder.lo |= (a.hi >> 63) & 1;
|
|
a = luint_shl(a, 1);
|
|
|
|
if (!(luint_lt(remainder, b))) {
|
|
remainder = luint_sub(remainder, b);
|
|
quotient.lo |= 1;
|
|
}
|
|
}
|
|
|
|
return quotient;
|
|
}
|
|
|
|
static inline sexp_luint_t luint_div_uint(sexp_luint_t a, sexp_uint_t b) {
|
|
return luint_div(a, luint_from_uint(b));
|
|
}
|
|
|
|
static inline sexp_luint_t luint_and(sexp_luint_t a, sexp_luint_t b) {
|
|
sexp_luint_t result;
|
|
result.hi = a.hi & b.hi;
|
|
result.lo = a.lo & b.lo;
|
|
return result;
|
|
}
|
|
|
|
static inline int luint_is_fixnum(sexp_luint_t x) {
|
|
return (x.hi == 0) && (x.lo <= SEXP_MAX_FIXNUM);
|
|
}
|
|
|
|
static inline int lsint_is_fixnum(sexp_lsint_t x) {
|
|
if (x.hi > 0)
|
|
return 0;
|
|
else if (x.hi == 0)
|
|
return x.lo <= SEXP_MAX_FIXNUM;
|
|
else if (x.hi == -1)
|
|
return SEXP_MIN_FIXNUM <= x.lo;
|
|
else return 0;
|
|
}
|
|
|
|
#endif
|
|
|
|
SEXP_API sexp_sint_t sexp_bignum_compare (sexp a, sexp b);
|
|
SEXP_API sexp sexp_compare (sexp ctx, sexp a, sexp b);
|
|
SEXP_API sexp sexp_make_bignum (sexp ctx, sexp_uint_t len);
|
|
SEXP_API sexp sexp_copy_bignum (sexp ctx, sexp dst, sexp a, sexp_uint_t len);
|
|
SEXP_API sexp sexp_bignum_normalize (sexp a);
|
|
SEXP_API sexp_uint_t sexp_bignum_hi (sexp a);
|
|
SEXP_API sexp sexp_fixnum_to_bignum (sexp ctx, sexp a);
|
|
SEXP_API double sexp_bignum_to_double (sexp a);
|
|
SEXP_API sexp sexp_double_to_bignum (sexp ctx, double f);
|
|
SEXP_API double sexp_to_double (sexp ctx, sexp x);
|
|
SEXP_API sexp sexp_bignum_fxadd (sexp ctx, sexp a, sexp_uint_t b);
|
|
SEXP_API sexp sexp_bignum_fxsub (sexp ctx, sexp a, sexp_uint_t b);
|
|
SEXP_API sexp sexp_bignum_fxmul (sexp ctx, sexp d, sexp a, sexp_uint_t b, int offset);
|
|
SEXP_API sexp_uint_t sexp_bignum_fxdiv (sexp ctx, sexp a, sexp_uint_t b, int offset);
|
|
SEXP_API sexp sexp_bignum_add (sexp ctx, sexp dst, sexp a, sexp b);
|
|
SEXP_API sexp sexp_bignum_sub (sexp ctx, sexp dst, sexp a, sexp b);
|
|
SEXP_API sexp sexp_bignum_mul (sexp ctx, sexp dst, sexp a, sexp b);
|
|
SEXP_API sexp sexp_bignum_div (sexp ctx, sexp dst, sexp a, sexp b);
|
|
SEXP_API sexp sexp_bignum_expt (sexp ctx, sexp n, sexp e);
|
|
SEXP_API sexp sexp_bignum_sqrt (sexp ctx, sexp a, sexp* rem);
|
|
SEXP_API sexp sexp_add (sexp ctx, sexp a, sexp b);
|
|
SEXP_API sexp sexp_sub (sexp ctx, sexp a, sexp b);
|
|
SEXP_API sexp sexp_mul (sexp ctx, sexp a, sexp b);
|
|
SEXP_API sexp sexp_div (sexp ctx, sexp a, sexp b);
|
|
SEXP_API sexp sexp_quotient (sexp ctx, sexp a, sexp b);
|
|
SEXP_API sexp sexp_remainder (sexp ctx, sexp a, sexp b);
|
|
#if SEXP_USE_RATIOS
|
|
SEXP_API sexp sexp_double_to_ratio (sexp ctx, double f);
|
|
SEXP_API sexp sexp_double_to_ratio_2 (sexp ctx, double f);
|
|
SEXP_API double sexp_ratio_to_double (sexp ctx, sexp rat);
|
|
SEXP_API sexp sexp_make_ratio (sexp ctx, sexp num, sexp den);
|
|
SEXP_API sexp sexp_ratio_normalize (sexp ctx, sexp rat, sexp in);
|
|
SEXP_API sexp sexp_ratio_round (sexp ctx, sexp a);
|
|
SEXP_API sexp sexp_ratio_trunc (sexp ctx, sexp a);
|
|
SEXP_API sexp sexp_ratio_floor (sexp ctx, sexp a);
|
|
SEXP_API sexp sexp_ratio_ceiling (sexp ctx, sexp a);
|
|
SEXP_API sexp sexp_ratio_compare (sexp ctx, sexp a, sexp b);
|
|
#endif
|
|
#if SEXP_USE_COMPLEX
|
|
SEXP_API sexp sexp_make_complex (sexp ctx, sexp real, sexp image);
|
|
SEXP_API sexp sexp_complex_normalize (sexp real);
|
|
SEXP_API sexp sexp_complex_math_error (sexp ctx, sexp z);
|
|
SEXP_API sexp sexp_complex_sqrt (sexp ctx, sexp z);
|
|
SEXP_API sexp sexp_complex_exp (sexp ctx, sexp z);
|
|
SEXP_API sexp sexp_complex_expt (sexp ctx, sexp a, sexp b);
|
|
SEXP_API sexp sexp_complex_log (sexp ctx, sexp z);
|
|
SEXP_API sexp sexp_complex_sin (sexp ctx, sexp z);
|
|
SEXP_API sexp sexp_complex_cos (sexp ctx, sexp z);
|
|
SEXP_API sexp sexp_complex_tan (sexp ctx, sexp z);
|
|
SEXP_API sexp sexp_complex_asin (sexp ctx, sexp z);
|
|
SEXP_API sexp sexp_complex_acos (sexp ctx, sexp z);
|
|
SEXP_API sexp sexp_complex_atan (sexp ctx, sexp z);
|
|
#endif
|
|
|
|
#endif /* ! SEXP_BIGNUM_H */
|
|
|