[tor-dev] prob_distr.c: LogLogistic fails stochastic tests on 32-bits mingw

[Taylor R Campbell]

> Date: Thu, 13 Dec 2018 03:25:30 +0000
> From: Taylor R Campbell <campbell+tor-dev at mumble.net>
> 
> Binary80 arithmetic tickled a problem in the lerp used for binning --
> [...]

Correction: while making the lerp less naive addresses this problem,
it also arises only when binary80 arithmetic and binary64 arithmetic
are mixed, which you get on 32-bit x86 in C double with the x87 unit
in binary80 mode so (a) double is binary64, but (b) intermediate
expressions are evaluated in binary80.

The attached program shows this, by putting either one (bad) or two
(good) intermediate expressions in double variables.  On x86, if you
compile it with -mfpmath=387 (default on 32-bit), you'll see the bad
result, a negative answer; if you compile it with -mfpmath=sse
(default on 64-bit), you'll see only good results, zero.  Convert
everything to use long double (and %Le) instead so that all the
arithmetic and intermediate quantities are binary80, and it's fine.

% cc -o loser -mfpmath=387 loser.c && ./loser
bad -2.45760000000000000e+04
good 0.00000000000000000e+00
% cc -o loser -mfpmath=sse loser.c && ./loser
bad 0.00000000000000000e+00
good 0.00000000000000000e+00

(This is why I don't like x87 and the automagic evaluation of
expressions in higher precision...)
-------------- next part --------------
#include <stdio.h>

#ifdef __NetBSD__

#include <ieeefp.h>

#else

#include <stdint.h>

#define	FP_PS	0
#define	FP_PRS	1
#define	FP_PD	2
#define	FP_PE	3

int
fpsetprec(int nprec)
{
	uint32_t ocw, ncw;
	int oprec;

	asm volatile("fnstcw %0" : "=m"(ocw));

	oprec = (ocw >> 8) & 3;
	ncw = (ocw & ~(uint32_t)(3 << 8)) | (uint32_t)((nprec & 3) << 8);

	asm volatile("fldcw %0" : : "m"(ncw));

	return oprec;
}

#endif

int
main(void)
{
	volatile double lo = 2.4608250784829636e-20;
	volatile double hi = 3.0026742508190853e+20;
	volatile double w, d;
	volatile size_t n = 100;

	fpsetprec(FP_PE);
	w = (hi - lo)/(n - 2);
	d = (n - 2)*w;
	printf("bad %.17e\n", hi - (n - 2)*w);
	printf("good %.17e\n", hi - d);

	fflush(stdout);
	return ferror(stdout);
}