At the moment I am seeing floating point exceptions FE_INEXACT errors
on 12.1-RELEASE-p6 amd64 and FreeBSD 13.0-CURRENT r362257 amd64 as 
well as FreeBSD 13.0-CURRENT r362037 ppc64.

The IEEE-754 2008 FP64 data format should allow for a single bit to be
set anywhere within the significand portion in order to represent a
perfect power of 2 which is added to the implied leading bit 1. Thus a
value of 1.5 is 1 + ( 2 ^ ( -1 ) ) which in binary is 1.1 exactly. The
exponent has an offset of 1023 and in this case the exponent will need
to be zero with a sign of positive represented by a zero. The offset for
the eleven bit exponent will be 1023 and thus the decimal 1.5 in FP64
format will be : 

        0   01111111111  1 0000 ... 0000 

Here I isolate the sign bit as the leading zero and then we see the
value 1023 as the exponent in the next eleven bits. The significand has
the single bit set in the uppermost position which is added to the 
now invisible implied leading one bit.  This is evaluated thus : 

      sign +1   *    2^( exponent 0 )  *   ( binary 1.1 )

wherein we get the decimal 1.5 precisely.

This would then imply that the decimal base ten value 1.5 should be
perfectly reasonable to convert to binary FP64 format by strtod() and
there will be no floating point exception FE_INEXACT nor any other of
the possible exception conditions. 

The same logic may be extended to any other perfect power of 2 added
to the value one wherein we use the possible bit positions from the left
most upper bit of the significand, that is to say the most significant
bit down to the unit of least precision in bit number 63 if we count 
with the sign bit at position zero.  Therefore we may use strtod() to
convert the value one plus two raised to the power negative fifty with
no loss of precision.  The data for 1 + 2 ^ ( -50 ) therefore looks like
so : 

        0   01111111111  1 0000 ... 0100 = 0x3f f0 00 00 00 00 00 04

This will evaluate to the precise decimal value : 

        1.000000000000000888178419700125232338905334472656250

The strtod() call should accept this string and be able to perfectly
encode the FP64 data without a floating point exception of any manner.

I therefore present the following test code and demonstrate results :


/*********************************************************************
 * The Open Group Base Specifications Issue 6
 * IEEE Std 1003.1, 2004 Edition
 *    An XSI-conforming application should ensure that the feature
 *    test macro _XOPEN_SOURCE is defined with the value 600 before
 *    inclusion of any header. This is needed to enable the
 *    functionality described in The _POSIX_C_SOURCE Feature Test
 *    Macro and in addition to enable the XSI extension.
 *********************************************************************/
#define _XOPEN_SOURCE 600

#include <errno.h>
#include <locale.h>
#include <math.h>
#include <stdint.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <sys/resource.h>
#include <sys/utsname.h>
#include <time.h>
#include <unistd.h>

#include <fenv.h>
#pragma STDC FENV_ACCESS ON

#define __STDC_FORMAT_MACROS
#include <inttypes.h>

/* Accept ascii data on the command line and attempt to convert
 * to IEEE 754-2008 floating point FP64 data.
 * 
 * FreeBSD 12 RELEASE triggers a floating point exception 
 * FE_INEXACT on some data that should be accepted perfectly. */

uint64_t system_memory(void);
int sysinfo(void);

int main ( int argc, char *argv[] )
{

    double candidate_double, num;
    int fpe_raised = 0;

    if ( argc < 2 ) {
        fprintf(stderr,"FAIL : provide a decimal number\n");
        return ( EXIT_FAILURE );
    }

    char *buf = malloc((size_t)32);
    if ( buf == NULL ) {
        perror ("malloc!");
        fprintf (stderr,"FAIL : malloc failed for buf\n");
        return(EXIT_FAILURE);
    }

    if ( argc > 3 ) {
        printf ("\nINFO : You suggest a locale of %s\n", argv[2]);
        buf = setlocale ( LC_ALL, argv[2] );
    } else {
        buf = setlocale ( LC_ALL, "POSIX" );
    }

    if ( buf == NULL ) {
        fprintf (stderr,"FAIL : setlocale fail\n");
        return(EXIT_FAILURE);
    }

    sysinfo();

    errno = 0;
    feclearexcept(FE_ALL_EXCEPT);

    candidate_double = strtod(argv[1], (char **)NULL);

    fpe_raised = fetestexcept(FE_ALL_EXCEPT);
    if (fpe_raised!=0){
        printf("INFO : FP Exception raised is");
        if ( fpe_raised & FE_INEXACT ) printf(" FE_INEXACT");
        if ( fpe_raised & FE_DIVBYZERO ) printf(" FE_DIVBYZERO");
        if ( fpe_raised & FE_UNDERFLOW ) printf(" FE_UNDERFLOW");
        if ( fpe_raised & FE_OVERFLOW ) printf(" FE_OVERFLOW");
        if ( fpe_raised & FE_INVALID ) printf(" FE_INVALID");
        printf("\n");
    }

    if ( fpe_raised & FE_INEXACT ) {
        printf("WARN : FE_INEXACT returned by strtod()\n");
    }

    if ( ( errno == ERANGE ) || ( errno == EINVAL ) ){
        fprintf(stderr,"FAIL : number not understood\n");
        perror("     ");
        return ( EXIT_FAILURE );
    }

    if ( !isnormal(candidate_double) && ( candidate_double != 0.0 ) ) {
        fprintf(stderr,"FAIL : number is not normal\n");
        fprintf(stderr,"     : looks like %-+22.16e\n", candidate_double);
        return ( EXIT_FAILURE );
    }

    feclearexcept(FE_ALL_EXCEPT);

    num = candidate_double;
    /* slightly wide format spec to see many digits which should
     * be well past the FP64 precision */
    printf ("INFO : seems like a decimal number %-+36.34g\n", num);
    return ( EXIT_SUCCESS );

}

