Quote:
Originally Posted by phx
With misconfigured I mean something else than round-to-nearest and extended precision.
|
Well, it certainly has to be set to "double precision" or "single precision", according to the library. As stated, a double-round (first extended precision round, then round to target precision) is in general not identical to rouding to the target precision (compute in double, round in double generates the same result).
Quote:
Originally Posted by phx
Not a problem for plain portable C programs, but AmigaOS programs which launch new tasks are at risk. This is good to keep in mind.  |
It's mentioned in the RKRM libraries that you have open each library separately from each task. Actually, the RKRMs state "all libraries", but the real problem are the math libraries.
Quote:
Originally Posted by phx
Ok. The C startup should load a NULL state. At least vbcc's does that, as soon as you link an mlib with FPU support. But not for mieee.
Does the OS do any FPU initialization, when a new task is launched?
|
Well, in the sense that if you create a new task, the stack frame is initialized such that the FPU is in NULL state, and hence, as soon as you try to use the FPU, it initializes itself to its default setup, which is extended precision, round to nearest.
Beware, however, that the shell does not create a new task for each program launched. Executables just overlay the shell, and thus "borrow" or "replace" the FPU state the last shell program left behind. IOWs, to use the FPU reliably, you either have to initialize it yourself, or you need to open the math libraries.
Quote:
Originally Posted by phx
Typo? Last time you wrote "round to zero"...?
|
Round to zero for "up to 3.1.4", round to nearest from 3.2 on.
Quote:
Originally Posted by phx
OMG! This is the big one! 
Which means I cannot use mathieeesingbas for float and mathieeedoubbas for double in the same task? This is quite common, I would say. And I always expected that to work.
|
Correct, it's really problematic. For 3.1.4, this did work, but this is due to the "round to zero" policy. That is, "round to zero in high precision" followed by "round to zero in low precision" is identical to "round to zero in low precision". This does not hold for "round to nearest", which can create different results (two times rounding != one time rounding).
3.1.4 configured the FPU to "round to zero, double precision" in both libraries, and this works due to the above property of rouding to zero (provable). Unfortunately, "round to zero" is not a good choice for numerics, so that was in a sense broken, too.
The problem is that the math libraries, on the 68030 and 68020 at least, have to use the extended precision opcodes (68040 onwards have "single precision" and "double precision" opcodes). The only other option would be to set the FPCR on each call to the right precision. I can already hear the speed geeks cry.
SAS/C got around it by using FFP for float, and otherwise "float = double" which is permissible by ANSI/C, but you really don't want to use mathffp. It's a rotten format, has no gradual underflow, and no INFs or NANs. Really bad for numerics.
The way how I got around it in DMandel is to save the FPCR after opening each library, and switch it manually before calling into the single or double precision, but that's also a bad hack.
Concerning ANSI/C, isn't it so that computing (float) (op) (float) should be implemented as (round to float) ((double) (op) (double)), that is, your code should rather perform the computation in double precision, and then convert the result to single precision upon storing it in the variable?
Quote:
Originally Posted by phx
Maybe I didn't configure it correctly, or there are problems in the target-specific support functions I have to provide. But this is what I get with Paranoia and -lmsoft (which statically links Hauser's Softfloat with vbcc) and a 68000 system
|
Interesting, thanks for your findings. pow() is indeed a rather troublesome function (and for that reason also new in 3.9). Most folks just implement it by a^b=exp(b*log(a)), but this is not a good solution as it is numerically not stable. I don't recall what Hauser's code does here.
I always had the plan to also provide a mathieeeextbas and mathieeeexttrans library, and extend the API to also include conversion to and from ASCII, but I afraid this will have to wait.