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 <nativehelper/JNIHelp.h>
22 #include <nativehelper/JniConstants.h>
23 #include <nativehelper/ScopedUtfChars.h>
24 #include "cbigint.h"
25 #include "JniException.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 y = D = D2 = NULL;
290
291 if (e >= 0 && k >= 0)
292 {
293 xLength = sizeOfTenToTheE (e) + length;
294 allocateU64 (x, xLength);
295 memset (x + length, 0, sizeof (uint64_t) * (xLength - length));
296 memcpy (x, f, sizeof (uint64_t) * length);
297 timesTenToTheEHighPrecision (x, xLength, e);
298
299 yLength = (k >> 6) + 2;
300 allocateU64 (y, yLength);
301 memset (y + 1, 0, sizeof (uint64_t) * (yLength - 1));
302 *y = m;
303 simpleShiftLeftHighPrecision (y, yLength, k);
304 }
305 else if (e >= 0)
306 {
307 xLength = sizeOfTenToTheE (e) + length + ((-k) >> 6) + 1;
308 allocateU64 (x, xLength);
309 memset (x + length, 0, sizeof (uint64_t) * (xLength - length));
310 memcpy (x, f, sizeof (uint64_t) * length);
311 timesTenToTheEHighPrecision (x, xLength, e);
312 simpleShiftLeftHighPrecision (x, xLength, -k);
313
314 yLength = 1;
315 allocateU64 (y, 1);
316 *y = m;
317 }
318 else if (k >= 0)
319 {
320 xLength = length;
321 x = f;
322
323 yLength = sizeOfTenToTheE (-e) + 2 + (k >> 6);
324 allocateU64 (y, yLength);
325 memset (y + 1, 0, sizeof (uint64_t) * (yLength - 1));
326 *y = m;
327 timesTenToTheEHighPrecision (y, yLength, -e);
328 simpleShiftLeftHighPrecision (y, yLength, k);
329 }
330 else
331 {
332 xLength = length + ((-k) >> 6) + 1;
333 allocateU64 (x, xLength);
334 memset (x + length, 0, sizeof (uint64_t) * (xLength - length));
335 memcpy (x, f, sizeof (uint64_t) * length);
336 simpleShiftLeftHighPrecision (x, xLength, -k);
337
338 yLength = sizeOfTenToTheE (-e) + 1;
339 allocateU64 (y, yLength);
340 memset (y + 1, 0, sizeof (uint64_t) * (yLength - 1));
341 *y = m;
342 timesTenToTheEHighPrecision (y, yLength, -e);
343 }
344
345 comparison = compareHighPrecision (x, xLength, y, yLength);
346 if (comparison > 0)
347 { /* x > y */
348 DLength = xLength;
349 allocateU64 (D, DLength);
350 memcpy (D, x, DLength * sizeof (uint64_t));
351 subtractHighPrecision (D, DLength, y, yLength);
352 }
353 else if (comparison)
354 { /* y > x */
355 DLength = yLength;
356 allocateU64 (D, DLength);
357 memcpy (D, y, DLength * sizeof (uint64_t));
358 subtractHighPrecision (D, DLength, x, xLength);
359 }
360 else
361 { /* y == x */
362 DLength = 1;
363 allocateU64 (D, 1);
364 *D = 0;
365 }
366
367 D2Length = DLength + 1;
368 allocateU64 (D2, D2Length);
369 m <<= 1;
370 multiplyHighPrecision (D, DLength, &m, 1, D2, D2Length);
371 m >>= 1;
372
373 comparison2 = compareHighPrecision (D2, D2Length, y, yLength);
374 if (comparison2 < 0)
375 {
376 if (comparison < 0 && m == DOUBLE_NORMAL_MASK
377 && DOUBLE_TO_LONGBITS(z) != DOUBLE_NORMAL_MASK)
378 {
379 simpleShiftLeftHighPrecision (D2, D2Length, 1);
380 if (compareHighPrecision (D2, D2Length, y, yLength) > 0)
381 {
382 --DOUBLE_TO_LONGBITS (z);
383 }
384 else
385 {
386 break;
387 }
388 }
389 else
390 {
391 break;
392 }
393 }
394 else if (comparison2 == 0)
395 {
396 if ((LOW_U32_FROM_VAR (m) & 1) == 0)
397 {
398 if (comparison < 0 && m == DOUBLE_NORMAL_MASK)
399 {
400 --DOUBLE_TO_LONGBITS (z);
401 }
402 else
403 {
404 break;
405 }
406 }
407 else if (comparison < 0)
408 {
409 --DOUBLE_TO_LONGBITS (z);
410 break;
411 }
412 else
413 {
414 ++DOUBLE_TO_LONGBITS (z);
415 break;
416 }
417 }
418 else if (comparison < 0)
419 {
420 --DOUBLE_TO_LONGBITS (z);
421 }
422 else
423 {
424 if (DOUBLE_TO_LONGBITS (z) == DOUBLE_INFINITE_LONGBITS)
425 break;
426 ++DOUBLE_TO_LONGBITS (z);
427 }
428 }
429 while (1);
430
431 if (x && x != f)
432 free(x);
433 free(y);
434 free(D);
435 free(D2);
436 return z;
437
438 OutOfMemory:
439 if (x && x != f)
440 free(x);
441 free(y);
442 free(D);
443 free(D2);
444 jniThrowOutOfMemoryError(env, NULL);
445 return z;
446 }
447
448
449
450 #define MAX_FLOAT_ACCURACY_WIDTH 8
451
452 #define DEFAULT_FLOAT_WIDTH MAX_FLOAT_ACCURACY_WIDTH
453
454 static jfloat createFloat1(JNIEnv* env, uint64_t* f, int32_t length, jint e);
455 static jfloat floatAlgorithm(JNIEnv* env, uint64_t* f, int32_t length, jint e, jfloat z);
456
457 static const uint32_t float_tens[] = {
458 0x3f800000,
459 0x41200000,
460 0x42c80000,
461 0x447a0000,
462 0x461c4000,
463 0x47c35000,
464 0x49742400,
465 0x4b189680,
466 0x4cbebc20,
467 0x4e6e6b28,
468 0x501502f9 /* 10 ^ 10 in float */
469 };
470
471 #define floatTenToTheE(e) (*reinterpret_cast<const jfloat*>(float_tens + (e)))
472 #define FLOAT_LOG5_OF_TWO_TO_THE_N 11
473
474 #define FLOAT_INFINITE_INTBITS (0x7F800000)
475 #define FLOAT_MINIMUM_INTBITS (1)
476
477 #define FLOAT_MANTISSA_MASK (0x007FFFFF)
478 #define FLOAT_EXPONENT_MASK (0x7F800000)
479 #define FLOAT_NORMAL_MASK (0x00800000)
480
createFloat(JNIEnv * env,const char * s,jint e)481 static jfloat createFloat(JNIEnv* env, const char* s, jint e) {
482 /* assumes s is a null terminated string with at least one
483 * character in it */
484 uint64_t def[DEFAULT_FLOAT_WIDTH];
485 uint64_t defBackup[DEFAULT_FLOAT_WIDTH];
486 uint64_t* f;
487 uint64_t* fNoOverflow;
488 uint64_t* g;
489 uint64_t* tempBackup;
490 uint32_t overflow;
491 jfloat result;
492 int32_t index = 1;
493 int unprocessedDigits = 0;
494
495 f = def;
496 fNoOverflow = defBackup;
497 *f = 0;
498 tempBackup = g = 0;
499 do
500 {
501 if (*s >= '0' && *s <= '9')
502 {
503 /* Make a back up of f before appending, so that we can
504 * back out of it if there is no more room, i.e. index >
505 * MAX_FLOAT_ACCURACY_WIDTH.
506 */
507 memcpy (fNoOverflow, f, sizeof (uint64_t) * index);
508 overflow =
509 simpleAppendDecimalDigitHighPrecision (f, index, *s - '0');
510 if (overflow)
511 {
512
513 f[index++] = overflow;
514 /* There is an overflow, but there is no more room
515 * to store the result. We really only need the top 52
516 * bits anyway, so we must back out of the overflow,
517 * and ignore the rest of the string.
518 */
519 if (index >= MAX_FLOAT_ACCURACY_WIDTH)
520 {
521 index--;
522 memcpy (f, fNoOverflow, sizeof (uint64_t) * index);
523 break;
524 }
525 if (tempBackup)
526 {
527 fNoOverflow = tempBackup;
528 }
529 }
530 }
531 else
532 index = -1;
533 }
534 while (index > 0 && *(++s) != '\0');
535
536 /* We've broken out of the parse loop either because we've reached
537 * the end of the string or we've overflowed the maximum accuracy
538 * limit of a double. If we still have unprocessed digits in the
539 * given string, then there are three possible results:
540 * 1. (unprocessed digits + e) == 0, in which case we simply
541 * convert the existing bits that are already parsed
542 * 2. (unprocessed digits + e) < 0, in which case we simply
543 * convert the existing bits that are already parsed along
544 * with the given e
545 * 3. (unprocessed digits + e) > 0 indicates that the value is
546 * simply too big to be stored as a double, so return Infinity
547 */
548 if ((unprocessedDigits = strlen (s)) > 0)
549 {
550 e += unprocessedDigits;
551 if (index > -1)
552 {
553 if (e <= 0)
554 {
555 result = createFloat1 (env, f, index, e);
556 }
557 else
558 {
559 FLOAT_TO_INTBITS (result) = FLOAT_INFINITE_INTBITS;
560 }
561 }
562 else
563 {
564 result = INTBITS_TO_FLOAT(index);
565 }
566 }
567 else
568 {
569 if (index > -1)
570 {
571 result = createFloat1 (env, f, index, e);
572 }
573 else
574 {
575 result = INTBITS_TO_FLOAT(index);
576 }
577 }
578
579 return result;
580 }
581
createFloat1(JNIEnv * env,uint64_t * f,int32_t length,jint e)582 static jfloat createFloat1 (JNIEnv* env, uint64_t* f, int32_t length, jint e) {
583 int32_t numBits;
584 jdouble dresult;
585 jfloat result;
586
587 numBits = highestSetBitHighPrecision (f, length) + 1;
588 if (numBits < 25 && e >= 0 && e < FLOAT_LOG5_OF_TWO_TO_THE_N)
589 {
590 return ((jfloat) LOW_I32_FROM_PTR (f)) * floatTenToTheE (e);
591 }
592 else if (numBits < 25 && e < 0 && (-e) < FLOAT_LOG5_OF_TWO_TO_THE_N)
593 {
594 return ((jfloat) LOW_I32_FROM_PTR (f)) / floatTenToTheE (-e);
595 }
596 else if (e >= 0 && e < 39)
597 {
598 result = (jfloat) (toDoubleHighPrecision (f, length) * pow (10.0, (double) e));
599 }
600 else if (e >= 39)
601 {
602 /* Convert the partial result to make sure that the
603 * non-exponential part is not zero. This check fixes the case
604 * where the user enters 0.0e309! */
605 result = (jfloat) toDoubleHighPrecision (f, length);
606
607 if (result == 0.0)
608
609 FLOAT_TO_INTBITS (result) = FLOAT_MINIMUM_INTBITS;
610 else
611 FLOAT_TO_INTBITS (result) = FLOAT_INFINITE_INTBITS;
612 }
613 else if (e > -309)
614 {
615 int dexp;
616 uint32_t fmant, fovfl;
617 uint64_t dmant;
618 dresult = toDoubleHighPrecision (f, length) / pow (10.0, (double) -e);
619 if (IS_DENORMAL_DBL (dresult))
620 {
621 FLOAT_TO_INTBITS (result) = 0;
622 return result;
623 }
624 dexp = doubleExponent (dresult) + 51;
625 dmant = doubleMantissa (dresult);
626 /* Is it too small to be represented by a single-precision
627 * float? */
628 if (dexp <= -155)
629 {
630 FLOAT_TO_INTBITS (result) = 0;
631 return result;
632 }
633 /* Is it a denormalized single-precision float? */
634 if ((dexp <= -127) && (dexp > -155))
635 {
636 /* Only interested in 24 msb bits of the 53-bit double mantissa */
637 fmant = (uint32_t) (dmant >> 29);
638 fovfl = ((uint32_t) (dmant & 0x1FFFFFFF)) << 3;
639 while ((dexp < -127) && ((fmant | fovfl) != 0))
640 {
641 if ((fmant & 1) != 0)
642 {
643 fovfl |= 0x80000000;
644 }
645 fovfl >>= 1;
646 fmant >>= 1;
647 dexp++;
648 }
649 if ((fovfl & 0x80000000) != 0)
650 {
651 if ((fovfl & 0x7FFFFFFC) != 0)
652 {
653 fmant++;
654 }
655 else if ((fmant & 1) != 0)
656 {
657 fmant++;
658 }
659 }
660 else if ((fovfl & 0x40000000) != 0)
661 {
662 if ((fovfl & 0x3FFFFFFC) != 0)
663 {
664 fmant++;
665 }
666 }
667 FLOAT_TO_INTBITS (result) = fmant;
668 }
669 else
670 {
671 result = (jfloat) dresult;
672 }
673 }
674
675 /* Don't go straight to zero as the fact that x*0 = 0 independent
676 * of x might cause the algorithm to produce an incorrect result.
677 * Instead try the min value first and let it fall to zero if need
678 * be.
679 */
680 if (e <= -309 || FLOAT_TO_INTBITS (result) == 0)
681 FLOAT_TO_INTBITS (result) = FLOAT_MINIMUM_INTBITS;
682
683 return floatAlgorithm (env, f, length, e, (jfloat) result);
684 }
685
686 /* The algorithm for the function floatAlgorithm() below can be found
687 * in:
688 *
689 * "How to Read Floating-Point Numbers Accurately", William D.
690 * Clinger, Proceedings of the ACM SIGPLAN '90 Conference on
691 * Programming Language Design and Implementation, June 20-22,
692 * 1990, pp. 92-101.
693 */
floatAlgorithm(JNIEnv * env,uint64_t * f,int32_t length,jint e,jfloat z)694 static jfloat floatAlgorithm(JNIEnv* env, uint64_t* f, int32_t length, jint e, jfloat z) {
695 uint64_t m;
696 int32_t k, comparison, comparison2;
697 uint64_t* x;
698 uint64_t* y;
699 uint64_t* D;
700 uint64_t* D2;
701 int32_t xLength, yLength, DLength, D2Length;
702
703 x = y = D = D2 = 0;
704 xLength = yLength = DLength = D2Length = 0;
705
706 do
707 {
708 m = floatMantissa (z);
709 k = floatExponent (z);
710
711 if (x && x != f)
712 free(x);
713
714 free(y);
715 free(D);
716 free(D2);
717 y = D = D2 = NULL;
718
719 if (e >= 0 && k >= 0)
720 {
721 xLength = sizeOfTenToTheE (e) + length;
722 allocateU64 (x, xLength);
723 memset (x + length, 0, sizeof (uint64_t) * (xLength - length));
724 memcpy (x, f, sizeof (uint64_t) * length);
725 timesTenToTheEHighPrecision (x, xLength, e);
726
727 yLength = (k >> 6) + 2;
728 allocateU64 (y, yLength);
729 memset (y + 1, 0, sizeof (uint64_t) * (yLength - 1));
730 *y = m;
731 simpleShiftLeftHighPrecision (y, yLength, k);
732 }
733 else if (e >= 0)
734 {
735 xLength = sizeOfTenToTheE (e) + length + ((-k) >> 6) + 1;
736 allocateU64 (x, xLength);
737 memset (x + length, 0, sizeof (uint64_t) * (xLength - length));
738 memcpy (x, f, sizeof (uint64_t) * length);
739 timesTenToTheEHighPrecision (x, xLength, e);
740 simpleShiftLeftHighPrecision (x, xLength, -k);
741
742 yLength = 1;
743 allocateU64 (y, 1);
744 *y = m;
745 }
746 else if (k >= 0)
747 {
748 xLength = length;
749 x = f;
750
751 yLength = sizeOfTenToTheE (-e) + 2 + (k >> 6);
752 allocateU64 (y, yLength);
753 memset (y + 1, 0, sizeof (uint64_t) * (yLength - 1));
754 *y = m;
755 timesTenToTheEHighPrecision (y, yLength, -e);
756 simpleShiftLeftHighPrecision (y, yLength, k);
757 }
758 else
759 {
760 xLength = length + ((-k) >> 6) + 1;
761 allocateU64 (x, xLength);
762 memset (x + length, 0, sizeof (uint64_t) * (xLength - length));
763 memcpy (x, f, sizeof (uint64_t) * length);
764 simpleShiftLeftHighPrecision (x, xLength, -k);
765
766 yLength = sizeOfTenToTheE (-e) + 1;
767 allocateU64 (y, yLength);
768 memset (y + 1, 0, sizeof (uint64_t) * (yLength - 1));
769 *y = m;
770 timesTenToTheEHighPrecision (y, yLength, -e);
771 }
772
773 comparison = compareHighPrecision (x, xLength, y, yLength);
774 if (comparison > 0)
775 { /* x > y */
776 DLength = xLength;
777 allocateU64 (D, DLength);
778 memcpy (D, x, DLength * sizeof (uint64_t));
779 subtractHighPrecision (D, DLength, y, yLength);
780 }
781 else if (comparison)
782 { /* y > x */
783 DLength = yLength;
784 allocateU64 (D, DLength);
785 memcpy (D, y, DLength * sizeof (uint64_t));
786 subtractHighPrecision (D, DLength, x, xLength);
787 }
788 else
789 { /* y == x */
790 DLength = 1;
791 allocateU64 (D, 1);
792 *D = 0;
793 }
794
795 D2Length = DLength + 1;
796 allocateU64 (D2, D2Length);
797 m <<= 1;
798 multiplyHighPrecision (D, DLength, &m, 1, D2, D2Length);
799 m >>= 1;
800
801 comparison2 = compareHighPrecision (D2, D2Length, y, yLength);
802 if (comparison2 < 0)
803 {
804 if (comparison < 0 && m == FLOAT_NORMAL_MASK
805 && FLOAT_TO_INTBITS(z) != FLOAT_NORMAL_MASK)
806 {
807 simpleShiftLeftHighPrecision (D2, D2Length, 1);
808 if (compareHighPrecision (D2, D2Length, y, yLength) > 0)
809 {
810 --FLOAT_TO_INTBITS (z);
811 }
812 else
813 {
814 break;
815 }
816 }
817 else
818 {
819 break;
820 }
821 }
822 else if (comparison2 == 0)
823 {
824 if ((m & 1) == 0)
825 {
826 if (comparison < 0 && m == FLOAT_NORMAL_MASK)
827 {
828 --FLOAT_TO_INTBITS (z);
829 }
830 else
831 {
832 break;
833 }
834 }
835 else if (comparison < 0)
836 {
837 --FLOAT_TO_INTBITS (z);
838 break;
839 }
840 else
841 {
842 ++FLOAT_TO_INTBITS (z);
843 break;
844 }
845 }
846 else if (comparison < 0)
847 {
848 --FLOAT_TO_INTBITS (z);
849 }
850 else
851 {
852 if (FLOAT_TO_INTBITS (z) == FLOAT_EXPONENT_MASK)
853 break;
854 ++FLOAT_TO_INTBITS (z);
855 }
856 }
857 while (1);
858
859 if (x && x != f)
860 free(x);
861 free(y);
862 free(D);
863 free(D2);
864 return z;
865
866 OutOfMemory:
867 if (x && x != f)
868 free(x);
869 free(y);
870 free(D);
871 free(D2);
872 jniThrowOutOfMemoryError(env, NULL);
873 return z;
874 }
875
StringToReal_parseFltImpl(JNIEnv * env,jclass,jstring s,jint e)876 static jfloat StringToReal_parseFltImpl(JNIEnv* env, jclass, jstring s, jint e) {
877 ScopedUtfChars str(env, s);
878 if (str.c_str() == NULL) {
879 return 0.0;
880 }
881 return createFloat(env, str.c_str(), e);
882 }
883
StringToReal_parseDblImpl(JNIEnv * env,jclass,jstring s,jint e)884 static jdouble StringToReal_parseDblImpl(JNIEnv* env, jclass, jstring s, jint e) {
885 ScopedUtfChars str(env, s);
886 if (str.c_str() == NULL) {
887 return 0.0;
888 }
889 return createDouble(env, str.c_str(), e);
890 }
891
892 static JNINativeMethod gMethods[] = {
893 NATIVE_METHOD(StringToReal, parseFltImpl, "(Ljava/lang/String;I)F"),
894 NATIVE_METHOD(StringToReal, parseDblImpl, "(Ljava/lang/String;I)D"),
895 };
register_java_lang_StringToReal(JNIEnv * env)896 void register_java_lang_StringToReal(JNIEnv* env) {
897 jniRegisterNativeMethods(env, "java/lang/StringToReal", gMethods, NELEM(gMethods));
898 }
899