/*
* Licensed to the Apache Software Foundation (ASF) under one or more
* contributor license agreements. See the NOTICE file distributed with
* this work for additional information regarding copyright ownership.
* The ASF licenses this file to You under the Apache License, Version 2.0
* (the "License"); you may not use this file except in compliance with
* the License. You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/
#include <string.h>
#include "cbigint.h"
#if defined(__linux__) || defined(__APPLE__)
#define USE_LL
#endif
#if __BYTE_ORDER == __LITTLE_ENDIAN
#define at(i) (i)
#else
#define at(i) ((i)^1)
/* the sequence for halfAt is -1, 2, 1, 4, 3, 6, 5, 8... */
/* and it should correspond to 0, 1, 2, 3, 4, 5, 6, 7... */
#define halfAt(i) (-((-(i)) ^ 1))
#endif
#define HIGH_IN_U64(u64) ((u64) >> 32)
#if defined(USE_LL)
#define LOW_IN_U64(u64) ((u64) & 0x00000000FFFFFFFFLL)
#else
#if defined(USE_L)
#define LOW_IN_U64(u64) ((u64) & 0x00000000FFFFFFFFL)
#else
#define LOW_IN_U64(u64) ((u64) & 0x00000000FFFFFFFF)
#endif /* USE_L */
#endif /* USE_LL */
#if defined(USE_LL)
#define TEN_E1 (0xALL)
#define TEN_E2 (0x64LL)
#define TEN_E3 (0x3E8LL)
#define TEN_E4 (0x2710LL)
#define TEN_E5 (0x186A0LL)
#define TEN_E6 (0xF4240LL)
#define TEN_E7 (0x989680LL)
#define TEN_E8 (0x5F5E100LL)
#define TEN_E9 (0x3B9ACA00LL)
#define TEN_E19 (0x8AC7230489E80000LL)
#else
#if defined(USE_L)
#define TEN_E1 (0xAL)
#define TEN_E2 (0x64L)
#define TEN_E3 (0x3E8L)
#define TEN_E4 (0x2710L)
#define TEN_E5 (0x186A0L)
#define TEN_E6 (0xF4240L)
#define TEN_E7 (0x989680L)
#define TEN_E8 (0x5F5E100L)
#define TEN_E9 (0x3B9ACA00L)
#define TEN_E19 (0x8AC7230489E80000L)
#else
#define TEN_E1 (0xA)
#define TEN_E2 (0x64)
#define TEN_E3 (0x3E8)
#define TEN_E4 (0x2710)
#define TEN_E5 (0x186A0)
#define TEN_E6 (0xF4240)
#define TEN_E7 (0x989680)
#define TEN_E8 (0x5F5E100)
#define TEN_E9 (0x3B9ACA00)
#define TEN_E19 (0x8AC7230489E80000)
#endif /* USE_L */
#endif /* USE_LL */
#define TIMES_TEN(x) (((x) << 3) + ((x) << 1))
#define bitSection(x, mask, shift) (((x) & (mask)) >> (shift))
#define CREATE_DOUBLE_BITS(normalizedM, e) (((normalizedM) & MANTISSA_MASK) | ((static_cast<uint64_t>((e) + E_OFFSET)) << 52))
#if defined(USE_LL)
#define MANTISSA_MASK (0x000FFFFFFFFFFFFFLL)
#define EXPONENT_MASK (0x7FF0000000000000LL)
#define NORMAL_MASK (0x0010000000000000LL)
#define SIGN_MASK (0x8000000000000000LL)
#else
#if defined(USE_L)
#define MANTISSA_MASK (0x000FFFFFFFFFFFFFL)
#define EXPONENT_MASK (0x7FF0000000000000L)
#define NORMAL_MASK (0x0010000000000000L)
#define SIGN_MASK (0x8000000000000000L)
#else
#define MANTISSA_MASK (0x000FFFFFFFFFFFFF)
#define EXPONENT_MASK (0x7FF0000000000000)
#define NORMAL_MASK (0x0010000000000000)
#define SIGN_MASK (0x8000000000000000)
#endif /* USE_L */
#endif /* USE_LL */
#define E_OFFSET (1075)
#define FLOAT_MANTISSA_MASK (0x007FFFFF)
#define FLOAT_EXPONENT_MASK (0x7F800000)
#define FLOAT_NORMAL_MASK (0x00800000)
#define FLOAT_E_OFFSET (150)
int32_t
simpleAddHighPrecision (uint64_t * arg1, int32_t length, uint64_t arg2)
{
/* assumes length > 0 */
int32_t index = 1;
*arg1 += arg2;
if (arg2 <= *arg1)
return 0;
else if (length == 1)
return 1;
while (++arg1[index] == 0 && ++index < length) {
}
return index == length;
}
int32_t
addHighPrecision (uint64_t * arg1, int32_t length1, uint64_t * arg2, int32_t length2)
{
/* addition is limited by length of arg1 as it this function is
* storing the result in arg1 */
/* fix for cc (GCC) 3.2 20020903 (Red Hat Linux 8.0 3.2-7): code generated does not
* do the temp1 + temp2 + carry addition correct. carry is 64 bit because gcc has
* subtle issues when you mix 64 / 32 bit maths. */
uint64_t temp1, temp2, temp3; /* temporary variables to help the SH-4, and gcc */
uint64_t carry;
int32_t index;
if (length1 == 0 || length2 == 0)
{
return 0;
}
else if (length1 < length2)
{
length2 = length1;
}
carry = 0;
index = 0;
do
{
temp1 = arg1[index];
temp2 = arg2[index];
temp3 = temp1 + temp2;
arg1[index] = temp3 + carry;
if (arg2[index] < arg1[index])
carry = 0;
else if (arg2[index] != arg1[index])
carry = 1;
}
while (++index < length2);
if (!carry)
return 0;
else if (index == length1)
return 1;
while (++arg1[index] == 0 && ++index < length1) {
}
return index == length1;
}
void
subtractHighPrecision (uint64_t * arg1, int32_t length1, uint64_t * arg2, int32_t length2)
{
/* assumes arg1 > arg2 */
int32_t index;
for (index = 0; index < length1; ++index)
arg1[index] = ~arg1[index];
simpleAddHighPrecision (arg1, length1, 1);
while (length2 > 0 && arg2[length2 - 1] == 0)
--length2;
addHighPrecision (arg1, length1, arg2, length2);
for (index = 0; index < length1; ++index)
arg1[index] = ~arg1[index];
simpleAddHighPrecision (arg1, length1, 1);
}
static uint32_t simpleMultiplyHighPrecision(uint64_t* arg1, int32_t length, uint64_t arg2) {
/* assumes arg2 only holds 32 bits of information */
uint64_t product;
int32_t index;
index = 0;
product = 0;
do
{
product =
HIGH_IN_U64 (product) + arg2 * LOW_U32_FROM_PTR (arg1 + index);
LOW_U32_FROM_PTR (arg1 + index) = LOW_U32_FROM_VAR (product);
product =
HIGH_IN_U64 (product) + arg2 * HIGH_U32_FROM_PTR (arg1 + index);
HIGH_U32_FROM_PTR (arg1 + index) = LOW_U32_FROM_VAR (product);
}
while (++index < length);
return HIGH_U32_FROM_VAR (product);
}
static void
simpleMultiplyAddHighPrecision (uint64_t * arg1, int32_t length, uint64_t arg2,
uint32_t * result)
{
/* Assumes result can hold the product and arg2 only holds 32 bits
of information */
uint64_t product;
int32_t index, resultIndex;
index = resultIndex = 0;
product = 0;
do
{
product =
HIGH_IN_U64 (product) + result[at (resultIndex)] +
arg2 * LOW_U32_FROM_PTR (arg1 + index);
result[at (resultIndex)] = LOW_U32_FROM_VAR (product);
++resultIndex;
product =
HIGH_IN_U64 (product) + result[at (resultIndex)] +
arg2 * HIGH_U32_FROM_PTR (arg1 + index);
result[at (resultIndex)] = LOW_U32_FROM_VAR (product);
++resultIndex;
}
while (++index < length);
result[at (resultIndex)] += HIGH_U32_FROM_VAR (product);
if (result[at (resultIndex)] < HIGH_U32_FROM_VAR (product))
{
/* must be careful with ++ operator and macro expansion */
++resultIndex;
while (++result[at (resultIndex)] == 0)
++resultIndex;
}
}
#if __BYTE_ORDER != __LITTLE_ENDIAN
void simpleMultiplyAddHighPrecisionBigEndianFix(uint64_t* arg1, int32_t length, uint64_t arg2, uint32_t* result) {
/* Assumes result can hold the product and arg2 only holds 32 bits of information */
int32_t index = 0;
int32_t resultIndex = 0;
uint64_t product = 0;
do {
product = HIGH_IN_U64(product) + result[halfAt(resultIndex)] + arg2 * LOW_U32_FROM_PTR(arg1 + index);
result[halfAt(resultIndex)] = LOW_U32_FROM_VAR(product);
++resultIndex;
product = HIGH_IN_U64(product) + result[halfAt(resultIndex)] + arg2 * HIGH_U32_FROM_PTR(arg1 + index);
result[halfAt(resultIndex)] = LOW_U32_FROM_VAR(product);
++resultIndex;
} while (++index < length);
result[halfAt(resultIndex)] += HIGH_U32_FROM_VAR(product);
if (result[halfAt(resultIndex)] < HIGH_U32_FROM_VAR(product)) {
/* must be careful with ++ operator and macro expansion */
++resultIndex;
while (++result[halfAt(resultIndex)] == 0) ++resultIndex;
}
}
#endif
void
multiplyHighPrecision (uint64_t * arg1, int32_t length1, uint64_t * arg2, int32_t length2,
uint64_t * result, int32_t length)
{
/* assumes result is large enough to hold product */
uint64_t* temp;
uint32_t* resultIn32;
int32_t count, index;
if (length1 < length2)
{
temp = arg1;
arg1 = arg2;
arg2 = temp;
count = length1;
length1 = length2;
length2 = count;
}
memset (result, 0, sizeof (uint64_t) * length);
/* length1 > length2 */
resultIn32 = reinterpret_cast<uint32_t*>(result);
index = -1;
for (count = 0; count < length2; ++count)
{
simpleMultiplyAddHighPrecision (arg1, length1, LOW_IN_U64 (arg2[count]),
resultIn32 + (++index));
#if __BYTE_ORDER == __LITTLE_ENDIAN
simpleMultiplyAddHighPrecision(arg1, length1, HIGH_IN_U64(arg2[count]), resultIn32 + (++index));
#else
simpleMultiplyAddHighPrecisionBigEndianFix(arg1, length1, HIGH_IN_U64(arg2[count]), resultIn32 + (++index));
#endif
}
}
uint32_t
simpleAppendDecimalDigitHighPrecision (uint64_t * arg1, int32_t length, uint64_t digit)
{
/* assumes digit is less than 32 bits */
uint64_t arg;
int32_t index = 0;
digit <<= 32;
do
{
arg = LOW_IN_U64 (arg1[index]);
digit = HIGH_IN_U64 (digit) + TIMES_TEN (arg);
LOW_U32_FROM_PTR (arg1 + index) = LOW_U32_FROM_VAR (digit);
arg = HIGH_IN_U64 (arg1[index]);
digit = HIGH_IN_U64 (digit) + TIMES_TEN (arg);
HIGH_U32_FROM_PTR (arg1 + index) = LOW_U32_FROM_VAR (digit);
}
while (++index < length);
return HIGH_U32_FROM_VAR (digit);
}
void
simpleShiftLeftHighPrecision (uint64_t * arg1, int32_t length, int32_t arg2)
{
/* assumes length > 0 */
int32_t index, offset;
if (arg2 >= 64)
{
offset = arg2 >> 6;
index = length;
while (--index - offset >= 0)
arg1[index] = arg1[index - offset];
do
{
arg1[index] = 0;
}
while (--index >= 0);
arg2 &= 0x3F;
}
if (arg2 == 0)
return;
while (--length > 0)
{
arg1[length] = arg1[length] << arg2 | arg1[length - 1] >> (64 - arg2);
}
*arg1 <<= arg2;
}
int32_t
highestSetBit (uint64_t * y)
{
uint32_t x;
int32_t result;
if (*y == 0)
return 0;
#if defined(USE_LL)
if (*y & 0xFFFFFFFF00000000LL)
{
x = HIGH_U32_FROM_PTR (y);
result = 32;
}
else
{
x = LOW_U32_FROM_PTR (y);
result = 0;
}
#else
#if defined(USE_L)
if (*y & 0xFFFFFFFF00000000L)
{
x = HIGH_U32_FROM_PTR (y);
result = 32;
}
else
{
x = LOW_U32_FROM_PTR (y);
result = 0;
}
#else
if (*y & 0xFFFFFFFF00000000)
{
x = HIGH_U32_FROM_PTR (y);
result = 32;
}
else
{
x = LOW_U32_FROM_PTR (y);
result = 0;
}
#endif /* USE_L */
#endif /* USE_LL */
if (x & 0xFFFF0000)
{
x = bitSection (x, 0xFFFF0000, 16);
result += 16;
}
if (x & 0xFF00)
{
x = bitSection (x, 0xFF00, 8);
result += 8;
}
if (x & 0xF0)
{
x = bitSection (x, 0xF0, 4);
result += 4;
}
if (x > 0x7)
return result + 4;
else if (x > 0x3)
return result + 3;
else if (x > 0x1)
return result + 2;
else
return result + 1;
}
int32_t
lowestSetBit (uint64_t * y)
{
uint32_t x;
int32_t result;
if (*y == 0)
return 0;
#if defined(USE_LL)
if (*y & 0x00000000FFFFFFFFLL)
{
x = LOW_U32_FROM_PTR (y);
result = 0;
}
else
{
x = HIGH_U32_FROM_PTR (y);
result = 32;
}
#else
#if defined(USE_L)
if (*y & 0x00000000FFFFFFFFL)
{
x = LOW_U32_FROM_PTR (y);
result = 0;
}
else
{
x = HIGH_U32_FROM_PTR (y);
result = 32;
}
#else
if (*y & 0x00000000FFFFFFFF)
{
x = LOW_U32_FROM_PTR (y);
result = 0;
}
else
{
x = HIGH_U32_FROM_PTR (y);
result = 32;
}
#endif /* USE_L */
#endif /* USE_LL */
if (!(x & 0xFFFF))
{
x = bitSection (x, 0xFFFF0000, 16);
result += 16;
}
if (!(x & 0xFF))
{
x = bitSection (x, 0xFF00, 8);
result += 8;
}
if (!(x & 0xF))
{
x = bitSection (x, 0xF0, 4);
result += 4;
}
if (x & 0x1)
return result + 1;
else if (x & 0x2)
return result + 2;
else if (x & 0x4)
return result + 3;
else
return result + 4;
}
int32_t
highestSetBitHighPrecision (uint64_t * arg, int32_t length)
{
int32_t highBit;
while (--length >= 0)
{
highBit = highestSetBit (arg + length);
if (highBit)
return highBit + 64 * length;
}
return 0;
}
int32_t
lowestSetBitHighPrecision (uint64_t * arg, int32_t length)
{
int32_t lowBit, index = -1;
while (++index < length)
{
lowBit = lowestSetBit (arg + index);
if (lowBit)
return lowBit + 64 * index;
}
return 0;
}
int32_t
compareHighPrecision (uint64_t * arg1, int32_t length1, uint64_t * arg2, int32_t length2)
{
while (--length1 >= 0 && arg1[length1] == 0) {
}
while (--length2 >= 0 && arg2[length2] == 0) {
}
if (length1 > length2)
return 1;
else if (length1 < length2)
return -1;
else if (length1 > -1)
{
do
{
if (arg1[length1] > arg2[length1])
return 1;
else if (arg1[length1] < arg2[length1])
return -1;
}
while (--length1 >= 0);
}
return 0;
}
jdouble
toDoubleHighPrecision (uint64_t * arg, int32_t length)
{
int32_t highBit;
uint64_t mantissa, test64;
uint32_t test;
jdouble result;
while (length > 0 && arg[length - 1] == 0)
--length;
if (length == 0)
result = 0.0;
else if (length > 16)
{
DOUBLE_TO_LONGBITS (result) = EXPONENT_MASK;
}
else if (length == 1)
{
highBit = highestSetBit (arg);
if (highBit <= 53)
{
highBit = 53 - highBit;
mantissa = *arg << highBit;
DOUBLE_TO_LONGBITS (result) =
CREATE_DOUBLE_BITS (mantissa, -highBit);
}
else
{
highBit -= 53;
mantissa = *arg >> highBit;
DOUBLE_TO_LONGBITS (result) =
CREATE_DOUBLE_BITS (mantissa, highBit);
/* perform rounding, round to even in case of tie */
test = (LOW_U32_FROM_PTR (arg) << (11 - highBit)) & 0x7FF;
if (test > 0x400 || ((test == 0x400) && (mantissa & 1)))
DOUBLE_TO_LONGBITS (result) = DOUBLE_TO_LONGBITS (result) + 1;
}
}
else
{
highBit = highestSetBit (arg + (--length));
if (highBit <= 53)
{
highBit = 53 - highBit;
if (highBit > 0)
{
mantissa =
(arg[length] << highBit) | (arg[length - 1] >>
(64 - highBit));
}
else
{
mantissa = arg[length];
}
DOUBLE_TO_LONGBITS (result) =
CREATE_DOUBLE_BITS (mantissa, length * 64 - highBit);
/* perform rounding, round to even in case of tie */
test64 = arg[--length] << highBit;
if (test64 > SIGN_MASK || ((test64 == SIGN_MASK) && (mantissa & 1)))
DOUBLE_TO_LONGBITS (result) = DOUBLE_TO_LONGBITS (result) + 1;
else if (test64 == SIGN_MASK)
{
while (--length >= 0)
{
if (arg[length] != 0)
{
DOUBLE_TO_LONGBITS (result) =
DOUBLE_TO_LONGBITS (result) + 1;
break;
}
}
}
}
else
{
highBit -= 53;
mantissa = arg[length] >> highBit;
DOUBLE_TO_LONGBITS (result) =
CREATE_DOUBLE_BITS (mantissa, length * 64 + highBit);
/* perform rounding, round to even in case of tie */
test = (LOW_U32_FROM_PTR (arg + length) << (11 - highBit)) & 0x7FF;
if (test > 0x400 || ((test == 0x400) && (mantissa & 1)))
DOUBLE_TO_LONGBITS (result) = DOUBLE_TO_LONGBITS (result) + 1;
else if (test == 0x400)
{
do
{
if (arg[--length] != 0)
{
DOUBLE_TO_LONGBITS (result) =
DOUBLE_TO_LONGBITS (result) + 1;
break;
}
}
while (length > 0);
}
}
}
return result;
}
static uint64_t simpleMultiplyHighPrecision64(uint64_t* arg1, int32_t length, uint64_t arg2);
int32_t
timesTenToTheEHighPrecision (uint64_t * result, int32_t length, jint e)
{
/* assumes result can hold value */
uint64_t overflow;
int exp10 = e;
if (e == 0)
return length;
/* bad O(n) way of doing it, but simple */
/*
do {
overflow = simpleAppendDecimalDigitHighPrecision(result, length, 0);
if (overflow)
result[length++] = overflow;
} while (--e);
*/
/* Replace the current implementaion which performs a
* "multiplication" by 10 e number of times with an actual
* multiplication. 10e19 is the largest exponent to the power of ten
* that will fit in a 64-bit integer, and 10e9 is the largest exponent to
* the power of ten that will fit in a 64-bit integer. Not sure where the
* break-even point is between an actual multiplication and a
* simpleAappendDecimalDigit() so just pick 10e3 as that point for
* now.
*/
while (exp10 >= 19)
{
overflow = simpleMultiplyHighPrecision64 (result, length, TEN_E19);
if (overflow)
result[length++] = overflow;
exp10 -= 19;
}
while (exp10 >= 9)
{
overflow = simpleMultiplyHighPrecision (result, length, TEN_E9);
if (overflow)
result[length++] = overflow;
exp10 -= 9;
}
if (exp10 == 0)
return length;
else if (exp10 == 1)
{
overflow = simpleAppendDecimalDigitHighPrecision (result, length, 0);
if (overflow)
result[length++] = overflow;
}
else if (exp10 == 2)
{
overflow = simpleAppendDecimalDigitHighPrecision (result, length, 0);
if (overflow)
result[length++] = overflow;
overflow = simpleAppendDecimalDigitHighPrecision (result, length, 0);
if (overflow)
result[length++] = overflow;
}
else if (exp10 == 3)
{
overflow = simpleMultiplyHighPrecision (result, length, TEN_E3);
if (overflow)
result[length++] = overflow;
}
else if (exp10 == 4)
{
overflow = simpleMultiplyHighPrecision (result, length, TEN_E4);
if (overflow)
result[length++] = overflow;
}
else if (exp10 == 5)
{
overflow = simpleMultiplyHighPrecision (result, length, TEN_E5);
if (overflow)
result[length++] = overflow;
}
else if (exp10 == 6)
{
overflow = simpleMultiplyHighPrecision (result, length, TEN_E6);
if (overflow)
result[length++] = overflow;
}
else if (exp10 == 7)
{
overflow = simpleMultiplyHighPrecision (result, length, TEN_E7);
if (overflow)
result[length++] = overflow;
}
else if (exp10 == 8)
{
overflow = simpleMultiplyHighPrecision (result, length, TEN_E8);
if (overflow)
result[length++] = overflow;
}
return length;
}
uint64_t
doubleMantissa (jdouble z)
{
uint64_t m = DOUBLE_TO_LONGBITS (z);
if ((m & EXPONENT_MASK) != 0)
m = (m & MANTISSA_MASK) | NORMAL_MASK;
else
m = (m & MANTISSA_MASK);
return m;
}
int32_t
doubleExponent (jdouble z)
{
/* assumes positive double */
int32_t k = HIGH_U32_FROM_VAR (z) >> 20;
if (k)
k -= E_OFFSET;
else
k = 1 - E_OFFSET;
return k;
}
uint32_t floatMantissa(jfloat z) {
uint32_t m = FLOAT_TO_INTBITS (z);
if ((m & FLOAT_EXPONENT_MASK) != 0)
m = (m & FLOAT_MANTISSA_MASK) | FLOAT_NORMAL_MASK;
else
m = (m & FLOAT_MANTISSA_MASK);
return m;
}
int32_t
floatExponent (jfloat z)
{
/* assumes positive float */
int32_t k = FLOAT_TO_INTBITS (z) >> 23;
if (k)
k -= FLOAT_E_OFFSET;
else
k = 1 - FLOAT_E_OFFSET;
return k;
}
/* Allow a 64-bit value in arg2 */
uint64_t
simpleMultiplyHighPrecision64 (uint64_t * arg1, int32_t length, uint64_t arg2)
{
uint64_t intermediate, carry1, carry2, prod1, prod2, sum;
uint64_t* pArg1;
int32_t index;
uint32_t buf32;
index = 0;
intermediate = 0;
pArg1 = arg1 + index;
carry1 = carry2 = 0;
do
{
if ((*pArg1 != 0) || (intermediate != 0))
{
prod1 =
static_cast<uint64_t>(LOW_U32_FROM_VAR (arg2)) * static_cast<uint64_t>(LOW_U32_FROM_PTR (pArg1));
sum = intermediate + prod1;
if ((sum < prod1) || (sum < intermediate))
{
carry1 = 1;
}
else
{
carry1 = 0;
}
prod1 =
static_cast<uint64_t>(LOW_U32_FROM_VAR (arg2)) * static_cast<uint64_t>(HIGH_U32_FROM_PTR (pArg1));
prod2 =
static_cast<uint64_t>(HIGH_U32_FROM_VAR (arg2)) * static_cast<uint64_t>(LOW_U32_FROM_PTR (pArg1));
intermediate = carry2 + HIGH_IN_U64 (sum) + prod1 + prod2;
if ((intermediate < prod1) || (intermediate < prod2))
{
carry2 = 1;
}
else
{
carry2 = 0;
}
LOW_U32_FROM_PTR (pArg1) = LOW_U32_FROM_VAR (sum);
buf32 = HIGH_U32_FROM_PTR (pArg1);
HIGH_U32_FROM_PTR (pArg1) = LOW_U32_FROM_VAR (intermediate);
intermediate = carry1 + HIGH_IN_U64 (intermediate)
+ static_cast<uint64_t>(HIGH_U32_FROM_VAR (arg2)) * static_cast<uint64_t>(buf32);
}
pArg1++;
}
while (++index < length);
return intermediate;
}