netinet/in.h: add missing IP protocol numbers from the linux uapi headers

support configurable page size on mips, powerpc and microblaze

PAGE_SIZE was hardcoded to 4096, which is historically what most
systems use, but on several archs it is a kernel config parameter,
user space can only know it at execution time from the aux vector.

PAGE_SIZE and PAGESIZE are not defined on archs where page size is
a runtime parameter, applications should use sysconf(_SC_PAGE_SIZE)
to query it. Internally libc code defines PAGE_SIZE to libc.page_size,
which is set to aux[AT_PAGESZ] in __init_libc and early in __dynlink
as well. (Note that libc.page_size can be accessed without GOT, ie.
before relocations are done)

Some fpathconf settings are hardcoded to 4096, these should be actually
queried from the filesystem using statfs.

fix overflow in sysconf for _SC_MQ_PRIO_MAX

the value of MQ_PRIO_MAX does not fit, so it needs to use OFLOW.

fix child stack alignment on mips clone

unlike other archs, the mips version of clone was not doing anything
to align the stack pointer. this seems to have been the cause for some
SIGBUS crashes that were observed in posix_spawn.

fix mips sysv ipc bits headers

msg.h was wrong for big-endian (wrong endiannness padding).
shm.h was just plain wrong (mips is not supposed to have padding).

both changes were tested using libc-test on qemu-system-mips.

fix x86_64 lrintl asm, again

the underlying problem was not incorrect sign extension (fixed in the
previous commit to this file by nsz) but that code that treats "long"
as 32-bit was copied blindly from i386 to x86_64.

now lrintl is identical to llrintl on x86_64, as it should be.

do not use default when dynamic linker fails to open existing path file

if fopen fails for a reason other than ENOENT, we must assume the
intent is that the path file be used. failure may be due to
misconfiguration or intentional resource-exhaustion attack (against
suid programs), in which case falling back to loading libraries from
an unintended path could be dangerous.

math: remove STRICT_ASSIGN from exp2f (see previous commit)

math: remove STRICT_ASSIGN macro

gcc did not always drop excess precision according to c99 at assignments
before version 4.5 even if -std=c99 was requested which caused badly
broken mathematical functions on i386 when FLT_EVAL_METHOD!=0

but STRICT_ASSIGN was not used consistently and it is worked around for
old compilers with -ffloat-store so it is no longer needed

the new convention is to get the compiler respect c99 semantics and when
excess precision is not harmful use float_t or double_t or to specialize
code using FLT_EVAL_METHOD

math: support invalid ld80 representations in fpclassify

apparently gnulib requires invalid long double representations
to be handled correctly in printf so we classify them according
to how the fpu treats them: bad inf is nan, bad nan is nan,
bad normal is nan and bad subnormal/zero is minimal normal

math: fix atanh (overflow and underflow issues)

in atanh exception handling was left to the called log functions,
but the argument to those functions could underflow or overflow.

use double_t and float_t to avoid some useless stores on x86

math: remove libc.h include from libm.h

libc.h is only for weak_alias so include it directly where it is used

math: fix acoshf on negative values

acosh(x) is invalid for x<1, acoshf tried to be clever using
signed comparisions to handle all x<2 the same way, but the
formula was wrong on large negative values.

math: fix expm1l on x86_64 (avoid underflow for large negative x)

copy the fix from i386: return -1 instead of exp2l(x)-1 when x <= -65

math: fix lrintl.s on x86_64 (use movslq to signextend the result)

math: fix exp2l asm on x86 (raise underflow correctly)

there were two problems:
* omitted underflow on subnormal results: exp2l(-16383.5) was calculated
as sqrt(2)*2^-16384, the last bits of sqrt(2) are zero so the down scaling
does not underflow eventhough the result is in subnormal range
* spurious underflow for subnormal inputs: exp2l(0x1p-16400) was evaluated
as f2xm1(x)+1 and f2xm1 raised underflow (because inexact subnormal result)

the first issue is fixed by raising underflow manually if x is in
(-32768,-16382] and not integer (x-0x1p63+0x1p63 != x)

the second issue is fixed by treating x in (-0x1p64,0x1p64) specially

for these fixes the special case handling was completely rewritten

math: cosmetic cleanup (use explicit union instead of fshape and dshape)

math: remove *_WORD64 macros from libm.h

