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