• Home
  • Line#
  • Scopes#
  • Navigate#
  • Raw
  • Download
1 /*
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
8  *
9  *     http://www.apache.org/licenses/LICENSE-2.0
10  *
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.
16  */
17 
18 #include <string.h>
19 #include "cbigint.h"
20 
21 #if defined(__linux__) || defined(FREEBSD) || defined(ZOS) || defined(MACOSX) || defined(AIX)
22 #define USE_LL
23 #endif
24 
25 #if __BYTE_ORDER == __LITTLE_ENDIAN
26 #define at(i) (i)
27 #else
28 #define at(i) ((i)^1)
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))
32 #endif
33 
34 #define HIGH_IN_U64(u64) ((u64) >> 32)
35 #if defined(USE_LL)
36 #define LOW_IN_U64(u64) ((u64) & 0x00000000FFFFFFFFLL)
37 #else
38 #if defined(USE_L)
39 #define LOW_IN_U64(u64) ((u64) & 0x00000000FFFFFFFFL)
40 #else
41 #define LOW_IN_U64(u64) ((u64) & 0x00000000FFFFFFFF)
42 #endif /* USE_L */
43 #endif /* USE_LL */
44 
45 #if defined(USE_LL)
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)
56 #else
57 #if defined(USE_L)
58 #define TEN_E1 (0xAL)
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)
68 #else
69 #define TEN_E1 (0xA)
70 #define TEN_E2 (0x64)
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)
79 #endif /* USE_L */
80 #endif /* USE_LL */
81 
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))
85 
86 #if defined(USE_LL)
87 #define MANTISSA_MASK (0x000FFFFFFFFFFFFFLL)
88 #define EXPONENT_MASK (0x7FF0000000000000LL)
89 #define NORMAL_MASK (0x0010000000000000LL)
90 #define SIGN_MASK (0x8000000000000000LL)
91 #else
92 #if defined(USE_L)
93 #define MANTISSA_MASK (0x000FFFFFFFFFFFFFL)
94 #define EXPONENT_MASK (0x7FF0000000000000L)
95 #define NORMAL_MASK (0x0010000000000000L)
96 #define SIGN_MASK (0x8000000000000000L)
97 #else
98 #define MANTISSA_MASK (0x000FFFFFFFFFFFFF)
99 #define EXPONENT_MASK (0x7FF0000000000000)
100 #define NORMAL_MASK (0x0010000000000000)
101 #define SIGN_MASK (0x8000000000000000)
102 #endif /* USE_L */
103 #endif /* USE_LL */
104 
105 #define E_OFFSET (1075)
106 
107 #define FLOAT_MANTISSA_MASK (0x007FFFFF)
108 #define FLOAT_EXPONENT_MASK (0x7F800000)
109 #define FLOAT_NORMAL_MASK   (0x00800000)
110 #define FLOAT_E_OFFSET (150)
111 
112 int32_t
simpleAddHighPrecision(uint64_t * arg1,int32_t length,uint64_t arg2)113 simpleAddHighPrecision (uint64_t * arg1, int32_t length, uint64_t arg2)
114 {
115   /* assumes length > 0 */
116   int32_t index = 1;
117 
118   *arg1 += arg2;
119   if (arg2 <= *arg1)
120     return 0;
121   else if (length == 1)
122     return 1;
123 
124   while (++arg1[index] == 0 && ++index < length);
125 
126   return index == length;
127 }
128 
129 int32_t
addHighPrecision(uint64_t * arg1,int32_t length1,uint64_t * arg2,int32_t length2)130 addHighPrecision (uint64_t * arg1, int32_t length1, uint64_t * arg2, int32_t length2)
131 {
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 */
138   uint64_t carry;
139   int32_t index;
140 
141   if (length1 == 0 || length2 == 0)
142     {
143       return 0;
144     }
145   else if (length1 < length2)
146     {
147       length2 = length1;
148     }
149 
150   carry = 0;
151   index = 0;
152   do
153     {
154       temp1 = arg1[index];
155       temp2 = arg2[index];
156       temp3 = temp1 + temp2;
157       arg1[index] = temp3 + carry;
158       if (arg2[index] < arg1[index])
159         carry = 0;
160       else if (arg2[index] != arg1[index])
161         carry = 1;
162     }
163   while (++index < length2);
164   if (!carry)
165     return 0;
166   else if (index == length1)
167     return 1;
168 
169   while (++arg1[index] == 0 && ++index < length1);
170 
171   return index == length1;
172 }
173 
174 void
subtractHighPrecision(uint64_t * arg1,int32_t length1,uint64_t * arg2,int32_t length2)175 subtractHighPrecision (uint64_t * arg1, int32_t length1, uint64_t * arg2, int32_t length2)
176 {
177   /* assumes arg1 > arg2 */
178   int32_t index;
179   for (index = 0; index < length1; ++index)
180     arg1[index] = ~arg1[index];
181   simpleAddHighPrecision (arg1, length1, 1);
182 
183   while (length2 > 0 && arg2[length2 - 1] == 0)
184     --length2;
185 
186   addHighPrecision (arg1, length1, arg2, length2);
187 
188   for (index = 0; index < length1; ++index)
189     arg1[index] = ~arg1[index];
190   simpleAddHighPrecision (arg1, length1, 1);
191 }
192 
simpleMultiplyHighPrecision(uint64_t * arg1,int32_t length,uint64_t arg2)193 static uint32_t simpleMultiplyHighPrecision(uint64_t* arg1, int32_t length, uint64_t arg2) {
194   /* assumes arg2 only holds 32 bits of information */
195   uint64_t product;
196   int32_t index;
197 
198   index = 0;
199   product = 0;
200 
201   do
202     {
203       product =
204         HIGH_IN_U64 (product) + arg2 * LOW_U32_FROM_PTR (arg1 + index);
205       LOW_U32_FROM_PTR (arg1 + index) = LOW_U32_FROM_VAR (product);
206       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);
209     }
210   while (++index < length);
211 
212   return HIGH_U32_FROM_VAR (product);
213 }
214 
215 static void
simpleMultiplyAddHighPrecision(uint64_t * arg1,int32_t length,uint64_t arg2,uint32_t * result)216 simpleMultiplyAddHighPrecision (uint64_t * arg1, int32_t length, uint64_t arg2,
217                                 uint32_t * result)
218 {
219   /* Assumes result can hold the product and arg2 only holds 32 bits
220      of information */
221   uint64_t product;
222   int32_t index, resultIndex;
223 
224   index = resultIndex = 0;
225   product = 0;
226 
227   do
228     {
229       product =
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);
233       ++resultIndex;
234       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);
238       ++resultIndex;
239     }
240   while (++index < length);
241 
242   result[at (resultIndex)] += HIGH_U32_FROM_VAR (product);
243   if (result[at (resultIndex)] < HIGH_U32_FROM_VAR (product))
244     {
245       /* must be careful with ++ operator and macro expansion */
246       ++resultIndex;
247       while (++result[at (resultIndex)] == 0)
248         ++resultIndex;
249     }
250 }
251 
252 #if __BYTE_ORDER != __LITTLE_ENDIAN
simpleMultiplyAddHighPrecisionBigEndianFix(uint64_t * arg1,int32_t length,uint64_t arg2,uint32_t * result)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 of information */
255     int32_t index = 0;
256     int32_t resultIndex = 0;
257     uint64_t product = 0;
258 
259     do {
260         product = HIGH_IN_U64(product) + result[halfAt(resultIndex)] + arg2 * LOW_U32_FROM_PTR(arg1 + index);
261         result[halfAt(resultIndex)] = LOW_U32_FROM_VAR(product);
262         ++resultIndex;
263         product = HIGH_IN_U64(product) + result[halfAt(resultIndex)] + arg2 * HIGH_U32_FROM_PTR(arg1 + index);
264         result[halfAt(resultIndex)] = LOW_U32_FROM_VAR(product);
265         ++resultIndex;
266     } while (++index < length);
267 
268     result[halfAt(resultIndex)] += HIGH_U32_FROM_VAR(product);
269     if (result[halfAt(resultIndex)] < HIGH_U32_FROM_VAR(product)) {
270         /* must be careful with ++ operator and macro expansion */
271         ++resultIndex;
272         while (++result[halfAt(resultIndex)] == 0) ++resultIndex;
273     }
274 }
275 #endif
276 
277 void
multiplyHighPrecision(uint64_t * arg1,int32_t length1,uint64_t * arg2,int32_t length2,uint64_t * result,int32_t length)278 multiplyHighPrecision (uint64_t * arg1, int32_t length1, uint64_t * arg2, int32_t length2,
279                        uint64_t * result, int32_t length)
280 {
281   /* assumes result is large enough to hold product */
282   uint64_t* temp;
283   uint32_t* resultIn32;
284   int32_t count, index;
285 
286   if (length1 < length2)
287     {
288       temp = arg1;
289       arg1 = arg2;
290       arg2 = temp;
291       count = length1;
292       length1 = length2;
293       length2 = count;
294     }
295 
296   memset (result, 0, sizeof (uint64_t) * length);
297 
298   /* length1 > length2 */
299   resultIn32 = reinterpret_cast<uint32_t*>(result);
300   index = -1;
301   for (count = 0; count < length2; ++count)
302     {
303       simpleMultiplyAddHighPrecision (arg1, length1, LOW_IN_U64 (arg2[count]),
304                                       resultIn32 + (++index));
305 #if __BYTE_ORDER == __LITTLE_ENDIAN
306       simpleMultiplyAddHighPrecision(arg1, length1, HIGH_IN_U64(arg2[count]), resultIn32 + (++index));
307 #else
308       simpleMultiplyAddHighPrecisionBigEndianFix(arg1, length1, HIGH_IN_U64(arg2[count]), resultIn32 + (++index));
309 #endif
310     }
311 }
312 
313 uint32_t
simpleAppendDecimalDigitHighPrecision(uint64_t * arg1,int32_t length,uint64_t digit)314 simpleAppendDecimalDigitHighPrecision (uint64_t * arg1, int32_t length, uint64_t digit)
315 {
316   /* assumes digit is less than 32 bits */
317   uint64_t arg;
318   int32_t index = 0;
319 
320   digit <<= 32;
321   do
322     {
323       arg = LOW_IN_U64 (arg1[index]);
324       digit = HIGH_IN_U64 (digit) + TIMES_TEN (arg);
325       LOW_U32_FROM_PTR (arg1 + index) = LOW_U32_FROM_VAR (digit);
326 
327       arg = HIGH_IN_U64 (arg1[index]);
328       digit = HIGH_IN_U64 (digit) + TIMES_TEN (arg);
329       HIGH_U32_FROM_PTR (arg1 + index) = LOW_U32_FROM_VAR (digit);
330     }
331   while (++index < length);
332 
333   return HIGH_U32_FROM_VAR (digit);
334 }
335 
336 void
simpleShiftLeftHighPrecision(uint64_t * arg1,int32_t length,int32_t arg2)337 simpleShiftLeftHighPrecision (uint64_t * arg1, int32_t length, int32_t arg2)
338 {
339   /* assumes length > 0 */
340   int32_t index, offset;
341   if (arg2 >= 64)
342     {
343       offset = arg2 >> 6;
344       index = length;
345 
346       while (--index - offset >= 0)
347         arg1[index] = arg1[index - offset];
348       do
349         {
350           arg1[index] = 0;
351         }
352       while (--index >= 0);
353 
354       arg2 &= 0x3F;
355     }
356 
357   if (arg2 == 0)
358     return;
359   while (--length > 0)
360     {
361       arg1[length] = arg1[length] << arg2 | arg1[length - 1] >> (64 - arg2);
362     }
363   *arg1 <<= arg2;
364 }
365 
366 int32_t
highestSetBit(uint64_t * y)367 highestSetBit (uint64_t * y)
368 {
369   uint32_t x;
370   int32_t result;
371 
372   if (*y == 0)
373     return 0;
374 
375 #if defined(USE_LL)
376   if (*y & 0xFFFFFFFF00000000LL)
377     {
378       x = HIGH_U32_FROM_PTR (y);
379       result = 32;
380     }
381   else
382     {
383       x = LOW_U32_FROM_PTR (y);
384       result = 0;
385     }
386 #else
387 #if defined(USE_L)
388   if (*y & 0xFFFFFFFF00000000L)
389     {
390       x = HIGH_U32_FROM_PTR (y);
391       result = 32;
392     }
393   else
394     {
395       x = LOW_U32_FROM_PTR (y);
396       result = 0;
397     }
398 #else
399   if (*y & 0xFFFFFFFF00000000)
400     {
401       x = HIGH_U32_FROM_PTR (y);
402       result = 32;
403     }
404   else
405     {
406       x = LOW_U32_FROM_PTR (y);
407       result = 0;
408     }
409 #endif /* USE_L */
410 #endif /* USE_LL */
411 
412   if (x & 0xFFFF0000)
413     {
414       x = bitSection (x, 0xFFFF0000, 16);
415       result += 16;
416     }
417   if (x & 0xFF00)
418     {
419       x = bitSection (x, 0xFF00, 8);
420       result += 8;
421     }
422   if (x & 0xF0)
423     {
424       x = bitSection (x, 0xF0, 4);
425       result += 4;
426     }
427   if (x > 0x7)
428     return result + 4;
429   else if (x > 0x3)
430     return result + 3;
431   else if (x > 0x1)
432     return result + 2;
433   else
434     return result + 1;
435 }
436 
437 int32_t
lowestSetBit(uint64_t * y)438 lowestSetBit (uint64_t * y)
439 {
440   uint32_t x;
441   int32_t result;
442 
443   if (*y == 0)
444     return 0;
445 
446 #if defined(USE_LL)
447   if (*y & 0x00000000FFFFFFFFLL)
448     {
449       x = LOW_U32_FROM_PTR (y);
450       result = 0;
451     }
452   else
453     {
454       x = HIGH_U32_FROM_PTR (y);
455       result = 32;
456     }
457 #else
458 #if defined(USE_L)
459   if (*y & 0x00000000FFFFFFFFL)
460     {
461       x = LOW_U32_FROM_PTR (y);
462       result = 0;
463     }
464   else
465     {
466       x = HIGH_U32_FROM_PTR (y);
467       result = 32;
468     }
469 #else
470   if (*y & 0x00000000FFFFFFFF)
471     {
472       x = LOW_U32_FROM_PTR (y);
473       result = 0;
474     }
475   else
476     {
477       x = HIGH_U32_FROM_PTR (y);
478       result = 32;
479     }
480 #endif /* USE_L */
481 #endif /* USE_LL */
482 
483   if (!(x & 0xFFFF))
484     {
485       x = bitSection (x, 0xFFFF0000, 16);
486       result += 16;
487     }
488   if (!(x & 0xFF))
489     {
490       x = bitSection (x, 0xFF00, 8);
491       result += 8;
492     }
493   if (!(x & 0xF))
494     {
495       x = bitSection (x, 0xF0, 4);
496       result += 4;
497     }
498 
499   if (x & 0x1)
500     return result + 1;
501   else if (x & 0x2)
502     return result + 2;
503   else if (x & 0x4)
504     return result + 3;
505   else
506     return result + 4;
507 }
508 
509 int32_t
highestSetBitHighPrecision(uint64_t * arg,int32_t length)510 highestSetBitHighPrecision (uint64_t * arg, int32_t length)
511 {
512   int32_t highBit;
513 
514   while (--length >= 0)
515     {
516       highBit = highestSetBit (arg + length);
517       if (highBit)
518         return highBit + 64 * length;
519     }
520 
521   return 0;
522 }
523 
524 int32_t
lowestSetBitHighPrecision(uint64_t * arg,int32_t length)525 lowestSetBitHighPrecision (uint64_t * arg, int32_t length)
526 {
527   int32_t lowBit, index = -1;
528 
529   while (++index < length)
530     {
531       lowBit = lowestSetBit (arg + index);
532       if (lowBit)
533         return lowBit + 64 * index;
534     }
535 
536   return 0;
537 }
538 
539 int32_t
compareHighPrecision(uint64_t * arg1,int32_t length1,uint64_t * arg2,int32_t length2)540 compareHighPrecision (uint64_t * arg1, int32_t length1, uint64_t * arg2, int32_t length2)
541 {
542   while (--length1 >= 0 && arg1[length1] == 0);
543   while (--length2 >= 0 && arg2[length2] == 0);
544 
545   if (length1 > length2)
546     return 1;
547   else if (length1 < length2)
548     return -1;
549   else if (length1 > -1)
550     {
551       do
552         {
553           if (arg1[length1] > arg2[length1])
554             return 1;
555           else if (arg1[length1] < arg2[length1])
556             return -1;
557         }
558       while (--length1 >= 0);
559     }
560 
561   return 0;
562 }
563 
564 jdouble
toDoubleHighPrecision(uint64_t * arg,int32_t length)565 toDoubleHighPrecision (uint64_t * arg, int32_t length)
566 {
567   int32_t highBit;
568   uint64_t mantissa, test64;
569   uint32_t test;
570   jdouble result;
571 
572   while (length > 0 && arg[length - 1] == 0)
573     --length;
574 
575   if (length == 0)
576     result = 0.0;
577   else if (length > 16)
578     {
579       DOUBLE_TO_LONGBITS (result) = EXPONENT_MASK;
580     }
581   else if (length == 1)
582     {
583       highBit = highestSetBit (arg);
584       if (highBit <= 53)
585         {
586           highBit = 53 - highBit;
587           mantissa = *arg << highBit;
588           DOUBLE_TO_LONGBITS (result) =
589             CREATE_DOUBLE_BITS (mantissa, -highBit);
590         }
591       else
592         {
593           highBit -= 53;
594           mantissa = *arg >> highBit;
595           DOUBLE_TO_LONGBITS (result) =
596             CREATE_DOUBLE_BITS (mantissa, highBit);
597 
598           /* perform rounding, round to even in case of tie */
599           test = (LOW_U32_FROM_PTR (arg) << (11 - highBit)) & 0x7FF;
600           if (test > 0x400 || ((test == 0x400) && (mantissa & 1)))
601             DOUBLE_TO_LONGBITS (result) = DOUBLE_TO_LONGBITS (result) + 1;
602         }
603     }
604   else
605     {
606       highBit = highestSetBit (arg + (--length));
607       if (highBit <= 53)
608         {
609           highBit = 53 - highBit;
610           if (highBit > 0)
611             {
612               mantissa =
613                 (arg[length] << highBit) | (arg[length - 1] >>
614                                             (64 - highBit));
615             }
616           else
617             {
618               mantissa = arg[length];
619             }
620           DOUBLE_TO_LONGBITS (result) =
621             CREATE_DOUBLE_BITS (mantissa, length * 64 - highBit);
622 
623           /* perform rounding, round to even in case of tie */
624           test64 = arg[--length] << highBit;
625           if (test64 > SIGN_MASK || ((test64 == SIGN_MASK) && (mantissa & 1)))
626             DOUBLE_TO_LONGBITS (result) = DOUBLE_TO_LONGBITS (result) + 1;
627           else if (test64 == SIGN_MASK)
628             {
629               while (--length >= 0)
630                 {
631                   if (arg[length] != 0)
632                     {
633                       DOUBLE_TO_LONGBITS (result) =
634                         DOUBLE_TO_LONGBITS (result) + 1;
635                       break;
636                     }
637                 }
638             }
639         }
640       else
641         {
642           highBit -= 53;
643           mantissa = arg[length] >> highBit;
644           DOUBLE_TO_LONGBITS (result) =
645             CREATE_DOUBLE_BITS (mantissa, length * 64 + highBit);
646 
647           /* perform rounding, round to even in case of tie */
648           test = (LOW_U32_FROM_PTR (arg + length) << (11 - highBit)) & 0x7FF;
649           if (test > 0x400 || ((test == 0x400) && (mantissa & 1)))
650             DOUBLE_TO_LONGBITS (result) = DOUBLE_TO_LONGBITS (result) + 1;
651           else if (test == 0x400)
652             {
653               do
654                 {
655                   if (arg[--length] != 0)
656                     {
657                       DOUBLE_TO_LONGBITS (result) =
658                         DOUBLE_TO_LONGBITS (result) + 1;
659                       break;
660                     }
661                 }
662               while (length > 0);
663             }
664         }
665     }
666 
667   return result;
668 }
669 
670 static uint64_t simpleMultiplyHighPrecision64(uint64_t* arg1, int32_t length, uint64_t arg2);
671 
672 int32_t
timesTenToTheEHighPrecision(uint64_t * result,int32_t length,jint e)673 timesTenToTheEHighPrecision (uint64_t * result, int32_t length, jint e)
674 {
675   /* assumes result can hold value */
676   uint64_t overflow;
677   int exp10 = e;
678 
679   if (e == 0)
680     return length;
681 
682   /* bad O(n) way of doing it, but simple */
683   /*
684      do {
685      overflow = simpleAppendDecimalDigitHighPrecision(result, length, 0);
686      if (overflow)
687      result[length++] = overflow;
688      } while (--e);
689    */
690   /* Replace the current implementaion which performs a
691    * "multiplication" by 10 e number of times with an actual
692    * multiplication. 10e19 is the largest exponent to the power of ten
693    * that will fit in a 64-bit integer, and 10e9 is the largest exponent to
694    * the power of ten that will fit in a 64-bit integer. Not sure where the
695    * break-even point is between an actual multiplication and a
696    * simpleAappendDecimalDigit() so just pick 10e3 as that point for
697    * now.
698    */
699   while (exp10 >= 19)
700     {
701       overflow = simpleMultiplyHighPrecision64 (result, length, TEN_E19);
702       if (overflow)
703         result[length++] = overflow;
704       exp10 -= 19;
705     }
706   while (exp10 >= 9)
707     {
708       overflow = simpleMultiplyHighPrecision (result, length, TEN_E9);
709       if (overflow)
710         result[length++] = overflow;
711       exp10 -= 9;
712     }
713   if (exp10 == 0)
714     return length;
715   else if (exp10 == 1)
716     {
717       overflow = simpleAppendDecimalDigitHighPrecision (result, length, 0);
718       if (overflow)
719         result[length++] = overflow;
720     }
721   else if (exp10 == 2)
722     {
723       overflow = simpleAppendDecimalDigitHighPrecision (result, length, 0);
724       if (overflow)
725         result[length++] = overflow;
726       overflow = simpleAppendDecimalDigitHighPrecision (result, length, 0);
727       if (overflow)
728         result[length++] = overflow;
729     }
730   else if (exp10 == 3)
731     {
732       overflow = simpleMultiplyHighPrecision (result, length, TEN_E3);
733       if (overflow)
734         result[length++] = overflow;
735     }
736   else if (exp10 == 4)
737     {
738       overflow = simpleMultiplyHighPrecision (result, length, TEN_E4);
739       if (overflow)
740         result[length++] = overflow;
741     }
742   else if (exp10 == 5)
743     {
744       overflow = simpleMultiplyHighPrecision (result, length, TEN_E5);
745       if (overflow)
746         result[length++] = overflow;
747     }
748   else if (exp10 == 6)
749     {
750       overflow = simpleMultiplyHighPrecision (result, length, TEN_E6);
751       if (overflow)
752         result[length++] = overflow;
753     }
754   else if (exp10 == 7)
755     {
756       overflow = simpleMultiplyHighPrecision (result, length, TEN_E7);
757       if (overflow)
758         result[length++] = overflow;
759     }
760   else if (exp10 == 8)
761     {
762       overflow = simpleMultiplyHighPrecision (result, length, TEN_E8);
763       if (overflow)
764         result[length++] = overflow;
765     }
766   return length;
767 }
768 
769 uint64_t
doubleMantissa(jdouble z)770 doubleMantissa (jdouble z)
771 {
772   uint64_t m = DOUBLE_TO_LONGBITS (z);
773 
774   if ((m & EXPONENT_MASK) != 0)
775     m = (m & MANTISSA_MASK) | NORMAL_MASK;
776   else
777     m = (m & MANTISSA_MASK);
778 
779   return m;
780 }
781 
782 int32_t
doubleExponent(jdouble z)783 doubleExponent (jdouble z)
784 {
785   /* assumes positive double */
786   int32_t k = HIGH_U32_FROM_VAR (z) >> 20;
787 
788   if (k)
789     k -= E_OFFSET;
790   else
791     k = 1 - E_OFFSET;
792 
793   return k;
794 }
795 
floatMantissa(jfloat z)796 uint32_t floatMantissa(jfloat z) {
797   uint32_t m = FLOAT_TO_INTBITS (z);
798 
799   if ((m & FLOAT_EXPONENT_MASK) != 0)
800     m = (m & FLOAT_MANTISSA_MASK) | FLOAT_NORMAL_MASK;
801   else
802     m = (m & FLOAT_MANTISSA_MASK);
803 
804   return m;
805 }
806 
807 int32_t
floatExponent(jfloat z)808 floatExponent (jfloat z)
809 {
810   /* assumes positive float */
811   int32_t k = FLOAT_TO_INTBITS (z) >> 23;
812   if (k)
813     k -= FLOAT_E_OFFSET;
814   else
815     k = 1 - FLOAT_E_OFFSET;
816 
817   return k;
818 }
819 
820 /* Allow a 64-bit value in arg2 */
821 uint64_t
simpleMultiplyHighPrecision64(uint64_t * arg1,int32_t length,uint64_t arg2)822 simpleMultiplyHighPrecision64 (uint64_t * arg1, int32_t length, uint64_t arg2)
823 {
824   uint64_t intermediate, carry1, carry2, prod1, prod2, sum;
825   uint64_t* pArg1;
826   int32_t index;
827   uint32_t buf32;
828 
829   index = 0;
830   intermediate = 0;
831   pArg1 = arg1 + index;
832   carry1 = carry2 = 0;
833 
834   do
835     {
836       if ((*pArg1 != 0) || (intermediate != 0))
837         {
838           prod1 =
839             static_cast<uint64_t>(LOW_U32_FROM_VAR (arg2)) * static_cast<uint64_t>(LOW_U32_FROM_PTR (pArg1));
840           sum = intermediate + prod1;
841           if ((sum < prod1) || (sum < intermediate))
842             {
843               carry1 = 1;
844             }
845           else
846             {
847               carry1 = 0;
848             }
849           prod1 =
850             static_cast<uint64_t>(LOW_U32_FROM_VAR (arg2)) * static_cast<uint64_t>(HIGH_U32_FROM_PTR (pArg1));
851           prod2 =
852             static_cast<uint64_t>(HIGH_U32_FROM_VAR (arg2)) * static_cast<uint64_t>(LOW_U32_FROM_PTR (pArg1));
853           intermediate = carry2 + HIGH_IN_U64 (sum) + prod1 + prod2;
854           if ((intermediate < prod1) || (intermediate < prod2))
855             {
856               carry2 = 1;
857             }
858           else
859             {
860               carry2 = 0;
861             }
862           LOW_U32_FROM_PTR (pArg1) = LOW_U32_FROM_VAR (sum);
863           buf32 = HIGH_U32_FROM_PTR (pArg1);
864           HIGH_U32_FROM_PTR (pArg1) = LOW_U32_FROM_VAR (intermediate);
865           intermediate = carry1 + HIGH_IN_U64 (intermediate)
866             + static_cast<uint64_t>(HIGH_U32_FROM_VAR (arg2)) * static_cast<uint64_t>(buf32);
867         }
868       pArg1++;
869     }
870   while (++index < length);
871   return intermediate;
872 }
873