Language designers do not like to add hair to their language to control precision.
The original C even specified that all floating arguments would be **double** and there was no standard way to call for a single precision square root.

Language designers came to feel obligated to define standard libraries thus to provide practical portability between platforms. It was indeed very convenient that the spelling of logarithm and square root routines became standardized across architectures. The desire for interoperability led many implementers to sacrifice performance for precision in mathematical routines. Those who made these tradeoffs were far removed from those in a position to know the relative benefits of precision over speed. Floating point performance of the hardware on standard arithmetic operations, was given high priority at an initially high cost in transistors, at least with the advent of IBM’s Power architecture. These different tradeoffs were not coordinated. Hardware technology was making such rapid strides that signals from the market back to those making the tradeoffs were never timely. I think that fast floating point on modern microprocessors is indeed a great utility but it has come about for its sex appeal, rather than its utility.

The above situation leads to math routines that in order to provide all of the precision of their argument carry intermediate values in some sort of special extra precision hardware. Adding such hardware support would not solve the problem for then the language designer would feel (properly) obligated to provide access to this greater precision and define math routines that supported this greater precision and then these routines would again be forced to do extra precision arithmetic, without help of hardware, for intermediate results. It seems we are stuck in a sub-optimal point regarding implementation of standard math routines. (See my fast square root for the special case of the PowerPC and sqrt.)

The HP calculator HP-35 has this problem internally. It provides 10 digit external precision but uses 12 digit arithmetic internally to provide 10 digit accuracy (most of the time) using simple calculations for transcendental functions. Quoting Why transcendentals and arbitrary precision?: “with comparable care in coding, correctly-rounded transcendental functions cost 2-4X more than conventional functions that aim for perhaps 0.53 ulps worst error.”.

Another sort of floating point dilemma arises in the compilation of an expression such as the harmonic mean:

1/(1/a + 1/b + 1/c + 1/d).

Except for issues of division by zero and exponent range this is the same as

(a*b)*(c*d)/((a+b)*(c*d) + (a*b)*(c+d)).

The second form is vastly superior for it replaces 4 divides with 5 multiplies when repeated subexpressions are noted.
Furthermore the divides are likely to be serial for lack of hardware resources while the multiplys are both much faster and likely to be overlapped.
Compilers don’t optimize thus for the floating exponent field of intermediate results may overflow unless
the inputs are in the middle quarter of the exponent range which, incidentally, they typically are.
Even when the programmer knows this there is no way in current languages to declare this to the compiler.
(Perhaps Fortran 90 provides such semantics!)
In practice the compiler would have to advise the programmer of the benefits of such a declaration.
Note that when the hardware exponent range does not suffice for the problem there is a run time warning available.
This is available in both schemes; it strikes for a narrower range of exponents in the optimized cases.
Regarding division by 0: in most cases for continuity it is desirable to define the harmonic mean of a set of numbers that include 0, as 0.
The optimized expression is usually superior in such cases.
When the input exponents are very small it is important to warn of exponent underflow.
This is a rare exception in my experience to the general rule of silent generation of zero for underflow.
An ancient precursor to this construct is Heron’s formula for the area of a triangle: √(s(s−a)(s−b)(s−c)).

There are a variety of similar issues handling fixed point numbers. The details differ but the solutions are probably related. Here are some suggestions.

The only principled solution that I see is to augment languages to support declarations of precision and range. I could declare that a variable or subexpression provided only 50 bits of precision and thus its square root could be done much more quickly. Alas I cannot argue that this extra language complexity is worth the cost. It will be several decades yet before the computing art becomes mature enough to warrant attending to this issue. Those employing FPGAs or building ASICs may address these issues first.