only fma used these macros and the explicit union is clearer

math: remove old longdbl.h

math: long double fix (use ldshape union)

* use new ldshape union consistently
* add ld128 support to frexpl
* simplify sqrtl comment (ld64 is not just arm)

math: use float_t and double_t in scalbnf and scalbn

remove STRICT_ASSIGN (c99 semantics is assumed) and use the conventional
union to prepare the scaling factor (so libm.h is no longer needed)

math: fix remaining old long double code (erfl, fmal, lgammal, scalbnl)

in lgammal don't handle 1 and 2 specially, in fma use the new ldshape
union instead of ld80 one.

math: cbrt cleanup and long double fix

* use float_t and double_t
* cleanup subnormal handling
* bithacks according to the new convention (ldshape for long double
and explicit unions for float and double)

math: fix underflow in exp*.c and long double handling in exp2l

* don't care about inexact flag
* use double_t and float_t (faster, smaller, more precise on x86)
* exp: underflow when result is zero or subnormal and not -inf
* exp2: underflow when result is zero or subnormal and not exact
* expm1: underflow when result is zero or subnormal
* expl: don't underflow on -inf
* exp2: fix incorrect comment
* expm1: simplify special case handling and overflow properly
* expm1: cleanup final scaling and fix negative left shift ub (twopk)

math: long double trigonometric cleanup (cosl, sinl, sincosl, tanl)

ld128 support was added to internal kernel functions (__cosl, __sinl,
__tanl, __rem_pio2l) from freebsd (not tested, but should be a good
start for when ld128 arch arrives)

__rem_pio2l had some code cleanup, the freebsd ld128 code seems to
gather the results of a large reduction with precision loss (fixed
the bug but a todo comment was added for later investigation)

the old copyright was removed from the non-kernel wrapper functions
(cosl, sinl, sincosl, tanl) since these are trivial and the interesting
parts and comments had been already rewritten.

math: long double inverse trigonometric cleanup (acosl, asinl, atanl, atan2l)

* added ld128 support from freebsd fdlibm (untested)
* using new ldshape union instead of IEEEl2bits
* inexact status flag is not supported

math: rewrite hypot

method: if there is a large difference between the scale of x and y
then the larger magnitude dominates, otherwise reduce x,y so the
argument of sqrt (x*x+y*y) does not overflow or underflow and calculate
the argument precisely using exact multiplication. If the argument
has less error than 1/sqrt(2) ~ 0.7 ulp, then the result has less error
than 1 ulp in nearest rounding mode.

the original fdlibm method was the same, except it used bit hacks
instead of dekker-veltkamp algorithm, which is problematic for long
double where different representations are supported. (the new hypot
and hypotl code should be smaller and faster on 32bit cpu archs with
fast fpu), the new code behaves differently in non-nearest rounding,
but the error should be still less than 2ulps.

ld80 and ld128 are supported

math: rewrite remainder functions (remainder, remquo, fmod, modf)

* results are exact
* modfl follows truncl (raises inexact flag spuriously now)
* modf and modff only had cosmetic cleanup
* remainder is just a wrapper around remquo now
* using iterative shift+subtract for remquo and fmod
* ld80 and ld128 are supported as well

math: rewrite rounding functions (ceil, floor, trunc, round, rint)

* faster, smaller, cleaner implementation than the bit hacks of fdlibm
* use arithmetics like y=(double)(x+0x1p52)-0x1p52, which is an integer
neighbor of x in all rounding modes (0<=x<0x1p52) and only use bithacks
when that's faster and smaller (for float it usually is)
* the code assumes standard excess precision handling for casts
* long double code supports both ld80 and ld128
* nearbyint is not changed (it is a wrapper around rint)

math: fix logb(-0.0) in downward rounding mode

use -1/(x*x) instead of -1/(x+0) to return -inf, -0+0 is -0 in
downward rounding mode

math: ilogb cleanup

* consistent code style
* explicit union instead of typedef for double and float bit access
* turn FENV_ACCESS ON to make 0/0.0f raise invalid flag
* (untested) ld128 version of ilogbl (used by logbl which has ld128 support)

long double cleanup, initial commit

new ldshape union, ld128 support is kept, code that used the old
ldshape union was rewritten (IEEEl2bits union of freebsd libm is
not touched yet)

ld80 __fpclassifyl no longer tries to handle invalid representation

