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