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