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#define HIGH_IN_U64(u64) ((u64) >> 32)
22#define LOW_IN_U64(u64) ((u64) & 0x00000000FFFFFFFFLL)
23
24#define TEN_E1 (0xALL)
25#define TEN_E2 (0x64LL)
26#define TEN_E3 (0x3E8LL)
27#define TEN_E4 (0x2710LL)
28#define TEN_E5 (0x186A0LL)
29#define TEN_E6 (0xF4240LL)
30#define TEN_E7 (0x989680LL)
31#define TEN_E8 (0x5F5E100LL)
32#define TEN_E9 (0x3B9ACA00LL)
33#define TEN_E19 (0x8AC7230489E80000LL)
34
35#define TIMES_TEN(x) (((x) << 3) + ((x) << 1))
36#define bitSection(x, mask, shift) (((x) & (mask)) >> (shift))
37#define CREATE_DOUBLE_BITS(normalizedM, e) (((normalizedM) & MANTISSA_MASK) | ((static_cast<uint64_t>((e) + E_OFFSET)) << 52))
38
39#define MANTISSA_MASK (0x000FFFFFFFFFFFFFLL)
40#define EXPONENT_MASK (0x7FF0000000000000LL)
41#define NORMAL_MASK (0x0010000000000000LL)
42#define SIGN_MASK (0x8000000000000000LL)
43
44#define E_OFFSET (1075)
45
46#define FLOAT_MANTISSA_MASK (0x007FFFFF)
47#define FLOAT_EXPONENT_MASK (0x7F800000)
48#define FLOAT_NORMAL_MASK   (0x00800000)
49#define FLOAT_E_OFFSET (150)
50
51int32_t
52simpleAddHighPrecision (uint64_t * arg1, int32_t length, uint64_t arg2)
53{
54  /* assumes length > 0 */
55  int32_t index = 1;
56
57  *arg1 += arg2;
58  if (arg2 <= *arg1)
59    return 0;
60  else if (length == 1)
61    return 1;
62
63  while (++arg1[index] == 0 && ++index < length) {
64  }
65
66  return index == length;
67}
68
69int32_t
70addHighPrecision (uint64_t * arg1, int32_t length1, uint64_t * arg2, int32_t length2)
71{
72  /* addition is limited by length of arg1 as it this function is
73   * storing the result in arg1 */
74  /* fix for cc (GCC) 3.2 20020903 (Red Hat Linux 8.0 3.2-7): code generated does not
75   * do the temp1 + temp2 + carry addition correct.  carry is 64 bit because gcc has
76   * subtle issues when you mix 64 / 32 bit maths. */
77  uint64_t temp1, temp2, temp3;     /* temporary variables to help the SH-4, and gcc */
78  uint64_t carry;
79  int32_t index;
80
81  if (length1 == 0 || length2 == 0)
82    {
83      return 0;
84    }
85  else if (length1 < length2)
86    {
87      length2 = length1;
88    }
89
90  carry = 0;
91  index = 0;
92  do
93    {
94      temp1 = arg1[index];
95      temp2 = arg2[index];
96      temp3 = temp1 + temp2;
97      arg1[index] = temp3 + carry;
98      if (arg2[index] < arg1[index])
99        carry = 0;
100      else if (arg2[index] != arg1[index])
101        carry = 1;
102    }
103  while (++index < length2);
104  if (!carry)
105    return 0;
106  else if (index == length1)
107    return 1;
108
109  while (++arg1[index] == 0 && ++index < length1) {
110  }
111
112  return index == length1;
113}
114
115void
116subtractHighPrecision (uint64_t * arg1, int32_t length1, uint64_t * arg2, int32_t length2)
117{
118  /* assumes arg1 > arg2 */
119  int32_t index;
120  for (index = 0; index < length1; ++index)
121    arg1[index] = ~arg1[index];
122  simpleAddHighPrecision (arg1, length1, 1);
123
124  while (length2 > 0 && arg2[length2 - 1] == 0)
125    --length2;
126
127  addHighPrecision (arg1, length1, arg2, length2);
128
129  for (index = 0; index < length1; ++index)
130    arg1[index] = ~arg1[index];
131  simpleAddHighPrecision (arg1, length1, 1);
132}
133
134static uint32_t simpleMultiplyHighPrecision(uint64_t* arg1, int32_t length, uint64_t arg2) {
135  /* assumes arg2 only holds 32 bits of information */
136  uint64_t product;
137  int32_t index;
138
139  index = 0;
140  product = 0;
141
142  do
143    {
144      product =
145        HIGH_IN_U64 (product) + arg2 * LOW_U32_FROM_PTR (arg1 + index);
146      LOW_U32_FROM_PTR (arg1 + index) = LOW_U32_FROM_VAR (product);
147      product =
148        HIGH_IN_U64 (product) + arg2 * HIGH_U32_FROM_PTR (arg1 + index);
149      HIGH_U32_FROM_PTR (arg1 + index) = LOW_U32_FROM_VAR (product);
150    }
151  while (++index < length);
152
153  return HIGH_U32_FROM_VAR (product);
154}
155
156static void
157simpleMultiplyAddHighPrecision (uint64_t * arg1, int32_t length, uint64_t arg2,
158                                uint32_t * result)
159{
160  /* Assumes result can hold the product and arg2 only holds 32 bits
161     of information */
162  uint64_t product;
163  int32_t index, resultIndex;
164
165  index = resultIndex = 0;
166  product = 0;
167
168  do
169    {
170      product =
171        HIGH_IN_U64 (product) + result[resultIndex] +
172        arg2 * LOW_U32_FROM_PTR (arg1 + index);
173      result[resultIndex] = LOW_U32_FROM_VAR (product);
174      ++resultIndex;
175      product =
176        HIGH_IN_U64 (product) + result[resultIndex] +
177        arg2 * HIGH_U32_FROM_PTR (arg1 + index);
178      result[resultIndex] = LOW_U32_FROM_VAR (product);
179      ++resultIndex;
180    }
181  while (++index < length);
182
183  result[resultIndex] += HIGH_U32_FROM_VAR (product);
184  if (result[resultIndex] < HIGH_U32_FROM_VAR (product))
185    {
186      /* must be careful with ++ operator and macro expansion */
187      ++resultIndex;
188      while (++result[resultIndex] == 0)
189        ++resultIndex;
190    }
191}
192
193void
194multiplyHighPrecision (uint64_t * arg1, int32_t length1, uint64_t * arg2, int32_t length2,
195                       uint64_t * result, int32_t length)
196{
197  /* assumes result is large enough to hold product */
198  uint64_t* temp;
199  uint32_t* resultIn32;
200  int32_t count, index;
201
202  if (length1 < length2)
203    {
204      temp = arg1;
205      arg1 = arg2;
206      arg2 = temp;
207      count = length1;
208      length1 = length2;
209      length2 = count;
210    }
211
212  memset (result, 0, sizeof (uint64_t) * length);
213
214  /* length1 > length2 */
215  resultIn32 = reinterpret_cast<uint32_t*>(result);
216  index = -1;
217  for (count = 0; count < length2; ++count)
218    {
219      simpleMultiplyAddHighPrecision (arg1, length1, LOW_IN_U64 (arg2[count]),
220                                      resultIn32 + (++index));
221      simpleMultiplyAddHighPrecision(arg1, length1, HIGH_IN_U64(arg2[count]), resultIn32 + (++index));
222    }
223}
224
225uint32_t
226simpleAppendDecimalDigitHighPrecision (uint64_t * arg1, int32_t length, uint64_t digit)
227{
228  /* assumes digit is less than 32 bits */
229  uint64_t arg;
230  int32_t index = 0;
231
232  digit <<= 32;
233  do
234    {
235      arg = LOW_IN_U64 (arg1[index]);
236      digit = HIGH_IN_U64 (digit) + TIMES_TEN (arg);
237      LOW_U32_FROM_PTR (arg1 + index) = LOW_U32_FROM_VAR (digit);
238
239      arg = HIGH_IN_U64 (arg1[index]);
240      digit = HIGH_IN_U64 (digit) + TIMES_TEN (arg);
241      HIGH_U32_FROM_PTR (arg1 + index) = LOW_U32_FROM_VAR (digit);
242    }
243  while (++index < length);
244
245  return HIGH_U32_FROM_VAR (digit);
246}
247
248void
249simpleShiftLeftHighPrecision (uint64_t * arg1, int32_t length, int32_t arg2)
250{
251  /* assumes length > 0 */
252  int32_t index, offset;
253  if (arg2 >= 64)
254    {
255      offset = arg2 >> 6;
256      index = length;
257
258      while (--index - offset >= 0)
259        arg1[index] = arg1[index - offset];
260      do
261        {
262          arg1[index] = 0;
263        }
264      while (--index >= 0);
265
266      arg2 &= 0x3F;
267    }
268
269  if (arg2 == 0)
270    return;
271  while (--length > 0)
272    {
273      arg1[length] = arg1[length] << arg2 | arg1[length - 1] >> (64 - arg2);
274    }
275  *arg1 <<= arg2;
276}
277
278int32_t
279highestSetBit (uint64_t * y)
280{
281  uint32_t x;
282  int32_t result;
283
284  if (*y == 0)
285    return 0;
286
287  if (*y & 0xFFFFFFFF00000000LL)
288    {
289      x = HIGH_U32_FROM_PTR (y);
290      result = 32;
291    }
292  else
293    {
294      x = LOW_U32_FROM_PTR (y);
295      result = 0;
296    }
297
298  if (x & 0xFFFF0000)
299    {
300      x = bitSection (x, 0xFFFF0000, 16);
301      result += 16;
302    }
303  if (x & 0xFF00)
304    {
305      x = bitSection (x, 0xFF00, 8);
306      result += 8;
307    }
308  if (x & 0xF0)
309    {
310      x = bitSection (x, 0xF0, 4);
311      result += 4;
312    }
313  if (x > 0x7)
314    return result + 4;
315  else if (x > 0x3)
316    return result + 3;
317  else if (x > 0x1)
318    return result + 2;
319  else
320    return result + 1;
321}
322
323int32_t
324lowestSetBit (uint64_t * y)
325{
326  uint32_t x;
327  int32_t result;
328
329  if (*y == 0)
330    return 0;
331
332  if (*y & 0x00000000FFFFFFFFLL)
333    {
334      x = LOW_U32_FROM_PTR (y);
335      result = 0;
336    }
337  else
338    {
339      x = HIGH_U32_FROM_PTR (y);
340      result = 32;
341    }
342
343  if (!(x & 0xFFFF))
344    {
345      x = bitSection (x, 0xFFFF0000, 16);
346      result += 16;
347    }
348  if (!(x & 0xFF))
349    {
350      x = bitSection (x, 0xFF00, 8);
351      result += 8;
352    }
353  if (!(x & 0xF))
354    {
355      x = bitSection (x, 0xF0, 4);
356      result += 4;
357    }
358
359  if (x & 0x1)
360    return result + 1;
361  else if (x & 0x2)
362    return result + 2;
363  else if (x & 0x4)
364    return result + 3;
365  else
366    return result + 4;
367}
368
369int32_t
370highestSetBitHighPrecision (uint64_t * arg, int32_t length)
371{
372  int32_t highBit;
373
374  while (--length >= 0)
375    {
376      highBit = highestSetBit (arg + length);
377      if (highBit)
378        return highBit + 64 * length;
379    }
380
381  return 0;
382}
383
384int32_t
385lowestSetBitHighPrecision (uint64_t * arg, int32_t length)
386{
387  int32_t lowBit, index = -1;
388
389  while (++index < length)
390    {
391      lowBit = lowestSetBit (arg + index);
392      if (lowBit)
393        return lowBit + 64 * index;
394    }
395
396  return 0;
397}
398
399int32_t
400compareHighPrecision (uint64_t * arg1, int32_t length1, uint64_t * arg2, int32_t length2)
401{
402  while (--length1 >= 0 && arg1[length1] == 0) {
403  }
404  while (--length2 >= 0 && arg2[length2] == 0) {
405  }
406
407  if (length1 > length2)
408    return 1;
409  else if (length1 < length2)
410    return -1;
411  else if (length1 > -1)
412    {
413      do
414        {
415          if (arg1[length1] > arg2[length1])
416            return 1;
417          else if (arg1[length1] < arg2[length1])
418            return -1;
419        }
420      while (--length1 >= 0);
421    }
422
423  return 0;
424}
425
426jdouble
427toDoubleHighPrecision (uint64_t * arg, int32_t length)
428{
429  int32_t highBit;
430  uint64_t mantissa, test64;
431  uint32_t test;
432  jdouble result;
433
434  while (length > 0 && arg[length - 1] == 0)
435    --length;
436
437  if (length == 0)
438    result = 0.0;
439  else if (length > 16)
440    {
441      DOUBLE_TO_LONGBITS (result) = EXPONENT_MASK;
442    }
443  else if (length == 1)
444    {
445      highBit = highestSetBit (arg);
446      if (highBit <= 53)
447        {
448          highBit = 53 - highBit;
449          mantissa = *arg << highBit;
450          DOUBLE_TO_LONGBITS (result) =
451            CREATE_DOUBLE_BITS (mantissa, -highBit);
452        }
453      else
454        {
455          highBit -= 53;
456          mantissa = *arg >> highBit;
457          DOUBLE_TO_LONGBITS (result) =
458            CREATE_DOUBLE_BITS (mantissa, highBit);
459
460          /* perform rounding, round to even in case of tie */
461          test = (LOW_U32_FROM_PTR (arg) << (11 - highBit)) & 0x7FF;
462          if (test > 0x400 || ((test == 0x400) && (mantissa & 1)))
463            DOUBLE_TO_LONGBITS (result) = DOUBLE_TO_LONGBITS (result) + 1;
464        }
465    }
466  else
467    {
468      highBit = highestSetBit (arg + (--length));
469      if (highBit <= 53)
470        {
471          highBit = 53 - highBit;
472          if (highBit > 0)
473            {
474              mantissa =
475                (arg[length] << highBit) | (arg[length - 1] >>
476                                            (64 - highBit));
477            }
478          else
479            {
480              mantissa = arg[length];
481            }
482          DOUBLE_TO_LONGBITS (result) =
483            CREATE_DOUBLE_BITS (mantissa, length * 64 - highBit);
484
485          /* perform rounding, round to even in case of tie */
486          test64 = arg[--length] << highBit;
487          if (test64 > SIGN_MASK || ((test64 == SIGN_MASK) && (mantissa & 1)))
488            DOUBLE_TO_LONGBITS (result) = DOUBLE_TO_LONGBITS (result) + 1;
489          else if (test64 == SIGN_MASK)
490            {
491              while (--length >= 0)
492                {
493                  if (arg[length] != 0)
494                    {
495                      DOUBLE_TO_LONGBITS (result) =
496                        DOUBLE_TO_LONGBITS (result) + 1;
497                      break;
498                    }
499                }
500            }
501        }
502      else
503        {
504          highBit -= 53;
505          mantissa = arg[length] >> highBit;
506          DOUBLE_TO_LONGBITS (result) =
507            CREATE_DOUBLE_BITS (mantissa, length * 64 + highBit);
508
509          /* perform rounding, round to even in case of tie */
510          test = (LOW_U32_FROM_PTR (arg + length) << (11 - highBit)) & 0x7FF;
511          if (test > 0x400 || ((test == 0x400) && (mantissa & 1)))
512            DOUBLE_TO_LONGBITS (result) = DOUBLE_TO_LONGBITS (result) + 1;
513          else if (test == 0x400)
514            {
515              do
516                {
517                  if (arg[--length] != 0)
518                    {
519                      DOUBLE_TO_LONGBITS (result) =
520                        DOUBLE_TO_LONGBITS (result) + 1;
521                      break;
522                    }
523                }
524              while (length > 0);
525            }
526        }
527    }
528
529  return result;
530}
531
532static uint64_t simpleMultiplyHighPrecision64(uint64_t* arg1, int32_t length, uint64_t arg2);
533
534int32_t
535timesTenToTheEHighPrecision (uint64_t * result, int32_t length, jint e)
536{
537  /* assumes result can hold value */
538  uint64_t overflow;
539  int exp10 = e;
540
541  if (e == 0)
542    return length;
543
544  /* bad O(n) way of doing it, but simple */
545  /*
546     do {
547     overflow = simpleAppendDecimalDigitHighPrecision(result, length, 0);
548     if (overflow)
549     result[length++] = overflow;
550     } while (--e);
551   */
552  /* Replace the current implementaion which performs a
553   * "multiplication" by 10 e number of times with an actual
554   * multiplication. 10e19 is the largest exponent to the power of ten
555   * that will fit in a 64-bit integer, and 10e9 is the largest exponent to
556   * the power of ten that will fit in a 64-bit integer. Not sure where the
557   * break-even point is between an actual multiplication and a
558   * simpleAappendDecimalDigit() so just pick 10e3 as that point for
559   * now.
560   */
561  while (exp10 >= 19)
562    {
563      overflow = simpleMultiplyHighPrecision64 (result, length, TEN_E19);
564      if (overflow)
565        result[length++] = overflow;
566      exp10 -= 19;
567    }
568  while (exp10 >= 9)
569    {
570      overflow = simpleMultiplyHighPrecision (result, length, TEN_E9);
571      if (overflow)
572        result[length++] = overflow;
573      exp10 -= 9;
574    }
575  if (exp10 == 0)
576    return length;
577  else if (exp10 == 1)
578    {
579      overflow = simpleAppendDecimalDigitHighPrecision (result, length, 0);
580      if (overflow)
581        result[length++] = overflow;
582    }
583  else if (exp10 == 2)
584    {
585      overflow = simpleAppendDecimalDigitHighPrecision (result, length, 0);
586      if (overflow)
587        result[length++] = overflow;
588      overflow = simpleAppendDecimalDigitHighPrecision (result, length, 0);
589      if (overflow)
590        result[length++] = overflow;
591    }
592  else if (exp10 == 3)
593    {
594      overflow = simpleMultiplyHighPrecision (result, length, TEN_E3);
595      if (overflow)
596        result[length++] = overflow;
597    }
598  else if (exp10 == 4)
599    {
600      overflow = simpleMultiplyHighPrecision (result, length, TEN_E4);
601      if (overflow)
602        result[length++] = overflow;
603    }
604  else if (exp10 == 5)
605    {
606      overflow = simpleMultiplyHighPrecision (result, length, TEN_E5);
607      if (overflow)
608        result[length++] = overflow;
609    }
610  else if (exp10 == 6)
611    {
612      overflow = simpleMultiplyHighPrecision (result, length, TEN_E6);
613      if (overflow)
614        result[length++] = overflow;
615    }
616  else if (exp10 == 7)
617    {
618      overflow = simpleMultiplyHighPrecision (result, length, TEN_E7);
619      if (overflow)
620        result[length++] = overflow;
621    }
622  else if (exp10 == 8)
623    {
624      overflow = simpleMultiplyHighPrecision (result, length, TEN_E8);
625      if (overflow)
626        result[length++] = overflow;
627    }
628  return length;
629}
630
631uint64_t
632doubleMantissa (jdouble z)
633{
634  uint64_t m = DOUBLE_TO_LONGBITS (z);
635
636  if ((m & EXPONENT_MASK) != 0)
637    m = (m & MANTISSA_MASK) | NORMAL_MASK;
638  else
639    m = (m & MANTISSA_MASK);
640
641  return m;
642}
643
644int32_t
645doubleExponent (jdouble z)
646{
647  /* assumes positive double */
648  int32_t k = HIGH_U32_FROM_VAR (z) >> 20;
649
650  if (k)
651    k -= E_OFFSET;
652  else
653    k = 1 - E_OFFSET;
654
655  return k;
656}
657
658uint32_t floatMantissa(jfloat z) {
659  uint32_t m = FLOAT_TO_INTBITS (z);
660
661  if ((m & FLOAT_EXPONENT_MASK) != 0)
662    m = (m & FLOAT_MANTISSA_MASK) | FLOAT_NORMAL_MASK;
663  else
664    m = (m & FLOAT_MANTISSA_MASK);
665
666  return m;
667}
668
669int32_t
670floatExponent (jfloat z)
671{
672  /* assumes positive float */
673  int32_t k = FLOAT_TO_INTBITS (z) >> 23;
674  if (k)
675    k -= FLOAT_E_OFFSET;
676  else
677    k = 1 - FLOAT_E_OFFSET;
678
679  return k;
680}
681
682/* Allow a 64-bit value in arg2 */
683uint64_t
684simpleMultiplyHighPrecision64 (uint64_t * arg1, int32_t length, uint64_t arg2)
685{
686  uint64_t intermediate, carry1, carry2, prod1, prod2, sum;
687  uint64_t* pArg1;
688  int32_t index;
689  uint32_t buf32;
690
691  index = 0;
692  intermediate = 0;
693  pArg1 = arg1 + index;
694  carry1 = carry2 = 0;
695
696  do
697    {
698      if ((*pArg1 != 0) || (intermediate != 0))
699        {
700          prod1 =
701            static_cast<uint64_t>(LOW_U32_FROM_VAR (arg2)) * static_cast<uint64_t>(LOW_U32_FROM_PTR (pArg1));
702          sum = intermediate + prod1;
703          if ((sum < prod1) || (sum < intermediate))
704            {
705              carry1 = 1;
706            }
707          else
708            {
709              carry1 = 0;
710            }
711          prod1 =
712            static_cast<uint64_t>(LOW_U32_FROM_VAR (arg2)) * static_cast<uint64_t>(HIGH_U32_FROM_PTR (pArg1));
713          prod2 =
714            static_cast<uint64_t>(HIGH_U32_FROM_VAR (arg2)) * static_cast<uint64_t>(LOW_U32_FROM_PTR (pArg1));
715          intermediate = carry2 + HIGH_IN_U64 (sum) + prod1 + prod2;
716          if ((intermediate < prod1) || (intermediate < prod2))
717            {
718              carry2 = 1;
719            }
720          else
721            {
722              carry2 = 0;
723            }
724          LOW_U32_FROM_PTR (pArg1) = LOW_U32_FROM_VAR (sum);
725          buf32 = HIGH_U32_FROM_PTR (pArg1);
726          HIGH_U32_FROM_PTR (pArg1) = LOW_U32_FROM_VAR (intermediate);
727          intermediate = carry1 + HIGH_IN_U64 (intermediate)
728            + static_cast<uint64_t>(HIGH_U32_FROM_VAR (arg2)) * static_cast<uint64_t>(buf32);
729        }
730      pArg1++;
731    }
732  while (++index < length);
733  return intermediate;
734}
735