uint64_t system_memory()
{
    /* should return the amount of memory available in bytes */
    long en;
    uint64_t pages, page_size;

    en = sysconf(_SC_PHYS_PAGES);
    if ( en < 0 ){
        perror("sysconf(_SC_PHYS_PAGES) : ");
        exit(EXIT_FAILURE);
    }
    pages = (uint64_t) en;

    page_size = (uint64_t)sysconf(_SC_PAGE_SIZE);
    return ( pages * page_size );
}

int sysinfo(void) {

    struct utsname uname_data;
    uint64_t sysmem = system_memory();
    uint64_t pagesize = (uint64_t)sysconf(_SC_PAGESIZE);
    int fp_round_mode;

    /* guess the architecture endianess? */
    int end_check = 1;
    int little_endian = (*(uint8_t*)&end_check == 1) ? 1 : 0;

    setlocale( LC_MESSAGES, "POSIX" );
    if ( uname( &uname_data ) < 0 ) {
        fprintf ( stderr,
                 "WARNING : Could not attain system uname data.\n" );
        perror ( "uname" );
    } else {
        printf ( "-------------------------------" );
        printf ( "------------------------------\n" );
        printf ( "           system name = %s\n", uname_data.sysname );
        printf ( "             node name = %s\n", uname_data.nodename );
        printf ( "               release = %s\n", uname_data.release );
        printf ( "               version = %s\n", uname_data.version );
        printf ( "               machine = %s\n", uname_data.machine );
        printf ( "             page size = %" PRIu64 "\n", pagesize );
        printf ( "          avail memory = %" PRIu64 "\n", sysmem );
        printf ( "                       = %" PRIu64 " kB\n", sysmem/1024 );
        printf ( "                       = %" PRIu64 " MB\n", sysmem/1048576 );
        /*  does not really work for memory size near GB boundaries
         *  if ( sysmem > ( 1024 * 1048576 ) )
         *      printf ( "                    = %" PRIu64 " GB\n",
         *              sysmem/( 1024 * 1048576 ) );
         */

        printf ( "                endian = ");
        if ( little_endian ) {
            printf ( "little");
        } else {
            printf ( "big");
        }
        printf ( " endian\n" );

        printf ( " sizeof(unsigned long) = %lu\n", sizeof(unsigned long) );
        printf ( "           sizeof(int) = %lu\n", sizeof(int) );
        printf ( "         sizeof(void*) = %lu\n", sizeof(void*) );

        fp_round_mode = fegetround();
        printf("     fp rounding mode is ");
        switch(fp_round_mode){
            case FE_TONEAREST:
                printf("FE_TONEAREST\n");
                break;
            case FE_TOWARDZERO:
                printf("FE_TOWARDZERO\n");
                break;
            case FE_UPWARD:
                printf("FE_UPWARD\n");
                break;
            case FE_DOWNWARD:
                printf("FE_DOWNWARD\n");
                break;
            default:
                printf("unknown!\n");
                break;
        }

        printf ( "-------------------------------" );
        printf ( "------------------------------" );
    }
    printf ("\n");

    return ( EXIT_SUCCESS );

}


If we compile and run the above on an Oracle Solaris UNIX server with 
the Fujitsu SPARC64-VII+ processors and use the Oracle Studio C99 tools
then we see : 

alpha$ 
alpha$ uname -a 
SunOS alpha 5.10 Generic_150400-65 sun4u sparc SUNW,SPARC-Enterprise
alpha$ echo '60k 1 2 _50 ^ p + pq' | dc 
.000000000000000888178419700125232338905334472656250000000000
1.000000000000000888178419700125232338905334472656250000000000
alpha$ ./fp64_inexact 1.000000000000000888178419700125232338905334472656250
-------------------------------------------------------------
           system name = SunOS
             node name = alpha
               release = 5.10
               version = Generic_150400-65
               machine = sun4u
             page size = 8192
          avail memory = 17179869184
                       = 16777216 kB
                       = 16384 MB
                endian = big endian
 sizeof(unsigned long) = 8
           sizeof(int) = 4
         sizeof(void*) = 8
     fp rounding mode is FE_TONEAREST
