SHARE
TWEET

Untitled

a guest Apr 21st, 2019 123 Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1. public static float linearCombination(float a1_, float b1_, float a2_, float b2_)
  2. {
  3.     double a1 = (double)a1_;
  4.     double a2 = (double)a2_;
  5.     double b1 = (double)b1_;
  6.     double b2 = (double)b2_;
  7.  
  8.     // the code below is split in many additions/subtractions that may
  9.     // appear redundant. However, they should NOT be simplified, as they
  10.     // use IEEE754 floating point arithmetic rounding properties.
  11.     // The variable naming conventions are that xyzHigh contains the most significant
  12.     // bits of xyz and xyzLow contains its least significant bits. So theoretically
  13.     // xyz is the sum xyzHigh + xyzLow, but in many cases below, this sum cannot
  14.     // be represented in only one double precision number so we preserve two numbers
  15.     // to hold it as long as we can, combining the high and low order bits together
  16.     // only at the end, after cancellation may have occurred on high order bits
  17.  
  18.     // split a1 and b1 as one 26 bits number and one 27 bits number
  19.     double a1High = System.BitConverter.Int64BitsToDouble(System.BitConverter.DoubleToInt64Bits(a1) & ((-1L) << 27));
  20.     double a1Low = a1 - a1High;
  21.     double b1High = System.BitConverter.Int64BitsToDouble(System.BitConverter.DoubleToInt64Bits(b1) & ((-1L) << 27));
  22.     double b1Low = b1 - b1High;
  23.  
  24.     // accurate multiplication a1 * b1
  25.     double prod1High = a1 * b1;
  26.     double prod1Low = a1Low * b1Low - (((prod1High - a1High * b1High) - a1Low * b1High) - a1High * b1Low);
  27.  
  28.     // split a2 and b2 as one 26 bits number and one 27 bits number
  29.     double a2High = System.BitConverter.Int64BitsToDouble(System.BitConverter.DoubleToInt64Bits(a2) & ((-1L) << 27));
  30.     double a2Low = a2 - a2High;
  31.     double b2High = System.BitConverter.Int64BitsToDouble(System.BitConverter.DoubleToInt64Bits(b2) & ((-1L) << 27));
  32.     double b2Low = b2 - b2High;
  33.  
  34.     // accurate multiplication a2 * b2
  35.     double prod2High = a2 * b2;
  36.     double prod2Low = a2Low * b2Low - (((prod2High - a2High * b2High) - a2Low * b2High) - a2High * b2Low);
  37.  
  38.     // accurate addition a1 * b1 + a2 * b2
  39.     double s12High = prod1High + prod2High;
  40.     double s12Prime = s12High - prod2High;
  41.     double s12Low = (prod2High - (s12High - s12Prime)) + (prod1High - s12Prime);
  42.  
  43.     // final rounding, s12 may have suffered many cancellations, we try
  44.     // to recover some bits from the extra words we have saved up to now
  45.     double result = s12High + (prod1Low + prod2Low + s12Low);
  46.  
  47.     if (double.IsNaN(result))
  48.     {
  49.         // either we have split infinite numbers or some coefficients were NaNs,
  50.         // just rely on the naive implementation and let IEEE754 handle this
  51.         result = a1 * b1 + a2 * b2;
  52.     }
  53.  
  54.     return (float)result;
  55. }
RAW Paste Data
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy. OK, I Understand
 
Top