2 * Licensed to the Apache Software Foundation (ASF) under one or more
3 * contributor license agreements. See the NOTICE file distributed with
4 * this work for additional information regarding copyright ownership.
5 * The ASF licenses this file to You under the Apache License, Version 2.0
6 * (the "License"); you may not use this file except in compliance with
7 * the License. You may obtain a copy of the License at
9 * http://www.apache.org/licenses/LICENSE-2.0
11 * Unless required by applicable law or agreed to in writing, software
12 * distributed under the License is distributed on an "AS IS" BASIS,
13 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14 * See the License for the specific language governing permissions and
15 * limitations under the License.
21 #if defined(__linux__) || defined(FREEBSD) || defined(ZOS) || defined(MACOSX) || defined(AIX)
25 #if __BYTE_ORDER == __LITTLE_ENDIAN
29 /* the sequence for halfAt is -1, 2, 1, 4, 3, 6, 5, 8... */
30 /* and it should correspond to 0, 1, 2, 3, 4, 5, 6, 7... */
31 #define halfAt(i) (-((-(i)) ^ 1))
34 #define HIGH_IN_U64(u64) ((u64) >> 32)
36 #define LOW_IN_U64(u64) ((u64) & 0x00000000FFFFFFFFLL)
39 #define LOW_IN_U64(u64) ((u64) & 0x00000000FFFFFFFFL)
41 #define LOW_IN_U64(u64) ((u64) & 0x00000000FFFFFFFF)
46 #define TEN_E1 (0xALL)
47 #define TEN_E2 (0x64LL)
48 #define TEN_E3 (0x3E8LL)
49 #define TEN_E4 (0x2710LL)
50 #define TEN_E5 (0x186A0LL)
51 #define TEN_E6 (0xF4240LL)
52 #define TEN_E7 (0x989680LL)
53 #define TEN_E8 (0x5F5E100LL)
54 #define TEN_E9 (0x3B9ACA00LL)
55 #define TEN_E19 (0x8AC7230489E80000LL)
59 #define TEN_E2 (0x64L)
60 #define TEN_E3 (0x3E8L)
61 #define TEN_E4 (0x2710L)
62 #define TEN_E5 (0x186A0L)
63 #define TEN_E6 (0xF4240L)
64 #define TEN_E7 (0x989680L)
65 #define TEN_E8 (0x5F5E100L)
66 #define TEN_E9 (0x3B9ACA00L)
67 #define TEN_E19 (0x8AC7230489E80000L)
71 #define TEN_E3 (0x3E8)
72 #define TEN_E4 (0x2710)
73 #define TEN_E5 (0x186A0)
74 #define TEN_E6 (0xF4240)
75 #define TEN_E7 (0x989680)
76 #define TEN_E8 (0x5F5E100)
77 #define TEN_E9 (0x3B9ACA00)
78 #define TEN_E19 (0x8AC7230489E80000)
82 #define TIMES_TEN(x) (((x) << 3) + ((x) << 1))
83 #define bitSection(x, mask, shift) (((x) & (mask)) >> (shift))
84 #define CREATE_DOUBLE_BITS(normalizedM, e) (((normalizedM) & MANTISSA_MASK) | ((static_cast<uint64_t>((e) + E_OFFSET)) << 52))
87 #define MANTISSA_MASK (0x000FFFFFFFFFFFFFLL)
88 #define EXPONENT_MASK (0x7FF0000000000000LL)
89 #define NORMAL_MASK (0x0010000000000000LL)
90 #define SIGN_MASK (0x8000000000000000LL)
93 #define MANTISSA_MASK (0x000FFFFFFFFFFFFFL)
94 #define EXPONENT_MASK (0x7FF0000000000000L)
95 #define NORMAL_MASK (0x0010000000000000L)
96 #define SIGN_MASK (0x8000000000000000L)
98 #define MANTISSA_MASK (0x000FFFFFFFFFFFFF)
99 #define EXPONENT_MASK (0x7FF0000000000000)
100 #define NORMAL_MASK (0x0010000000000000)
101 #define SIGN_MASK (0x8000000000000000)
105 #define E_OFFSET (1075)
107 #define FLOAT_MANTISSA_MASK (0x007FFFFF)
108 #define FLOAT_EXPONENT_MASK (0x7F800000)
109 #define FLOAT_NORMAL_MASK (0x00800000)
110 #define FLOAT_E_OFFSET (150)
113 simpleAddHighPrecision (uint64_t * arg1, int32_t length, uint64_t arg2)
115 /* assumes length > 0 */
121 else if (length == 1)
124 while (++arg1[index] == 0 && ++index < length);
126 return index == length;
130 addHighPrecision (uint64_t * arg1, int32_t length1, uint64_t * arg2, int32_t length2)
132 /* addition is limited by length of arg1 as it this function is
133 * storing the result in arg1 */
134 /* fix for cc (GCC) 3.2 20020903 (Red Hat Linux 8.0 3.2-7): code generated does not
135 * do the temp1 + temp2 + carry addition correct. carry is 64 bit because gcc has
136 * subtle issues when you mix 64 / 32 bit maths. */
137 uint64_t temp1, temp2, temp3; /* temporary variables to help the SH-4, and gcc */
141 if (length1 == 0 || length2 == 0)
145 else if (length1 < length2)
156 temp3 = temp1 + temp2;
157 arg1[index] = temp3 + carry;
158 if (arg2[index] < arg1[index])
160 else if (arg2[index] != arg1[index])
163 while (++index < length2);
166 else if (index == length1)
169 while (++arg1[index] == 0 && ++index < length1);
171 return index == length1;
175 subtractHighPrecision (uint64_t * arg1, int32_t length1, uint64_t * arg2, int32_t length2)
177 /* assumes arg1 > arg2 */
179 for (index = 0; index < length1; ++index)
180 arg1[index] = ~arg1[index];
181 simpleAddHighPrecision (arg1, length1, 1);
183 while (length2 > 0 && arg2[length2 - 1] == 0)
186 addHighPrecision (arg1, length1, arg2, length2);
188 for (index = 0; index < length1; ++index)
189 arg1[index] = ~arg1[index];
190 simpleAddHighPrecision (arg1, length1, 1);
193 static uint32_t simpleMultiplyHighPrecision(uint64_t* arg1, int32_t length, uint64_t arg2) {
194 /* assumes arg2 only holds 32 bits of information */
204 HIGH_IN_U64 (product) + arg2 * LOW_U32_FROM_PTR (arg1 + index);
205 LOW_U32_FROM_PTR (arg1 + index) = LOW_U32_FROM_VAR (product);
207 HIGH_IN_U64 (product) + arg2 * HIGH_U32_FROM_PTR (arg1 + index);
208 HIGH_U32_FROM_PTR (arg1 + index) = LOW_U32_FROM_VAR (product);
210 while (++index < length);
212 return HIGH_U32_FROM_VAR (product);
216 simpleMultiplyAddHighPrecision (uint64_t * arg1, int32_t length, uint64_t arg2,
219 /* Assumes result can hold the product and arg2 only holds 32 bits
222 int32_t index, resultIndex;
224 index = resultIndex = 0;
230 HIGH_IN_U64 (product) + result[at (resultIndex)] +
231 arg2 * LOW_U32_FROM_PTR (arg1 + index);
232 result[at (resultIndex)] = LOW_U32_FROM_VAR (product);
235 HIGH_IN_U64 (product) + result[at (resultIndex)] +
236 arg2 * HIGH_U32_FROM_PTR (arg1 + index);
237 result[at (resultIndex)] = LOW_U32_FROM_VAR (product);
240 while (++index < length);
242 result[at (resultIndex)] += HIGH_U32_FROM_VAR (product);
243 if (result[at (resultIndex)] < HIGH_U32_FROM_VAR (product))
245 /* must be careful with ++ operator and macro expansion */
247 while (++result[at (resultIndex)] == 0)
252 #if __BYTE_ORDER != __LITTLE_ENDIAN
253 void simpleMultiplyAddHighPrecisionBigEndianFix(uint64_t* arg1, int32_t length, uint64_t arg2, uint32_t* result) {
254 /* Assumes result can hold the product and arg2 only holds 32 bits
257 int32_t index, resultIndex;
259 index = resultIndex = 0;
263 product = HIGH_IN_U64(product) + result[halfAt(resultIndex)] + arg2 * LOW_U32_FROM_PTR(arg1 + index);
264 result[halfAt(resultIndex)] = LOW_U32_FROM_VAR(product);
266 product = HIGH_IN_U64(product) + result[halfAt(resultIndex)] + arg2 * HIGH_U32_FROM_PTR(arg1 + index);
267 result[halfAt(resultIndex)] = LOW_U32_FROM_VAR(product);
269 } while (++index < length);
271 result[halfAt(resultIndex)] += HIGH_U32_FROM_VAR(product);
272 if (result[halfAt(resultIndex)] < HIGH_U32_FROM_VAR(product)) {
273 /* must be careful with ++ operator and macro expansion */
275 while (++result[halfAt(resultIndex)] == 0) ++resultIndex;
281 multiplyHighPrecision (uint64_t * arg1, int32_t length1, uint64_t * arg2, int32_t length2,
282 uint64_t * result, int32_t length)
284 /* assumes result is large enough to hold product */
286 uint32_t* resultIn32;
287 int32_t count, index;
289 if (length1 < length2)
299 memset (result, 0, sizeof (uint64_t) * length);
301 /* length1 > length2 */
302 resultIn32 = reinterpret_cast<uint32_t*>(result);
304 for (count = 0; count < length2; ++count)
306 simpleMultiplyAddHighPrecision (arg1, length1, LOW_IN_U64 (arg2[count]),
307 resultIn32 + (++index));
308 #if __BYTE_ORDER == __LITTLE_ENDIAN
309 simpleMultiplyAddHighPrecision(arg1, length1, HIGH_IN_U64(arg2[count]), resultIn32 + (++index));
311 simpleMultiplyAddHighPrecisionBigEndianFix(arg1, length1, HIGH_IN_U64(arg2[count]), resultIn32 + (++index));
317 simpleAppendDecimalDigitHighPrecision (uint64_t * arg1, int32_t length, uint64_t digit)
319 /* assumes digit is less than 32 bits */
326 arg = LOW_IN_U64 (arg1[index]);
327 digit = HIGH_IN_U64 (digit) + TIMES_TEN (arg);
328 LOW_U32_FROM_PTR (arg1 + index) = LOW_U32_FROM_VAR (digit);
330 arg = HIGH_IN_U64 (arg1[index]);
331 digit = HIGH_IN_U64 (digit) + TIMES_TEN (arg);
332 HIGH_U32_FROM_PTR (arg1 + index) = LOW_U32_FROM_VAR (digit);
334 while (++index < length);
336 return HIGH_U32_FROM_VAR (digit);
340 simpleShiftLeftHighPrecision (uint64_t * arg1, int32_t length, int32_t arg2)
342 /* assumes length > 0 */
343 int32_t index, offset;
349 while (--index - offset >= 0)
350 arg1[index] = arg1[index - offset];
355 while (--index >= 0);
364 arg1[length] = arg1[length] << arg2 | arg1[length - 1] >> (64 - arg2);
370 highestSetBit (uint64_t * y)
379 if (*y & 0xFFFFFFFF00000000LL)
381 x = HIGH_U32_FROM_PTR (y);
386 x = LOW_U32_FROM_PTR (y);
391 if (*y & 0xFFFFFFFF00000000L)
393 x = HIGH_U32_FROM_PTR (y);
398 x = LOW_U32_FROM_PTR (y);
402 if (*y & 0xFFFFFFFF00000000)
404 x = HIGH_U32_FROM_PTR (y);
409 x = LOW_U32_FROM_PTR (y);
417 x = bitSection (x, 0xFFFF0000, 16);
422 x = bitSection (x, 0xFF00, 8);
427 x = bitSection (x, 0xF0, 4);
441 lowestSetBit (uint64_t * y)
450 if (*y & 0x00000000FFFFFFFFLL)
452 x = LOW_U32_FROM_PTR (y);
457 x = HIGH_U32_FROM_PTR (y);
462 if (*y & 0x00000000FFFFFFFFL)
464 x = LOW_U32_FROM_PTR (y);
469 x = HIGH_U32_FROM_PTR (y);
473 if (*y & 0x00000000FFFFFFFF)
475 x = LOW_U32_FROM_PTR (y);
480 x = HIGH_U32_FROM_PTR (y);
488 x = bitSection (x, 0xFFFF0000, 16);
493 x = bitSection (x, 0xFF00, 8);
498 x = bitSection (x, 0xF0, 4);
513 highestSetBitHighPrecision (uint64_t * arg, int32_t length)
517 while (--length >= 0)
519 highBit = highestSetBit (arg + length);
521 return highBit + 64 * length;
528 lowestSetBitHighPrecision (uint64_t * arg, int32_t length)
530 int32_t lowBit, index = -1;
532 while (++index < length)
534 lowBit = lowestSetBit (arg + index);
536 return lowBit + 64 * index;
543 compareHighPrecision (uint64_t * arg1, int32_t length1, uint64_t * arg2, int32_t length2)
545 while (--length1 >= 0 && arg1[length1] == 0);
546 while (--length2 >= 0 && arg2[length2] == 0);
548 if (length1 > length2)
550 else if (length1 < length2)
552 else if (length1 > -1)
556 if (arg1[length1] > arg2[length1])
558 else if (arg1[length1] < arg2[length1])
561 while (--length1 >= 0);
568 toDoubleHighPrecision (uint64_t * arg, int32_t length)
571 uint64_t mantissa, test64;
575 while (length > 0 && arg[length - 1] == 0)
580 else if (length > 16)
582 DOUBLE_TO_LONGBITS (result) = EXPONENT_MASK;
584 else if (length == 1)
586 highBit = highestSetBit (arg);
589 highBit = 53 - highBit;
590 mantissa = *arg << highBit;
591 DOUBLE_TO_LONGBITS (result) =
592 CREATE_DOUBLE_BITS (mantissa, -highBit);
597 mantissa = *arg >> highBit;
598 DOUBLE_TO_LONGBITS (result) =
599 CREATE_DOUBLE_BITS (mantissa, highBit);
601 /* perform rounding, round to even in case of tie */
602 test = (LOW_U32_FROM_PTR (arg) << (11 - highBit)) & 0x7FF;
603 if (test > 0x400 || ((test == 0x400) && (mantissa & 1)))
604 DOUBLE_TO_LONGBITS (result) = DOUBLE_TO_LONGBITS (result) + 1;
609 highBit = highestSetBit (arg + (--length));
612 highBit = 53 - highBit;
616 (arg[length] << highBit) | (arg[length - 1] >>
621 mantissa = arg[length];
623 DOUBLE_TO_LONGBITS (result) =
624 CREATE_DOUBLE_BITS (mantissa, length * 64 - highBit);
626 /* perform rounding, round to even in case of tie */
627 test64 = arg[--length] << highBit;
628 if (test64 > SIGN_MASK || ((test64 == SIGN_MASK) && (mantissa & 1)))
629 DOUBLE_TO_LONGBITS (result) = DOUBLE_TO_LONGBITS (result) + 1;
630 else if (test64 == SIGN_MASK)
632 while (--length >= 0)
634 if (arg[length] != 0)
636 DOUBLE_TO_LONGBITS (result) =
637 DOUBLE_TO_LONGBITS (result) + 1;
646 mantissa = arg[length] >> highBit;
647 DOUBLE_TO_LONGBITS (result) =
648 CREATE_DOUBLE_BITS (mantissa, length * 64 + highBit);
650 /* perform rounding, round to even in case of tie */
651 test = (LOW_U32_FROM_PTR (arg + length) << (11 - highBit)) & 0x7FF;
652 if (test > 0x400 || ((test == 0x400) && (mantissa & 1)))
653 DOUBLE_TO_LONGBITS (result) = DOUBLE_TO_LONGBITS (result) + 1;
654 else if (test == 0x400)
658 if (arg[--length] != 0)
660 DOUBLE_TO_LONGBITS (result) =
661 DOUBLE_TO_LONGBITS (result) + 1;
673 static uint64_t simpleMultiplyHighPrecision64(uint64_t* arg1, int32_t length, uint64_t arg2);
676 timesTenToTheEHighPrecision (uint64_t * result, int32_t length, jint e)
678 /* assumes result can hold value */
685 /* bad O(n) way of doing it, but simple */
688 overflow = simpleAppendDecimalDigitHighPrecision(result, length, 0);
690 result[length++] = overflow;
693 /* Replace the current implementaion which performs a
694 * "multiplication" by 10 e number of times with an actual
695 * multiplication. 10e19 is the largest exponent to the power of ten
696 * that will fit in a 64-bit integer, and 10e9 is the largest exponent to
697 * the power of ten that will fit in a 64-bit integer. Not sure where the
698 * break-even point is between an actual multiplication and a
699 * simpleAappendDecimalDigit() so just pick 10e3 as that point for
704 overflow = simpleMultiplyHighPrecision64 (result, length, TEN_E19);
706 result[length++] = overflow;
711 overflow = simpleMultiplyHighPrecision (result, length, TEN_E9);
713 result[length++] = overflow;
720 overflow = simpleAppendDecimalDigitHighPrecision (result, length, 0);
722 result[length++] = overflow;
726 overflow = simpleAppendDecimalDigitHighPrecision (result, length, 0);
728 result[length++] = overflow;
729 overflow = simpleAppendDecimalDigitHighPrecision (result, length, 0);
731 result[length++] = overflow;
735 overflow = simpleMultiplyHighPrecision (result, length, TEN_E3);
737 result[length++] = overflow;
741 overflow = simpleMultiplyHighPrecision (result, length, TEN_E4);
743 result[length++] = overflow;
747 overflow = simpleMultiplyHighPrecision (result, length, TEN_E5);
749 result[length++] = overflow;
753 overflow = simpleMultiplyHighPrecision (result, length, TEN_E6);
755 result[length++] = overflow;
759 overflow = simpleMultiplyHighPrecision (result, length, TEN_E7);
761 result[length++] = overflow;
765 overflow = simpleMultiplyHighPrecision (result, length, TEN_E8);
767 result[length++] = overflow;
773 doubleMantissa (jdouble z)
775 uint64_t m = DOUBLE_TO_LONGBITS (z);
777 if ((m & EXPONENT_MASK) != 0)
778 m = (m & MANTISSA_MASK) | NORMAL_MASK;
780 m = (m & MANTISSA_MASK);
786 doubleExponent (jdouble z)
788 /* assumes positive double */
789 int32_t k = HIGH_U32_FROM_VAR (z) >> 20;
799 uint32_t floatMantissa(jfloat z) {
800 uint32_t m = FLOAT_TO_INTBITS (z);
802 if ((m & FLOAT_EXPONENT_MASK) != 0)
803 m = (m & FLOAT_MANTISSA_MASK) | FLOAT_NORMAL_MASK;
805 m = (m & FLOAT_MANTISSA_MASK);
811 floatExponent (jfloat z)
813 /* assumes positive float */
814 int32_t k = FLOAT_TO_INTBITS (z) >> 23;
818 k = 1 - FLOAT_E_OFFSET;
823 /* Allow a 64-bit value in arg2 */
825 simpleMultiplyHighPrecision64 (uint64_t * arg1, int32_t length, uint64_t arg2)
827 uint64_t intermediate, carry1, carry2, prod1, prod2, sum;
834 pArg1 = arg1 + index;
839 if ((*pArg1 != 0) || (intermediate != 0))
842 static_cast<uint64_t>(LOW_U32_FROM_VAR (arg2)) * static_cast<uint64_t>(LOW_U32_FROM_PTR (pArg1));
843 sum = intermediate + prod1;
844 if ((sum < prod1) || (sum < intermediate))
853 static_cast<uint64_t>(LOW_U32_FROM_VAR (arg2)) * static_cast<uint64_t>(HIGH_U32_FROM_PTR (pArg1));
855 static_cast<uint64_t>(HIGH_U32_FROM_VAR (arg2)) * static_cast<uint64_t>(LOW_U32_FROM_PTR (pArg1));
856 intermediate = carry2 + HIGH_IN_U64 (sum) + prod1 + prod2;
857 if ((intermediate < prod1) || (intermediate < prod2))
865 LOW_U32_FROM_PTR (pArg1) = LOW_U32_FROM_VAR (sum);
866 buf32 = HIGH_U32_FROM_PTR (pArg1);
867 HIGH_U32_FROM_PTR (pArg1) = LOW_U32_FROM_VAR (intermediate);
868 intermediate = carry1 + HIGH_IN_U64 (intermediate)
869 + static_cast<uint64_t>(HIGH_U32_FROM_VAR (arg2)) * static_cast<uint64_t>(buf32);
873 while (++index < length);