1 | /*
|
---|
2 | * Copyright (c) 2011 Jiri Zarevucky
|
---|
3 | * All rights reserved.
|
---|
4 | *
|
---|
5 | * Redistribution and use in source and binary forms, with or without
|
---|
6 | * modification, are permitted provided that the following conditions
|
---|
7 | * are met:
|
---|
8 | *
|
---|
9 | * - Redistributions of source code must retain the above copyright
|
---|
10 | * notice, this list of conditions and the following disclaimer.
|
---|
11 | * - Redistributions in binary form must reproduce the above copyright
|
---|
12 | * notice, this list of conditions and the following disclaimer in the
|
---|
13 | * documentation and/or other materials provided with the distribution.
|
---|
14 | * - The name of the author may not be used to endorse or promote products
|
---|
15 | * derived from this software without specific prior written permission.
|
---|
16 | *
|
---|
17 | * THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR
|
---|
18 | * IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
|
---|
19 | * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
|
---|
20 | * IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT,
|
---|
21 | * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT
|
---|
22 | * NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
|
---|
23 | * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
|
---|
24 | * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
|
---|
25 | * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
|
---|
26 | * THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
|
---|
27 | */
|
---|
28 |
|
---|
29 | /** @addtogroup libposix
|
---|
30 | * @{
|
---|
31 | */
|
---|
32 | /** @file Backend for floating point conversions.
|
---|
33 | */
|
---|
34 |
|
---|
35 | #define LIBPOSIX_INTERNAL
|
---|
36 | #define __POSIX_DEF__(x) posix_##x
|
---|
37 |
|
---|
38 | #include "../internal/common.h"
|
---|
39 | #include "libc/stdbool.h"
|
---|
40 | #include "posix/stdlib.h"
|
---|
41 |
|
---|
42 | #include "posix/assert.h"
|
---|
43 | #include "posix/ctype.h"
|
---|
44 | #include "posix/stdint.h"
|
---|
45 | #include "posix/strings.h"
|
---|
46 | #include "posix/errno.h"
|
---|
47 | #include "posix/limits.h"
|
---|
48 |
|
---|
49 | #include "posix/float.h"
|
---|
50 |
|
---|
51 | #ifndef HUGE_VALL
|
---|
52 | #define HUGE_VALL (+1.0l / +0.0l)
|
---|
53 | #endif
|
---|
54 |
|
---|
55 | #ifndef abs
|
---|
56 | #define abs(x) (((x) < 0) ? -(x) : (x))
|
---|
57 | #endif
|
---|
58 |
|
---|
59 | /* If the constants are not defined, use double precision as default. */
|
---|
60 | #ifndef LDBL_MANT_DIG
|
---|
61 | #define LDBL_MANT_DIG 53
|
---|
62 | #endif
|
---|
63 | #ifndef LDBL_MAX_EXP
|
---|
64 | #define LDBL_MAX_EXP 1024
|
---|
65 | #endif
|
---|
66 | #ifndef LDBL_MIN_EXP
|
---|
67 | #define LDBL_MIN_EXP (-1021)
|
---|
68 | #endif
|
---|
69 | #ifndef LDBL_DIG
|
---|
70 | #define LDBL_DIG 15
|
---|
71 | #endif
|
---|
72 | #ifndef LDBL_MIN
|
---|
73 | #define LDBL_MIN 2.2250738585072014E-308
|
---|
74 | #endif
|
---|
75 |
|
---|
76 | /* power functions ************************************************************/
|
---|
77 |
|
---|
78 | #if LDBL_MAX_EXP >= 16384
|
---|
79 | const int MAX_POW5 = 12;
|
---|
80 | #else
|
---|
81 | const int MAX_POW5 = 8;
|
---|
82 | #endif
|
---|
83 |
|
---|
84 | /* The value at index i is approximately 5**(2**i). */
|
---|
85 | long double pow5[] = {
|
---|
86 | 0x5p0l,
|
---|
87 | 0x19p0l,
|
---|
88 | 0x271p0l,
|
---|
89 | 0x5F5E1p0l,
|
---|
90 | 0x2386F26FC1p0l,
|
---|
91 | 0x4EE2D6D415B85ACEF81p0l,
|
---|
92 | 0x184F03E93FF9F4DAA797ED6E38ED6p36l,
|
---|
93 | 0x127748F9301D319BF8CDE66D86D62p185l,
|
---|
94 | 0x154FDD7F73BF3BD1BBB77203731FDp482l,
|
---|
95 | #if LDBL_MAX_EXP >= 16384
|
---|
96 | 0x1C633415D4C1D238D98CAB8A978A0p1076l,
|
---|
97 | 0x192ECEB0D02EA182ECA1A7A51E316p2265l,
|
---|
98 | 0x13D1676BB8A7ABBC94E9A519C6535p4643l,
|
---|
99 | 0x188C0A40514412F3592982A7F0094p9398l,
|
---|
100 | #endif
|
---|
101 | };
|
---|
102 |
|
---|
103 | #if LDBL_MAX_EXP >= 16384
|
---|
104 | const int MAX_POW2 = 15;
|
---|
105 | #else
|
---|
106 | const int MAX_POW2 = 9;
|
---|
107 | #endif
|
---|
108 |
|
---|
109 | /* Powers of two. */
|
---|
110 | long double pow2[] = {
|
---|
111 | 0x1P1l,
|
---|
112 | 0x1P2l,
|
---|
113 | 0x1P4l,
|
---|
114 | 0x1P8l,
|
---|
115 | 0x1P16l,
|
---|
116 | 0x1P32l,
|
---|
117 | 0x1P64l,
|
---|
118 | 0x1P128l,
|
---|
119 | 0x1P256l,
|
---|
120 | 0x1P512l,
|
---|
121 | #if LDBL_MAX_EXP >= 16384
|
---|
122 | 0x1P1024l,
|
---|
123 | 0x1P2048l,
|
---|
124 | 0x1P4096l,
|
---|
125 | 0x1P8192l,
|
---|
126 | #endif
|
---|
127 | };
|
---|
128 |
|
---|
129 | /**
|
---|
130 | * Multiplies a number by a power of five.
|
---|
131 | * The result may be inexact and may not be the best possible approximation.
|
---|
132 | *
|
---|
133 | * @param mant Number to be multiplied.
|
---|
134 | * @param exp Base 5 exponent.
|
---|
135 | * @return mant multiplied by 5**exp
|
---|
136 | */
|
---|
137 | static long double mul_pow5(long double mant, int exp)
|
---|
138 | {
|
---|
139 | if (mant == 0.0l || mant == HUGE_VALL) {
|
---|
140 | return mant;
|
---|
141 | }
|
---|
142 |
|
---|
143 | if (abs(exp) >> (MAX_POW5 + 1) != 0) {
|
---|
144 | /* Too large exponent. */
|
---|
145 | errno = ERANGE;
|
---|
146 | return exp < 0 ? LDBL_MIN : HUGE_VALL;
|
---|
147 | }
|
---|
148 |
|
---|
149 | if (exp < 0) {
|
---|
150 | exp = abs(exp);
|
---|
151 | for (int bit = 0; bit <= MAX_POW5; ++bit) {
|
---|
152 | /* Multiply by powers of five bit-by-bit. */
|
---|
153 | if (((exp >> bit) & 1) != 0) {
|
---|
154 | mant /= pow5[bit];
|
---|
155 | if (mant == 0.0l) {
|
---|
156 | /* Underflow. */
|
---|
157 | mant = LDBL_MIN;
|
---|
158 | errno = ERANGE;
|
---|
159 | break;
|
---|
160 | }
|
---|
161 | }
|
---|
162 | }
|
---|
163 | } else {
|
---|
164 | for (int bit = 0; bit <= MAX_POW5; ++bit) {
|
---|
165 | /* Multiply by powers of five bit-by-bit. */
|
---|
166 | if (((exp >> bit) & 1) != 0) {
|
---|
167 | mant *= pow5[bit];
|
---|
168 | if (mant == HUGE_VALL) {
|
---|
169 | /* Overflow. */
|
---|
170 | errno = ERANGE;
|
---|
171 | break;
|
---|
172 | }
|
---|
173 | }
|
---|
174 | }
|
---|
175 | }
|
---|
176 |
|
---|
177 | return mant;
|
---|
178 | }
|
---|
179 |
|
---|
180 | /**
|
---|
181 | * Multiplies a number by a power of two. This is always exact.
|
---|
182 | *
|
---|
183 | * @param mant Number to be multiplied.
|
---|
184 | * @param exp Base 2 exponent.
|
---|
185 | * @return mant multiplied by 2**exp.
|
---|
186 | */
|
---|
187 | static long double mul_pow2(long double mant, int exp)
|
---|
188 | {
|
---|
189 | if (mant == 0.0l || mant == HUGE_VALL) {
|
---|
190 | return mant;
|
---|
191 | }
|
---|
192 |
|
---|
193 | if (exp > LDBL_MAX_EXP || exp < LDBL_MIN_EXP) {
|
---|
194 | errno = ERANGE;
|
---|
195 | return exp < 0 ? LDBL_MIN : HUGE_VALL;
|
---|
196 | }
|
---|
197 |
|
---|
198 | if (exp < 0) {
|
---|
199 | exp = abs(exp);
|
---|
200 | for (int i = 0; i <= MAX_POW2; ++i) {
|
---|
201 | if (((exp >> i) & 1) != 0) {
|
---|
202 | mant /= pow2[i];
|
---|
203 | if (mant == 0.0l) {
|
---|
204 | mant = LDBL_MIN;
|
---|
205 | errno = ERANGE;
|
---|
206 | break;
|
---|
207 | }
|
---|
208 | }
|
---|
209 | }
|
---|
210 | } else {
|
---|
211 | for (int i = 0; i <= MAX_POW2; ++i) {
|
---|
212 | if (((exp >> i) & 1) != 0) {
|
---|
213 | mant *= pow2[i];
|
---|
214 | if (mant == HUGE_VALL) {
|
---|
215 | errno = ERANGE;
|
---|
216 | break;
|
---|
217 | }
|
---|
218 | }
|
---|
219 | }
|
---|
220 | }
|
---|
221 |
|
---|
222 | return mant;
|
---|
223 | }
|
---|
224 |
|
---|
225 | /* end power functions ********************************************************/
|
---|
226 |
|
---|
227 |
|
---|
228 |
|
---|
229 | /**
|
---|
230 | * Convert decimal string representation of the floating point number.
|
---|
231 | * Function expects the string pointer to be already pointed at the first
|
---|
232 | * digit (i.e. leading optional sign was already consumed by the caller).
|
---|
233 | *
|
---|
234 | * @param sptr Pointer to the storage of the string pointer. Upon successful
|
---|
235 | * conversion, the string pointer is updated to point to the first
|
---|
236 | * unrecognized character.
|
---|
237 | * @return An approximate representation of the input floating-point number.
|
---|
238 | */
|
---|
239 | static long double parse_decimal(const char **sptr)
|
---|
240 | {
|
---|
241 | assert(sptr != NULL);
|
---|
242 | assert (*sptr != NULL);
|
---|
243 |
|
---|
244 | const int DEC_BASE = 10;
|
---|
245 | const char DECIMAL_POINT = '.';
|
---|
246 | const char EXPONENT_MARK = 'e';
|
---|
247 |
|
---|
248 | const char *str = *sptr;
|
---|
249 | long double significand = 0;
|
---|
250 | long exponent = 0;
|
---|
251 |
|
---|
252 | /* number of digits parsed so far */
|
---|
253 | int parsed_digits = 0;
|
---|
254 | bool after_decimal = false;
|
---|
255 |
|
---|
256 | while (isdigit(*str) || (!after_decimal && *str == DECIMAL_POINT)) {
|
---|
257 | if (*str == DECIMAL_POINT) {
|
---|
258 | after_decimal = true;
|
---|
259 | str++;
|
---|
260 | continue;
|
---|
261 | }
|
---|
262 |
|
---|
263 | if (parsed_digits == 0 && *str == '0') {
|
---|
264 | /* Nothing, just skip leading zeros. */
|
---|
265 | } else if (parsed_digits < LDBL_DIG) {
|
---|
266 | significand = significand * DEC_BASE + (*str - '0');
|
---|
267 | parsed_digits++;
|
---|
268 | } else {
|
---|
269 | exponent++;
|
---|
270 | }
|
---|
271 |
|
---|
272 | if (after_decimal) {
|
---|
273 | /* Decrement exponent if we are parsing the fractional part. */
|
---|
274 | exponent--;
|
---|
275 | }
|
---|
276 |
|
---|
277 | str++;
|
---|
278 | }
|
---|
279 |
|
---|
280 | /* exponent */
|
---|
281 | if (tolower(*str) == EXPONENT_MARK) {
|
---|
282 | str++;
|
---|
283 |
|
---|
284 | /* Returns MIN/MAX value on error, which is ok. */
|
---|
285 | long exp = strtol(str, (char **) &str, DEC_BASE);
|
---|
286 |
|
---|
287 | if (exponent > 0 && exp > LONG_MAX - exponent) {
|
---|
288 | exponent = LONG_MAX;
|
---|
289 | } else if (exponent < 0 && exp < LONG_MIN - exponent) {
|
---|
290 | exponent = LONG_MIN;
|
---|
291 | } else {
|
---|
292 | exponent += exp;
|
---|
293 | }
|
---|
294 | }
|
---|
295 |
|
---|
296 | *sptr = str;
|
---|
297 |
|
---|
298 | /* Return multiplied by a power of ten. */
|
---|
299 | return mul_pow2(mul_pow5(significand, exponent), exponent);
|
---|
300 | }
|
---|
301 |
|
---|
302 | /**
|
---|
303 | * Derive a hexadecimal digit from its character representation.
|
---|
304 | *
|
---|
305 | * @param ch Character representation of the hexadecimal digit.
|
---|
306 | * @return Digit value represented by an integer.
|
---|
307 | */
|
---|
308 | static inline int hex_value(char ch)
|
---|
309 | {
|
---|
310 | if (ch <= '9') {
|
---|
311 | return ch - '0';
|
---|
312 | } else {
|
---|
313 | return 10 + tolower(ch) - 'a';
|
---|
314 | }
|
---|
315 | }
|
---|
316 |
|
---|
317 | /**
|
---|
318 | * Convert hexadecimal string representation of the floating point number.
|
---|
319 | * Function expects the string pointer to be already pointed at the first
|
---|
320 | * digit (i.e. leading optional sign and 0x prefix were already consumed
|
---|
321 | * by the caller).
|
---|
322 | *
|
---|
323 | * @param sptr Pointer to the storage of the string pointer. Upon successful
|
---|
324 | * conversion, the string pointer is updated to point to the first
|
---|
325 | * unrecognized character.
|
---|
326 | * @return Representation of the input floating-point number.
|
---|
327 | */
|
---|
328 | static long double parse_hexadecimal(const char **sptr)
|
---|
329 | {
|
---|
330 | assert(sptr != NULL && *sptr != NULL);
|
---|
331 |
|
---|
332 | const int DEC_BASE = 10;
|
---|
333 | const int HEX_BASE = 16;
|
---|
334 | const char DECIMAL_POINT = '.';
|
---|
335 | const char EXPONENT_MARK = 'p';
|
---|
336 |
|
---|
337 | const char *str = *sptr;
|
---|
338 | long double significand = 0;
|
---|
339 | long exponent = 0;
|
---|
340 |
|
---|
341 | /* number of bits parsed so far */
|
---|
342 | int parsed_bits = 0;
|
---|
343 | bool after_decimal = false;
|
---|
344 |
|
---|
345 | while (posix_isxdigit(*str) || (!after_decimal && *str == DECIMAL_POINT)) {
|
---|
346 | if (*str == DECIMAL_POINT) {
|
---|
347 | after_decimal = true;
|
---|
348 | str++;
|
---|
349 | continue;
|
---|
350 | }
|
---|
351 |
|
---|
352 | if (parsed_bits == 0 && *str == '0') {
|
---|
353 | /* Nothing, just skip leading zeros. */
|
---|
354 | } else if (parsed_bits <= LDBL_MANT_DIG) {
|
---|
355 | significand = significand * HEX_BASE + hex_value(*str);
|
---|
356 | parsed_bits += 4;
|
---|
357 | } else {
|
---|
358 | exponent += 4;
|
---|
359 | }
|
---|
360 |
|
---|
361 | if (after_decimal) {
|
---|
362 | exponent -= 4;
|
---|
363 | }
|
---|
364 |
|
---|
365 | str++;
|
---|
366 | }
|
---|
367 |
|
---|
368 | /* exponent */
|
---|
369 | if (tolower(*str) == EXPONENT_MARK) {
|
---|
370 | str++;
|
---|
371 |
|
---|
372 | /* Returns MIN/MAX value on error, which is ok. */
|
---|
373 | long exp = strtol(str, (char **) &str, DEC_BASE);
|
---|
374 |
|
---|
375 | if (exponent > 0 && exp > LONG_MAX - exponent) {
|
---|
376 | exponent = LONG_MAX;
|
---|
377 | } else if (exponent < 0 && exp < LONG_MIN - exponent) {
|
---|
378 | exponent = LONG_MIN;
|
---|
379 | } else {
|
---|
380 | exponent += exp;
|
---|
381 | }
|
---|
382 | }
|
---|
383 |
|
---|
384 | *sptr = str;
|
---|
385 |
|
---|
386 | /* Return multiplied by a power of two. */
|
---|
387 | return mul_pow2(significand, exponent);
|
---|
388 | }
|
---|
389 |
|
---|
390 | /**
|
---|
391 | * Converts a string representation of a floating-point number to
|
---|
392 | * its native representation. Largely POSIX compliant, except for
|
---|
393 | * locale differences (always uses '.' at the moment) and rounding.
|
---|
394 | * Decimal strings are NOT guaranteed to be correctly rounded. This function
|
---|
395 | * should return a good enough approximation for most purposes but if you
|
---|
396 | * depend on a precise conversion, use hexadecimal representation.
|
---|
397 | * Hexadecimal strings are currently always rounded towards zero, regardless
|
---|
398 | * of the current rounding mode.
|
---|
399 | *
|
---|
400 | * @param nptr Input string.
|
---|
401 | * @param endptr If non-NULL, *endptr is set to the position of the first
|
---|
402 | * unrecognized character.
|
---|
403 | * @return An approximate representation of the input floating-point number.
|
---|
404 | */
|
---|
405 | long double posix_strtold(const char *restrict nptr, char **restrict endptr)
|
---|
406 | {
|
---|
407 | assert(nptr != NULL);
|
---|
408 |
|
---|
409 | const int RADIX = '.';
|
---|
410 |
|
---|
411 | /* minus sign */
|
---|
412 | bool negative = false;
|
---|
413 | /* current position in the string */
|
---|
414 | int i = 0;
|
---|
415 |
|
---|
416 | /* skip whitespace */
|
---|
417 | while (isspace(nptr[i])) {
|
---|
418 | i++;
|
---|
419 | }
|
---|
420 |
|
---|
421 | /* parse sign */
|
---|
422 | switch (nptr[i]) {
|
---|
423 | case '-':
|
---|
424 | negative = true;
|
---|
425 | /* fallthrough */
|
---|
426 | case '+':
|
---|
427 | i++;
|
---|
428 | }
|
---|
429 |
|
---|
430 | /* check for NaN */
|
---|
431 | if (posix_strncasecmp(&nptr[i], "nan", 3) == 0) {
|
---|
432 | // FIXME: return NaN
|
---|
433 | // TODO: handle the parenthesised case
|
---|
434 |
|
---|
435 | if (endptr != NULL) {
|
---|
436 | *endptr = (char *) nptr;
|
---|
437 | }
|
---|
438 | errno = EINVAL;
|
---|
439 | return 0;
|
---|
440 | }
|
---|
441 |
|
---|
442 | /* check for Infinity */
|
---|
443 | if (posix_strncasecmp(&nptr[i], "inf", 3) == 0) {
|
---|
444 | i += 3;
|
---|
445 | if (posix_strncasecmp(&nptr[i], "inity", 5) == 0) {
|
---|
446 | i += 5;
|
---|
447 | }
|
---|
448 |
|
---|
449 | if (endptr != NULL) {
|
---|
450 | *endptr = (char *) &nptr[i];
|
---|
451 | }
|
---|
452 | return negative ? -HUGE_VALL : +HUGE_VALL;
|
---|
453 | }
|
---|
454 |
|
---|
455 | /* check for a hex number */
|
---|
456 | if (nptr[i] == '0' && tolower(nptr[i + 1]) == 'x' &&
|
---|
457 | (posix_isxdigit(nptr[i + 2]) ||
|
---|
458 | (nptr[i + 2] == RADIX && posix_isxdigit(nptr[i + 3])))) {
|
---|
459 | i += 2;
|
---|
460 |
|
---|
461 | const char *ptr = &nptr[i];
|
---|
462 | /* this call sets errno if appropriate. */
|
---|
463 | long double result = parse_hexadecimal(&ptr);
|
---|
464 | if (endptr != NULL) {
|
---|
465 | *endptr = (char *) ptr;
|
---|
466 | }
|
---|
467 | return negative ? -result : result;
|
---|
468 | }
|
---|
469 |
|
---|
470 | /* check for a decimal number */
|
---|
471 | if (isdigit(nptr[i]) || (nptr[i] == RADIX && isdigit(nptr[i + 1]))) {
|
---|
472 | const char *ptr = &nptr[i];
|
---|
473 | /* this call sets errno if appropriate. */
|
---|
474 | long double result = parse_decimal(&ptr);
|
---|
475 | if (endptr != NULL) {
|
---|
476 | *endptr = (char *) ptr;
|
---|
477 | }
|
---|
478 | return negative ? -result : result;
|
---|
479 | }
|
---|
480 |
|
---|
481 | /* nothing to parse */
|
---|
482 | if (endptr != NULL) {
|
---|
483 | *endptr = (char *) nptr;
|
---|
484 | }
|
---|
485 | errno = EINVAL;
|
---|
486 | return 0;
|
---|
487 | }
|
---|
488 |
|
---|
489 | /** @}
|
---|
490 | */
|
---|