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 #include <stdlib.h>
18 #include <math.h>
19
20 #include "commonDblParce.h"
21 #include "cbigint.h"
22
23
24 /* ************************* Defines ************************* */
25 #if defined(LINUX) || defined(FREEBSD)
26 #define USE_LL
27 #endif
28
29 #define LOW_I32_FROM_VAR(u64) LOW_I32_FROM_LONG64(u64)
30 #define LOW_I32_FROM_PTR(u64ptr) LOW_I32_FROM_LONG64_PTR(u64ptr)
31 #define HIGH_I32_FROM_VAR(u64) HIGH_I32_FROM_LONG64(u64)
32 #define HIGH_I32_FROM_PTR(u64ptr) HIGH_I32_FROM_LONG64_PTR(u64ptr)
33
34 #define MAX_ACCURACY_WIDTH 17
35
36 #define DEFAULT_WIDTH MAX_ACCURACY_WIDTH
37
38 #if defined(USE_LL)
39 #define INFINITE_LONGBITS (0x7FF0000000000000LL)
40 #else
41 #if defined(USE_L)
42 #define INFINITE_LONGBITS (0x7FF0000000000000L)
43 #else
44 #define INFINITE_LONGBITS (0x7FF0000000000000)
45 #endif /* USE_L */
46 #endif /* USE_LL */
47
48 #define MINIMUM_LONGBITS (0x1)
49
50 #if defined(USE_LL)
51 #define MANTISSA_MASK (0x000FFFFFFFFFFFFFLL)
52 #define EXPONENT_MASK (0x7FF0000000000000LL)
53 #define NORMAL_MASK (0x0010000000000000LL)
54 #else
55 #if defined(USE_L)
56 #define MANTISSA_MASK (0x000FFFFFFFFFFFFFL)
57 #define EXPONENT_MASK (0x7FF0000000000000L)
58 #define NORMAL_MASK (0x0010000000000000L)
59 #else
60 #define MANTISSA_MASK (0x000FFFFFFFFFFFFF)
61 #define EXPONENT_MASK (0x7FF0000000000000)
62 #define NORMAL_MASK (0x0010000000000000)
63 #endif /* USE_L */
64 #endif /* USE_LL */
65
66 #define DOUBLE_TO_LONGBITS(dbl) (*((U_64 *)(&dbl)))
67
68 /* Keep a count of the number of times we decrement and increment to
69 * approximate the double, and attempt to detect the case where we
70 * could potentially toggle back and forth between decrementing and
71 * incrementing. It is possible for us to be stuck in the loop when
72 * incrementing by one or decrementing by one may exceed or stay below
73 * the value that we are looking for. In this case, just break out of
74 * the loop if we toggle between incrementing and decrementing for more
75 * than twice.
76 */
77 #define INCREMENT_DOUBLE(_x, _decCount, _incCount) \
78 { \
79 ++DOUBLE_TO_LONGBITS(_x); \
80 _incCount++; \
81 if( (_incCount > 2) && (_decCount > 2) ) { \
82 if( _decCount > _incCount ) { \
83 DOUBLE_TO_LONGBITS(_x) += _decCount - _incCount; \
84 } else if( _incCount > _decCount ) { \
85 DOUBLE_TO_LONGBITS(_x) -= _incCount - _decCount; \
86 } \
87 break; \
88 } \
89 }
90 #define DECREMENT_DOUBLE(_x, _decCount, _incCount) \
91 { \
92 --DOUBLE_TO_LONGBITS(_x); \
93 _decCount++; \
94 if( (_incCount > 2) && (_decCount > 2) ) { \
95 if( _decCount > _incCount ) { \
96 DOUBLE_TO_LONGBITS(_x) += _decCount - _incCount; \
97 } else if( _incCount > _decCount ) { \
98 DOUBLE_TO_LONGBITS(_x) -= _incCount - _decCount; \
99 } \
100 break; \
101 } \
102 }
103
104 #define allocateU64(x, n) if (!((x) = (U_64*) malloc((n) * sizeof(U_64)))) goto OutOfMemory;
105 #define release(r) if ((r)) free((r));
106
107 /* *********************************************************** */
108
109 /* ************************ local data ************************ */
110 static const jdouble tens[] = {
111 1.0,
112 1.0e1,
113 1.0e2,
114 1.0e3,
115 1.0e4,
116 1.0e5,
117 1.0e6,
118 1.0e7,
119 1.0e8,
120 1.0e9,
121 1.0e10,
122 1.0e11,
123 1.0e12,
124 1.0e13,
125 1.0e14,
126 1.0e15,
127 1.0e16,
128 1.0e17,
129 1.0e18,
130 1.0e19,
131 1.0e20,
132 1.0e21,
133 1.0e22
134 };
135 /* *********************************************************** */
136
137 /* ************** private function declarations ************** */
138 static U_64 dblparse_shiftRight64 (U_64 * lp, volatile int mbe);
139
140 static jdouble createDouble1 (JNIEnv * env, U_64 * f, IDATA length, jint e);
141 static jdouble doubleAlgorithm (JNIEnv * env, U_64 * f, IDATA length, jint e,
142 jdouble z);
143 /* *********************************************************** */
144
145 #define tenToTheE(e) (*(tens + (e)))
146 #define LOG5_OF_TWO_TO_THE_N 23
147
148 #define sizeOfTenToTheE(e) (((e) / 19) + 1)
149
150 jdouble
createDouble(JNIEnv * env,const char * s,jint e)151 createDouble (JNIEnv * env, const char *s, jint e)
152 {
153 /* assumes s is a null terminated string with at least one
154 * character in it */
155 U_64 def[DEFAULT_WIDTH];
156 U_64 defBackup[DEFAULT_WIDTH];
157 U_64 *f, *fNoOverflow, *g, *tempBackup;
158 U_32 overflow;
159 jdouble result;
160 IDATA index = 1;
161 int unprocessedDigits = 0;
162
163 f = def;
164 fNoOverflow = defBackup;
165 *f = 0;
166 tempBackup = g = 0;
167 do
168 {
169 if (*s >= '0' && *s <= '9')
170 {
171 /* Make a back up of f before appending, so that we can
172 * back out of it if there is no more room, i.e. index >
173 * MAX_ACCURACY_WIDTH.
174 */
175 memcpy (fNoOverflow, f, sizeof (U_64) * index);
176 overflow =
177 simpleAppendDecimalDigitHighPrecision (f, index, *s - '0');
178 if (overflow)
179 {
180 f[index++] = overflow;
181 /* There is an overflow, but there is no more room
182 * to store the result. We really only need the top 52
183 * bits anyway, so we must back out of the overflow,
184 * and ignore the rest of the string.
185 */
186 if (index >= MAX_ACCURACY_WIDTH)
187 {
188 index--;
189 memcpy (f, fNoOverflow, sizeof (U_64) * index);
190 break;
191 }
192 if (tempBackup)
193 {
194 fNoOverflow = tempBackup;
195 }
196 }
197 }
198 else
199 index = -1;
200 }
201 while (index > 0 && *(++s) != '\0');
202
203 /* We've broken out of the parse loop either because we've reached
204 * the end of the string or we've overflowed the maximum accuracy
205 * limit of a double. If we still have unprocessed digits in the
206 * given string, then there are three possible results:
207 * 1. (unprocessed digits + e) == 0, in which case we simply
208 * convert the existing bits that are already parsed
209 * 2. (unprocessed digits + e) < 0, in which case we simply
210 * convert the existing bits that are already parsed along
211 * with the given e
212 * 3. (unprocessed digits + e) > 0 indicates that the value is
213 * simply too big to be stored as a double, so return Infinity
214 */
215 if ((unprocessedDigits = strlen (s)) > 0)
216 {
217 e += unprocessedDigits;
218 if (index > -1)
219 {
220 if (e == 0)
221 result = toDoubleHighPrecision (f, index);
222 else if (e < 0)
223 result = createDouble1 (env, f, index, e);
224 else
225 {
226 DOUBLE_TO_LONGBITS (result) = INFINITE_LONGBITS;
227 }
228 }
229 else
230 {
231 LOW_I32_FROM_VAR (result) = -1;
232 HIGH_I32_FROM_VAR (result) = -1;
233 }
234 }
235 else
236 {
237 if (index > -1)
238 {
239 if (e == 0)
240 result = toDoubleHighPrecision (f, index);
241 else
242 result = createDouble1 (env, f, index, e);
243 }
244 else
245 {
246 LOW_I32_FROM_VAR (result) = -1;
247 HIGH_I32_FROM_VAR (result) = -1;
248 }
249 }
250
251 return result;
252
253 }
254
255 jdouble
createDouble1(JNIEnv * env,U_64 * f,IDATA length,jint e)256 createDouble1 (JNIEnv * env, U_64 * f, IDATA length, jint e)
257 {
258 IDATA numBits;
259 jdouble result;
260
261 #define APPROX_MIN_MAGNITUDE -309
262
263 #define APPROX_MAX_MAGNITUDE 309
264
265 numBits = highestSetBitHighPrecision (f, length) + 1;
266 numBits -= lowestSetBitHighPrecision (f, length);
267 if (numBits < 54 && e >= 0 && e < LOG5_OF_TWO_TO_THE_N)
268 {
269 return toDoubleHighPrecision (f, length) * tenToTheE (e);
270 }
271 else if (numBits < 54 && e < 0 && (-e) < LOG5_OF_TWO_TO_THE_N)
272 {
273 return toDoubleHighPrecision (f, length) / tenToTheE (-e);
274 }
275 else if (e >= 0 && e < APPROX_MAX_MAGNITUDE)
276 {
277 result = toDoubleHighPrecision (f, length) * pow (10.0, e);
278 }
279 else if (e >= APPROX_MAX_MAGNITUDE)
280 {
281 /* Convert the partial result to make sure that the
282 * non-exponential part is not zero. This check fixes the case
283 * where the user enters 0.0e309! */
284 result = toDoubleHighPrecision (f, length);
285 /* Don't go straight to zero as the fact that x*0 = 0 independent of x might
286 cause the algorithm to produce an incorrect result. Instead try the min value
287 first and let it fall to zero if need be. */
288
289 if (result == 0.0)
290 {
291 DOUBLE_TO_LONGBITS (result) = MINIMUM_LONGBITS;
292 }
293 else
294 {
295 DOUBLE_TO_LONGBITS (result) = INFINITE_LONGBITS;
296 return result;
297 }
298 }
299 else if (e > APPROX_MIN_MAGNITUDE)
300 {
301 result = toDoubleHighPrecision (f, length) / pow (10.0, -e);
302 }
303
304 if (e <= APPROX_MIN_MAGNITUDE)
305 {
306
307 result = toDoubleHighPrecision (f, length) * pow (10.0, e + 52);
308 result = result * pow (10.0, -52);
309
310 }
311
312 /* Don't go straight to zero as the fact that x*0 = 0 independent of x might
313 cause the algorithm to produce an incorrect result. Instead try the min value
314 first and let it fall to zero if need be. */
315
316 if (result == 0.0)
317
318 DOUBLE_TO_LONGBITS (result) = MINIMUM_LONGBITS;
319
320 return doubleAlgorithm (env, f, length, e, result);
321 }
322
323 static U_64
dblparse_shiftRight64(U_64 * lp,volatile int mbe)324 dblparse_shiftRight64 (U_64 * lp, volatile int mbe)
325 {
326 U_64 b1Value = 0;
327 U_32 hi = HIGH_U32_FROM_LONG64_PTR (lp);
328 U_32 lo = LOW_U32_FROM_LONG64_PTR (lp);
329 int srAmt;
330
331 if (mbe == 0)
332 return 0;
333 if (mbe >= 128)
334 {
335 HIGH_U32_FROM_LONG64_PTR (lp) = 0;
336 LOW_U32_FROM_LONG64_PTR (lp) = 0;
337 return 0;
338 }
339
340 /* Certain platforms do not handle de-referencing a 64-bit value
341 * from a pointer on the stack correctly (e.g. MVL-hh/XScale)
342 * because the pointer may not be properly aligned, so we'll have
343 * to handle two 32-bit chunks. */
344 if (mbe < 32)
345 {
346 LOW_U32_FROM_LONG64 (b1Value) = 0;
347 HIGH_U32_FROM_LONG64 (b1Value) = lo << (32 - mbe);
348 LOW_U32_FROM_LONG64_PTR (lp) = (hi << (32 - mbe)) | (lo >> mbe);
349 HIGH_U32_FROM_LONG64_PTR (lp) = hi >> mbe;
350 }
351 else if (mbe == 32)
352 {
353 LOW_U32_FROM_LONG64 (b1Value) = 0;
354 HIGH_U32_FROM_LONG64 (b1Value) = lo;
355 LOW_U32_FROM_LONG64_PTR (lp) = hi;
356 HIGH_U32_FROM_LONG64_PTR (lp) = 0;
357 }
358 else if (mbe < 64)
359 {
360 srAmt = mbe - 32;
361 LOW_U32_FROM_LONG64 (b1Value) = lo << (32 - srAmt);
362 HIGH_U32_FROM_LONG64 (b1Value) = (hi << (32 - srAmt)) | (lo >> srAmt);
363 LOW_U32_FROM_LONG64_PTR (lp) = hi >> srAmt;
364 HIGH_U32_FROM_LONG64_PTR (lp) = 0;
365 }
366 else if (mbe == 64)
367 {
368 LOW_U32_FROM_LONG64 (b1Value) = lo;
369 HIGH_U32_FROM_LONG64 (b1Value) = hi;
370 LOW_U32_FROM_LONG64_PTR (lp) = 0;
371 HIGH_U32_FROM_LONG64_PTR (lp) = 0;
372 }
373 else if (mbe < 96)
374 {
375 srAmt = mbe - 64;
376 b1Value = *lp;
377 HIGH_U32_FROM_LONG64_PTR (lp) = 0;
378 LOW_U32_FROM_LONG64_PTR (lp) = 0;
379 LOW_U32_FROM_LONG64 (b1Value) >>= srAmt;
380 LOW_U32_FROM_LONG64 (b1Value) |= (hi << (32 - srAmt));
381 HIGH_U32_FROM_LONG64 (b1Value) >>= srAmt;
382 }
383 else if (mbe == 96)
384 {
385 LOW_U32_FROM_LONG64 (b1Value) = hi;
386 HIGH_U32_FROM_LONG64 (b1Value) = 0;
387 HIGH_U32_FROM_LONG64_PTR (lp) = 0;
388 LOW_U32_FROM_LONG64_PTR (lp) = 0;
389 }
390 else
391 {
392 LOW_U32_FROM_LONG64 (b1Value) = hi >> (mbe - 96);
393 HIGH_U32_FROM_LONG64 (b1Value) = 0;
394 HIGH_U32_FROM_LONG64_PTR (lp) = 0;
395 LOW_U32_FROM_LONG64_PTR (lp) = 0;
396 }
397
398 return b1Value;
399 }
400
401 #if defined(WIN32)
402 /* disable global optimizations on the microsoft compiler for the
403 * doubleAlgorithm function otherwise it won't compile */
404 #pragma optimize("g",off)
405 #endif
406
407
408 /* The algorithm for the function doubleAlgorithm() below can be found
409 * in:
410 *
411 * "How to Read Floating-Point Numbers Accurately", William D.
412 * Clinger, Proceedings of the ACM SIGPLAN '90 Conference on
413 * Programming Language Design and Implementation, June 20-22,
414 * 1990, pp. 92-101.
415 *
416 * There is a possibility that the function will end up in an endless
417 * loop if the given approximating floating-point number (a very small
418 * floating-point whose value is very close to zero) straddles between
419 * two approximating integer values. We modified the algorithm slightly
420 * to detect the case where it oscillates back and forth between
421 * incrementing and decrementing the floating-point approximation. It
422 * is currently set such that if the oscillation occurs more than twice
423 * then return the original approximation.
424 */
425 static jdouble
doubleAlgorithm(JNIEnv * env,U_64 * f,IDATA length,jint e,jdouble z)426 doubleAlgorithm (JNIEnv * env, U_64 * f, IDATA length, jint e, jdouble z)
427 {
428 U_64 m;
429 IDATA k, comparison, comparison2;
430 U_64 *x, *y, *D, *D2;
431 IDATA xLength, yLength, DLength, D2Length, decApproxCount, incApproxCount;
432
433 x = y = D = D2 = 0;
434 xLength = yLength = DLength = D2Length = 0;
435 decApproxCount = incApproxCount = 0;
436
437 do
438 {
439 m = doubleMantissa (z);
440 k = doubleExponent (z);
441
442 if (x && x != f)
443 free(x);
444
445 release (y);
446 release (D);
447 release (D2);
448
449 if (e >= 0 && k >= 0)
450 {
451 xLength = sizeOfTenToTheE (e) + length;
452 allocateU64 (x, xLength);
453 memset (x + length, 0, sizeof (U_64) * (xLength - length));
454 memcpy (x, f, sizeof (U_64) * length);
455 timesTenToTheEHighPrecision (x, xLength, e);
456
457 yLength = (k >> 6) + 2;
458 allocateU64 (y, yLength);
459 memset (y + 1, 0, sizeof (U_64) * (yLength - 1));
460 *y = m;
461 simpleShiftLeftHighPrecision (y, yLength, k);
462 }
463 else if (e >= 0)
464 {
465 xLength = sizeOfTenToTheE (e) + length + ((-k) >> 6) + 1;
466 allocateU64 (x, xLength);
467 memset (x + length, 0, sizeof (U_64) * (xLength - length));
468 memcpy (x, f, sizeof (U_64) * length);
469 timesTenToTheEHighPrecision (x, xLength, e);
470 simpleShiftLeftHighPrecision (x, xLength, -k);
471
472 yLength = 1;
473 allocateU64 (y, 1);
474 *y = m;
475 }
476 else if (k >= 0)
477 {
478 xLength = length;
479 x = f;
480
481 yLength = sizeOfTenToTheE (-e) + 2 + (k >> 6);
482 allocateU64 (y, yLength);
483 memset (y + 1, 0, sizeof (U_64) * (yLength - 1));
484 *y = m;
485 timesTenToTheEHighPrecision (y, yLength, -e);
486 simpleShiftLeftHighPrecision (y, yLength, k);
487 }
488 else
489 {
490 xLength = length + ((-k) >> 6) + 1;
491 allocateU64 (x, xLength);
492 memset (x + length, 0, sizeof (U_64) * (xLength - length));
493 memcpy (x, f, sizeof (U_64) * length);
494 simpleShiftLeftHighPrecision (x, xLength, -k);
495
496 yLength = sizeOfTenToTheE (-e) + 1;
497 allocateU64 (y, yLength);
498 memset (y + 1, 0, sizeof (U_64) * (yLength - 1));
499 *y = m;
500 timesTenToTheEHighPrecision (y, yLength, -e);
501 }
502
503 comparison = compareHighPrecision (x, xLength, y, yLength);
504 if (comparison > 0)
505 { /* x > y */
506 DLength = xLength;
507 allocateU64 (D, DLength);
508 memcpy (D, x, DLength * sizeof (U_64));
509 subtractHighPrecision (D, DLength, y, yLength);
510 }
511 else if (comparison)
512 { /* y > x */
513 DLength = yLength;
514 allocateU64 (D, DLength);
515 memcpy (D, y, DLength * sizeof (U_64));
516 subtractHighPrecision (D, DLength, x, xLength);
517 }
518 else
519 { /* y == x */
520 DLength = 1;
521 allocateU64 (D, 1);
522 *D = 0;
523 }
524
525 D2Length = DLength + 1;
526 allocateU64 (D2, D2Length);
527 m <<= 1;
528 multiplyHighPrecision (D, DLength, &m, 1, D2, D2Length);
529 m >>= 1;
530
531 comparison2 = compareHighPrecision (D2, D2Length, y, yLength);
532 if (comparison2 < 0)
533 {
534 if (comparison < 0 && m == NORMAL_MASK)
535 {
536 simpleShiftLeftHighPrecision (D2, D2Length, 1);
537 if (compareHighPrecision (D2, D2Length, y, yLength) > 0)
538 {
539 DECREMENT_DOUBLE (z, decApproxCount, incApproxCount);
540 }
541 else
542 {
543 break;
544 }
545 }
546 else
547 {
548 break;
549 }
550 }
551 else if (comparison2 == 0)
552 {
553 if ((LOW_U32_FROM_VAR (m) & 1) == 0)
554 {
555 if (comparison < 0 && m == NORMAL_MASK)
556 {
557 DECREMENT_DOUBLE (z, decApproxCount, incApproxCount);
558 }
559 else
560 {
561 break;
562 }
563 }
564 else if (comparison < 0)
565 {
566 DECREMENT_DOUBLE (z, decApproxCount, incApproxCount);
567 break;
568 }
569 else
570 {
571 INCREMENT_DOUBLE (z, decApproxCount, incApproxCount);
572 break;
573 }
574 }
575 else if (comparison < 0)
576 {
577 DECREMENT_DOUBLE (z, decApproxCount, incApproxCount);
578 }
579 else
580 {
581 if (DOUBLE_TO_LONGBITS (z) == INFINITE_LONGBITS)
582 break;
583 INCREMENT_DOUBLE (z, decApproxCount, incApproxCount);
584 }
585 }
586 while (1);
587
588 if (x && x != f)
589 free(x);
590 release (y);
591 release (D);
592 release (D2);
593 return z;
594
595 OutOfMemory:
596 if (x && x != f)
597 free(x);
598 release (y);
599 release (y);
600 release (D);
601 release (D2);
602
603 DOUBLE_TO_LONGBITS (z) = -2;
604
605 return z;
606 }
607