FYI: Pentium FPATAN can lose 36 bits (out of 64)

[Zhishun Alex Liu]

/*
 * On a Pentium with the FDIV bug, atanl loses about 36 bits (out of a total
 * of 64) in this case
 *
 * x                 0x3FFE EEEEE008 2B40C010 (+9.333324451392881479853e-01)
 * Correct  atanl(x) 0x3FFE C03CDB1D D839EA5B (+7.509285877097010891072e-01)
 * Computed atanl(x) 0x3FFE C03CDB26 60BFEF2E (+7.509285896965136483617e-01)
 *
 * On a 486 (with the same binary):
 *
 * x                 0x3FFE EEEEE008 2B40C010 (+9.333324451392881479853e-01)
 * Correct  atanl(x) 0x3FFE C03CDB1D D839EA5B (+7.509285877097010891072e-01)
 * Computed atanl(x) 0x3FFE C03CDB1D D839EA5B (+7.509285877097010891072e-01)
 *
 * A little more investigation by masking off trailing bits of the above x
 * one by one shows that the shortest bit pattern whose atanl still loses
 * about 36 bits is
 * x                 0x3FFE EEEEE008 2A000000 (+9.333324451381486142054e-01)
 *
 * The following 29-bit number whose atanl loses only 34 bits
 * x                 0x3FFE EEEEE008 00000000 (+9.333324450999498367310e-01)
 *
 * BTW, our implementation of atanl(x) is the trivial 4-liner:
 *
 *	fldt	4(%esp)			/ push x
 *	fld1				/ push 1.0L
 *	fpatan				/ atanl(x / 1.0L)
 *	ret
 */
#include <stdio.h>
#include <stdlib.h>

#define	xdump(str, x)	printf("%s 0x%04.4X %08.8X %08.8X (%+.21Le)\n", str,\
	((int *)&x)[2] & 0xffff, ((int *)&x)[1], ((int *)&x)[0], x)
/*
 * BeEF hit this case with parameters 1024 and 25000 - in less than 10 CPU
 * hours on our Solaris 2.4 FDIV-challenged Pentium box.
 *
 * 9.3331146e-01  -0.547   9.3341501832660104762e-01 0
 * 9.3342590e-01   4e+10   9.3333244513928814799e-01
 */
extern long double atanl(long double);

void
main(void) {
	int i;
	long double x = 9.3333244513928814798525701590392600338e-01L, y;
	long double ref = 7.5092858770970108910720527961579762177e-01L;

	y = atanl(x);
	xdump("x                ", x);
	xdump("Correct  atanl(x)", ref);
	xdump("Computed atanl(x)", y);
	exit(0);
}