如果没有 FMA，您可以使用双产品算法

只要指数不溢出或下溢，两个积就会产生精确的乘法，形式为两个浮点数的和：

a x b = round(a x b) + residue

如果余数为正，则应向上舍入，那么正确答案是 round(axb) 之后朝正无穷方向的下一个浮点数。

对于除法，我们也可以通过乘法得到余数：

round(a/b) x b = round( round(a/b) x b) + residue_1

我们寻找的残留物有点不同：

(a/b) x b = round(a/b) x b + residue

因此，感兴趣的残基为：

residue = a - round(a/b) x b
residue = (a - round( round(a/b) x b)) - residue_1

类似地，对于 sqrt：

round(sqrt(a)) x round(sqrt(a)) = round(round(sqrt(a)) x round(sqrt(a))) + residue_1

a = round(sqrt(a)) x round(sqrt(a)) + residue

residue = (a - round(round(sqrt(a)) x round(sqrt(a)))) - residue_1

由于我们感兴趣的是余数的符号而不是真值，因此我们可以对操作数应用 2^n 比例，以避免上溢/下溢情况

这是使用二和二积来查找残差舍入误差的解决方案的草图。

该解决方案需要使用标准 nextafter() 函数的某些实现，可以从多个地方找到，例如 https://github.com/scijs/nextafter/blob/master/README.md

请原谅我的 javascript 技能可能非常近似，但您会明白的。

const TO_NEAREST = 0;
const TO_NEGATIVE_INFINITY: = 1;
const TO_POSITIVE_INFINITY: = 2;
const TO_ZERO: = 3;
const ROUNDING_DIRECTION=[Number.NaN , Number.Number.NEGATIVE_INFINITY , Number.POSITIVE_INFINITY , 0.0];

function split(a)
// helper function for two product, split a number intwo parts a=x+y
// each part having no more than 26 consecutive bits set
{
    let splitter=1+Math.pow(2,26); // splitter of two-product algorithm for double precision
    let c = factor*a; // BEWARE: might overflow for big a, but this case has been excluded a priori
    let x = c - (c-a);
    let y = a - x;
    return [x,y];
}

function mul_residue(a,b,c)
// answer the residue of exact operation a*b-c
// using the two product algorithm
// could be expressed as fma(a,b,-c)
// limitations: might fail if operations overflow/underflow
{
    let [a1,a2] = split(a);
    let [b1,b2] = split(b);
    let r1 = a*b - (a1*b1);
    let r2 = r1 - (a1*b2);
    let r3 = r2 - (a2*b1);
    let res1 = a2*b2 - r3;
    let res2 = a*b - c;
    return res1 + res2;
}

function mul_residue_scaled(a,b,c)
// answer the residue of exact operation a*b-c
// after applying some scaling to avoid overflow/underflow
{
    let [a_scaled , expa] = frexp(a);
    let [b_scaled , expb] = frexp(b);
    let c_scaled = ldexp(c , -expa-expb );
    return mul_residue(a_scaled,b_scaled,c_scaled)
}

function sum_residue(a,b,c)
// answer the residue of operation c=a+b using the two sum algorithm
{
    let a1=c-b;
    let b1=c-a;
    let delta_a=a-a1;
    let delta_b=b-b1;
    let res=delta_a+delta_b;
    return res;
}

function round_toward(c,res,mode)
// round the result c to upper or lower, depending on the residue and rounding mode
{
    let direction=ROUNDING_DIRECTION[ mode ];
    if(res > 0.0 && c < direction) return nextafter(c,direction);
    if(res < 0.0 && c > direction) return nextafter(c,direction);
    return c;
}

function round_sum(a,b,mode)
{
    let c=a+b;
    // directed rounding might avoid case of overflow
    if ( (! isFinite(c)) && isFinite(a) && isFinite(b) && (b != 0)) return round_toward(c,-c,mode);
    res=sum_residue(a,b,c);
    return round_toward(c,res,mode);
}

function round_sub(a,b,mode)
{
    return round_sum(a,-b,mode);
}

function round_mul(a,b,mode)
{
    let c=a*b;
    if ( (! isFinite(c)) && isFinite(a) && isFinite(b) && (b != 0)) return round_toward(c,-c,mode);
    let res=mul_residue_scaled(a,b,c);
    return round_toward(c,res,mode);
}

function round_div(a,b,mode)
{
    let c=a/b;
    if ( (! isFinite(c)) && isFinite(a) && isFinite(b) && (b != 0)) return round_toward(c,-c,mode);
    let res= mul_residue_scaled(c,b,a);
    return round_toward(c,-res * Math.sign(b),mode);
}

function round_sqrt(a,mode)
{
    let c=Math.sqrt(a);
    let res=mul_residue_scaled(c,c,a);
    return round_toward(c,-res,mode);
}

这是模拟 frexp 和 ldexp 的候选实现。
请注意，乘以 2 的幂是一个精确运算，除非结果下溢，此 ldexp 模拟会处理这种情况。

function exponent_bits(x)
// extract the 11 bits of float representation encoding the exponent part
{
    let float = new Float64Array(1),
    let bytes = new Uint8Array(float.buffer);
    float[0] = x;
    return ((bytes[7] & 0x7f) << 4 | bytes[6] >> 4);
}

function exponent(x)
// answer the exponent of x
// where non zero finite x can be represented as
//    sign * significand * 2^exponent
// with significand in the interval (0.5,1(
// undefined behaviour for special values infinity and NaN
{
    if(x === 0) return 0;
    let e = exponent_bits(x);
    if(e === 0) return exponent( x*Math.pow(2.0,53) ) - 53;
    return e - 1022;
}

function frexp(x)
// emulate stdlib frexp function
// answer the number scaled by a power of two so that it's magnitude fall in interval (0.5,1(
// and the associated power
{
    if( ! isFinite(x)) return [x,0];
    if( x == 0 ) return [x,0];
    let e = exponent(x);
    return [ldexp(x,-e) , e];
}

function ldexp(x,n)
// emulate stdlib ldexp function
// scale a number by a power of two : x * 2^n
{
    if( ! isFinite(x)) return x;
    if( x == 0 ) return x;
    if( n > 1023 ) return ldexp(ldexp(x,1023),n-1023);
    if( n >= -1074) return x * Math.pow(2.0,n);
    // care to avoid inexact result in case of underflow
    let delta_to_underflow = Math.max( exponent(x) - 1022 , -1074 );
    // handle case of total cancellation
    if (delta_to_underflow >= 0) return ldexp(ldexp(x,-1022),n+1022);
    ldexp(ldexp(x,delta_to_underflow),n-delta_to_underflow);
}