fix typo in comment in __randname

fix mips-specific bug in synccall (too little space for signal mask)

switch to the new __block_all_sigs/__restore_sigs internal API to
clean up the code too.

in synccall, ignore the signal before any threads' signal handlers return

this protects against deadlock from spurious signals (e.g. sent by
another process) arriving after the controlling thread releases the
other threads from the sync operation.

fix invalid pointer in synccall (multithread setuid, etc.)

the head pointer was not being reset between calls to synccall, so any
use of this interface more than once would build the linked list
incorrectly, keeping the (now invalid) list nodes from the previous
call.

fix special-case breakage in popen due to reversed argument order

add workaround header for wait.h remapping to sys/wait.h

fix missing return value warning in faccessat, minor cleanup

clone will pass the return value of the start function to SYS_exit
anyway; there's no need to call the syscall directly.

fix invalid %m format crash in wide scanf variants

the wide variant was missed in the previous commit.

avoid crash in scanf when invalid %m format is encountered

invalid format strings invoke undefined behavior, so this is not a
conformance issue, but it's nicer for scanf to report the error safely
instead of calling free on a potentially-uninitialized pointer or a
pointer to memory belonging to the caller.

remove incorrect cancellation points from realpath

debloat realpath's allocation strategy

rather than allocating a PATH_MAX-sized buffer when the caller does
not provide an output buffer, work first with a PATH_MAX-sized temp
buffer with automatic storage, and either copy it to the caller's
buffer or strdup it on success. this not only avoids massive memory
waste, but also avoids pulling in free (and thus the full malloc
implementation) unnecessarily in static programs.

make realpath use O_PATH when opening the file

this avoids failure if the file is not readable and avoids odd
behavior for device nodes, etc. on old kernels that lack O_PATH, the
old behavior (O_RDONLY) will naturally happen as the fallback.

fix regression in creation of ldso symlink

DESTDIR was wrongly included in the symlink contents.

fix breakage in synccall due to incorrect signal restoration in sigqueue

commit 07827d1a82fb33262f686eda959857f0d28cd8fa seems to have
introduced this issue. sigqueue is called from the synccall core, at
which time, even implementation-internal signals are blocked. however,
pthread_sigmask removes the implementation-internal signals from the
old mask before returning, so that a process which began life with
them blocked will not be able to save a signal mask that has them
blocked, possibly causing them to become re-blocked later. however,
this was causing sigqueue to unblock the implementation-internal
signals during synccall, leading to deadlock.

fix typo in release notes

release notes for 0.9.13

only expose struct tcphdr under _GNU_SOURCE

the BSD and GNU versions of this structure differ, so exposing it in
the default _BSD_SOURCE profile is possibly problematic. both versions
could be simultaneously supported with anonymous unions if needed in
the future, but for now, just omitting it except under _GNU_SOURCE
should be safe.

add struct tcphdr in netinet/tcp.h

remove -Wcast-align from --enable-warnings

I originally added this warning option based on a misunderstanding of
how it works. it does not warn whenever the destination of the cast
has stricter alignment; it only warns in cases where misaligned
dereference could lead to a fault. thus, it's essentially a no-op for
i386, which had me wrongly believing the code was clean for this
warning level. on other archs, numerous diagnostic messages are
produced, and all of them are false-positives, so it's better just not
to use it.

optimized C memcpy

unlike the old C memcpy, this version handles word-at-a-time reads and
writes even for misaligned copies. it does not require that the cpu
support misaligned accesses; instead, it performs bit shifts to
realign the bytes for the destination.

essentially, this is the C version of the ARM assembly language
memcpy. the ideas are all the same, and it should perform well on any
arch with a decent number of general-purpose registers that has a
barrel shift operation. since the barrel shifter is an optional cpu
feature on microblaze, it may be desirable to provide an alternate asm
implementation on microblaze, but otherwise the C code provides a
competitive implementation for "generic risc-y" cpu archs that should
alleviate the urgent need for arch-specific memcpy asm.

stdbool.h should define __bool_true_false_are_defined even for C++

while the incorporation of this requirement from C99 into C++11 was
likely an accident, some software expects it to be defined, and it
doesn't hurt. if the requirement is removed, then presumably
__bool_true_false_are_defined would just be in the implementation
namespace and thus defining it would still be legal.

