fork download
copy 
  1. #include <float.h>
  2. #include <math.h>
  3. #include <stdlib.h>
  4. #include <stdio.h>
  5.  
  6. /* returns the correctly rounded to float result of 
  7.    the mathematical expression a*b+c. 
  8.    Assumes round-to-nearest, a≥0, b≥0, c≥0, all finite. */
  9. float myfma(float a, float b, float c) 
  10. {
  11.   double p1 = a * (double) b; 
  12.   double p2 = c;
  13.   double s1 = p1 + p2;
  14.   if (p1 < p2)
  15.     {
  16.       double tmp = p1;
  17.       p1 = p2;
  18.       p2 = tmp;
  19.     }
  20.   double r1 = s1 - p1 - p2; /* fma = s1 + r1 */
  21.  
  22.   float f1 = s1;            
  23.  
  24.   if (r1 == 0) return f1;
  25.  
  26.   double t = s1 - f1;
  27.  
  28.   if (t == 0) return f1; /* definitely not a halfway point */
  29.  
  30.   double dir = copysign(1.0 / 0.0, t);
  31.   float f2 = nextafterf(f1, dir);
  32.   double ulp = f2 - (double) f1;
  33.  
  34. #if 0
  35.   printf("f1:%a\nf2:%a\np1:%a\np2:%a\ns1:%a\nt :%a\nr1:%a\nul:%a\n", 
  36.           f1,    f2,    p1,    p2,    s1,    t,     r1,    ulp);
  37. #endif
  38.  
  39.   if (fabs(t) > fabs(0.75 * ulp))
  40.     return f2;
  41.  
  42.   if (fabs(t) < fabs(0.25 * ulp))
  43.     return f1;
  44.  
  45.   double r = t - 0.5 * ulp + r1;
  46.  
  47.   if (r > 0 && t > 0 || r < 0 && t < 0)
  48.     return f2;
  49.  
  50.   if (r > 0 && t < 0 || r < 0 && t > 0)
  51.     return f1;
  52.  
  53.   /* Round to even. */
  54.   return f1 + 0.5 * ulp;
  55. }
  56.  
  57. int main()
  58. {
  59.   while (1)
  60.     {
  61.       int a = (rand() & 0xFFFFF) << (rand() & 7);
  62.       int b = (rand() & 0xFFFFF) << (rand() & 7);
  63.       long long c = (long long)(rand() & 0xFFFFF) << (rand() & 31);
  64.  
  65.       float truefma = (long long) a * b + c;
  66.       float r = myfma(a, b, c);
  67.       if (r != truefma)
  68.         {
  69.           printf("bad: %d %d %lld\ntr:%a\nmy:%a\n\n", a, b, c, truefma, r);
  70.         }
  71.     }
  72. }
  73.  
Time limit exceeded #stdin #stdout 5s 1784KB
comments (?)
stdin
copy 
Standard input is empty
stdout
copy 
Standard output is empty
https://ideone.com/kx7MXE
language:
C (gcc 8.3)
created:
10 years ago
visibility:
public

Share or Embed source code

Discover > Sphere Engine API

The brand new service which powers Ideone!

Discover > IDE Widget

Widget for compiling and running the source code in a web browser!