-------------------------------------------------------------
INFO : seems like a decimal number +1.000000000000000888178419700125232
alpha$ 

Note that there is no floating point exception.

The same code on Debian Linux with IBM Power9 ppc64le reveals : 

C_titan$ 
C_titan$ uname -a 
Linux titan 5.5.0-1-powerpc64le #1 SMP Debian 5.5.13-2 (2020-03-30) ppc64le GNU/Linux
C_titan$ echo '60k 1 2 _50 ^ p + pq' | dc 
.000000000000000888178419700125232338905334472656250000000000
1.000000000000000888178419700125232338905334472656250000000000
C_titan$ ./fp64_inexact 1.000000000000000888178419700125232338905334472656250
-------------------------------------------------------------
           system name = Linux
             node name = titan
               release = 5.5.0-1-powerpc64le
               version = #1 SMP Debian 5.5.13-2 (2020-03-30)
               machine = ppc64le
             page size = 65536
          avail memory = 4255776768
                       = 4156032 kB
                       = 4058 MB
                endian = little endian
 sizeof(unsigned long) = 8
           sizeof(int) = 4
         sizeof(void*) = 8
     fp rounding mode is FE_TONEAREST
-------------------------------------------------------------
INFO : seems like a decimal number +1.000000000000000888178419700125232
C_titan$ 


Again we see perfect conversion with no exceptions. 

Here is the same test with FreeBSD 13.0 CURRENT on amd64 :

inferno$ uname -apKU
FreeBSD inferno 13.0-CURRENT FreeBSD 13.0-CURRENT #1 r362257: Wed Jun 17 17:24:45 UTC 2020     root@inferno:/opt/obj/usr/src/amd64.amd64/sys/GENERIC  amd64 amd64 1300098 1300097
inferno$ echo '60k 1 2 _50 ^ p + pq' | dc 
.000000000000000888178419700125232338905334472656250000000000
1.000000000000000888178419700125232338905334472656250000000000
inferno$ ./fp64_inexact 1.000000000000000888178419700125232338905334472656250
-------------------------------------------------------------
           system name = FreeBSD
             node name = inferno
               release = 13.0-CURRENT
               version = FreeBSD 13.0-CURRENT #1 r362257: Wed Jun 17 17:24:45 UTC 2020     root@inferno:/opt/obj/usr/src/amd64.amd64/sys/GENERIC
               machine = amd64
             page size = 4096
          avail memory = 17140547584
                       = 16738816 kB
                       = 16346 MB
                endian = little endian
 sizeof(unsigned long) = 8
           sizeof(int) = 4
         sizeof(void*) = 8
     fp rounding mode is FE_TONEAREST
-------------------------------------------------------------
INFO : seems like a decimal number +1.000000000000000888178419700125232
inferno$ 

This seems perfect.

However if we shift one bit we see : 

inferno$ 
inferno$ echo '60k 1 2 _49 ^ p + pq' | dc
.000000000000001776356839400250464677810668945312500000000000
1.000000000000001776356839400250464677810668945312500000000000
inferno$ ./fp64_inexact 1.00000000000000177635683940025046467781066894531250
-------------------------------------------------------------
           system name = FreeBSD
             node name = inferno
               release = 13.0-CURRENT
               version = FreeBSD 13.0-CURRENT #1 r362257: Wed Jun 17 17:24:45 UTC 2020     root@inferno:/opt/obj/usr/src/amd64.amd64/sys/GENERIC
               machine = amd64
             page size = 4096
          avail memory = 17140547584
                       = 16738816 kB
                       = 16346 MB
                endian = little endian
 sizeof(unsigned long) = 8
           sizeof(int) = 4
         sizeof(void*) = 8
     fp rounding mode is FE_TONEAREST
-------------------------------------------------------------
INFO : FP Exception raised is FE_INEXACT
WARN : FE_INEXACT returned by strtod()
INFO : seems like a decimal number +1.000000000000001776356839400250465
inferno$ 

This should not happen. 

The same test on various other systems and architectures operates with
no exceptions. 

I therefore wrote a script to test the situation on various machines and
various operating systems.  See tarball attached.  The strtod() code on
FreeBSD signals exception FE_INEXACT across a wide range of perfectly
reasonable input data.  Other systems accept the data flawlessly. Thus I
suggest that there may be a bug in the David M. Gay code for strtod.c as
well as possibly other places. I have not yet begun to single step into
this code but felt I should do the research first and confirm data with
multiple systems and compilers. 

Please see tarball attached of various systems and results as well as a
trivial script which will run on just about any reasonable machine.

Dennis M. Clarke
dclarke@blastwave.org