fix invalid instruction mnemonics in powerpc fenv asm

there is no non-dot version of the andis instruction, but there's no
harm in updating the flags anyway, so just use the dot version.

optimized C memset

this version of memset is optimized both for small and large values of
n, and makes no misaligned writes, so it is usable (and near-optimal)
on all archs. it is capable of filling up to 52 or 56 bytes without
entering a loop and with at most 7 branches, all of which can be fully
predicted if memset is called multiple times with the same size.

it also uses the attribute extension to inform the compiler that it is
violating the aliasing rules, unlike the previous code which simply
assumed it was safe to violate the aliasing rules since translation
unit boundaries hide the violations from the compiler. for non-GNUC
compilers, 100% portable fallback code in the form of a naive loop is
provided. I intend to eventually apply this approach to all of the
string/memory functions which are doing word-at-a-time accesses.

add attribute((may_alias)) checking in configure

this will be needed for upcoming commits to the string/mem functions
to correct their unannounced use of aliasing violations for
word-at-a-time search, fill, and copy operations.

add the %s (seconds since the epoch) format to strftime

this is a nonstandard extension but will be required in the next
version of POSIX, and it's widely used/useful in shell scripts
utilizing the date utility.

fix strftime regression in %e format

%e pads with spaces instead of zeros.

properly fill in tzname[] for old (pre-64-bit-format) zoneinfo files

in this case, the first standard-time and first daylight-time rules
should be taken as the "default" ones to expose.

minor fix to tz name checking

if a zoneinfo file is not (or is no longer) in use, don't check the
abbrevs pointers, which may be invalid.

fix strftime handling of time zone data

this may need further revision in the future, since POSIX is rather
unclear on the requirements, and is designed around the assumption of
POSIX TZ specifiers which are not sufficiently powerful to represent
real-world timezones (this is why zoneinfo support was added).

the basic issue is that strftime gets the string and numeric offset
for the timezone from the extra fields in struct tm, which are
initialized when calling localtime/gmtime/etc. however, a conforming
application might have created its own struct tm without initializing
these fields, in which case using __tm_zone (a pointer) could crash.
other zoneinfo-based implementations simply check for a null pointer,
but otherwise can still crash of the field contains junk.

simply ignoring __tm_zone and using tzname[] would "work" but would
give incorrect results in time zones with more complex rules. I feel
like this would lower the quality of implementation.

instead, simply validate __tm_zone: unless it points to one of the
zone name strings managed by the timezone system, assume it's invalid.

this commit also fixes several other minor bugs with formatting:
tm_isdst being negative is required to suppress printing of the zone
formats, and %z was using the wrong format specifiers since the type
of val was changed, resulting in bogus output.

make dlopen honor the rpath of the main program

this seems to match what other systems do, and seems useful for
programs that have their libraries and plugins stored relative to the
executable.

fix mishandling of empty or blank TZ environment variable

the empty TZ string was matching equal to the initial value of the
cached TZ name, thus causing do_tzset never to run and never to
initialize the time zone data.

fix regression in dn_expand/reverse dns

off-by-one error copying the name components was yielding junk at the
beginning and truncating one character at the end (of every
component).

fix bugs in $ORIGIN handling

1. an occurrence of ${ORIGIN} before $ORIGIN would be ignored due to
the strstr logic. (note that rpath contains multiple :-delimited paths
to be searched.)

2. data read by readlink was not null-terminated.

use AT_EXECFN, if available, for dynamic linker to identify main program

fallback to argv[0] as before. unlike argv[0], AT_EXECFN was a valid
(but possibly relative) pathname for the new program image at the time
the execve syscall was made.

as a special case, ignore AT_EXECFN if it begins with "/proc/", in
order not to give bogus (and possibly harmful) results when fexecve
was used.

add rpath $ORIGIN processing to dynamic linker

add recursive rpath support to dynamic linker

previously, rpath was only honored for direct dependencies. in other
words, if A depends on B and B depends on C, only B's rpath (if any),
not A's rpath, was being searched for C. this limitation made
rpath-based deployment difficult in the presence of multiple levels of
library dependency.

at present, $ORIGIN processing in rpath is still unsupported.

fix missing string.h in strftime.c (needed by new strftime code)

this bug was masked by local experimental CFLAGS in my config.mak.

fix some documentation typos

add strftime and wcsftime field widths

