Bug in egcc optimization (f point ops give wrong answers)

[Enrique Zanardi <ezanardi at ull dot es>]

Hello!

I have been trying to fix a bug in g77, and have found it is caused by bad
behaviour in egcc optimizations.

My test suite are the following three files:

 ... testing.c, a tiny program that computes log(1+j), where j is the
 ... imaginary unit

#include <stdio.h>
typedef struct { float r, i; } complex;
extern void c_log(complex*, complex*);
void main(void) {
  complex x,temp;
  temp.i = 1.0e0;
  temp.r = 1.0e0;
  c_log(&x,&temp);
  printf(" (%11.9f,%11.9f)\n",x.r,x.i);
}

 ... c_log.c, adapted from egcs-1.0.2 sources
 ... (egcs-1.0.2/gcc/f/runtime/libF77/c_log.c)
 ... c_log(&zout,&zin) computes zout=log(zin), where zin and zout are complex

#include <math.h>
typedef struct { float r, i; } complex;
extern double f__cabs(double, double);
void c_log(complex *resx, complex *z) {
  complex res;
  res.i = atan2(z->i, z->r);
  res.r = log( f__cabs(z->r, z->i) );
  resx->r = res.r;
  resx->i = res.i;
}

 ... cabs.c, adapted from egcs-1.0.2 sources 
 ... (egcs-1.0.2/gcc/f/runtime/libF77/cabs.c)
 ... f__cabs(x,y) returns sqrt(x*x+y*y)

#include <math.h>
double f__cabs(double real, double imag) {
  double temp;
  if(real < 0)
	real = -real;
  if(imag < 0)
	imag = -imag;
  if(imag > real){
	temp = real;
	real = imag;
	imag = temp;
  }
  if((real+imag) == real)
	return(real);
  temp = imag/real;
  temp = real*sqrt(1.0 + temp*temp);  /*overflow!!*/
  return(temp);
}

 ...

When compiled with -O, the program computes the right answer
(log[sqrt(2)],PI/2). 
$ egcc -Wall -O -o testing testing.c c_log.c cabs.c
testing.c:14: warning: return type of `main' is not `int'
$ ./testing 
 (0.346573591,0.785398185)

But when compiled with -O2, the program gives a wrong value for the real
part. 
$ egcc -Wall -O2 -o testing testing.c c_log.c cabs.c    
testing.c:14: warning: return type of `main' is not `int'
$ ./testing 
 (0.240239993,0.785398185)

By looking at the disassembled program one can find the error is in
c_log:

0x8048190 <c_log>:	pushl  %ebp
0x8048191 <c_log+1>:	movl   %esp,%ebp
0x8048193 <c_log+3>:	subl   $0xc,%esp
0x8048196 <c_log+6>:	pushl  %ebx
0x8048197 <c_log+7>:	movl   0xc(%ebp),%eax	# eax=&z;
0x804819a <c_log+10>:	movl   0x8(%ebp),%ebx  	# ebx=&resx;
0x804819d <c_log+13>:	flds   (%eax)		# push z.r on the FPU stack
0x804819f <c_log+15>:	flds   0x4(%eax)	# push z.i on the FPU stack
0x80481a2 <c_log+18>:	fld    %st(1)		# F-push z.r again 
0x80481a4 <c_log+20>:	fpatan 			# F-pop z.r and z.i, compute
						# atan(z.i/z.r) and F-push
						# the result (let's call
						# it 'a')
0x80481a6 <c_log+22>:	fld    %st(0)		# F-push 'a' again
0x80481a8 <c_log+24>:	fstps  0xfffffff4(%ebp)	# F-pop 'a', put it on ebp-12
0x80481ab <c_log+27>:	flds   0xfffffff4(%ebp)	# F-push 'a' again
0x80481ae <c_log+30>:	fstps  0xfffffffc(%ebp)	# F-pop 'a', put it on ebp-4
0x80481b1 <c_log+33>:	subl   $0x8,%esp	# leave room for f__cabs
						# second argument
0x80481b4 <c_log+36>:	fstpl  (%esp,1)		#*** F-pop 'a' and use it as
						#*** f__cabs second argument!
0x80481b7 <c_log+39>:	subl   $0x8,%esp	# leave room for f__cabs
						# first  argument
0x80481ba <c_log+42>:	fstpl  (%esp,1)		# F-pop z.r and use it as
						# f__cabs first argument
0x80481bd <c_log+45>:	call   0x80481e0 <f__cabs>
[...]

Indeed, a=PI/2, and log[sqrt(1+a*a)]=0.240239991, the wrong answer we saw
before.

Also, that bug is related to force-mem optimization. If I disable that
optimization explicitly the bug doesn't appear:

$ egcc -Wall -O2 -fno-force-mem -o testing testing.c c_log.c cabs.c
testing.c:4: warning: return type of `main' is not `int'
$ ./testing 
 (0.346573591,0.785398185)

OS: Debian GNU/Linux 2.0 (with glibc 2.0.7pre1-4)
egcc version: egcs-2.90.28 980423 (egcs-1.0.3 prerelease)

I hope that helps.

	Thanks,
-- 
Enrique Zanardi					ezanardi@ull.es
						ezanard@debian.org
-- 
Enrique Zanardi						ezanardi@ull.es