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 <stdlib.h>
19 #include <string.h>
20 #include <math.h>
21 #include "JNIHelp.h"
22 #include "JniConstants.h"
23 #include "JniException.h"
24 #include "ScopedUtfChars.h"
25 #include "cbigint.h"
26
27 /* ************************* Defines ************************* */
28 #if defined(__linux__) || defined(__APPLE__)
29 #define USE_LL
30 #endif
31
32 #define LOW_I32_FROM_VAR(u64) LOW_I32_FROM_LONG64(u64)
33 #define LOW_I32_FROM_PTR(u64ptr) LOW_I32_FROM_LONG64_PTR(u64ptr)
34 #define HIGH_I32_FROM_VAR(u64) HIGH_I32_FROM_LONG64(u64)
35 #define HIGH_I32_FROM_PTR(u64ptr) HIGH_I32_FROM_LONG64_PTR(u64ptr)
36
37 #define MAX_DOUBLE_ACCURACY_WIDTH 17
38
39 #define DEFAULT_DOUBLE_WIDTH MAX_DOUBLE_ACCURACY_WIDTH
40
41 #if defined(USE_LL)
42 #define DOUBLE_INFINITE_LONGBITS (0x7FF0000000000000LL)
43 #else
44 #if defined(USE_L)
45 #define DOUBLE_INFINITE_LONGBITS (0x7FF0000000000000L)
46 #else
47 #define DOUBLE_INFINITE_LONGBITS (0x7FF0000000000000)
48 #endif /* USE_L */
49 #endif /* USE_LL */
50
51 #define DOUBLE_MINIMUM_LONGBITS (0x1)
52
53 #if defined(USE_LL)
54 #define DOUBLE_MANTISSA_MASK (0x000FFFFFFFFFFFFFLL)
55 #define DOUBLE_EXPONENT_MASK (0x7FF0000000000000LL)
56 #define DOUBLE_NORMAL_MASK (0x0010000000000000LL)
57 #else
58 #if defined(USE_L)
59 #define DOUBLE_MANTISSA_MASK (0x000FFFFFFFFFFFFFL)
60 #define DOUBLE_EXPONENT_MASK (0x7FF0000000000000L)
61 #define DOUBLE_NORMAL_MASK (0x0010000000000000L)
62 #else
63 #define DOUBLE_MANTISSA_MASK (0x000FFFFFFFFFFFFF)
64 #define DOUBLE_EXPONENT_MASK (0x7FF0000000000000)
65 #define DOUBLE_NORMAL_MASK (0x0010000000000000)
66 #endif /* USE_L */
67 #endif /* USE_LL */
68
69 /* Keep a count of the number of times we decrement and increment to
70 * approximate the double, and attempt to detect the case where we
71 * could potentially toggle back and forth between decrementing and
72 * incrementing. It is possible for us to be stuck in the loop when
73 * incrementing by one or decrementing by one may exceed or stay below
74 * the value that we are looking for. In this case, just break out of
75 * the loop if we toggle between incrementing and decrementing for more
76 * than twice.
77 */
78 #define INCREMENT_DOUBLE(_x, _decCount, _incCount) \
79 { \
80 ++DOUBLE_TO_LONGBITS(_x); \
81 _incCount++; \
82 if( (_incCount > 2) && (_decCount > 2) ) { \
83 if( _decCount > _incCount ) { \
84 DOUBLE_TO_LONGBITS(_x) += _decCount - _incCount; \
85 } else if( _incCount > _decCount ) { \
86 DOUBLE_TO_LONGBITS(_x) -= _incCount - _decCount; \
87 } \
88 break; \
89 } \
90 }
91 #define DECREMENT_DOUBLE(_x, _decCount, _incCount) \
92 { \
93 --DOUBLE_TO_LONGBITS(_x); \
94 _decCount++; \
95 if( (_incCount > 2) && (_decCount > 2) ) { \
96 if( _decCount > _incCount ) { \
97 DOUBLE_TO_LONGBITS(_x) += _decCount - _incCount; \
98 } else if( _incCount > _decCount ) { \
99 DOUBLE_TO_LONGBITS(_x) -= _incCount - _decCount; \
100 } \
101 break; \
102 } \
103 }
104
105 #define allocateU64(x, n) if (!((x) = reinterpret_cast<uint64_t*>(malloc((n) * sizeof(uint64_t))))) goto OutOfMemory;
106
107 /* *********************************************************** */
108
109 /* ************************ local data ************************ */
110 static const jdouble double_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 jdouble createDouble1 (JNIEnv* env, uint64_t * f, int32_t length, jint e);
139 static jdouble doubleAlgorithm (JNIEnv* env, uint64_t * f, int32_t length, jint e,
140 jdouble z);
141 /* *********************************************************** */
142
143 #define doubleTenToTheE(e) (*(double_tens + (e)))
144 #define DOUBLE_LOG5_OF_TWO_TO_THE_N 23
145
146 #define sizeOfTenToTheE(e) (((e) / 19) + 1)
147
createDouble(JNIEnv * env,const char * s,jint e)148 static jdouble createDouble(JNIEnv* env, const char* s, jint e) {
149 /* assumes s is a null terminated string with at least one
150 * character in it */
151 uint64_t def[DEFAULT_DOUBLE_WIDTH];
152 uint64_t defBackup[DEFAULT_DOUBLE_WIDTH];
153 uint64_t* f;
154 uint64_t* fNoOverflow;
155 uint64_t* g;
156 uint64_t* tempBackup;
157 uint32_t overflow;
158 jdouble result;
159 int32_t index = 1;
160 int unprocessedDigits = 0;
161
162 f = def;
163 fNoOverflow = defBackup;
164 *f = 0;
165 tempBackup = g = 0;
166 do
167 {
168 if (*s >= '0' && *s <= '9')
169 {
170 /* Make a back up of f before appending, so that we can
171 * back out of it if there is no more room, i.e. index >
172 * MAX_DOUBLE_ACCURACY_WIDTH.
173 */
174 memcpy (fNoOverflow, f, sizeof (uint64_t) * index);
175 overflow =
176 simpleAppendDecimalDigitHighPrecision (f, index, *s - '0');
177 if (overflow)
178 {
179 f[index++] = overflow;
180 /* There is an overflow, but there is no more room
181 * to store the result. We really only need the top 52
182 * bits anyway, so we must back out of the overflow,
183 * and ignore the rest of the string.
184 */
185 if (index >= MAX_DOUBLE_ACCURACY_WIDTH)
186 {
187 index--;
188 memcpy (f, fNoOverflow, sizeof (uint64_t) * index);
189 break;
190 }
191 if (tempBackup)
192 {
193 fNoOverflow = tempBackup;
194 }
195 }
196 }
197 else
198 index = -1;
199 }
200 while (index > 0 && *(++s) != '\0');
201
202 /* We've broken out of the parse loop either because we've reached
203 * the end of the string or we've overflowed the maximum accuracy
204 * limit of a double. If we still have unprocessed digits in the
205 * given string, then there are three possible results:
206 * 1. (unprocessed digits + e) == 0, in which case we simply
207 * convert the existing bits that are already parsed
208 * 2. (unprocessed digits + e) < 0, in which case we simply
209 * convert the existing bits that are already parsed along
210 * with the given e
211 * 3. (unprocessed digits + e) > 0 indicates that the value is
212 * simply too big to be stored as a double, so return Infinity
213 */
214 if ((unprocessedDigits = strlen (s)) > 0)
215 {
216 e += unprocessedDigits;
217 if (index > -1)
218 {
219 if (e == 0)
220 result = toDoubleHighPrecision (f, index);
221 else if (e < 0)
222 result = createDouble1 (env, f, index, e);
223 else
224 {
225 DOUBLE_TO_LONGBITS (result) = DOUBLE_INFINITE_LONGBITS;
226 }
227 }
228 else
229 {
230 LOW_I32_FROM_VAR (result) = -1;
231 HIGH_I32_FROM_VAR (result) = -1;
232 }
233 }
234 else
235 {
236 if (index > -1)
237 {
238 if (e == 0)
239 result = toDoubleHighPrecision (f, index);
240 else
241 result = createDouble1 (env, f, index, e);
242 }
243 else
244 {
245 LOW_I32_FROM_VAR (result) = -1;
246 HIGH_I32_FROM_VAR (result) = -1;
247 }
248 }
249
250 return result;
251
252 }
253
createDouble1(JNIEnv * env,uint64_t * f,int32_t length,jint e)254 static jdouble createDouble1(JNIEnv* env, uint64_t* f, int32_t length, jint e) {
255 int32_t numBits;
256 jdouble result;
257
258 static const jint APPROX_MIN_MAGNITUDE = -309;
259 static const jint APPROX_MAX_MAGNITUDE = 309;
260
261 numBits = highestSetBitHighPrecision (f, length) + 1;
262 numBits -= lowestSetBitHighPrecision (f, length);
263 if (numBits < 54 && e >= 0 && e < DOUBLE_LOG5_OF_TWO_TO_THE_N)
264 {
265 return toDoubleHighPrecision (f, length) * doubleTenToTheE (e);
266 }
267 else if (numBits < 54 && e < 0 && (-e) < DOUBLE_LOG5_OF_TWO_TO_THE_N)
268 {
269 return toDoubleHighPrecision (f, length) / doubleTenToTheE (-e);
270 }
271 else if (e >= 0 && e < APPROX_MAX_MAGNITUDE)
272 {
273 result = toDoubleHighPrecision (f, length) * pow (10.0, e);
274 }
275 else if (e >= APPROX_MAX_MAGNITUDE)
276 {
277 /* Convert the partial result to make sure that the
278 * non-exponential part is not zero. This check fixes the case
279 * where the user enters 0.0e309! */
280 result = toDoubleHighPrecision (f, length);
281 /* Don't go straight to zero as the fact that x*0 = 0 independent of x might
282 cause the algorithm to produce an incorrect result. Instead try the min value
283 first and let it fall to zero if need be. */
284
285 if (result == 0.0)
286 {
287 DOUBLE_TO_LONGBITS (result) = DOUBLE_MINIMUM_LONGBITS;
288 }
289 else
290 {
291 DOUBLE_TO_LONGBITS (result) = DOUBLE_INFINITE_LONGBITS;
292 return result;
293 }
294 }
295 else if (e > APPROX_MIN_MAGNITUDE)
296 {
297 result = toDoubleHighPrecision (f, length) / pow (10.0, -e);
298 }
299
300 if (e <= APPROX_MIN_MAGNITUDE)
301 {
302
303 result = toDoubleHighPrecision (f, length) * pow (10.0, e + 52);
304 result = result * pow (10.0, -52);
305
306 }
307
308 /* Don't go straight to zero as the fact that x*0 = 0 independent of x might
309 cause the algorithm to produce an incorrect result. Instead try the min value
310 first and let it fall to zero if need be. */
311
312 if (result == 0.0)
313
314 DOUBLE_TO_LONGBITS (result) = DOUBLE_MINIMUM_LONGBITS;
315
316 return doubleAlgorithm (env, f, length, e, result);
317 }
318
319 /* The algorithm for the function doubleAlgorithm() below can be found
320 * in:
321 *
322 * "How to Read Floating-Point Numbers Accurately", William D.
323 * Clinger, Proceedings of the ACM SIGPLAN '90 Conference on
324 * Programming Language Design and Implementation, June 20-22,
325 * 1990, pp. 92-101.
326 *
327 * There is a possibility that the function will end up in an endless
328 * loop if the given approximating floating-point number (a very small
329 * floating-point whose value is very close to zero) straddles between
330 * two approximating integer values. We modified the algorithm slightly
331 * to detect the case where it oscillates back and forth between
332 * incrementing and decrementing the floating-point approximation. It
333 * is currently set such that if the oscillation occurs more than twice
334 * then return the original approximation.
335 */
doubleAlgorithm(JNIEnv * env,uint64_t * f,int32_t length,jint e,jdouble z)336 static jdouble doubleAlgorithm(JNIEnv* env, uint64_t* f, int32_t length, jint e, jdouble z) {
337 uint64_t m;
338 int32_t k, comparison, comparison2;
339 uint64_t* x;
340 uint64_t* y;
341 uint64_t* D;
342 uint64_t* D2;
343 int32_t xLength, yLength, DLength, D2Length, decApproxCount, incApproxCount;
344
345 x = y = D = D2 = 0;
346 xLength = yLength = DLength = D2Length = 0;
347 decApproxCount = incApproxCount = 0;
348
349 do
350 {
351 m = doubleMantissa (z);
352 k = doubleExponent (z);
353
354 if (x && x != f)
355 free(x);
356
357 free(y);
358 free(D);
359 free(D2);
360
361 if (e >= 0 && k >= 0)
362 {
363 xLength = sizeOfTenToTheE (e) + length;
364 allocateU64 (x, xLength);
365 memset (x + length, 0, sizeof (uint64_t) * (xLength - length));
366 memcpy (x, f, sizeof (uint64_t) * length);
367 timesTenToTheEHighPrecision (x, xLength, e);
368
369 yLength = (k >> 6) + 2;
370 allocateU64 (y, yLength);
371 memset (y + 1, 0, sizeof (uint64_t) * (yLength - 1));
372 *y = m;
373 simpleShiftLeftHighPrecision (y, yLength, k);
374 }
375 else if (e >= 0)
376 {
377 xLength = sizeOfTenToTheE (e) + length + ((-k) >> 6) + 1;
378 allocateU64 (x, xLength);
379 memset (x + length, 0, sizeof (uint64_t) * (xLength - length));
380 memcpy (x, f, sizeof (uint64_t) * length);
381 timesTenToTheEHighPrecision (x, xLength, e);
382 simpleShiftLeftHighPrecision (x, xLength, -k);
383
384 yLength = 1;
385 allocateU64 (y, 1);
386 *y = m;
387 }
388 else if (k >= 0)
389 {
390 xLength = length;
391 x = f;
392
393 yLength = sizeOfTenToTheE (-e) + 2 + (k >> 6);
394 allocateU64 (y, yLength);
395 memset (y + 1, 0, sizeof (uint64_t) * (yLength - 1));
396 *y = m;
397 timesTenToTheEHighPrecision (y, yLength, -e);
398 simpleShiftLeftHighPrecision (y, yLength, k);
399 }
400 else
401 {
402 xLength = length + ((-k) >> 6) + 1;
403 allocateU64 (x, xLength);
404 memset (x + length, 0, sizeof (uint64_t) * (xLength - length));
405 memcpy (x, f, sizeof (uint64_t) * length);
406 simpleShiftLeftHighPrecision (x, xLength, -k);
407
408 yLength = sizeOfTenToTheE (-e) + 1;
409 allocateU64 (y, yLength);
410 memset (y + 1, 0, sizeof (uint64_t) * (yLength - 1));
411 *y = m;
412 timesTenToTheEHighPrecision (y, yLength, -e);
413 }
414
415 comparison = compareHighPrecision (x, xLength, y, yLength);
416 if (comparison > 0)
417 { /* x > y */
418 DLength = xLength;
419 allocateU64 (D, DLength);
420 memcpy (D, x, DLength * sizeof (uint64_t));
421 subtractHighPrecision (D, DLength, y, yLength);
422 }
423 else if (comparison)
424 { /* y > x */
425 DLength = yLength;
426 allocateU64 (D, DLength);
427 memcpy (D, y, DLength * sizeof (uint64_t));
428 subtractHighPrecision (D, DLength, x, xLength);
429 }
430 else
431 { /* y == x */
432 DLength = 1;
433 allocateU64 (D, 1);
434 *D = 0;
435 }
436
437 D2Length = DLength + 1;
438 allocateU64 (D2, D2Length);
439 m <<= 1;
440 multiplyHighPrecision (D, DLength, &m, 1, D2, D2Length);
441 m >>= 1;
442
443 comparison2 = compareHighPrecision (D2, D2Length, y, yLength);
444 if (comparison2 < 0)
445 {
446 if (comparison < 0 && m == DOUBLE_NORMAL_MASK)
447 {
448 simpleShiftLeftHighPrecision (D2, D2Length, 1);
449 if (compareHighPrecision (D2, D2Length, y, yLength) > 0)
450 {
451 DECREMENT_DOUBLE (z, decApproxCount, incApproxCount);
452 }
453 else
454 {
455 break;
456 }
457 }
458 else
459 {
460 break;
461 }
462 }
463 else if (comparison2 == 0)
464 {
465 if ((LOW_U32_FROM_VAR (m) & 1) == 0)
466 {
467 if (comparison < 0 && m == DOUBLE_NORMAL_MASK)
468 {
469 DECREMENT_DOUBLE (z, decApproxCount, incApproxCount);
470 }
471 else
472 {
473 break;
474 }
475 }
476 else if (comparison < 0)
477 {
478 DECREMENT_DOUBLE (z, decApproxCount, incApproxCount);
479 break;
480 }
481 else
482 {
483 INCREMENT_DOUBLE (z, decApproxCount, incApproxCount);
484 break;
485 }
486 }
487 else if (comparison < 0)
488 {
489 DECREMENT_DOUBLE (z, decApproxCount, incApproxCount);
490 }
491 else
492 {
493 if (DOUBLE_TO_LONGBITS (z) == DOUBLE_INFINITE_LONGBITS)
494 break;
495 INCREMENT_DOUBLE (z, decApproxCount, incApproxCount);
496 }
497 }
498 while (1);
499
500 if (x && x != f)
501 free(x);
502 free(y);
503 free(D);
504 free(D2);
505 return z;
506
507 OutOfMemory:
508 if (x && x != f)
509 free(x);
510 free(y);
511 free(D);
512 free(D2);
513 jniThrowOutOfMemoryError(env, NULL);
514 return z;
515 }
516
517
518
519 #define MAX_FLOAT_ACCURACY_WIDTH 8
520
521 #define DEFAULT_FLOAT_WIDTH MAX_FLOAT_ACCURACY_WIDTH
522
523 static jfloat createFloat1(JNIEnv* env, uint64_t* f, int32_t length, jint e);
524 static jfloat floatAlgorithm(JNIEnv* env, uint64_t* f, int32_t length, jint e, jfloat z);
525
526 static const uint32_t float_tens[] = {
527 0x3f800000,
528 0x41200000,
529 0x42c80000,
530 0x447a0000,
531 0x461c4000,
532 0x47c35000,
533 0x49742400,
534 0x4b189680,
535 0x4cbebc20,
536 0x4e6e6b28,
537 0x501502f9 /* 10 ^ 10 in float */
538 };
539
540 #define floatTenToTheE(e) (*reinterpret_cast<const jfloat*>(float_tens + (e)))
541 #define FLOAT_LOG5_OF_TWO_TO_THE_N 11
542
543 #define FLOAT_INFINITE_INTBITS (0x7F800000)
544 #define FLOAT_MINIMUM_INTBITS (1)
545
546 #define FLOAT_MANTISSA_MASK (0x007FFFFF)
547 #define FLOAT_EXPONENT_MASK (0x7F800000)
548 #define FLOAT_NORMAL_MASK (0x00800000)
549
550 /* Keep a count of the number of times we decrement and increment to
551 * approximate the double, and attempt to detect the case where we
552 * could potentially toggle back and forth between decrementing and
553 * incrementing. It is possible for us to be stuck in the loop when
554 * incrementing by one or decrementing by one may exceed or stay below
555 * the value that we are looking for. In this case, just break out of
556 * the loop if we toggle between incrementing and decrementing for more
557 * than twice.
558 */
559 #define INCREMENT_FLOAT(_x, _decCount, _incCount) \
560 { \
561 ++FLOAT_TO_INTBITS(_x); \
562 _incCount++; \
563 if( (_incCount > 2) && (_decCount > 2) ) { \
564 if( _decCount > _incCount ) { \
565 FLOAT_TO_INTBITS(_x) += _decCount - _incCount; \
566 } else if( _incCount > _decCount ) { \
567 FLOAT_TO_INTBITS(_x) -= _incCount - _decCount; \
568 } \
569 break; \
570 } \
571 }
572 #define DECREMENT_FLOAT(_x, _decCount, _incCount) \
573 { \
574 --FLOAT_TO_INTBITS(_x); \
575 _decCount++; \
576 if( (_incCount > 2) && (_decCount > 2) ) { \
577 if( _decCount > _incCount ) { \
578 FLOAT_TO_INTBITS(_x) += _decCount - _incCount; \
579 } else if( _incCount > _decCount ) { \
580 FLOAT_TO_INTBITS(_x) -= _incCount - _decCount; \
581 } \
582 break; \
583 } \
584 }
585
createFloat(JNIEnv * env,const char * s,jint e)586 static jfloat createFloat(JNIEnv* env, const char* s, jint e) {
587 /* assumes s is a null terminated string with at least one
588 * character in it */
589 uint64_t def[DEFAULT_FLOAT_WIDTH];
590 uint64_t defBackup[DEFAULT_FLOAT_WIDTH];
591 uint64_t* f;
592 uint64_t* fNoOverflow;
593 uint64_t* g;
594 uint64_t* tempBackup;
595 uint32_t overflow;
596 jfloat result;
597 int32_t index = 1;
598 int unprocessedDigits = 0;
599
600 f = def;
601 fNoOverflow = defBackup;
602 *f = 0;
603 tempBackup = g = 0;
604 do
605 {
606 if (*s >= '0' && *s <= '9')
607 {
608 /* Make a back up of f before appending, so that we can
609 * back out of it if there is no more room, i.e. index >
610 * MAX_FLOAT_ACCURACY_WIDTH.
611 */
612 memcpy (fNoOverflow, f, sizeof (uint64_t) * index);
613 overflow =
614 simpleAppendDecimalDigitHighPrecision (f, index, *s - '0');
615 if (overflow)
616 {
617
618 f[index++] = overflow;
619 /* There is an overflow, but there is no more room
620 * to store the result. We really only need the top 52
621 * bits anyway, so we must back out of the overflow,
622 * and ignore the rest of the string.
623 */
624 if (index >= MAX_FLOAT_ACCURACY_WIDTH)
625 {
626 index--;
627 memcpy (f, fNoOverflow, sizeof (uint64_t) * index);
628 break;
629 }
630 if (tempBackup)
631 {
632 fNoOverflow = tempBackup;
633 }
634 }
635 }
636 else
637 index = -1;
638 }
639 while (index > 0 && *(++s) != '\0');
640
641 /* We've broken out of the parse loop either because we've reached
642 * the end of the string or we've overflowed the maximum accuracy
643 * limit of a double. If we still have unprocessed digits in the
644 * given string, then there are three possible results:
645 * 1. (unprocessed digits + e) == 0, in which case we simply
646 * convert the existing bits that are already parsed
647 * 2. (unprocessed digits + e) < 0, in which case we simply
648 * convert the existing bits that are already parsed along
649 * with the given e
650 * 3. (unprocessed digits + e) > 0 indicates that the value is
651 * simply too big to be stored as a double, so return Infinity
652 */
653 if ((unprocessedDigits = strlen (s)) > 0)
654 {
655 e += unprocessedDigits;
656 if (index > -1)
657 {
658 if (e <= 0)
659 {
660 result = createFloat1 (env, f, index, e);
661 }
662 else
663 {
664 FLOAT_TO_INTBITS (result) = FLOAT_INFINITE_INTBITS;
665 }
666 }
667 else
668 {
669 result = INTBITS_TO_FLOAT(index);
670 }
671 }
672 else
673 {
674 if (index > -1)
675 {
676 result = createFloat1 (env, f, index, e);
677 }
678 else
679 {
680 result = INTBITS_TO_FLOAT(index);
681 }
682 }
683
684 return result;
685
686 }
687
createFloat1(JNIEnv * env,uint64_t * f,int32_t length,jint e)688 static jfloat createFloat1 (JNIEnv* env, uint64_t* f, int32_t length, jint e) {
689 int32_t numBits;
690 jdouble dresult;
691 jfloat result;
692
693 numBits = highestSetBitHighPrecision (f, length) + 1;
694 if (numBits < 25 && e >= 0 && e < FLOAT_LOG5_OF_TWO_TO_THE_N)
695 {
696 return ((jfloat) LOW_I32_FROM_PTR (f)) * floatTenToTheE (e);
697 }
698 else if (numBits < 25 && e < 0 && (-e) < FLOAT_LOG5_OF_TWO_TO_THE_N)
699 {
700 return ((jfloat) LOW_I32_FROM_PTR (f)) / floatTenToTheE (-e);
701 }
702 else if (e >= 0 && e < 39)
703 {
704 result = (jfloat) (toDoubleHighPrecision (f, length) * pow (10.0, (double) e));
705 }
706 else if (e >= 39)
707 {
708 /* Convert the partial result to make sure that the
709 * non-exponential part is not zero. This check fixes the case
710 * where the user enters 0.0e309! */
711 result = (jfloat) toDoubleHighPrecision (f, length);
712
713 if (result == 0.0)
714
715 FLOAT_TO_INTBITS (result) = FLOAT_MINIMUM_INTBITS;
716 else
717 FLOAT_TO_INTBITS (result) = FLOAT_INFINITE_INTBITS;
718 }
719 else if (e > -309)
720 {
721 int dexp;
722 uint32_t fmant, fovfl;
723 uint64_t dmant;
724 dresult = toDoubleHighPrecision (f, length) / pow (10.0, (double) -e);
725 if (IS_DENORMAL_DBL (dresult))
726 {
727 FLOAT_TO_INTBITS (result) = 0;
728 return result;
729 }
730 dexp = doubleExponent (dresult) + 51;
731 dmant = doubleMantissa (dresult);
732 /* Is it too small to be represented by a single-precision
733 * float? */
734 if (dexp <= -155)
735 {
736 FLOAT_TO_INTBITS (result) = 0;
737 return result;
738 }
739 /* Is it a denormalized single-precision float? */
740 if ((dexp <= -127) && (dexp > -155))
741 {
742 /* Only interested in 24 msb bits of the 53-bit double mantissa */
743 fmant = (uint32_t) (dmant >> 29);
744 fovfl = ((uint32_t) (dmant & 0x1FFFFFFF)) << 3;
745 while ((dexp < -127) && ((fmant | fovfl) != 0))
746 {
747 if ((fmant & 1) != 0)
748 {
749 fovfl |= 0x80000000;
750 }
751 fovfl >>= 1;
752 fmant >>= 1;
753 dexp++;
754 }
755 if ((fovfl & 0x80000000) != 0)
756 {
757 if ((fovfl & 0x7FFFFFFC) != 0)
758 {
759 fmant++;
760 }
761 else if ((fmant & 1) != 0)
762 {
763 fmant++;
764 }
765 }
766 else if ((fovfl & 0x40000000) != 0)
767 {
768 if ((fovfl & 0x3FFFFFFC) != 0)
769 {
770 fmant++;
771 }
772 }
773 FLOAT_TO_INTBITS (result) = fmant;
774 }
775 else
776 {
777 result = (jfloat) dresult;
778 }
779 }
780
781 /* Don't go straight to zero as the fact that x*0 = 0 independent
782 * of x might cause the algorithm to produce an incorrect result.
783 * Instead try the min value first and let it fall to zero if need
784 * be.
785 */
786 if (e <= -309 || FLOAT_TO_INTBITS (result) == 0)
787 FLOAT_TO_INTBITS (result) = FLOAT_MINIMUM_INTBITS;
788
789 return floatAlgorithm (env, f, length, e, (jfloat) result);
790 }
791
792 /* The algorithm for the function floatAlgorithm() below can be found
793 * in:
794 *
795 * "How to Read Floating-Point Numbers Accurately", William D.
796 * Clinger, Proceedings of the ACM SIGPLAN '90 Conference on
797 * Programming Language Design and Implementation, June 20-22,
798 * 1990, pp. 92-101.
799 *
800 * There is a possibility that the function will end up in an endless
801 * loop if the given approximating floating-point number (a very small
802 * floating-point whose value is very close to zero) straddles between
803 * two approximating integer values. We modified the algorithm slightly
804 * to detect the case where it oscillates back and forth between
805 * incrementing and decrementing the floating-point approximation. It
806 * is currently set such that if the oscillation occurs more than twice
807 * then return the original approximation.
808 */
floatAlgorithm(JNIEnv * env,uint64_t * f,int32_t length,jint e,jfloat z)809 static jfloat floatAlgorithm(JNIEnv* env, uint64_t* f, int32_t length, jint e, jfloat z) {
810 uint64_t m;
811 int32_t k, comparison, comparison2;
812 uint64_t* x;
813 uint64_t* y;
814 uint64_t* D;
815 uint64_t* D2;
816 int32_t xLength, yLength, DLength, D2Length;
817 int32_t decApproxCount, incApproxCount;
818
819 x = y = D = D2 = 0;
820 xLength = yLength = DLength = D2Length = 0;
821 decApproxCount = incApproxCount = 0;
822
823 do
824 {
825 m = floatMantissa (z);
826 k = floatExponent (z);
827
828 if (x && x != f)
829 free(x);
830
831 free(y);
832 free(D);
833 free(D2);
834
835 if (e >= 0 && k >= 0)
836 {
837 xLength = sizeOfTenToTheE (e) + length;
838 allocateU64 (x, xLength);
839 memset (x + length, 0, sizeof (uint64_t) * (xLength - length));
840 memcpy (x, f, sizeof (uint64_t) * length);
841 timesTenToTheEHighPrecision (x, xLength, e);
842
843 yLength = (k >> 6) + 2;
844 allocateU64 (y, yLength);
845 memset (y + 1, 0, sizeof (uint64_t) * (yLength - 1));
846 *y = m;
847 simpleShiftLeftHighPrecision (y, yLength, k);
848 }
849 else if (e >= 0)
850 {
851 xLength = sizeOfTenToTheE (e) + length + ((-k) >> 6) + 1;
852 allocateU64 (x, xLength);
853 memset (x + length, 0, sizeof (uint64_t) * (xLength - length));
854 memcpy (x, f, sizeof (uint64_t) * length);
855 timesTenToTheEHighPrecision (x, xLength, e);
856 simpleShiftLeftHighPrecision (x, xLength, -k);
857
858 yLength = 1;
859 allocateU64 (y, 1);
860 *y = m;
861 }
862 else if (k >= 0)
863 {
864 xLength = length;
865 x = f;
866
867 yLength = sizeOfTenToTheE (-e) + 2 + (k >> 6);
868 allocateU64 (y, yLength);
869 memset (y + 1, 0, sizeof (uint64_t) * (yLength - 1));
870 *y = m;
871 timesTenToTheEHighPrecision (y, yLength, -e);
872 simpleShiftLeftHighPrecision (y, yLength, k);
873 }
874 else
875 {
876 xLength = length + ((-k) >> 6) + 1;
877 allocateU64 (x, xLength);
878 memset (x + length, 0, sizeof (uint64_t) * (xLength - length));
879 memcpy (x, f, sizeof (uint64_t) * length);
880 simpleShiftLeftHighPrecision (x, xLength, -k);
881
882 yLength = sizeOfTenToTheE (-e) + 1;
883 allocateU64 (y, yLength);
884 memset (y + 1, 0, sizeof (uint64_t) * (yLength - 1));
885 *y = m;
886 timesTenToTheEHighPrecision (y, yLength, -e);
887 }
888
889 comparison = compareHighPrecision (x, xLength, y, yLength);
890 if (comparison > 0)
891 { /* x > y */
892 DLength = xLength;
893 allocateU64 (D, DLength);
894 memcpy (D, x, DLength * sizeof (uint64_t));
895 subtractHighPrecision (D, DLength, y, yLength);
896 }
897 else if (comparison)
898 { /* y > x */
899 DLength = yLength;
900 allocateU64 (D, DLength);
901 memcpy (D, y, DLength * sizeof (uint64_t));
902 subtractHighPrecision (D, DLength, x, xLength);
903 }
904 else
905 { /* y == x */
906 DLength = 1;
907 allocateU64 (D, 1);
908 *D = 0;
909 }
910
911 D2Length = DLength + 1;
912 allocateU64 (D2, D2Length);
913 m <<= 1;
914 multiplyHighPrecision (D, DLength, &m, 1, D2, D2Length);
915 m >>= 1;
916
917 comparison2 = compareHighPrecision (D2, D2Length, y, yLength);
918 if (comparison2 < 0)
919 {
920 if (comparison < 0 && m == FLOAT_NORMAL_MASK)
921 {
922 simpleShiftLeftHighPrecision (D2, D2Length, 1);
923 if (compareHighPrecision (D2, D2Length, y, yLength) > 0)
924 {
925 DECREMENT_FLOAT (z, decApproxCount, incApproxCount);
926 }
927 else
928 {
929 break;
930 }
931 }
932 else
933 {
934 break;
935 }
936 }
937 else if (comparison2 == 0)
938 {
939 if ((m & 1) == 0)
940 {
941 if (comparison < 0 && m == FLOAT_NORMAL_MASK)
942 {
943 DECREMENT_FLOAT (z, decApproxCount, incApproxCount);
944 }
945 else
946 {
947 break;
948 }
949 }
950 else if (comparison < 0)
951 {
952 DECREMENT_FLOAT (z, decApproxCount, incApproxCount);
953 break;
954 }
955 else
956 {
957 INCREMENT_FLOAT (z, decApproxCount, incApproxCount);
958 break;
959 }
960 }
961 else if (comparison < 0)
962 {
963 DECREMENT_FLOAT (z, decApproxCount, incApproxCount);
964 }
965 else
966 {
967 if (FLOAT_TO_INTBITS (z) == FLOAT_EXPONENT_MASK)
968 break;
969 INCREMENT_FLOAT (z, decApproxCount, incApproxCount);
970 }
971 }
972 while (1);
973
974 if (x && x != f)
975 free(x);
976 free(y);
977 free(D);
978 free(D2);
979 return z;
980
981 OutOfMemory:
982 if (x && x != f)
983 free(x);
984 free(y);
985 free(D);
986 free(D2);
987 jniThrowOutOfMemoryError(env, NULL);
988 return z;
989 }
990
StringToReal_parseFltImpl(JNIEnv * env,jclass,jstring s,jint e)991 static jfloat StringToReal_parseFltImpl(JNIEnv* env, jclass, jstring s, jint e) {
992 ScopedUtfChars str(env, s);
993 if (str.c_str() == NULL) {
994 return 0.0;
995 }
996 return createFloat(env, str.c_str(), e);
997 }
998
StringToReal_parseDblImpl(JNIEnv * env,jclass,jstring s,jint e)999 static jdouble StringToReal_parseDblImpl(JNIEnv* env, jclass, jstring s, jint e) {
1000 ScopedUtfChars str(env, s);
1001 if (str.c_str() == NULL) {
1002 return 0.0;
1003 }
1004 return createDouble(env, str.c_str(), e);
1005 }
1006
1007 static JNINativeMethod gMethods[] = {
1008 NATIVE_METHOD(StringToReal, parseFltImpl, "(Ljava/lang/String;I)F"),
1009 NATIVE_METHOD(StringToReal, parseDblImpl, "(Ljava/lang/String;I)D"),
1010 };
register_java_lang_StringToReal(JNIEnv * env)1011 void register_java_lang_StringToReal(JNIEnv* env) {
1012 jniRegisterNativeMethods(env, "java/lang/StringToReal", gMethods, NELEM(gMethods));
1013 }
1014