at present, since POSIX requires %F to behave as %+4Y-%m-%d and ISO C
requires %F to behave as %Y-%m-%d, the default behavior for %Y has
been changed to match %+4Y. this seems to be the only way to conform
to the requirements of both standards, and it does not affect years
prior to the year 10000. depending on the outcome of interpretations
from the standards bodies, this may be adjusted at some point.

simplify strftime and fix integer overflows

use a long long value so that even with offsets, values cannot
overflow. instead of using different format strings for different
numeric formats, simply use a per-format width and %0*lld for all of
them.

this width specifier is not for use with strftime field widths; that
will be a separate step in the caller.

strftime cleanup: avoid recomputing strlen when it's known

more strftime refactoring

make __strftime_fmt_1 return a string (possibly in the caller-provided
temp buffer) rather than writing into the output buffer. this approach
makes more sense when padding to a minimum field width might be
required, and it's also closer to what wcsftime wants.

begin refactoring strftime to make adding field widths easier

add SUN_LEN macro to sys/un.h under appropriate feature tests

this is ugly and useless, but it seems to be the least-ugly way to
provide it...

unbreak vwarn: print ": " before errno message

patch by Strake. this seems to be a regression caused by fixing the
behavior of perror("") to match perror(0) at some point in the past.

fix two bugs in sed code configure uses to save command line

one place where semicolon (non-portable) was still used in place of
separate -e options (copied over from an old version of this code),
and use of a literal slash in the bracket expression for the final
command, despite slash being used as the delimiter for the s command.

re-add logic for ignoring failure of ld.so symlink installation

this was inadvertently removed when switching to the new install.sh.

fix fenv exception functions to mask their argument

fesetround.c is a wrapper to do the arch independent argument
check (on archs where rounding mode is not stored in 2 bits
__fesetround still has to check its arguments)

on powerpc fe*except functions do not accept the extra invalid
flags of its fpscr register

the useless FENV_ACCESS pragma was removed from feupdateenv

optimize x86 feclearexcept: only use save/restore x87 fenv if needed

the x87 exception summary (ES) and stack fault (SF) flags may be
spuriously cleared by feclearexcept using the fnclex instruction,
but these flags are not observable through libc hence maintaining
their state is not critical.

remove the __mxcsr member from fenv_t on i386 to follow the glibc abi

in the previous commit sse fenv support was added, but there is no
need to save mxcsr (sse fenv register) so fix the abi incompatibility
with glibc.

add sse fenv support on i386 through hwcap

the sse and x87 rounding modes should be always the same,
the visible exception flags are the bitwise or of the two
fenv states (so it's enough to query the rounding mode or
raise exceptions on one fenv)

fix i386 fesetenv: FE_DFL_ENV is (fenv_t*)-1 not 0

remove spurious tmp file present since initial git check-in

replace system's install command with a shell script

the historical (non-standardized) install command is really
inappropriate for installing binaries/libraries on a system that
utilizes memory-mapped executable files. rather than replacing an
existing file atomically, it overwrites the existing file. this can
cause running programs to see a partially-modified version of the
file, resulting in unpredictable behavior, or SIGBUS. a MAP_COPY mode
for mmap would get around this problem, but Linux lacks MAP_COPY.

the shell script added with this commit works around the problem by
writing temporary files and moving them into place. unlike the
historical install utility, it also support a -l option for installing
a symbolic link atomically, via the same method.

add hkscs/big5-2003/eten extensions to iconv big5

with these changes, the character set implemented as "big5" in musl is
a pure superset of cp950, the canonical "big5", and agrees with the
normative parts of Unicode. this means it has minor differences from
both hkscs and big5-2003:

- the range A2CC-A2CE maps to CJK ideographs rather than numerals,
  contrary to changes made in big5-2003.

- C6CD maps to a CJK ideograph rather than its corresponding Kangxi
  radical character, contrary to changes made in hkscs.

- F9FE maps to U+2593 rather than U+FFED.

of these differences, none but the last are visually distinct, and the
last is a character used purely for text-based graphics, not to convey
linguistic content.

should there be future demand for strict conformance to big5-2003 or
hkscs mappings, the present charset aliases can be replaced with
distinct variants.

reportedly there are other non-standard big5 extensions in common use
in Taiwan and perhaps elsewhere, which could also be added as layers
on top of the existing big5 support.

