For now, solve the precision problem with brute force,

by using long double where available.
Unfortunately, it won't be available everywhere, so a better solution
would still be nice.
Also, sometimes rounding of smaller sizes doesn't work right yet.
This commit is contained in:
Olaf Seibert 2021-01-18 23:11:00 +01:00
parent 37abe35427
commit 8aea1498c5
3 changed files with 91 additions and 4 deletions

41
parse.c
View File

@ -3,6 +3,7 @@
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <float.h>
#include <ctype.h>
#include "parse.h" /* my own definitions */
@ -318,6 +319,38 @@ int get_mode(
return TRUE;
}
/*
* We need 56 bits of mantissa.
*
* Try to detect if it is needed, possible and useful to use
* long double instead of double, when parsing floating point numbers.
*/
#if DBL_MANT_DIG >= 56
/* plain double seems big enough */
# define USE_LONG_DOUBLE 0
/* long double exists and seems big enough */
#elif LDBL_MANT_DIG >= 56
# define USE_LONG_DOUBLE 1
#elif defined(LDBL_MANT_DIG)
/* long double exists but is probably still too small */
# define USE_LONG_DOUBLE 1
#else
/* long double does not exist and plain double is too small */
# define USE_LONG_DOUBLE 0
#endif
#if USE_LONG_DOUBLE
# define DOUBLE long double
# define SCANF_FMT "%Lf"
# define FREXP frexpl
#else
# define DOUBLE double
# define SCANF_FMT "%lf"
# define FREXP frexp
#endif
/* Parse PDP-11 64-bit floating point format. */
/* Give a pointer to "size" words to receive the result. */
/* Note: there are probably degenerate cases that store incorrect
@ -334,8 +367,8 @@ int parse_float(
int size,
unsigned *flt)
{
double d; /* value */
double frac; /* fractional value */
DOUBLE d; /* value */
DOUBLE frac; /* fractional value */
ulong64 ufrac; /* fraction converted to 49 bit
unsigned integer */
int i; /* Number of fields converted by sscanf */
@ -344,7 +377,7 @@ int parse_float(
unsigned uexp; /* Unsigned excess-128 exponent */
unsigned sign = 0; /* Sign mask */
i = sscanf(cp, "%lf%n", &d, &n);
i = sscanf(cp, SCANF_FMT "%n", &d, &n);
if (i == 0)
return 0; /* Wasn't able to convert */
@ -359,7 +392,7 @@ int parse_float(
return 1; /* Good job. */
}
frac = frexp(d, &sexp); /* Separate into exponent and mantissa */
frac = FREXP(d, &sexp); /* Separate into exponent and mantissa */
if (sexp < -127 || sexp > 127)
return 0; /* Exponent out of range. */

50
tests/test-float.lst.ok Normal file
View File

@ -0,0 +1,50 @@
1 ;;;;;
2 ;
3 ; Test floating point numbers
4
5 000000 040200 .word ^F 1.0 ; 040200
6 000002 140200 .word ^F-1.0 ; 140200
7 000004 137600 .word -^F 1.0 ; 137600
8 000006 037600 .word -^F-1.0 ; 037600
9
10 000010 040706 .word ^F6.2 ; 040706
11 000012 137071 .word ^C^F6.2 ; 137071
12 000014 137071 .word ^C<^F6.2> ; 137071
13
14 000016 040706 063146 .flt2 6.2 ; 040706 063146
15 000022 040706 063146 063146 .flt4 6.2 ; 040706 063146 063146 063146
000030 063146
16
17 000032 040300 000000 .flt2 1.5 ; 040300 000000
18 000036 040300 000000 000000 .flt4 1.5 ; 040300 000000 000000 000000
000044 000000
19
20 000046 056200 .word ^F 72057594037927935 ; 056200 (rounded!)
21 000050 056200 000000 .flt2 72057594037927935 ; 056200 000000 (rounded!)
22 000054 056177 177777 177777 .flt4 72057594037927935 ; 056177 177777 177777 177777 (exact!)
000062 177777
23
24 000064 056200 .word ^F 72057594037927936 ; 056200
25 000066 056200 000000 .flt2 72057594037927936 ; 056200 000000
26 000072 056200 000000 000000 .flt4 72057594037927936 ; 056200 000000 000000 000000
000100 000000
27
28 000102 056200 000000 000000 .flt4 72057594037927937 ; 056200 000000 000000 000001
000110 000000
29
30 000112 056400 000000 000000 .flt4 144115188075855873 ; 1 << 57 +1
000120 000000
31 ; 056400 000000 000000 000000 (rounded!)
32
32
Symbol table
. ******R 001
Program sections:
. ABS. 000000 000 (RW,I,GBL,ABS,OVR,NOSAV)
000122 001 (RW,I,LCL,REL,CON,NOSAV)

View File

@ -26,3 +26,7 @@
.flt4 72057594037927936 ; 056200 000000 000000 000000
.flt4 72057594037927937 ; 056200 000000 000000 000001
.flt4 144115188075855873 ; 1 << 57 +1
; 056400 000000 000000 000000 (rounded!)