> Date: Thu, 13 Dec 2018 03:25:30 +0000 > From: Taylor R Campbell <campbell+tor-...@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...)
#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); }
_______________________________________________ tor-dev mailing list tor-dev@lists.torproject.org https://lists.torproject.org/cgi-bin/mailman/listinfo/tor-dev