|
|
|
/*
|
|
|
|
r128.h: 128-bit (64.64) signed fixed-point arithmetic. Version 1.6.0
|
|
|
|
|
|
|
|
Modified by tcsullivan to be used as (4.124) fixed-point.
|
|
|
|
Unused operators also removed.
|
|
|
|
|
|
|
|
COMPILATION
|
|
|
|
-----------
|
|
|
|
Drop this header file somewhere in your project and include it wherever it is
|
|
|
|
needed. There is no separate .c file for this library. To get the code, in ONE
|
|
|
|
file in your project, put:
|
|
|
|
|
|
|
|
#define R128_IMPLEMENTATION
|
|
|
|
|
|
|
|
before you include this file. You may also provide a definition for R128_ASSERT
|
|
|
|
to force the library to use a custom assert macro.
|
|
|
|
|
|
|
|
COMPILER/LIBRARY SUPPORT
|
|
|
|
------------------------
|
|
|
|
This library requires a C89 compiler with support for 64-bit integers. If your
|
|
|
|
compiler does not support the long long data type, the R128_U64, etc. macros
|
|
|
|
must be set appropriately. On x86 and x64 targets, Intel intrinsics are used
|
|
|
|
for speed. If your compiler does not support these intrinsics, you can add
|
|
|
|
#define R128_STDC_ONLY
|
|
|
|
in your implementation file before including r128.h.
|
|
|
|
|
|
|
|
The only C runtime library functionality used by this library is <assert.h>.
|
|
|
|
This can be avoided by defining an R128_ASSERT macro in your implementation
|
|
|
|
file. Since this library uses 64-bit arithmetic, this may implicitly add a
|
|
|
|
runtime library dependency on 32-bit platforms.
|
|
|
|
|
|
|
|
C++ SUPPORT
|
|
|
|
-----------
|
|
|
|
Operator overloads are supplied for C++ files that include this file. Since all
|
|
|
|
C++ functions are declared inline (or static inline), the R128_IMPLEMENTATION
|
|
|
|
file can be either C++ or C.
|
|
|
|
|
|
|
|
LICENSE
|
|
|
|
-------
|
|
|
|
This is free and unencumbered software released into the public domain.
|
|
|
|
|
|
|
|
Anyone is free to copy, modify, publish, use, compile, sell, or
|
|
|
|
distribute this software, either in source code form or as a compiled
|
|
|
|
binary, for any purpose, commercial or non-commercial, and by any
|
|
|
|
means.
|
|
|
|
|
|
|
|
In jurisdictions that recognize copyright laws, the author or authors
|
|
|
|
of this software dedicate any and all copyright interest in the
|
|
|
|
software to the public domain. We make this dedication for the benefit
|
|
|
|
of the public at large and to the detriment of our heirs and
|
|
|
|
successors. We intend this dedication to be an overt act of
|
|
|
|
relinquishment in perpetuity of all present and future rights to this
|
|
|
|
software under copyright law.
|
|
|
|
|
|
|
|
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 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.
|
|
|
|
*/
|
|
|
|
|
|
|
|
#ifndef H_R128_H
|
|
|
|
#define H_R128_H
|
|
|
|
|
|
|
|
#include <stddef.h>
|
|
|
|
#include <stdint.h>
|
|
|
|
|
|
|
|
#define R128_S32 int32_t
|
|
|
|
#define R128_U32 uint32_t
|
|
|
|
#define R128_S64 int64_t
|
|
|
|
#define R128_U64 uint64_t
|
|
|
|
|
|
|
|
#ifdef __cplusplus
|
|
|
|
extern "C" {
|
|
|
|
#endif
|
|
|
|
|
|
|
|
typedef struct R128 {
|
|
|
|
R128_U64 lo;
|
|
|
|
R128_U64 hi;
|
|
|
|
|
|
|
|
#ifdef __cplusplus
|
|
|
|
R128();
|
|
|
|
R128(double);
|
|
|
|
R128(int);
|
|
|
|
R128(R128_S64);
|
|
|
|
R128(R128_U64 low, R128_U64 high);
|
|
|
|
|
|
|
|
operator double() const;
|
|
|
|
|
|
|
|
R128 operator-() const;
|
|
|
|
R128 &operator+=(const R128 &rhs);
|
|
|
|
R128 &operator-=(const R128 &rhs);
|
|
|
|
R128 &operator*=(const R128 &rhs);
|
|
|
|
#endif //__cplusplus
|
|
|
|
} __attribute__ ((packed)) R128;
|
|
|
|
|
|
|
|
// Type conversion
|
|
|
|
extern void r128FromInt(R128 *dst, R128_S64 v);
|
|
|
|
extern void r128FromFloat(R128 *dst, double v);
|
|
|
|
extern double r128ToFloat(const R128 *v);
|
|
|
|
|
|
|
|
// Copy
|
|
|
|
extern void r128Copy(R128 *dst, const R128 *src);
|
|
|
|
|
|
|
|
// Sign manipulation
|
|
|
|
extern void r128Neg(R128 *dst, const R128 *v); // -v
|
|
|
|
|
|
|
|
// Bitwise operations
|
|
|
|
extern void r128Shl(R128 *dst, const R128 *src, int amount); // shift left by amount mod 128
|
|
|
|
extern void r128Shr(R128 *dst, const R128 *src, int amount); // shift right logical by amount mod 128
|
|
|
|
|
|
|
|
// Arithmetic
|
|
|
|
extern void r128Add(R128 *dst, const R128 *a, const R128 *b); // a + b
|
|
|
|
extern void r128Sub(R128 *dst, const R128 *a, const R128 *b); // a - b
|
|
|
|
extern void r128Mul(R128 *dst, const R128 *a, const R128 *b); // a * b
|
|
|
|
|
|
|
|
// Comparison
|
|
|
|
extern int r128Cmp(const R128 *a, const R128 *b); // sign of a-b
|
|
|
|
extern int r128IsNeg(const R128 *v); // quick check for < 0
|
|
|
|
|
|
|
|
// Constants
|
|
|
|
extern const R128 R128_min; // minimum (most negative) value
|
|
|
|
extern const R128 R128_max; // maximum (most positive) value
|
|
|
|
extern const R128 R128_smallest; // smallest positive value
|
|
|
|
extern const R128 R128_zero; // zero
|
|
|
|
extern const R128 R128_one; // 1.0
|
|
|
|
|
|
|
|
extern char R128_decimal; // decimal point character used by r128From/ToString. defaults to '.'
|
|
|
|
|
|
|
|
#ifdef __cplusplus
|
|
|
|
}
|
|
|
|
|
|
|
|
inline R128::R128() {}
|
|
|
|
|
|
|
|
inline R128::R128(double v)
|
|
|
|
{
|
|
|
|
r128FromFloat(this, v);
|
|
|
|
}
|
|
|
|
|
|
|
|
inline R128::R128(int v)
|
|
|
|
{
|
|
|
|
r128FromInt(this, v);
|
|
|
|
}
|
|
|
|
|
|
|
|
inline R128::R128(R128_S64 v)
|
|
|
|
{
|
|
|
|
r128FromInt(this, v);
|
|
|
|
}
|
|
|
|
|
|
|
|
inline R128::R128(R128_U64 low, R128_U64 high)
|
|
|
|
{
|
|
|
|
lo = low;
|
|
|
|
hi = high;
|
|
|
|
}
|
|
|
|
|
|
|
|
inline R128::operator double() const
|
|
|
|
{
|
|
|
|
return r128ToFloat(this);
|
|
|
|
}
|
|
|
|
|
|
|
|
inline R128 R128::operator-() const
|
|
|
|
{
|
|
|
|
R128 r;
|
|
|
|
r128Neg(&r, this);
|
|
|
|
return r;
|
|
|
|
}
|
|
|
|
|
|
|
|
inline R128 &R128::operator+=(const R128 &rhs)
|
|
|
|
{
|
|
|
|
r128Add(this, this, &rhs);
|
|
|
|
return *this;
|
|
|
|
}
|
|
|
|
|
|
|
|
inline R128 &R128::operator-=(const R128 &rhs)
|
|
|
|
{
|
|
|
|
r128Sub(this, this, &rhs);
|
|
|
|
return *this;
|
|
|
|
}
|
|
|
|
|
|
|
|
inline R128 &R128::operator*=(const R128 &rhs)
|
|
|
|
{
|
|
|
|
r128Mul(this, this, &rhs);
|
|
|
|
return *this;
|
|
|
|
}
|
|
|
|
|
|
|
|
static inline R128 operator+(const R128 &lhs, const R128 &rhs)
|
|
|
|
{
|
|
|
|
R128 r(lhs);
|
|
|
|
return r += rhs;
|
|
|
|
}
|
|
|
|
|
|
|
|
static inline R128 operator-(const R128 &lhs, const R128 &rhs)
|
|
|
|
{
|
|
|
|
R128 r(lhs);
|
|
|
|
return r -= rhs;
|
|
|
|
}
|
|
|
|
|
|
|
|
static inline R128 operator*(const R128 &lhs, const R128 &rhs)
|
|
|
|
{
|
|
|
|
R128 r(lhs);
|
|
|
|
return r *= rhs;
|
|
|
|
}
|
|
|
|
|
|
|
|
static inline bool operator<(const R128 &lhs, const R128 &rhs)
|
|
|
|
{
|
|
|
|
return r128Cmp(&lhs, &rhs) < 0;
|
|
|
|
}
|
|
|
|
|
|
|
|
static inline bool operator>(const R128 &lhs, const R128 &rhs)
|
|
|
|
{
|
|
|
|
return r128Cmp(&lhs, &rhs) > 0;
|
|
|
|
}
|
|
|
|
|
|
|
|
static inline bool operator<=(const R128 &lhs, const R128 &rhs)
|
|
|
|
{
|
|
|
|
return r128Cmp(&lhs, &rhs) <= 0;
|
|
|
|
}
|
|
|
|
|
|
|
|
static inline bool operator>=(const R128 &lhs, const R128 &rhs)
|
|
|
|
{
|
|
|
|
return r128Cmp(&lhs, &rhs) >= 0;
|
|
|
|
}
|
|
|
|
|
|
|
|
static inline bool operator==(const R128 &lhs, const R128 &rhs)
|
|
|
|
{
|
|
|
|
return lhs.lo == rhs.lo && lhs.hi == rhs.hi;
|
|
|
|
}
|
|
|
|
|
|
|
|
static inline bool operator!=(const R128 &lhs, const R128 &rhs)
|
|
|
|
{
|
|
|
|
return lhs.lo != rhs.lo || lhs.hi != rhs.hi;
|
|
|
|
}
|
|
|
|
|
|
|
|
#endif //__cplusplus
|
|
|
|
#endif //H_R128_H
|
|
|
|
|
|
|
|
#ifdef R128_IMPLEMENTATION
|
|
|
|
|
|
|
|
#define R128_SET2(x, l, h) { (x)->lo = (R128_U64)(l); (x)->hi = (R128_U64)(h); }
|
|
|
|
|
|
|
|
#include <stdlib.h> // for NULL
|
|
|
|
|
|
|
|
const R128 R128_min = { 0, 0x8000000000000000 };
|
|
|
|
const R128 R128_max = { 0xffffffffffffffff, 0x7fffffffffffffff };
|
|
|
|
const R128 R128_smallest = { 1, 0 };
|
|
|
|
const R128 R128_zero = { 0, 0 };
|
|
|
|
const R128 R128_one = { 0, 1 };
|
|
|
|
char R128_decimal = '.';
|
|
|
|
|
|
|
|
static void r128__neg(R128 *dst, const R128 *src)
|
|
|
|
{
|
|
|
|
if (src->lo) {
|
|
|
|
dst->lo = ~src->lo + 1;
|
|
|
|
dst->hi = ~src->hi;
|
|
|
|
} else {
|
|
|
|
dst->lo = 0;
|
|
|
|
dst->hi = ~src->hi + 1;
|
|
|
|
}
|
|
|
|
}
|
|
|
|
|
|
|
|
// 64*64->128
|
|
|
|
static void r128__umul128(R128 *dst, R128_U64 a, R128_U64 b)
|
|
|
|
{
|
|
|
|
R128_U32 alo = (R128_U32)a;
|
|
|
|
R128_U32 ahi = (R128_U32)(a >> 32);
|
|
|
|
R128_U32 blo = (R128_U32)b;
|
|
|
|
R128_U32 bhi = (R128_U32)(b >> 32);
|
|
|
|
R128_U64 p0, p1, p2, p3;
|
|
|
|
|
|
|
|
p0 = alo * (uint64_t)blo;
|
|
|
|
p1 = alo * (uint64_t)bhi;
|
|
|
|
p2 = ahi * (uint64_t)blo;
|
|
|
|
p3 = ahi * (uint64_t)bhi;
|
|
|
|
|
|
|
|
{
|
|
|
|
R128_U64 carry, lo, hi;
|
|
|
|
carry = ((R128_U64)(R128_U32)p1 + (R128_U64)(R128_U32)p2 + (p0 >> 32)) >> 32;
|
|
|
|
|
|
|
|
lo = p0 + ((p1 + p2) << 32);
|
|
|
|
hi = p3 + ((R128_U32)(p1 >> 32) + (R128_U32)(p2 >> 32)) + carry;
|
|
|
|
|
|
|
|
R128_SET2(dst, lo, hi);
|
|
|
|
}
|
|
|
|
}
|
|
|
|
|
|
|
|
static int r128__ucmp(const R128 *a, const R128 *b)
|
|
|
|
{
|
|
|
|
if (a->hi != b->hi) {
|
|
|
|
if (a->hi > b->hi) {
|
|
|
|
return 1;
|
|
|
|
} else {
|
|
|
|
return -1;
|
|
|
|
}
|
|
|
|
} else {
|
|
|
|
if (a->lo == b->lo) {
|
|
|
|
return 0;
|
|
|
|
} else if (a->lo > b->lo) {
|
|
|
|
return 1;
|
|
|
|
} else {
|
|
|
|
return -1;
|
|
|
|
}
|
|
|
|
}
|
|
|
|
}
|
|
|
|
|
|
|
|
static void r128__umul(R128 *dst, const R128 *a, const R128 *b)
|
|
|
|
{
|
|
|
|
R128 albl, ahbl, albh, ahbh;
|
|
|
|
|
|
|
|
r128__umul128(&albl, a->lo, b->lo);
|
|
|
|
r128__umul128(&ahbl, a->hi, b->lo);
|
|
|
|
r128__umul128(&albh, a->lo, b->hi);
|
|
|
|
r128__umul128(&ahbh, a->hi, b->hi);
|
|
|
|
|
|
|
|
R128 sum;
|
|
|
|
|
|
|
|
r128Shr(&ahbh, &ahbh, 60);
|
|
|
|
sum.lo = 0;
|
|
|
|
sum.hi = ahbh.lo;
|
|
|
|
|
|
|
|
albl.lo = albl.hi;
|
|
|
|
albl.hi = 0;
|
|
|
|
r128Add(&sum, &sum, &albl);
|
|
|
|
r128Add(&sum, &sum, &ahbl);
|
|
|
|
r128Add(&sum, &sum, &albh);
|
|
|
|
|
|
|
|
R128_SET2(dst, sum.lo, sum.hi);
|
|
|
|
}
|
|
|
|
|
|
|
|
void r128FromInt(R128 *dst, R128_S64 v)
|
|
|
|
{
|
|
|
|
dst->lo = 0;
|
|
|
|
dst->hi = (R128_U64)v << 60;
|
|
|
|
}
|
|
|
|
|
|
|
|
void r128FromFloat(R128 *dst, double v)
|
|
|
|
{
|
|
|
|
if (v < -9223372036854775808.0) {
|
|
|
|
r128Copy(dst, &R128_min);
|
|
|
|
} else if (v >= 9223372036854775808.0) {
|
|
|
|
r128Copy(dst, &R128_max);
|
|
|
|
} else {
|
|
|
|
R128 r;
|
|
|
|
int sign = 0;
|
|
|
|
|
|
|
|
if (v < 0) {
|
|
|
|
v = -v;
|
|
|
|
sign = 1;
|
|
|
|
}
|
|
|
|
|
|
|
|
r.hi = ((R128_U64)(R128_S32)v) << 60;
|
|
|
|
v -= (R128_S32)v;
|
|
|
|
v *= (double)((R128_U64)1 << 60);
|
|
|
|
r.hi = (r.hi & 0xF000000000000000) + (R128_U64)v;
|
|
|
|
r.lo = (R128_U64)(v * 18446744073709551616.0);
|
|
|
|
|
|
|
|
if (sign) {
|
|
|
|
r128__neg(&r, &r);
|
|
|
|
}
|
|
|
|
|
|
|
|
r128Copy(dst, &r);
|
|
|
|
}
|
|
|
|
}
|
|
|
|
|
|
|
|
R128_S64 r128ToInt(const R128 *v)
|
|
|
|
{
|
|
|
|
if ((R128_S64)v->hi < 0) {
|
|
|
|
return (R128_S64)v->hi + (v->lo != 0);
|
|
|
|
} else {
|
|
|
|
return (R128_S64)v->hi;
|
|
|
|
}
|
|
|
|
}
|
|
|
|
|
|
|
|
double r128ToFloat(const R128 *v)
|
|
|
|
{
|
|
|
|
R128 tmp;
|
|
|
|
int sign = 0;
|
|
|
|
double d;
|
|
|
|
|
|
|
|
R128_SET2(&tmp, v->lo, v->hi);
|
|
|
|
if (r128IsNeg(&tmp)) {
|
|
|
|
r128__neg(&tmp, &tmp);
|
|
|
|
sign = 1;
|
|
|
|
}
|
|
|
|
|
|
|
|
d = (tmp.hi >> 60);
|
|
|
|
d += (tmp.hi & 0xFFFFFFFFFFFFFF) / (double)((uint64_t)1 << 60);
|
|
|
|
d += tmp.lo / (double)((uint64_t)1 << 60) / 18446744073709551616.0;
|
|
|
|
|
|
|
|
if (sign) {
|
|
|
|
d = -d;
|
|
|
|
}
|
|
|
|
|
|
|
|
return d;
|
|
|
|
}
|
|
|
|
|
|
|
|
void r128Copy(R128 *dst, const R128 *src)
|
|
|
|
{
|
|
|
|
dst->lo = src->lo;
|
|
|
|
dst->hi = src->hi;
|
|
|
|
}
|
|
|
|
|
|
|
|
void r128Neg(R128 *dst, const R128 *v)
|
|
|
|
{
|
|
|
|
r128__neg(dst, v);
|
|
|
|
}
|
|
|
|
|
|
|
|
void r128Shl(R128 *dst, const R128 *src, int amount)
|
|
|
|
{
|
|
|
|
R128_U64 r[4];
|
|
|
|
|
|
|
|
r[0] = src->lo;
|
|
|
|
r[1] = src->hi;
|
|
|
|
|
|
|
|
amount &= 127;
|
|
|
|
if (amount >= 64) {
|
|
|
|
r[1] = r[0] << (amount - 64);
|
|
|
|
r[0] = 0;
|
|
|
|
} else if (amount) {
|
|
|
|
r[1] = (r[1] << amount) | (r[0] >> (64 - amount));
|
|
|
|
r[0] = r[0] << amount;
|
|
|
|
}
|
|
|
|
|
|
|
|
dst->lo = r[0];
|
|
|
|
dst->hi = r[1];
|
|
|
|
}
|
|
|
|
|
|
|
|
void r128Shr(R128 *dst, const R128 *src, int amount)
|
|
|
|
{
|
|
|
|
R128_U64 r[4];
|
|
|
|
|
|
|
|
r[2] = src->lo;
|
|
|
|
r[3] = src->hi;
|
|
|
|
|
|
|
|
amount &= 127;
|
|
|
|
if (amount >= 64) {
|
|
|
|
r[2] = r[3] >> (amount - 64);
|
|
|
|
r[3] = 0;
|
|
|
|
} else if (amount) {
|
|
|
|
r[2] = (r[2] >> amount) | (r[3] << (64 - amount));
|
|
|
|
r[3] = r[3] >> amount;
|
|
|
|
}
|
|
|
|
|
|
|
|
dst->lo = r[2];
|
|
|
|
dst->hi = r[3];
|
|
|
|
}
|
|
|
|
|
|
|
|
void r128Add(R128 *dst, const R128 *a, const R128 *b)
|
|
|
|
{
|
|
|
|
unsigned char carry = 0;
|
|
|
|
|
|
|
|
{
|
|
|
|
R128_U64 r = a->lo + b->lo;
|
|
|
|
carry = r < a->lo;
|
|
|
|
dst->lo = r;
|
|
|
|
dst->hi = a->hi + b->hi + carry;
|
|
|
|
}
|
|
|
|
}
|
|
|
|
|
|
|
|
void r128Sub(R128 *dst, const R128 *a, const R128 *b)
|
|
|
|
{
|
|
|
|
unsigned char borrow = 0;
|
|
|
|
|
|
|
|
{
|
|
|
|
R128_U64 r = a->lo - b->lo;
|
|
|
|
borrow = r > a->lo;
|
|
|
|
dst->lo = r;
|
|
|
|
dst->hi = a->hi - b->hi - borrow;
|
|
|
|
}
|
|
|
|
}
|
|
|
|
|
|
|
|
void r128Mul(R128 *dst, const R128 *a, const R128 *b)
|
|
|
|
{
|
|
|
|
int sign = 0;
|
|
|
|
R128 ta, tb, tc;
|
|
|
|
|
|
|
|
R128_SET2(&ta, a->lo, a->hi);
|
|
|
|
R128_SET2(&tb, b->lo, b->hi);
|
|
|
|
|
|
|
|
if (r128IsNeg(&ta)) {
|
|
|
|
r128__neg(&ta, &ta);
|
|
|
|
sign = !sign;
|
|
|
|
}
|
|
|
|
if (r128IsNeg(&tb)) {
|
|
|
|
r128__neg(&tb, &tb);
|
|
|
|
sign = !sign;
|
|
|
|
}
|
|
|
|
|
|
|
|
r128__umul(&tc, &ta, &tb);
|
|
|
|
if (sign) {
|
|
|
|
r128__neg(&tc, &tc);
|
|
|
|
}
|
|
|
|
|
|
|
|
r128Copy(dst, &tc);
|
|
|
|
}
|
|
|
|
|
|
|
|
int r128Cmp(const R128 *a, const R128 *b)
|
|
|
|
{
|
|
|
|
if (a->hi == b->hi) {
|
|
|
|
if (a->lo == b->lo) {
|
|
|
|
return 0;
|
|
|
|
} else if (a->lo > b->lo) {
|
|
|
|
return 1;
|
|
|
|
} else {
|
|
|
|
return -1;
|
|
|
|
}
|
|
|
|
} else if ((R128_S64)a->hi > (R128_S64)b->hi) {
|
|
|
|
return 1;
|
|
|
|
} else {
|
|
|
|
return -1;
|
|
|
|
}
|
|
|
|
}
|
|
|
|
|
|
|
|
int r128IsNeg(const R128 *v)
|
|
|
|
{
|
|
|
|
return (R128_S64)v->hi < 0;
|
|
|
|
}
|
|
|
|
|
|
|
|
#endif //R128_IMPLEMENTATION
|
|
|
|
|