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