1 /********************************************************************/
2 /* */
3 /* flt_rtl.c Primitive actions for the float type. */
4 /* Copyright (C) 1989 - 2021 Thomas Mertes */
5 /* */
6 /* This file is part of the Seed7 Runtime Library. */
7 /* */
8 /* The Seed7 Runtime Library is free software; you can */
9 /* redistribute it and/or modify it under the terms of the GNU */
10 /* Lesser General Public License as published by the Free Software */
11 /* Foundation; either version 2.1 of the License, or (at your */
12 /* option) any later version. */
13 /* */
14 /* The Seed7 Runtime Library is distributed in the hope that it */
15 /* will be useful, but WITHOUT ANY WARRANTY; without even the */
16 /* implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR */
17 /* PURPOSE. See the GNU Lesser General Public License for more */
18 /* details. */
19 /* */
20 /* You should have received a copy of the GNU Lesser General */
21 /* Public License along with this program; if not, write to the */
22 /* Free Software Foundation, Inc., 51 Franklin Street, */
23 /* Fifth Floor, Boston, MA 02110-1301, USA. */
24 /* */
25 /* Module: Seed7 Runtime Library */
26 /* File: seed7/src/flt_rtl.c */
27 /* Changes: 1993, 1994, 2005, 2009 - 2021 Thomas Mertes */
28 /* Content: Primitive actions for the float type. */
29 /* */
30 /********************************************************************/
31
32 #define LOG_FUNCTIONS 0
33 #define VERBOSE_EXCEPTIONS 0
34
35 #include "version.h"
36
37 #include "stdlib.h"
38 #include "stdio.h"
39 #include "string.h"
40 #include "limits.h"
41 #include "math.h"
42 #include "float.h"
43 #include "ctype.h"
44 #include "errno.h"
45
46 #include "common.h"
47 #include "data_rtl.h"
48 #include "heaputl.h"
49 #include "striutl.h"
50 #include "rtl_err.h"
51 #include "int_rtl.h"
52
53 #undef EXTERN
54 #define EXTERN
55 #include "flt_rtl.h"
56
57
58 #define USE_STRTOD 1
59 #define MAX_CSTRI_BUFFER_LEN 25
60 #define IPOW_EXPONENTIATION_BY_SQUARING 1
61 #define PRECISION_BUFFER_LEN 1000
62 /* The 3 additional chars below are for: "-1.". */
63 #define FLT_DGTS_ADDITIONAL_CHARS STRLEN("-1.")
64 #define FLT_DGTS_LEN (FLT_DGTS_ADDITIONAL_CHARS + DOUBLE_MAX_EXP10)
65 /* The 5 additional chars below are for: "-1.e+". */
66 #define FLT_SCI_ADDITIONAL_CHARS STRLEN("-1.e+")
67 #define FLT_SCI_LEN (FLT_SCI_ADDITIONAL_CHARS + MAX_PRINTED_EXPONENT_DIGITS)
68
69 #if FLOATTYPE_DOUBLE
70 #define MIN_NEGATIVE_SUBNORMAL -2.2250738585072009e-308
71 #define MAX_POSITIVE_SUBNORMAL 2.2250738585072009e-308
72 #define FREXP_SUBNORMAL_FACTOR 1.8446744073709552e+19 /* == 0x1p64 */
73 #else
74 #define MIN_NEGATIVE_SUBNORMAL -1.1754942106924411e-38
75 #define MAX_POSITIVE_SUBNORMAL 1.1754942106924411e-38
76 #define FREXP_SUBNORMAL_FACTOR 4.294967296e+09 /* == 0x1p32 */
77 #endif
78
79 /* Natural logarithm of 2: */
80 #define LN2 0.693147180559945309417232121458176568075500134360255254120680009493393
81
82 #if FLOAT_ZERO_DIV_ERROR
83 const rtlValueUnion f_const[] =
84 #if FLOATTYPE_DOUBLE
85 {{GENERIC_SUFFIX(0xfff8000000000000)},
86 {GENERIC_SUFFIX(0x7ff0000000000000)},
87 {GENERIC_SUFFIX(0xfff0000000000000)}};
88 #else
89 {{0xffc00000}, {0x7f800000}, {0xff800000}};
90 #endif
91 #endif
92
93 #if USE_NEGATIVE_ZERO_BITPATTERN
94 static const rtlValueUnion neg_zero_const =
95 #if FLOATTYPE_DOUBLE
96 {GENERIC_SUFFIX(0x8000000000000000)};
97 #else
98 {0x80000000};
99 #endif
100 #endif
101
102 floatType negativeZero;
103
104 #if !PRINTF_SUPPORTS_VARIABLE_FORMATS
105 static const const_cstriType fmt_e[] = {
106 "%1.0e", "%1.1e", "%1.2e", "%1.3e", "%1.4e", "%1.5e",
107 "%1.6e", "%1.7e", "%1.8e", "%1.9e", "%1.10e", "%1.11e",
108 "%1.12e", "%1.13e", "%1.14e", "%1.15e", "%1.16e", "%1.17e",
109 "%1.18e", "%1.19e", "%1.20e", "%1.21e", "%1.22e", "%1.23e",
110 "%1.24e", "%1.25e", "%1.26e", "%1.27e", "%1.28e"};
111
112 static const const_cstriType fmt_f[] = {
113 "%1.0f", "%1.1f", "%1.2f", "%1.3f", "%1.4f", "%1.5f",
114 "%1.6f", "%1.7f", "%1.8f", "%1.9f", "%1.10f", "%1.11f",
115 "%1.12f", "%1.13f", "%1.14f", "%1.15f", "%1.16f", "%1.17f",
116 "%1.18f", "%1.19f", "%1.20f", "%1.21f", "%1.22f", "%1.23f",
117 "%1.24f", "%1.25f", "%1.26f", "%1.27f", "%1.28f"};
118 #endif
119
120
121
122 /**
123 * Setup floating point computations and initialize constants.
124 * This function must be called before doing any floating
125 * point computations.
126 */
setupFloat(void)127 void setupFloat (void)
128
129 {
130 #if !USE_NEGATIVE_ZERO_BITPATTERN
131 floatType zero = 0.0;
132 floatType plusInf;
133 #endif
134
135 /* setupFloat */
136 #ifdef TURN_OFF_FP_EXCEPTIONS
137 _control87(MCW_EM, MCW_EM);
138 #endif
139 #if USE_NEGATIVE_ZERO_BITPATTERN
140 negativeZero = neg_zero_const.floatValue;
141 #else
142 plusInf = 1.0 / zero;
143 negativeZero = -1.0 / plusInf;
144 #endif
145 } /* setupFloat */
146
147
148
149 #ifdef DEFINE_MATHERR_FUNCTION
matherr(struct _exception * a)150 int matherr (struct _exception *a)
151
152 { /* matherr */
153 a->retval = a->retval;
154 return 1;
155 } /* matherr */
156 #endif
157
158
159
160 #ifdef DEFINE__MATHERR_FUNCTION
_matherr(struct _exception * a)161 int _matherr (struct _exception *a)
162
163 { /* _matherr */
164 a->retval = a->retval;
165 return 1;
166 } /* _matherr */
167 #endif
168
169
170
171 /**
172 * Get the integer mantissa and a binary exponent from a double.
173 * @param doubleValue Value that should be split into mantissa
174 * and exponent (NaN, Infinity and -Infinity are not allowed).
175 * @param binaryExponent Destination for the binary exponent.
176 * @return the integer mantissa.
177 */
getMantissaAndExponent(double doubleValue,int * binaryExponent)178 int64Type getMantissaAndExponent (double doubleValue, int *binaryExponent)
179
180 {
181 double mantissa;
182 #if FREXP_FUNCTION_OKAY
183 int exponent;
184 #else
185 intType exponent;
186 #endif
187 int64Type intMantissa;
188
189 /* getMantissaAndExponent */
190 #if FREXP_FUNCTION_OKAY
191 mantissa = frexp(doubleValue, &exponent);
192 #else
193 mantissa = fltDecompose(doubleValue, &exponent);
194 #endif
195 intMantissa = (int64Type) (mantissa * DOUBLE_MANTISSA_FACTOR);
196 *binaryExponent = (int) (exponent - DOUBLE_MANTISSA_SHIFT);
197 return intMantissa;
198 } /* getMantissaAndExponent */
199
200
201
202 /**
203 * Create double value from an integer mantissa and a binary exponent.
204 * @param intMantissa Integer mantissa of double value.
205 * @param binaryExponent Binary exponent of double value.
206 * @return the result of flt(intMantissa) * 2.0 ** binaryExponent.
207 */
setMantissaAndExponent(int64Type intMantissa,int binaryExponent)208 double setMantissaAndExponent (int64Type intMantissa, int binaryExponent)
209
210 {
211 double mantissa;
212 int exponent;
213 double doubleValue;
214
215 /* setMantissaAndExponent */
216 mantissa = (double) intMantissa;
217 exponent = binaryExponent;
218 doubleValue = fltLdexp(mantissa, exponent);
219 return doubleValue;
220 } /* setMantissaAndExponent */
221
222
223
224 /**
225 * Write the decimal representation of a double to a buffer.
226 * The result in buffer uses the style [-]ddd.ddd where there is at least
227 * one digit before and after the decimal point. The number of digits
228 * after the decimal point is determined automatically. Except for the
229 * case if there is only one zero digit after the decimal point,
230 * the last digit is never zero. Negative zero (-0.0) and positive
231 * zero (+0.0) are both converted to "0.0".
232 * @param doubleValue Number to be converted (NaN, Infinity and
233 * -Infinity are not allowed).
234 * @param largeNumber If abs(doubleValue) > largeValue holds
235 * the format %1.1f is used (largeNumber is either
236 * DOUBLE_STR_LARGE_NUMBER or FLOAT_STR_LARGE_NUMBER).
237 * @param format Format to be used if abs(doubleValue) <= largeValue
238 * holds (format is eitner FMT_E_DBL or FMT_E_FLT).
239 * @param buffer Destination buffer for the decimal representation.
240 * @return the number of characters in the destination buffer.
241 */
doubleToCharBuffer(const double doubleValue,const double largeNumber,const char * format,char * buffer)242 memSizeType doubleToCharBuffer (const double doubleValue,
243 const double largeNumber, const char *format, char *buffer)
244
245 {
246 int decimalExponent;
247 memSizeType pos;
248 memSizeType start;
249 memSizeType scale;
250 memSizeType len;
251
252 /* doubleToCharBuffer */
253 logFunction(printf("doubleToCharBuffer(" FMT_E_DBL ", " FMT_E_DBL
254 ", \"%s\", *)\n",
255 doubleValue, largeNumber, format););
256 #if FLOAT_ZERO_COMPARISON_OKAY
257 if (doubleValue == 0.0) {
258 #else
259 if (doubleValue == 0.0 || fltIsNegativeZero(doubleValue)) {
260 #endif
261 memcpy(buffer, "0.0", 3);
262 len = 3;
263 } else if (doubleValue < -largeNumber || doubleValue > largeNumber) {
264 len = (memSizeType) sprintf(buffer, "%1.1f", doubleValue);
265 } else {
266 len = (memSizeType) sprintf(buffer, format, doubleValue);
267 /* printf("buffer: \"%s\"\n", buffer); */
268 /* Subtract two more chars for sign and letter 'e': */
269 len -= MIN_PRINTED_EXPONENT_DIGITS + 2;
270 while (buffer[len] != 'e') {
271 len--;
272 } /* while */
273 pos = len + 2; /* skip the exponent sign */
274 decimalExponent = (int) buffer[pos] - '0';
275 pos++;
276 while (buffer[pos] >= '0') {
277 decimalExponent *= 10;
278 decimalExponent += (int) buffer[pos] - '0';
279 pos++;
280 } /* while */
281 if (buffer[len + 1] == '-') {
282 decimalExponent = -decimalExponent;
283 } /* if */
284 /* printf("decimalExponent: %ld\n", decimalExponent); */
285 do {
286 len--;
287 } while (buffer[len] == '0');
288 len++;
289 /* printf("len: " FMT_U_MEM "\n", len); */
290 start = buffer[0] == '-';
291 if (decimalExponent > 0) {
292 scale = (memSizeType) decimalExponent;
293 if (scale >= len - start - 2) {
294 /* No digits <> '0' after the decimal point */
295 memmove(&buffer[start + 1], &buffer[start + 2], len - start - 2);
296 memset(&buffer[len - 1], '0', scale + 2 - len + start);
297 memcpy(&buffer[start + scale + 1], ".0", 2);
298 len = start + scale + 3;
299 } else {
300 memmove(&buffer[start + 1], &buffer[start + 2], scale);
301 buffer[start + scale + 1] = '.';
302 } /* if */
303 } else if (decimalExponent < 0) {
304 /* No digits <> '0' before the decimal point */
305 scale = (memSizeType) -decimalExponent;
306 memmove(&buffer[start + 2 + scale], &buffer[start + 2], len - start - 2);
307 buffer[start + 1 + scale] = buffer[start];
308 memset(&buffer[start + 2], '0', scale - 1);
309 memcpy(&buffer[start], "0.", 2);
310 len += scale;
311 } else if (buffer[len - 1] == '.') {
312 /* Make sure that there is one digit after the decimal point */
313 len++;
314 } /* if */
315 /* printf("len: " FMT_U_MEM "\n", len);
316 buffer[len] = '\0';
317 printf("buffer: \"%s\"\n", buffer); */
318 } /* if */
319 logFunction(printf("doubleToCharBuffer(" FMT_E_DBL ", " FMT_E_DBL
320 ", \"%s\", \"%.*s\") --> " FMT_U_MEM "\n",
321 doubleValue, largeNumber, format, (int) len,
322 buffer, len););
323 return len;
324 } /* doubleToCharBuffer */
325
326
327
328 /**
329 * Compare two float numbers.
330 * Because fltCmp is used to sort float values, a unique
331 * sort sequence of all values is needed. Therefore fltCmp
332 * considers NaN as equal to itself and greater than Infinity.
333 * Negative zero (-0.0) is considered by fltCmp to be equal to
334 * positive zero (+0.0). This conforms to the behavior of all
335 * other float comparisons with zero.
336 * @return -1, 0 or 1 if the first argument is considered to be
337 * respectively less than, equal to, or greater than the
338 * second.
339 */
340 intType fltCmp (floatType number1, floatType number2)
341
342 {
343 intType signumValue;
344
345 /* fltCmp */
346 logFunction(printf("fltCmp(" FMT_E ", " FMT_E ")\n",
347 number1, number2););
348 #if FLOAT_COMPARISON_OKAY
349 if (number1 < number2) {
350 signumValue = -1;
351 } else if (number1 > number2) {
352 signumValue = 1;
353 } else {
354 /* The expression isnan(NaN) can return any value except 0. */
355 signumValue = (os_isnan(number1) != 0) - (os_isnan(number2) != 0);
356 } /* if */
357 #else
358 if (unlikely(os_isnan(number1))) {
359 /* The expression isnan(NaN) can return any value except 0. */
360 signumValue = os_isnan(number2) == 0;
361 } else if (unlikely(os_isnan(number2))) {
362 signumValue = -1;
363 #if !FLOAT_ZERO_COMPARISON_OKAY
364 } else if (fltLt(number1, number2)) {
365 signumValue = -1;
366 } else {
367 signumValue = fltGt(number1, number2);
368 } /* if */
369 #else
370 } else if (number1 < number2) {
371 signumValue = -1;
372 } else {
373 signumValue = number1 > number2;
374 } /* if */
375 #endif
376 #endif
377 logFunction(printf("fltCmp --> " FMT_D "\n", signumValue););
378 return signumValue;
379 } /* fltCmp */
380
381
382
383 /**
384 * Reinterpret the generic parameters as floatType and call fltCmp.
385 * Function pointers in C programs generated by the Seed7 compiler
386 * may point to this function. This assures correct behaviour even
387 * if sizeof(genericType) != sizeof(floatType).
388 * @return -1, 0 or 1 if the first argument is considered to be
389 * respectively less than, equal to, or greater than the
390 * second.
391 */
fltCmpGeneric(const genericType value1,const genericType value2)392 intType fltCmpGeneric (const genericType value1, const genericType value2)
393
394 { /* fltCmpGeneric */
395 return fltCmp(((const_rtlObjectType *) &value1)->value.floatValue,
396 ((const_rtlObjectType *) &value2)->value.floatValue);
397 } /* fltCmpGeneric */
398
399
400
401 /**
402 * Reinterpret the generic parameters as floatType and assign source to dest.
403 * Function pointers in C programs generated by the Seed7 compiler
404 * may point to this function. This assures correct behaviour even
405 * if sizeof(genericType) != sizeof(floatType).
406 */
fltCpyGeneric(genericType * const dest,const genericType source)407 void fltCpyGeneric (genericType *const dest, const genericType source)
408
409 { /* fltCpyGeneric */
410 ((rtlObjectType *) dest)->value.floatValue =
411 ((const_rtlObjectType *) &source)->value.floatValue;
412 } /* fltCpyGeneric */
413
414
415
416 #if !FREXP_FUNCTION_OKAY
417 /**
418 * Decompose float into normalized fraction and integral exponent for 2.
419 * If the argument (number) is 0.0, -0.0, Infinity, -Infinity or NaN the
420 * fraction is set to the argument and the exponent is set to 0.
421 * For all other arguments the fraction is set to an absolute value
422 * between 0.5(included) and 1.0(excluded) and the exponent is set such
423 * that number = fraction * 2.0 ** exponent holds.
424 * @param exponent Destination for the exponent of the decomposed number.
425 * @param number Number to be decomposed into fraction and exponent.
426 * @return the normalized fraction.
427 */
fltDecompose(const floatType number,intType * const exponent)428 floatType fltDecompose (const floatType number, intType *const exponent)
429
430 {
431 int exponent_var = 0;
432 floatType fraction;
433
434 /* fltDecompose */
435 logFunction(printf("fltDecompose(" FMT_E ", *)\n", number););
436 #if !FREXP_INFINITY_NAN_OKAY
437 if (unlikely(os_isnan(number) ||
438 number == POSITIVE_INFINITY ||
439 number == NEGATIVE_INFINITY)) {
440 fraction = number;
441 *exponent = 0;
442 } else
443 #endif
444 #if !FREXP_SUBNORMAL_OKAY
445 if (number >= MIN_NEGATIVE_SUBNORMAL &&
446 number <= MAX_POSITIVE_SUBNORMAL && number != 0.0) {
447 /* Subnormal number */
448 fraction = frexp(number * FREXP_SUBNORMAL_FACTOR, &exponent_var);
449 *exponent = (intType) exponent_var - FLOATTYPE_SIZE;
450 } else
451 #endif
452 {
453 fraction = frexp(number, &exponent_var);
454 *exponent = (intType) exponent_var;
455 } /* if */
456 logFunction(printf("fltDecompose(" FMT_E ", " FMT_D ") --> " FMT_E "\n",
457 number, *exponent, fraction););
458 return fraction;
459 } /* fltDecompose */
460 #endif
461
462
463
464 /**
465 * Convert a float to a string in decimal fixed point notation.
466 * The number is rounded to the specified number of digits ('precision').
467 * Halfway cases are rounded away from zero. Except for a 'precision' of
468 * zero the representation has a decimal point and at least one digit
469 * before and after the decimal point. Negative numbers are preceded by
470 * a minus sign (e.g.: "-1.25"). If all digits in the result are 0 a
471 * possible negative sign is omitted.
472 * @param precision Number of digits after the decimal point.
473 * If the 'precision' is zero the decimal point is omitted.
474 * @return the string result of the conversion.
475 * @exception RANGE_ERROR If the 'precision' is negative.
476 * @exception MEMORY_ERROR Not enough memory to represent the result.
477 */
fltDgts(floatType number,intType precision)478 striType fltDgts (floatType number, intType precision)
479
480 {
481 char buffer_1[PRECISION_BUFFER_LEN + FLT_DGTS_LEN + NULL_TERMINATION_LEN];
482 char *buffer = buffer_1;
483 const_cstriType buffer_ptr;
484 /* The 4 additional chars below are for "%1.f". */
485 char form_buffer[INTTYPE_DECIMAL_SIZE + STRLEN("%1.f") + NULL_TERMINATION_LEN];
486 memSizeType pos;
487 memSizeType len;
488 striType result;
489
490 /* fltDgts */
491 logFunction(printf("fltDgts(" FMT_E ", " FMT_D ")\n", number, precision););
492 if (unlikely(precision < 0)) {
493 logError(printf("fltDgts(" FMT_E ", " FMT_D "): Precision negative.\n",
494 number, precision););
495 raise_error(RANGE_ERROR);
496 result = NULL;
497 } else if (unlikely(precision > PRECISION_BUFFER_LEN &&
498 ((uintType) precision > MAX_CSTRI_LEN - FLT_DGTS_LEN ||
499 !ALLOC_CSTRI(buffer, (memSizeType) precision + FLT_DGTS_LEN)))) {
500 raise_error(MEMORY_ERROR);
501 result = NULL;
502 } else {
503 if (os_isnan(number)) {
504 buffer_ptr = "NaN";
505 len = STRLEN("NaN");
506 } else if (number == POSITIVE_INFINITY) {
507 buffer_ptr = "Infinity";
508 len = STRLEN("Infinity");
509 } else if (number == NEGATIVE_INFINITY) {
510 buffer_ptr = "-Infinity";
511 len = STRLEN("-Infinity");
512 } else {
513 #ifdef LIMIT_FMT_F_MAXIMUM_FLOAT_PRECISION
514 if (unlikely(precision > PRINTF_FMT_F_MAXIMUM_FLOAT_PRECISION)) {
515 len = (memSizeType) sprintf(buffer, "%1."
516 LIMIT_FMT_F_MAXIMUM_FLOAT_PRECISION "f", number);
517 } else
518 #endif
519 #if PRINTF_SUPPORTS_VARIABLE_FORMATS
520 if (unlikely(precision > INT_MAX)) {
521 sprintf(form_buffer, "%%1." FMT_D "f", precision);
522 len = (memSizeType) sprintf(buffer, form_buffer, number);
523 } else {
524 len = (memSizeType) sprintf(buffer, "%1.*f", (int) precision, number);
525 } /* if */
526 #else
527 if (unlikely(precision >= sizeof(fmt_f) / sizeof(char *))) {
528 sprintf(form_buffer, "%%1." FMT_D "f", precision);
529 len = (memSizeType) sprintf(buffer, form_buffer, number);
530 } else {
531 len = (memSizeType) sprintf(buffer, fmt_f[precision], number);
532 } /* if */
533 #endif
534 #ifdef PRINTF_FMT_F_MAXIMUM_FLOAT_PRECISION
535 /* Some C run-time libraries do not have a fixed maximum */
536 /* for the float precision of printf(). Instead the actual */
537 /* precision varies with each call of printf(). Up to a */
538 /* precision of PRINTF_FMT_F_MAXIMUM_FLOAT_PRECISION */
539 /* printf() will always work ok. */
540 if (precision > PRINTF_FMT_F_MAXIMUM_FLOAT_PRECISION) {
541 pos = (memSizeType) (strchr(buffer, '.') - buffer);
542 if ((memSizeType) precision >= len - pos) {
543 memSizeType numZeros = (memSizeType) precision - (len - pos) + 1;
544 memset(&buffer[len], '0', numZeros);
545 len += numZeros;
546 buffer[len] = '\0';
547 } /* if */
548 } /* if */
549 #endif
550 buffer_ptr = buffer;
551 if (buffer[0] == '-' && buffer[1] == '0') {
552 /* All forms of -0 are converted to 0 */
553 if (buffer[2] == '.') {
554 pos = 3;
555 while (buffer[pos] == '0') {
556 pos++;
557 } /* while */
558 if (buffer[pos] == '\0') {
559 buffer_ptr++;
560 len--;
561 } /* if */
562 } else if (buffer[2] == '\0') {
563 buffer_ptr++;
564 len--;
565 } /* if */
566 } /* if */
567 } /* if */
568 result = cstri_buf_to_stri(buffer_ptr, len);
569 if (buffer != buffer_1) {
570 UNALLOC_CSTRI(buffer, (memSizeType) precision + FLT_DGTS_LEN);
571 } /* if */
572 if (unlikely(result == NULL)) {
573 raise_error(MEMORY_ERROR);
574 } /* if */
575 } /* if */
576 logFunction(printf("fltDgts(" FMT_E ", " FMT_D ") --> \"%s\"\n",
577 number, precision, striAsUnquotedCStri(result)););
578 return result;
579 } /* fltDgts */
580
581
582
583 #if !FLOAT_COMPARISON_OKAY
584 /**
585 * Check if two float numbers are equal.
586 * According to IEEE 754 a NaN is not equal to any float value.
587 * Therefore 'NaN = any_value' and 'any_value = NaN'
588 * always return FALSE. Even 'NaN = NaN' returns FALSE.
589 * @return TRUE if both numbers are equal, FALSE otherwise.
590 */
fltEq(floatType number1,floatType number2)591 boolType fltEq (floatType number1, floatType number2)
592
593 {
594 boolType isEqual;
595
596 /* fltEq */
597 logFunction(printf("fltEq(" FMT_E ", " FMT_E ")\n",
598 number1, number2););
599 #if !FLOAT_NAN_COMPARISON_OKAY
600 if (os_isnan(number1) || os_isnan(number2)) {
601 isEqual = FALSE;
602 } else
603 #endif
604 #if !FLOAT_ZERO_COMPARISON_OKAY
605 if (fltIsNegativeZero(number1)) {
606 if (fltIsNegativeZero(number2)) {
607 isEqual = TRUE;
608 } else {
609 isEqual = 0.0 == number2;
610 } /* if */
611 } else if (fltIsNegativeZero(number2)) {
612 isEqual = number1 == 0.0;
613 } else
614 #endif
615 {
616 isEqual = number1 == number2;
617 } /* if */
618 logFunction(printf("fltEq --> %d\n", isEqual););
619 return isEqual;
620 } /* fltEq */
621 #endif
622
623
624
625 #if !EXP_FUNCTION_OKAY
626 /**
627 * Compute Euler's number e raised to the power of x.
628 * @return e raised to the power of x.
629 */
fltExp(floatType exponent)630 floatType fltExp (floatType exponent)
631
632 {
633 floatType power;
634
635 /* fltExp */
636 logFunction(printf("fltExp(" FMT_E ", *)\n", exponent););
637 #if !EXP_OF_NAN_OKAY
638 if (unlikely(os_isnan(exponent))) {
639 power = NOT_A_NUMBER;
640 } else
641 #endif
642 {
643 power = exp(exponent);
644 } /* if */
645 logFunction(printf("fltExp --> " FMT_E "\n", power););
646 return power;
647 } /* fltExp */
648 #endif
649
650
651
652 #if !HAS_EXPM1
653 /**
654 * Compute exp(x) - 1.0 (subtract one from e raised to the power of x).
655 * This function is used, if the function expm1() is missing from
656 * the C math library. This simple implementation is not accurate if
657 * the value of x is near zero.
658 * @return exp(x) - 1.0
659 */
fltExpM1(floatType exponent)660 floatType fltExpM1 (floatType exponent)
661
662 {
663 floatType power;
664
665 /* fltExpM1 */
666 logFunction(printf("fltExpM1(" FMT_E ", *)\n", exponent););
667 power = fltExp(exponent) - 1;
668 logFunction(printf("fltExpM1 --> " FMT_E "\n", power););
669 return power;
670 } /* fltExpM1 */
671 #endif
672
673
674
675 #if !FLOAT_COMPARISON_OKAY
676 /**
677 * Check if 'number1' is greater than or equal to 'number2'.
678 * According to IEEE 754 a NaN is neither less than,
679 * equal to, nor greater than any value, including itself.
680 * If 'number1' or 'number2' is NaN, the result is FALSE.
681 * @return TRUE if 'number1' is greater than or equal to 'number2',
682 * FALSE otherwise.
683 */
fltGe(floatType number1,floatType number2)684 boolType fltGe (floatType number1, floatType number2)
685
686 {
687 boolType isGreaterEqual;
688
689 /* fltGe */
690 logFunction(printf("fltGe(" FMT_E ", " FMT_E ")\n",
691 number1, number2););
692 #if !FLOAT_NAN_COMPARISON_OKAY
693 if (os_isnan(number1) || os_isnan(number2)) {
694 isGreaterEqual = FALSE;
695 } else
696 #endif
697 #if !FLOAT_ZERO_COMPARISON_OKAY
698 if (fltIsNegativeZero(number1)) {
699 if (fltIsNegativeZero(number2)) {
700 isGreaterEqual = TRUE;
701 } else {
702 isGreaterEqual = 0.0 >= number2;
703 } /* if */
704 } else if (fltIsNegativeZero(number2)) {
705 isGreaterEqual = number1 >= 0.0;
706 } else
707 #endif
708 {
709 isGreaterEqual = number1 >= number2;
710 } /* if */
711 logFunction(printf("fltGe --> %d\n", isGreaterEqual););
712 return isGreaterEqual;
713 } /* fltGe */
714 #endif
715
716
717
718 #if !FLOAT_COMPARISON_OKAY
719 /**
720 * Check if 'number1' is greater than 'number2'.
721 * According to IEEE 754 a NaN is neither less than,
722 * equal to, nor greater than any value, including itself.
723 * If 'number1' or 'number2' is NaN, the result is FALSE.
724 * @return TRUE if 'number1' is greater than 'number2',
725 * FALSE otherwise.
726 */
fltGt(floatType number1,floatType number2)727 boolType fltGt (floatType number1, floatType number2)
728
729 {
730 boolType isGreaterThan;
731
732 /* fltGt */
733 logFunction(printf("fltGt(" FMT_E ", " FMT_E ")\n",
734 number1, number2););
735 #if !FLOAT_NAN_COMPARISON_OKAY
736 if (os_isnan(number1) || os_isnan(number2)) {
737 isGreaterThan = FALSE;
738 } else
739 #endif
740 #if !FLOAT_ZERO_COMPARISON_OKAY
741 if (fltIsNegativeZero(number1)) {
742 if (fltIsNegativeZero(number2)) {
743 isGreaterThan = FALSE;
744 } else {
745 isGreaterThan = 0.0 > number2;
746 } /* if */
747 } else if (fltIsNegativeZero(number2)) {
748 isGreaterThan = number1 > 0.0;
749 } else
750 #endif
751 {
752 isGreaterThan = number1 > number2;
753 } /* if */
754 logFunction(printf("fltGt --> %d\n", isGreaterThan););
755 return isGreaterThan;
756 } /* fltGt */
757 #endif
758
759
760
761 /**
762 * Compute the exponentiation of a float 'base' with an integer 'exponent'.
763 * A ** 0 returns 1.0
764 * NaN ** 0 returns 1.0
765 * NaN ** B returns NaN for B <> 0
766 * 0.0 ** B returns 0.0 for B > 0
767 * 0.0 ** 0 returns 1.0
768 * 0.0 ** B returns Infinity for B < 0
769 * (-0.0) ** B returns -Infinity for B < 0 and odd(B)
770 * A ** B returns 1.0 / A ** (-B) for B < 0
771 * @return the result of the exponentiation.
772 */
fltIPow(floatType base,intType exponent)773 floatType fltIPow (floatType base, intType exponent)
774
775 {
776 #if IPOW_EXPONENTIATION_BY_SQUARING
777 uintType unsignedExponent;
778 boolType neg_exp;
779 #endif
780 floatType power;
781
782 /* fltIPow */
783 logFunction(printf("fltIPow(" FMT_E ", " FMT_D ")\n", base, exponent););
784 #if !FLOAT_NAN_COMPARISON_OKAY
785 /* This is checked first on purpose. NaN should not be equal */
786 /* to any value. E.g.: NaN == x should always return FALSE. */
787 /* Beyond that NaN should not be equal to itself also. Some */
788 /* C compilers do not compute comparisons with NaN correctly. */
789 /* As a consequence the NaN check is done first. */
790 if (unlikely(os_isnan(base))) {
791 if (unlikely(exponent == 0)) {
792 power = 1.0;
793 } else {
794 power = base;
795 } /* if */
796 } else
797 #endif
798 if (unlikely(base == 0.0)) {
799 if (exponent < 0) {
800 if (fltIsNegativeZero(base) && ((-(uintType) exponent) & 1)) {
801 power = NEGATIVE_INFINITY;
802 } else {
803 power = POSITIVE_INFINITY;
804 } /* if */
805 } else if (exponent == 0) {
806 power = 1.0;
807 } else { /* exponent > 0 */
808 if (exponent & 1) {
809 power = base; /* +0.0 respectively -0.0 are left as is. */
810 } else {
811 power = 0.0;
812 } /* if */
813 } /* if */
814 } else
815 #if IPOW_EXPONENTIATION_BY_SQUARING
816 {
817 if (exponent < 0) {
818 /* The unsigned value is negated to avoid a signed integer */
819 /* overflow if the smallest signed integer is negated. */
820 unsignedExponent = -(uintType) exponent;
821 neg_exp = TRUE;
822 } else {
823 unsignedExponent = (uintType) exponent;
824 neg_exp = FALSE;
825 } /* if */
826 if (unsignedExponent & 1) {
827 power = base;
828 } else {
829 power = 1.0;
830 } /* if */
831 unsignedExponent >>= 1;
832 while (unsignedExponent != 0) {
833 base *= base;
834 if (unsignedExponent & 1) {
835 power *= base;
836 } /* if */
837 unsignedExponent >>= 1;
838 } /* while */
839 if (neg_exp) {
840 #if CHECK_FLOAT_DIV_BY_ZERO
841 if (power == 0.0) {
842 if (fltIsNegativeZero(power)) {
843 power = NEGATIVE_INFINITY;
844 } else {
845 power = POSITIVE_INFINITY;
846 } /* if */
847 } else {
848 power = 1.0 / power;
849 } /* if */
850 #else
851 power = 1.0 / power;
852 #endif
853 } /* if */
854 }
855 #else
856 if (base < 0.0) {
857 if (exponent & 1) {
858 power = -pow(-base, (floatType) exponent);
859 } else {
860 power = pow(-base, (floatType) exponent);
861 } /* if */
862 } else { /* base > 0.0 */
863 power = pow(base, (floatType) exponent);
864 } /* if */
865 #endif
866 logFunction(printf("fltIPow --> " FMT_E "\n", power););
867 return power;
868 } /* fltIPow */
869
870
871
872 /**
873 * Determine if a number is -0.0.
874 * This function is the only possibility to determine if a number
875 * is -0.0. The comparison operators (=, <>, <, >, <=, >=) and
876 * the function 'compare' treat 0.0 and -0.0 as equal. The
877 * operators ''digits'' and ''sci'' and the function ''str''
878 * return the same [[string]] for -0.0 and +0.0.
879 * @return TRUE if the number is -0.0,
880 * FALSE otherwise.
881 */
fltIsNegativeZero(floatType number)882 boolType fltIsNegativeZero (floatType number)
883
884 {
885 boolType isNegativeZero;
886
887 /* fltIsNegativeZero */
888 isNegativeZero = memcmp(&number, &negativeZero, sizeof(floatType)) == 0;
889 logFunction(printf("fltIsNegativeZero(" FMT_E ") --> %d\n",
890 number, isNegativeZero););
891 return isNegativeZero;
892 } /* fltIsNegativeZero */
893
894
895
896 #if !LDEXP_FUNCTION_OKAY
897 /**
898 * Multiply 'number' by 2 raised to the power 'exponent'.
899 * @return number * 2 ** exponent, or
900 * NaN if 'number' is NaN.
901 */
fltLdexp(floatType number,int exponent)902 floatType fltLdexp (floatType number, int exponent)
903
904 {
905 floatType product;
906
907 /* fltLdexp */
908 logFunction(printf("fltLdexp(" FMT_E ", %d)\n", number, exponent););
909 #if !LDEXP_OF_NAN_OKAY
910 if (unlikely(os_isnan(number))) {
911 product = NOT_A_NUMBER;
912 } else
913 #endif
914 {
915 product = ldexp(number, exponent);
916 } /* if */
917 logFunction(printf("fltLdexp --> " FMT_E "\n", product););
918 return product;
919 } /* fltLdexp */
920 #endif
921
922
923
924 #if !FLOAT_COMPARISON_OKAY
925 /**
926 * Check if 'number1' is less than or equal to 'number2'.
927 * According to IEEE 754 a NaN is neither less than,
928 * equal to, nor greater than any value, including itself.
929 * If 'number1' or 'number2' is NaN, the result is FALSE.
930 * @return TRUE if 'number1' is less than or equal to 'number2',
931 * FALSE otherwise.
932 */
fltLe(floatType number1,floatType number2)933 boolType fltLe (floatType number1, floatType number2)
934
935 {
936 boolType isLessEqual;
937
938 /* fltLe */
939 logFunction(printf("fltLe(" FMT_E ", " FMT_E ")\n",
940 number1, number2););
941 #if !FLOAT_NAN_COMPARISON_OKAY
942 if (os_isnan(number1) || os_isnan(number2)) {
943 isLessEqual = FALSE;
944 } else
945 #endif
946 #if !FLOAT_ZERO_COMPARISON_OKAY
947 if (fltIsNegativeZero(number1)) {
948 if (fltIsNegativeZero(number2)) {
949 isLessEqual = TRUE;
950 } else {
951 isLessEqual = 0.0 <= number2;
952 } /* if */
953 } else if (fltIsNegativeZero(number2)) {
954 isLessEqual = number1 <= 0.0;
955 } else
956 #endif
957 {
958 isLessEqual = number1 <= number2;
959 } /* if */
960 logFunction(printf("fltLe --> %d\n", isLessEqual););
961 return isLessEqual;
962 } /* fltLe */
963 #endif
964
965
966
967 #if !LOG_FUNCTION_OKAY
968 /**
969 * Return the natural logarithm (base e) of x.
970 * log(NaN) returns NaN
971 * log(1.0) returns 0.0
972 * log(Infinity) returns Infinity
973 * log(0.0) returns -Infinity
974 * log(X) returns NaN for X < 0.0
975 * @return the natural logarithm of x.
976 */
fltLog(floatType number)977 floatType fltLog (floatType number)
978
979 {
980 floatType logarithm;
981
982 /* fltLog */
983 logFunction(printf("fltLog(" FMT_E ")\n", number););
984 #if !LOG_OF_NAN_OKAY || !FLOAT_NAN_COMPARISON_OKAY
985 /* This is checked first on purpose. NaN should not be equal */
986 /* to any value. E.g.: NaN == x should always return FALSE. */
987 /* Beyond that NaN should not be equal to itself also. Some */
988 /* C compilers do not compute comparisons with NaN correctly. */
989 /* As a consequence the NaN check is done first. */
990 if (unlikely(os_isnan(number))) {
991 logarithm = number;
992 } else
993 #endif
994 #if !LOG_OF_ZERO_OKAY
995 if (unlikely(number == 0.0)) {
996 logarithm = NEGATIVE_INFINITY;
997 } else
998 #endif
999 #if !LOG_OF_NEGATIVE_OKAY
1000 if (unlikely(number < 0.0)) {
1001 logarithm = NOT_A_NUMBER;
1002 } else
1003 #endif
1004 {
1005 logarithm = log(number);
1006 } /* if */
1007 logFunction(printf("fltLog(" FMT_E ") --> " FMT_E "\n",
1008 number, logarithm););
1009 return logarithm;
1010 } /* fltLog */
1011 #endif
1012
1013
1014
1015 #if !LOG10_FUNCTION_OKAY
1016 /**
1017 * Returns the base 10 logarithm of x.
1018 * log10(NaN) returns NaN
1019 * log10(1.0) returns 0.0
1020 * log10(Infinity) returns Infinity
1021 * log10(0.0) returns -Infinity
1022 * log10(X) returns NaN for X < 0.0
1023 * @return the base 10 logarithm of x.
1024 */
fltLog10(floatType number)1025 floatType fltLog10 (floatType number)
1026
1027 {
1028 floatType logarithm;
1029
1030 /* fltLog10 */
1031 logFunction(printf("fltLog10(" FMT_E ")\n", number););
1032 #if !LOG10_OF_NAN_OKAY || !FLOAT_NAN_COMPARISON_OKAY
1033 /* This is checked first on purpose. NaN should not be equal */
1034 /* to any value. E.g.: NaN == x should always return FALSE. */
1035 /* Beyond that NaN should not be equal to itself also. Some */
1036 /* C compilers do not compute comparisons with NaN correctly. */
1037 /* As a consequence the NaN check is done first. */
1038 if (unlikely(os_isnan(number))) {
1039 logarithm = number;
1040 } else
1041 #endif
1042 #if !LOG10_OF_ZERO_OKAY
1043 if (unlikely(number == 0.0)) {
1044 logarithm = NEGATIVE_INFINITY;
1045 } else
1046 #endif
1047 #if !LOG10_OF_NEGATIVE_OKAY
1048 if (unlikely(number < 0.0)) {
1049 logarithm = NOT_A_NUMBER;
1050 } else
1051 #endif
1052 {
1053 logarithm = log10(number);
1054 } /* if */
1055 logFunction(printf("fltLog10(" FMT_E ") --> " FMT_E "\n",
1056 number, logarithm););
1057 return logarithm;
1058 } /* fltLog10 */
1059 #endif
1060
1061
1062
1063 #if !HAS_LOG1P
1064 /**
1065 * Compute log(1.0 + x) (natural logarithm of the sum of 1 and x).
1066 * This function is used, if the function log1p() is missing from
1067 * the C math library. This simple implementation is not accurate if
1068 * the value of x is near zero.
1069 * @return log(1.0 + x)
1070 */
fltLog1p(floatType number)1071 floatType fltLog1p (floatType number)
1072
1073 {
1074 floatType logarithm;
1075
1076 /* fltLog1p */
1077 logFunction(printf("fltLog1p(" FMT_E ")\n", number););
1078 logarithm = log(1.0 + number);
1079 logFunction(printf("fltLog1p(" FMT_E ") --> " FMT_E "\n",
1080 number, logarithm););
1081 return logarithm;
1082 } /* fltLog1p */
1083 #endif
1084
1085
1086
1087 #if !LOG2_FUNCTION_OKAY
1088 /**
1089 * Returns the base 2 logarithm of x.
1090 * log2(NaN) returns NaN
1091 * log2(1.0) returns 0.0
1092 * log2(Infinity) returns Infinity
1093 * log2(0.0) returns -Infinity
1094 * log2(X) returns NaN for X < 0.0
1095 * @return the base 2 logarithm of x.
1096 */
fltLog2(floatType number)1097 floatType fltLog2 (floatType number)
1098
1099 {
1100 floatType logarithm;
1101
1102 /* fltLog2 */
1103 logFunction(printf("fltLog2(" FMT_E ")\n", number););
1104 #if HAS_LOG2
1105 #if !LOG2_OF_NAN_OKAY || !FLOAT_NAN_COMPARISON_OKAY
1106 /* This is checked first on purpose. NaN should not be equal */
1107 /* to any value. E.g.: NaN == x should always return FALSE. */
1108 /* Beyond that NaN should not be equal to itself also. Some */
1109 /* C compilers do not compute comparisons with NaN correctly. */
1110 /* As a consequence the NaN check is done first. */
1111 if (unlikely(os_isnan(number))) {
1112 logarithm = number;
1113 } else
1114 #endif
1115 #if !LOG2_OF_ZERO_OKAY
1116 if (unlikely(number == 0.0)) {
1117 logarithm = NEGATIVE_INFINITY;
1118 } else
1119 #endif
1120 #if !LOG2_OF_NEGATIVE_OKAY
1121 if (unlikely(number < 0.0)) {
1122 logarithm = NOT_A_NUMBER;
1123 } else
1124 #endif
1125 {
1126 logarithm = log2(number);
1127 } /* if */
1128 #else
1129 logarithm = fltLog(number) / LN2;
1130 #endif
1131 logFunction(printf("fltLog2(" FMT_E ") --> " FMT_E "\n",
1132 number, logarithm););
1133 return logarithm;
1134 } /* fltLog2 */
1135 #endif
1136
1137
1138
1139 #if !FLOAT_COMPARISON_OKAY
1140 /**
1141 * Check if 'number1' is less than 'number2'.
1142 * According to IEEE 754 a NaN is neither less than,
1143 * equal to, nor greater than any value, including itself.
1144 * If 'number1' or 'number2' is NaN, the result is FALSE.
1145 * @return TRUE if 'number1' is less than 'number2',
1146 * FALSE otherwise.
1147 */
fltLt(floatType number1,floatType number2)1148 boolType fltLt (floatType number1, floatType number2)
1149
1150 {
1151 boolType isLessThan;
1152
1153 /* fltLt */
1154 logFunction(printf("fltLt(" FMT_E ", " FMT_E ")\n",
1155 number1, number2););
1156 #if !FLOAT_NAN_COMPARISON_OKAY
1157 if (os_isnan(number1) || os_isnan(number2)) {
1158 isLessThan = FALSE;
1159 } else
1160 #endif
1161 #if !FLOAT_ZERO_COMPARISON_OKAY
1162 if (fltIsNegativeZero(number1)) {
1163 if (fltIsNegativeZero(number2)) {
1164 isLessThan = FALSE;
1165 } else {
1166 isLessThan = 0.0 < number2;
1167 } /* if */
1168 } else if (fltIsNegativeZero(number2)) {
1169 isLessThan = number1 < 0.0;
1170 } else
1171 #endif
1172 {
1173 isLessThan = number1 < number2;
1174 } /* if */
1175 logFunction(printf("fltLt --> %d\n", isLessThan););
1176 return isLessThan;
1177 } /* fltLt */
1178 #endif
1179
1180
1181
1182 /**
1183 * Compute the floating-point modulo of a division.
1184 * The modulo has the same sign as the divisor.
1185 * The modulo is dividend - floor(dividend / divisor) * divisor
1186 * A mod NaN returns NaN
1187 * NaN mod B returns NaN
1188 * A mod 0.0 returns NaN
1189 * Infinity mod B returns NaN
1190 * -Infinity mod B returns NaN
1191 * 0.0 mod B returns 0.0 for B <> 0.0
1192 * A mod Infinity returns A for A > 0
1193 * A mod Infinity returns Infinity for A < 0
1194 * A mod -Infinity returns A for A < 0
1195 * A mod -Infinity returns -Infinity for A > 0
1196 * @return the floating-point modulo of the division.
1197 */
fltMod(floatType dividend,floatType divisor)1198 floatType fltMod (floatType dividend, floatType divisor)
1199
1200 {
1201 floatType modulo;
1202
1203 /* fltMod */
1204 logFunction(printf("fltMod(" FMT_E ", " FMT_E ")\n", dividend, divisor););
1205 modulo = fltRem(dividend, divisor);
1206 #if FLOAT_COMPARISON_OKAY
1207 if ((dividend < 0.0) ^ (divisor < 0.0) && modulo != 0.0) {
1208 #else
1209 if (fltLt(dividend, 0.0) ^ fltLt(divisor, 0.0) && !fltEq(modulo, 0.0)) {
1210 #endif
1211 modulo += divisor;
1212 } /* if */
1213 logFunction(printf("fltMod --> " FMT_E "\n", modulo););
1214 return modulo;
1215 } /* fltMod */
1216
1217
1218
1219 /**
1220 * Convert a string to a float number.
1221 * @return the float result of the conversion.
1222 * @exception RANGE_ERROR If the string contains not a float literal.
1223 */
1224 floatType fltParse (const const_striType stri)
1225
1226 {
1227 #if USE_STRTOD
1228 char buffer[MAX_CSTRI_BUFFER_LEN + NULL_TERMINATION_LEN];
1229 const_cstriType buffer_ptr;
1230 const_cstriType cstri;
1231 char *next_ch;
1232 #else
1233 memSizeType position;
1234 floatType digitval;
1235 #endif
1236 errInfoType err_info = OKAY_NO_ERROR;
1237 floatType result;
1238
1239 /* fltParse */
1240 logFunction(printf("fltParse(\"%s\")\n", striAsUnquotedCStri(stri)););
1241 #if USE_STRTOD
1242 if (likely(stri->size <= MAX_CSTRI_BUFFER_LEN)) {
1243 cstri = NULL;
1244 buffer_ptr = conv_to_cstri(buffer, stri);
1245 if (unlikely(buffer_ptr == NULL)) {
1246 logError(printf("fltParse(\"%s\"): conv_to_cstri() failed.\n",
1247 striAsUnquotedCStri(stri)););
1248 err_info = RANGE_ERROR;
1249 result = 0.0;
1250 } /* if */
1251 } else {
1252 cstri = stri_to_cstri(stri, &err_info);
1253 buffer_ptr = cstri;
1254 } /* if */
1255 if (likely(buffer_ptr != NULL)) {
1256 if (isspace(buffer_ptr[0])) {
1257 logError(printf("fltParse(\"%s\"): String starts with whitespace.\n",
1258 striAsUnquotedCStri(stri)););
1259 err_info = RANGE_ERROR;
1260 result = 0.0;
1261 } else {
1262 result = (floatType) strtod(buffer_ptr, &next_ch);
1263 if (next_ch == buffer_ptr) {
1264 if (strcmp(buffer_ptr, "NaN") == 0) {
1265 result = NOT_A_NUMBER;
1266 } else if (strcmp(buffer_ptr, "Infinity") == 0) {
1267 result = POSITIVE_INFINITY;
1268 } else if (strcmp(buffer_ptr, "-Infinity") == 0) {
1269 result = NEGATIVE_INFINITY;
1270 } else {
1271 logError(printf("fltParse(\"%s\"): No digit or sign found.\n",
1272 striAsUnquotedCStri(stri)););
1273 err_info = RANGE_ERROR;
1274 } /* if */
1275 } else if (next_ch != &buffer_ptr[stri->size]) {
1276 logError(printf("fltParse(\"%s\"): Superfluous characters after float literal.\n",
1277 striAsUnquotedCStri(stri)););
1278 err_info = RANGE_ERROR;
1279 #if STRTOD_ACCEPTS_HEX_NUMBERS
1280 } else if (buffer_ptr[0] != '\0' &&
1281 (buffer_ptr[1] == 'x' || buffer_ptr[1] == 'X' ||
1282 (buffer_ptr[1] == '0' &&
1283 (buffer_ptr[2] == 'x' || buffer_ptr[2] == 'X')))) {
1284 logError(printf("fltParse(\"%s\"): Hex float literals are not supported.\n",
1285 striAsUnquotedCStri(stri)););
1286 err_info = RANGE_ERROR;
1287 #endif
1288 #if !STRTOD_ACCEPTS_DENORMAL_NUMBERS && ATOF_ACCEPTS_DENORMAL_NUMBERS
1289 } else if (result == 0.0 && errno == ERANGE) {
1290 result = (floatType) atof(buffer_ptr);
1291 #endif
1292 } /* if */
1293 } /* if */
1294 if (cstri != NULL) {
1295 free_cstri(cstri, stri);
1296 } /* if */
1297 } /* if */
1298 #else
1299 position = 0;
1300 result = 0.0;
1301 while (position < stri->size &&
1302 stri->mem[position] >= ((strElemType) '0') &&
1303 stri->mem[position] <= ((strElemType) '9')) {
1304 digitval = ((intType) stri->mem[position]) - ((intType) '0');
1305 if (result < MAX_DIV_10 ||
1306 (result == MAX_DIV_10 &&
1307 digitval <= MAX_REM_10)) {
1308 result = ((floatType) 10.0) * result + digitval;
1309 } else {
1310 err_info = RANGE_ERROR;
1311 } /* if */
1312 position++;
1313 } /* while */
1314 if (position == 0 || position < stri->size) {
1315 err_info = RANGE_ERROR;
1316 } /* if */
1317 #endif
1318 if (unlikely(err_info != OKAY_NO_ERROR)) {
1319 raise_error(err_info);
1320 result = 0.0;
1321 } /* if */
1322 logFunction(printf("fltParse(\"%s\") --> " FMT_E "\n",
1323 striAsUnquotedCStri(stri), result););
1324 return result;
1325 } /* fltParse */
1326
1327
1328
1329 #if !POW_FUNCTION_OKAY
1330 /**
1331 * Compute the exponentiation of a float 'base' with a float 'exponent'.
1332 * This function corrects errors of the C function pow().
1333 * A ** B returns NaN for A < 0.0 and B is not integer
1334 * A ** 0.0 returns 1.0
1335 * NaN ** 0.0 returns 1.0
1336 * NaN ** B returns NaN for B <> 0.0
1337 * 0.0 ** B returns 0.0 for B > 0.0
1338 * 0.0 ** 0.0 returns 1.0
1339 * 0.0 ** B returns Infinity for B < 0.0
1340 * (-0.0) ** B returns -Infinity for B < 0.0 and odd(B)
1341 * 1.0 ** B returns 1.0
1342 * 1.0 ** NaN returns 1.0
1343 * A ** NaN returns NaN for A <> 1.0
1344 * @return the result of the exponentiation.
1345 */
1346 floatType fltPow (floatType base, floatType exponent)
1347
1348 {
1349 floatType power;
1350 #if !POW_OF_NEGATIVE_OKAY
1351 floatType intPart;
1352 #endif
1353
1354 /* fltPow */
1355 logFunction(printf("fltPow(" FMT_E ", " FMT_E ")\n", base, exponent););
1356 #if !POW_OF_NAN_OKAY || !FLOAT_NAN_COMPARISON_OKAY
1357 /* This is checked first on purpose. NaN should not be equal */
1358 /* to any value. E.g.: NaN == x should always return FALSE. */
1359 /* Beyond that NaN should not be equal to itself also. Some */
1360 /* C compilers do not compute comparisons with NaN correctly. */
1361 /* As a consequence the NaN check is done first. */
1362 if (unlikely(os_isnan(base))) {
1363 #if FLOAT_NAN_COMPARISON_OKAY
1364 if (unlikely(exponent == 0.0)) {
1365 #else
1366 if (unlikely(!os_isnan(exponent) && exponent == 0.0)) {
1367 #endif
1368 power = 1.0;
1369 } else {
1370 power = base;
1371 } /* if */
1372 } else
1373 #endif
1374 #if !POW_OF_ZERO_OKAY
1375 if (unlikely(base == 0.0)) {
1376 if (unlikely(os_isnan(exponent))) {
1377 power = exponent;
1378 } else if (exponent < 0.0) {
1379 if (unlikely(fltIsNegativeZero(base) &&
1380 exponent >= -FLOATTYPE_TO_INT_CONVERSION_LIMIT &&
1381 floor(exponent) == exponent && /* exponent is an integer */
1382 (((uint64Type) -exponent) & 1))) { /* exponent is odd */
1383 power = NEGATIVE_INFINITY;
1384 } else {
1385 power = POSITIVE_INFINITY;
1386 } /* if */
1387 } else if (exponent == 0.0) {
1388 power = 1.0;
1389 } else { /* exponent > 0.0 */
1390 if (unlikely(exponent <= FLOATTYPE_TO_INT_CONVERSION_LIMIT &&
1391 floor(exponent) == exponent && /* exponent is an integer */
1392 (((uint64Type) exponent) & 1))) { /* exponent is odd */
1393 power = base; /* +0.0 respectively -0.0 are left as is. */
1394 } else {
1395 power = 0.0;
1396 } /* if */
1397 } /* if */
1398 } else
1399 #endif
1400 #if !POW_OF_NEGATIVE_OKAY
1401 if (unlikely(base < 0.0 && modf(exponent, &intPart) != 0.0)) {
1402 power = NOT_A_NUMBER;
1403 } else
1404 #endif
1405 #if !POW_OF_ONE_OKAY
1406 if (unlikely(base == 1.0)) {
1407 power = 1.0;
1408 } else
1409 #endif
1410 #if !POW_EXP_NAN_OKAY
1411 /* This is checked before checking for negative infinity on purpose. */
1412 if (unlikely(os_isnan(exponent))) {
1413 if (unlikely(base == 1.0)) {
1414 power = 1.0;
1415 } else {
1416 power = exponent;
1417 } /* if */
1418 } else
1419 #endif
1420 #if !POW_EXP_MINUS_INFINITY_OKAY
1421 if (unlikely(exponent == NEGATIVE_INFINITY)) {
1422 if (base < -1.0 || base > 1.0) {
1423 power = 0.0;
1424 } else if (base == 1.0) {
1425 power = 1.0;
1426 } else {
1427 power = POSITIVE_INFINITY;
1428 } /* if */
1429 } else
1430 #endif
1431 {
1432 power = pow(base, exponent);
1433 } /* if */
1434 logFunction(printf("fltPow --> " FMT_E "\n", power););
1435 return power;
1436 } /* fltPow */
1437 #endif
1438
1439
1440
1441 /**
1442 * Compute pseudo-random number in the range [low, high).
1443 * The random values are uniform distributed.
1444 * @return the computed pseudo-random number.
1445 * @exception RANGE_ERROR The range is empty (low >= high holds).
1446 */
1447 floatType fltRand (floatType low, floatType high)
1448
1449 {
1450 double delta;
1451 rtlValueUnion conv;
1452 floatType randomNumber;
1453
1454 /* fltRand */
1455 logFunction(printf("fltRand(" FMT_E ", " FMT_E ")\n", low, high););
1456 if (unlikely(low >= high)) {
1457 logError(printf("fltRand(" FMT_E ", " FMT_E "): "
1458 "The range is empty (low > high holds).\n",
1459 low, high););
1460 raise_error(RANGE_ERROR);
1461 randomNumber = 0.0;
1462 } else {
1463 delta = high - low;
1464 if (unlikely(os_isnan(delta))) {
1465 logError(printf("fltRand(" FMT_E ", " FMT_E "): "
1466 "NaN as upper or lower bound.\n",
1467 low, high););
1468 raise_error(RANGE_ERROR);
1469 randomNumber = 0.0;
1470 } else if (delta == POSITIVE_INFINITY) {
1471 do {
1472 conv.binaryValue = uintRand();
1473 randomNumber = conv.floatValue;
1474 } while (os_isnan(randomNumber) || randomNumber < low || randomNumber >= high);
1475 } else {
1476 do {
1477 conv.binaryValue = uintRandMantissa() |
1478 ((uintType) FLOATTYPE_EXPONENT_OFFSET << FLOATTYPE_MANTISSA_BITS);
1479 randomNumber = (conv.floatValue - 1.0) * delta + low;
1480 } while (randomNumber >= high);
1481 } /* if */
1482 } /* if */
1483 logFunction(printf("fltRand --> " FMT_E "\n", randomNumber););
1484 return randomNumber;
1485 } /* fltRand */
1486
1487
1488
1489 /**
1490 * Compute pseudo-random number in the range [low, high).
1491 * The random values are uniform distributed.
1492 * @return the computed pseudo-random number.
1493 */
1494 floatType fltRandNoChk (floatType low, floatType high)
1495
1496 {
1497 double delta;
1498 rtlValueUnion conv;
1499 floatType randomNumber;
1500
1501 /* fltRandNoChk */
1502 logFunction(printf("fltRandNoChk(" FMT_E ", " FMT_E ")\n", low, high););
1503 delta = high - low;
1504 do {
1505 conv.binaryValue = uintRandMantissa() |
1506 ((uintType) FLOATTYPE_EXPONENT_OFFSET << FLOATTYPE_MANTISSA_BITS);
1507 randomNumber = (conv.floatValue - 1.0) * delta + low;
1508 } while (randomNumber >= high);
1509 logFunction(printf("fltRandNoChk --> " FMT_E "\n", randomNumber););
1510 return randomNumber;
1511 } /* fltRandNoChk */
1512
1513
1514
1515 #if !FMOD_FUNCTION_OKAY
1516 /**
1517 * Compute the floating-point remainder of a division.
1518 * The remainder has the same sign as the dividend.
1519 * The remainder is dividend - flt(trunc(dividend / divisor)) * divisor
1520 * The remainder is computed without a conversion to integer.
1521 * A rem NaN returns NaN
1522 * NaN rem B returns NaN
1523 * A rem 0.0 returns NaN
1524 * Infinity rem B returns NaN
1525 * -Infinity rem B returns NaN
1526 * 0.0 rem B returns 0.0 for B <> 0.0
1527 * A rem Infinity returns A
1528 * @return the floating-point remainder of the division.
1529 */
1530 floatType fltRem (floatType dividend, floatType divisor)
1531
1532 {
1533 floatType remainder;
1534
1535 /* fltRem */
1536 logFunction(printf("fltRem(" FMT_E ", " FMT_E ")\n", dividend, divisor););
1537 #if !FMOD_DIVIDEND_NAN_OKAY
1538 if (unlikely(os_isnan(dividend))) {
1539 remainder = NOT_A_NUMBER;
1540 } else
1541 #endif
1542 #if !FMOD_DIVISOR_NAN_OKAY
1543 if (unlikely(os_isnan(divisor))) {
1544 remainder = NOT_A_NUMBER;
1545 } else
1546 #endif
1547 #if !FMOD_DIVIDEND_INFINITY_OKAY
1548 if (dividend == POSITIVE_INFINITY || dividend == NEGATIVE_INFINITY) {
1549 remainder = NOT_A_NUMBER;
1550 } else
1551 #endif
1552 #if !FMOD_DIVISOR_INFINITY_OKAY
1553 if (divisor == POSITIVE_INFINITY || divisor == NEGATIVE_INFINITY) {
1554 if (dividend == POSITIVE_INFINITY || dividend == NEGATIVE_INFINITY) {
1555 remainder = NOT_A_NUMBER;
1556 } else {
1557 remainder = dividend;
1558 } /* if */
1559 } else
1560 #endif
1561 #if !FMOD_DIVISOR_ZERO_OKAY
1562 if (divisor == 0.0) {
1563 remainder = NOT_A_NUMBER;
1564 } else
1565 #endif
1566 {
1567 remainder = fmod(dividend, divisor);
1568 } /* if */
1569 logFunction(printf("fltRem --> " FMT_E "\n", remainder););
1570 return remainder;
1571 } /* fltRem */
1572 #endif
1573
1574
1575
1576 /**
1577 * Convert a 'float' number to a [[string]] in scientific notation.
1578 * Scientific notation uses a decimal significand and a decimal exponent.
1579 * The significand has an optional sign and exactly one digit before the
1580 * decimal point. The fractional part of the significand is rounded
1581 * to the specified number of digits ('precision'). Halfway cases are
1582 * rounded away from zero. The fractional part is followed by the
1583 * letter e and an exponent, which is always signed. The value zero is
1584 * never written with a negative sign.
1585 * @param precision Number of digits after the decimal point.
1586 * If the 'precision' is zero the decimal point is omitted.
1587 * @return the string result of the conversion.
1588 * @exception RANGE_ERROR If the 'precision' is negative.
1589 * @exception MEMORY_ERROR Not enough memory to represent the result.
1590 */
1591 striType fltSci (floatType number, intType precision)
1592
1593 {
1594 char buffer_1[PRECISION_BUFFER_LEN + FLT_SCI_LEN + NULL_TERMINATION_LEN];
1595 char *buffer = buffer_1;
1596 const_cstriType buffer_ptr;
1597 /* The 4 additional chars below are for "%1.e". */
1598 char form_buffer[INTTYPE_DECIMAL_SIZE + STRLEN("%1.e") + NULL_TERMINATION_LEN];
1599 memSizeType startPos;
1600 memSizeType pos;
1601 memSizeType len;
1602 memSizeType after_zeros;
1603 striType result;
1604
1605 /* fltSci */
1606 logFunction(printf("fltSci(" FMT_E ", " FMT_D ")\n", number, precision););
1607 if (unlikely(precision < 0)) {
1608 logError(printf("fltSci(" FMT_E ", " FMT_D "): Precision negative.\n",
1609 number, precision););
1610 raise_error(RANGE_ERROR);
1611 result = NULL;
1612 } else if (unlikely(precision > PRECISION_BUFFER_LEN &&
1613 ((uintType) precision > MAX_CSTRI_LEN - FLT_SCI_LEN ||
1614 !ALLOC_CSTRI(buffer, (memSizeType) precision + FLT_SCI_LEN)))) {
1615 raise_error(MEMORY_ERROR);
1616 result = NULL;
1617 } else {
1618 if (os_isnan(number)) {
1619 buffer_ptr = "NaN";
1620 len = STRLEN("NaN");
1621 } else if (number == POSITIVE_INFINITY) {
1622 buffer_ptr = "Infinity";
1623 len = STRLEN("Infinity");
1624 } else if (number == NEGATIVE_INFINITY) {
1625 buffer_ptr = "-Infinity";
1626 len = STRLEN("-Infinity");
1627 } else {
1628 #ifdef LIMIT_FMT_E_MAXIMUM_FLOAT_PRECISION
1629 if (unlikely(precision > PRINTF_FMT_E_MAXIMUM_FLOAT_PRECISION)) {
1630 len = (memSizeType) sprintf(buffer, "%1."
1631 LIMIT_FMT_E_MAXIMUM_FLOAT_PRECISION "e", number);
1632 } else
1633 #endif
1634 #if PRINTF_SUPPORTS_VARIABLE_FORMATS
1635 if (unlikely(precision > INT_MAX)) {
1636 sprintf(form_buffer, "%%1." FMT_D "e", precision);
1637 len = (memSizeType) sprintf(buffer, form_buffer, number);
1638 } else {
1639 len = (memSizeType) sprintf(buffer, "%1.*e", (int) precision, number);
1640 } /* if */
1641 #else
1642 if (unlikely(precision >= sizeof(fmt_e) / sizeof(char *))) {
1643 sprintf(form_buffer, "%%1." FMT_D "e", precision);
1644 len = (memSizeType) sprintf(buffer, form_buffer, number);
1645 } else {
1646 len = (memSizeType) sprintf(buffer, fmt_e[precision], number);
1647 } /* if */
1648 #endif
1649 #ifdef PRINTF_FMT_E_MAXIMUM_FLOAT_PRECISION
1650 /* Some C run-time libraries do not have a fixed maximum */
1651 /* for the float precision of printf(). Instead the actual */
1652 /* precision varies with each call of printf(). Up to a */
1653 /* precision of PRINTF_FMT_E_MAXIMUM_FLOAT_PRECISION */
1654 /* printf() will always work ok. */
1655 if (precision > PRINTF_FMT_E_MAXIMUM_FLOAT_PRECISION) {
1656 pos = (memSizeType) (strchr(buffer, 'e') - buffer);
1657 if ((memSizeType) precision > pos - STRLEN("1.") -
1658 (buffer[0] == '-')){
1659 memSizeType numZeros = (memSizeType) precision -
1660 (pos - STRLEN("1.") - (buffer[0] == '-'));
1661 memmove(&buffer[pos + numZeros], &buffer[pos],
1662 len - pos + NULL_TERMINATION_LEN);
1663 memset(&buffer[pos], '0', numZeros);
1664 len += numZeros;
1665 } /* if */
1666 } /* if */
1667 #endif
1668 startPos = 0;
1669 if (buffer[0] == '-' && buffer[1] == '0') {
1670 /* All forms of -0 are converted to 0 */
1671 if (buffer[2] == '.') {
1672 pos = 3;
1673 while (buffer[pos] == '0') {
1674 pos++;
1675 } /* while */
1676 if (buffer[pos] == 'e' && buffer[pos + 2] == '0') {
1677 pos += 3;
1678 while (buffer[pos] == '0') {
1679 pos++;
1680 } /* while */
1681 if (buffer[pos] == '\0') {
1682 startPos++;
1683 len--;
1684 } /* if */
1685 } /* if */
1686 } else if (buffer[2] == 'e' && buffer[4] == '0') {
1687 pos = 5;
1688 while (buffer[pos] == '0') {
1689 pos++;
1690 } /* while */
1691 if (buffer[pos] == '\0') {
1692 startPos++;
1693 len--;
1694 } /* if */
1695 } /* if */
1696 } /* if */
1697 if (len > startPos) {
1698 pos = len;
1699 do {
1700 pos--;
1701 } while (pos > startPos && buffer[pos] != 'e');
1702 pos += 2;
1703 after_zeros = pos;
1704 while (buffer[after_zeros] == '0') {
1705 after_zeros++;
1706 } /* while */
1707 if (buffer[after_zeros] == '\0') {
1708 after_zeros--;
1709 } /* if */
1710 memmove(&buffer[pos], &buffer[after_zeros],
1711 sizeof(char) * (len - after_zeros + 1));
1712 len -= after_zeros - pos;
1713 } /* if */
1714 buffer_ptr = &buffer[startPos];
1715 } /* if */
1716 result = cstri_buf_to_stri(buffer_ptr, len);
1717 if (buffer != buffer_1) {
1718 UNALLOC_CSTRI(buffer, (memSizeType) precision + FLT_SCI_LEN);
1719 } /* if */
1720 if (unlikely(result == NULL)) {
1721 raise_error(MEMORY_ERROR);
1722 } /* if */
1723 } /* if */
1724 logFunction(printf("fltSci(" FMT_E ", " FMT_D ") --> \"%s\"\n",
1725 number, precision, striAsUnquotedCStri(result)););
1726 return result;
1727 } /* fltSci */
1728
1729
1730
1731 #if !SQRT_FUNCTION_OKAY
1732 /**
1733 * Returns the non-negative square root of x.
1734 * Used if sqrt() does not return NaN for negative numbers.
1735 * @return the square root of x.
1736 */
1737 floatType fltSqrt (floatType radicand)
1738
1739 {
1740 floatType squareRoot;
1741
1742 /* fltSqrt */
1743 logFunction(printf("fltSqrt(" FMT_E ")\n", radicand););
1744 #if !SQRT_OF_NAN_OKAY
1745 if (unlikely(os_isnan(radicand))) {
1746 squareRoot = NOT_A_NUMBER;
1747 } else
1748 #endif
1749 #if !SQRT_OF_NEGATIVE_OKAY
1750 if (radicand < 0.0) {
1751 squareRoot = NOT_A_NUMBER;
1752 } else
1753 #endif
1754 {
1755 squareRoot = sqrt(radicand);
1756 } /* if */
1757 logFunction(printf("fltSqrt(" FMT_E ") --> " FMT_E "\n",
1758 radicand, squareRoot););
1759 return squareRoot;
1760 } /* fltSqrt */
1761 #endif
1762
1763
1764
1765 /**
1766 * Convert a float number to a string.
1767 * The number is converted to a string with decimal representation.
1768 * The result string has the style [-]ddd.ddd where there is at least
1769 * one digit before and after the decimal point. The number of digits
1770 * after the decimal point is determined automatically. Except for the
1771 * case if there is only one zero digit after the decimal point,
1772 * the last digit is never zero. Negative zero (-0.0) and positive
1773 * zero (+0.0) are both converted to "0.0".
1774 * str(16.125) returns "16.125"
1775 * str(-0.0) returns "0.0"
1776 * str(Infinity) returns "Infinity"
1777 * str(-Infinity) returns "-Infinity"
1778 * str(NaN) returns "NaN"
1779 * @return the string result of the conversion.
1780 * @exception MEMORY_ERROR Not enough memory to represent the result.
1781 */
1782 striType fltStr (floatType number)
1783
1784 {
1785 char buffer[DOUBLE_TO_CHAR_BUFFER_SIZE];
1786 const_cstriType buffer_ptr;
1787 memSizeType len;
1788 striType result;
1789
1790 /* fltStr */
1791 logFunction(printf("fltStr(" FMT_E ")\n", number););
1792 if (os_isnan(number)) {
1793 buffer_ptr = "NaN";
1794 len = STRLEN("NaN");
1795 } else if (number == POSITIVE_INFINITY) {
1796 buffer_ptr = "Infinity";
1797 len = STRLEN("Infinity");
1798 } else if (number == NEGATIVE_INFINITY) {
1799 buffer_ptr = "-Infinity";
1800 len = STRLEN("-Infinity");
1801 } else {
1802 buffer_ptr = buffer;
1803 #if FLOATTYPE_DOUBLE
1804 len = doubleToCharBuffer(number, DOUBLE_STR_LARGE_NUMBER,
1805 FMT_E_DBL, buffer);
1806 #else
1807 len = doubleToCharBuffer(number, FLOAT_STR_LARGE_NUMBER,
1808 FMT_E_FLT, buffer);
1809 #endif
1810 } /* if */
1811 if (unlikely(!ALLOC_STRI_SIZE_OK(result, len))) {
1812 raise_error(MEMORY_ERROR);
1813 } else {
1814 result->size = len;
1815 memcpy_to_strelem(result->mem, (const_ustriType) buffer_ptr, len);
1816 } /* if */
1817 logFunction(printf("fltStr(" FMT_E ") --> \"%s\"\n",
1818 number, striAsUnquotedCStri(result)););
1819 return result;
1820 } /* fltStr */
1821