public static float linearCombination(float a1_, float b1_, float a2_, float b2_)
{
    double a1 = (double)a1_;
    double a2 = (double)a2_;
    double b1 = (double)b1_;
    double b2 = (double)b2_;

    // the code below is split in many additions/subtractions that may
    // appear redundant. However, they should NOT be simplified, as they
    // use IEEE754 floating point arithmetic rounding properties.
    // The variable naming conventions are that xyzHigh contains the most significant
    // bits of xyz and xyzLow contains its least significant bits. So theoretically
    // xyz is the sum xyzHigh + xyzLow, but in many cases below, this sum cannot
    // be represented in only one double precision number so we preserve two numbers
    // to hold it as long as we can, combining the high and low order bits together
    // only at the end, after cancellation may have occurred on high order bits

    // split a1 and b1 as one 26 bits number and one 27 bits number
    double a1High = System.BitConverter.Int64BitsToDouble(System.BitConverter.DoubleToInt64Bits(a1) & ((-1L) << 27));
    double a1Low = a1 - a1High;
    double b1High = System.BitConverter.Int64BitsToDouble(System.BitConverter.DoubleToInt64Bits(b1) & ((-1L) << 27));
    double b1Low = b1 - b1High;

    // accurate multiplication a1 * b1
    double prod1High = a1 * b1;
    double prod1Low = a1Low * b1Low - (((prod1High - a1High * b1High) - a1Low * b1High) - a1High * b1Low);

    // split a2 and b2 as one 26 bits number and one 27 bits number
    double a2High = System.BitConverter.Int64BitsToDouble(System.BitConverter.DoubleToInt64Bits(a2) & ((-1L) << 27));
    double a2Low = a2 - a2High;
    double b2High = System.BitConverter.Int64BitsToDouble(System.BitConverter.DoubleToInt64Bits(b2) & ((-1L) << 27));
    double b2Low = b2 - b2High;

    // accurate multiplication a2 * b2
    double prod2High = a2 * b2;
    double prod2Low = a2Low * b2Low - (((prod2High - a2High * b2High) - a2Low * b2High) - a2High * b2Low);

    // accurate addition a1 * b1 + a2 * b2
    double s12High = prod1High + prod2High;
    double s12Prime = s12High - prod2High;
    double s12Low = (prod2High - (s12High - s12Prime)) + (prod1High - s12Prime);

    // final rounding, s12 may have suffered many cancellations, we try
    // to recover some bits from the extra words we have saved up to now
    double result = s12High + (prod1Low + prod2Low + s12Low);

    if (double.IsNaN(result))
    {
        // either we have split infinite numbers or some coefficients were NaNs,
        // just rely on the naive implementation and let IEEE754 handle this
        result = a1 * b1 + a2 * b2;
    }

    return (float)result;
}