there may be additional characters which should be added to the hkscs
table: the whatwg standard for big5 defines what appears to be a
superset of hkscs.

make configure store its command line in config.mak for easy re-run

proper shell quoting and pretty-printing (avoiding ugly gratuitous
quoting and bad quoting style) is included.

fix atomicity and other issues installing dynamic linker symlink

ln -sf is non-atomic; it unlinks the destination first. instead, make
a temporary link and rename it into place.

this commit also fixes some of the dependency tracking behavior for
the link. depending on the directory it's to be installed in is not
reasonable; it causes a new link to be attempted if the library
directory has been modified, but does not attempt to make a new link
just because libc has been updated. instead, depend on the target to
be linked to. this will ensure that, if prefix has changed but
syslibdir has not, the link will be updated to point to the new
prefix.

some initial math asm for armhf (fabs[f] and sqrt[f])

fix detection of arm hardfloat

it turns out that __SOFTFP__ does not indicate the ABI in use but
rather that fpu instructions are not to be used at all. this is
specified in ARM's documentation so I'm unclear on how I previously
got the wrong idea. unfortunately, this resulted in the 0.9.12 release
producing a dynamic linker with the wrong name. fortunately, there do
not yet seem to be any public toolchain builds using the wrong name.

the __ARM_PCS_VFP macro does not seem to be official from ARM, and in
fact it was missing from the very earliest gcc versions (around 4.5.x)
that added -mfloat-abi=hard. it would be possible on such versions to
perform some ugly linker-based tests instead in hopes that the linker
will reject ABI-mismatching object files, if there is demand for
supporting such versions. I would probably prefer to document which
versions are broken and warn users to manually add -D__ARM_PCS_VFP if
using such a version.

there's definitely an argument to be made that the fenv macros should
be exposed even in -mfloat-abi=softfp mode. for now, I have chosen not
to expose them in this case, since the math library will not
necessarily have the capability to raise exceptions (it depends on the
CFLAGS used to compile it), and since exceptions are officially
excluded from the ARM EABI, which the plain "arm" arch aims to
follow.

support floating point environment (fenv) on armhf (hard float) subarchs

patch by nsz. I've tested it on an armhf machine and it seems to be
working correctly.

fix build of x86_64 expl assembly

apparently this label change was not carried over when adapting the
changes from the i386 version.

add function types to arm crt assembly

without these, calls may be resolved incorrectly if the calling code
has been compiled to thumb instead of arm. it's not clear to me at
this point whether crt_arch.h is even working if crt1.c is built as
thumb; this needs testing. but the _init and _fini issues were known
to cause crashes in static-linked apps when libc was built as thumb,
and this commit should fix that issue.

math: fix pow(x,-1) to raise underflow properly

if FLT_EVAL_METHOD!=0 check if (double)(1/x) is subnormal and not a
power of 2 (if 1/x is power of 2 then either it is exact or the
long double to double rounding already raised inexact and underflow)

math: fix i386 atan2.s to raise underflow for subnormal results

math: clean up atan2.c

* remove volatile hacks
* don't care about inexact flag for now (removed all the +-tiny)
* fix atanl to raise underflow properly
* remove signed int arithmetics
* use pi/2 instead of pi_o_2 (gcc generates the same code, which is not
correct, but it does not matter: we mainly care about nearest rounding)

math: fix x86 asin, atan, exp, log1p to raise underflow

underflow is raised by an inexact subnormal float store,
since subnormal operations are slow, check the underflow
flag and skip the store if it's already raised

math: fix x86 expl.s to raise underflow and clean up special case handling

math: fix asin, atan, log1p, tanh to raise underflow on subnormal

for these functions f(x)=x for small inputs, because f(0)=0 and
f'(0)=1, but for subnormal values they should raise the underflow
flag (required by annex F), if they are approximated by a polynomial
around 0 then spurious underflow should be avoided (not required by
annex F)

all these functions should raise inexact flag for small x if x!=0,
but it's not required by the standard and it does not seem a worthy
goal, so support for it is removed in some cases.

raising underflow:
- x*x may not raise underflow for subnormal x if FLT_EVAL_METHOD!=0
- x*x may raise spurious underflow for normal x if FLT_EVAL_METHOD==0
- in case of double subnormal x, store x as float
- in case of float subnormal x, store x*x as float