On 2010-Jan-01 11:56:02 +0000, "Dr. David Kirkby" <david.kir...@onetel.net> wrote: >But if FreeBSD's libm for the 387 chips is based on Sun's library for a SPARC >processor one needs to be cautious before making that assumption. If FreeBSD's >implementation was based on Sun's library for the x86 processor, then I'd say >FreeBSD got something wrong, since exp(1.0) is computed accurately with Open >Solaris running on an Xeon processor.
Sun's libm probably started on M68K. There are two distinct sets of code: The i386 code uses i387 instructions to do: x1 = x * log2(e) exp(x) = 2^int(x1) * 2^fract(x1) in 80-bit precision. The amd64 and SPARC code uses a semi-numerical expansion in C, assuming a 64-bit FPU. >Virtually all computations on reals are going to introduce some >errors. The fact the 387 can do some of those calculations in 80 bits >means there is a reason for having a more accurate library than on >the SPARC processor, which is only 64-bits. AFAIK, the original justification for including the 80-bit FP format was to reduce the need for multi-precision or semi-numerical code to provide accurate double-precision results in floating point libraries. >I do not know whether the SPARC processor has an instruction which computes >exp(x) directly. It doesn't - the SPARC only provide the basic IEEE operations. > Finding such a simple bit of information in manuals which are >100's of pages long, is no easy task. You want a SPARC Assembler Language Manual - googling turned up suitable extracts. >Based on the 64/80 bit difference, I can see a justification for >using a more accurate implementation of exp(x) on a x86 chip, than I >can on a SPARC chip. The core of i386 exp() is 11 x87 instructions (and another 22 instructions to handle special cases and ensure the correct FPU mode). The FreeBSD C version is 56 lines of C and 17 magic constants and claims less than 1ULP error. The code in glibc-2.9 comprises 100 lines of C with a fallback to multi-precision code. > If it going to take you long to compute something to an extra bit of >accuracy, only to lose that bit as soon as you do any non-trivial >computation with it, then one would have to question whether it was >worth spending the extra time in the first place. Calculating perfectly rounded results for transcendental functions is between very hard and impractical without relying on (very slow) multi-precision arithmetic. Most implementations allow a few ULP error. Whether those few bits are critical is a numerical analysis job. > Because the x86 processor works internally at 80 bits, it >should be able to preserve all 64 bits when it writes to RAM. There are some gotchas here. On an i387 in default (80-bit) mode, a local double variable may be 80-bits if kept in a register but 64-bits if written to memory. This means the precision of working variables in a function can alternate between 64 and 80 bits depending on the code, optimisation level and whim of the compiler. More precision is not always good. This extra precision (and particularly the arbitrariness of its existence) can wreak havoc in a function that expects 'double' to be 64-bits. >result more accurately, I'm inclined to think this is not a bug. I tend to agree. -- Peter Jeremy
pgpeNrBnb4TY9.pgp
Description: PGP signature
