i386/133583: fma does not respect rounding mode using extended precision

[Abramo Bagnara]

>Number:         133583
>Category:       i386
>Synopsis:       fma does not respect rounding mode using extended precision
>Confidential:   no
>Severity:       serious
>Priority:       medium
>Responsible:    freebsd-i386
>State:          open
>Quarter:        
>Keywords:       
>Date-Required:
>Class:          sw-bug
>Submitter-Id:   current-users
>Arrival-Date:   Fri Apr 10 22:40:01 UTC 2009
>Closed-Date:
>Last-Modified:
>Originator:     Abramo Bagnara
>Release:        7.1
>Organization:
>Environment:
FreeBSD freebsd.homenet.telecomitalia.it 7.1-RELEASE FreeBSD 7.1-RELEASE #0: Thu Jan  1 14:37:25 UTC 2009     root at logan.cse.buffalo.edu:/usr/obj/usr/src/sys/GENERIC  i386
>Description:
After fpsetprec(FP_PE) libc fma does not respect rounding mode.
>How-To-Repeat:
$ gcc bug.c -lm
$ ./a.out
Exact result: 9973859.79298831405913006165064871311187744140625
Default precision, DOWNWARD  +* : 9973859.79298831336200237274169921875
Default precision, UPWARD    +* : 9973859.79298831522464752197265625
Default precision, DOWNWARD  fma: 9973859.79298831336200237274169921875
Default precision, UPWARD    fma: 9973859.79298831522464752197265625
Extended precision, DOWNWARD  +* : 9973859.79298831336200237274169921875
Extended precision, UPWARD    +* : 9973859.79298831522464752197265625
Extended precision, DOWNWARD  fma: 9973859.7929883114993572235107421875
Extended precision, UPWARD    fma: 9973859.79298831336200237274169921875
>Fix:


Patch attached with submission follows:

#include <math.h>
#include <fenv.h>
#include <stdio.h>

#ifdef __FreeBSD__
#include <ieeefp.h>
#endif

double to = 8248384;
double x = 2871;
double y = 601.00166944908187360852025449275970458984375;


// Exact result is: 9973859.79298831405913006165064871311187744140625
int main() {
  printf("Exact result: 9973859.79298831405913006165064871311187744140625\n");

  fesetround(FE_DOWNWARD);
  printf("Default precision, DOWNWARD  +* : %.1000g\n", to + x * y);
  fesetround(FE_UPWARD);
  printf("Default precision, UPWARD    +* : %.1000g\n", to + x * y);

  fesetround(FE_DOWNWARD);
  printf("Default precision, DOWNWARD  fma: %.1000g\n", fma(x, y, to));
  fesetround(FE_UPWARD);
  printf("Default precision, UPWARD    fma: %.1000g\n", fma(x, y, to));

#ifdef __FreeBSD__
  fpsetprec(FP_PE);

  fesetround(FE_DOWNWARD);
  printf("Extended precision, DOWNWARD  +* : %.1000g\n", to + x * y);
  fesetround(FE_UPWARD);
  printf("Extended precision, UPWARD    +* : %.1000g\n", to + x * y);

  fesetround(FE_DOWNWARD);
  printf("Extended precision, DOWNWARD  fma: %.1000g\n", fma(x, y, to));
  fesetround(FE_UPWARD);
  printf("Extended precision, UPWARD    fma: %.1000g\n", fma(x, y, to));
#endif
}


>Release-Note:
>Audit-Trail:
>Unformatted: