Menu
▾
▴
[Mingw-notify] [ mingw-Bugs-2117590 ] correction of function asinh - inverse hyperbolic sine
|
From: SourceForge.net <no...@so...> - 2009-01-08 14:10:50
|
Bugs item #2117590, was opened at 2008-09-18 09:41 Message generated for change (Comment added) made by keithmarshall You can respond by visiting: https://sourceforge.net/tracker/?func=detail&atid=102435&aid=2117590&group_id=2435 Please note that this message will contain a full copy of the comment thread, including the initial issue submission, for this request, not just the latest update. Category: None Group: None Status: Open Resolution: None Priority: 5 Private: No Submitted By: Roumen Petrov (rumen) >Assigned to: Keith Marshall (keithmarshall) Summary: correction of function asinh - inverse hyperbolic sine Initial Comment: The current implementation of mingwex function with declaration "double asinh(double x)" return unexpected value for some arguments. Please to replace implementation with new one that pass at least following tests(argument->result): 1) 0.0->0.0 2) -0.0->-0.0 3) -1.0000000000000002e+299->-689.16608998577965 The current code return -0.0 for argument 0.0. Uncommenting code for "/* Avoid setting FPU underflow exception flag in x * x. */" isn't enough since the case 3) will return NaN. ---------------------------------------------------------------------- >Comment By: Keith Marshall (keithmarshall) Date: 2009-01-08 14:10 Message: Please test the attached generic implementation. Unpack it in a clean directory, then configure [--build=whatever --host=mingw32] [CFLAGS='<your choice>'] && make will create a local copy of libmingwex.a, with the three asinh() functions replaced; `make check', (or you may need `make EMULATOR=wine check', if you are cross-hosting without misc-binaries support), will run a test case, for each of the three functions, with each of the input values you've reported as contentious. Feel free to add additional cases in testcase.c; if this fits your requirements, I'll fold it in for mingwrt-3.16. File Added: asinh_generic.tar.gz ---------------------------------------------------------------------- Comment By: Roumen Petrov (rumen) Date: 2009-01-06 21:50 Message: How to processed further with this case ? What about to create more cases - first one for big numbers and second one for small numbers ? Initially I open the issue only for 0.0/-0.0 case, i.e to use "copysign(z,x)" instead expression "( x > 0.0 ? z : -z)" in all tree asinh functions. ---------------------------------------------------------------------- Comment By: Roumen Petrov (rumen) Date: 2008-11-11 00:08 Message: About overflow It seems to me that 1./sqrt(EPSILON) isn't so good. If Z is near to this value Z*Z and Z*Z+1 are different: 1./sqrt(DBL_EPSILON) = 67108864 Z*Z-1=4503599627370495 Z*Z =4503599627370496 Z*Z+1=4503599627370497 :( But if Z number is greater then ceil(sqrt(2^DBL_MANT_DIG)) the multiplication (Z*Z) is greater then 2^DBL_MANT_DIG. In this case adding 1.0 is irrelevant to the result. About underflow: Why to compute log1p(Z)? May be it has to be log(Z). If (Z*Z) < (2^-DBL_MANT_DIG) (note minus sign) then 1+Z*Z is equal to 1 and the result has to be copysign(log(Z),X). But the tests with DBL_MIN and 2^-28 return: asinh(2.2250738585072014e-308) = 2.2250738585072014e-308 asinh(3.7252902984619140625e-009) = 3.7252902984619141e-009 , so may be don't need to change code. Also updated patch attached. File Added: mingwex-asinhX.diff ---------------------------------------------------------------------- Comment By: Roumen Petrov (rumen) Date: 2008-09-26 21:26 Message: About big numbers - ok for 1/sqrt(x) as limit in asinhl. About implementation of similar limit in asinh/asinhf - I don't know gnu assembler (I forgot all other assemblers) and all __fast_log functions looks identical :) . No idea for small numbers. ----- [ 2009559 ] @CFLAGS@ not substituted in Makefiles" -> http://sourceforge.net/tracker/index.php?func=detail&aid=2009559&group_id=2435&atid=302435 ----- My repository(<repo>) is from CVSROOT=:pserver:an...@cy...:/cvs/src . In a separate directory(<SRCCOPY>) I keep a copy of subdirectories mingw and win32 from <repo>/src/winsup/. I run configure script and build only in <SRCCOPY>/mingw/mingwex/ . ----- flag -nostdinc and as example xxx_EPSILON: $ cd /empty_directory 1) $ printf "%s\n%s\n" "#include <float.h>" "long double a = LDBL_EPSILON" | i386-mingw32...-gcc -E -nostdinc - # 1 "<stdin>" # 1 "<built-in>" # 1 "<command line>" # 1 "<stdin>" <stdin>:1:19: no include path in which to search for float.h long double a = LDBL_EPSILON 2) $ printf .... | i386-mingw32...-/gcc -E -nostdinc -I/path_to_mingwrt/include - In file included from <stdin>:1: /path_to_mingwrt/include/float.h:19:23: no include path in which to search for float.h .... # 2 "<stdin>" 2 long double a = LDBL_EPSILON 3) $ printf .... | i386-mingw32...-/gcc -E - # 1 "<stdin>" # 1 "<built-in>" # 1 "<command line>" # 1 "<stdin>" # 1 "..../lib/gcc/i386-mingw32.../3.4.5/include/float.h" 1 3 4 # 2 "<stdin>" 2 long double a = 1.08420217248550443401e-19L As I understand -nostdinc flag is problem(not so important) in my build environment. Since #include <float.h> will work in you and Chris environments - then it is fine code to use generic standard defines. ---------------------------------------------------------------------- Comment By: Keith Marshall (keithmarshall) Date: 2008-09-26 10:07 Message: Roumen, >> Also for long double I would like to propose an additional >> modification that deal with big long double. > > Ok, thanks. That's basically a useful refinement of the technique > I proposed; I've no issue with adopting it... However, I am curious as to why you propose if (z > 0.5/__LDBL_EPSILON__) as the cut-off point, for switching from the analytically exact computation to the rounded approximation? Surely a more analytically appropriate expression is if( z > 1.0L / sqrtl( LDBL_EPSILON ) ) (and BTW, we should `#include <float.h>', and use the generic standard LDBL_EPSILON, rather than rely on a compiler specific definition such as __LDBL_EPSILON__). Granted, with this expression, the compiler is unlikely to recognise that the RHS is a constant, so will not be able to optimise away the sqrtl() call -- we'd need to help it out, by providing a suitable manifest constant definition, if we want that optimisation. Is this the basis of your proposal? You know `0.5 / EPSILON' > `1 / sqrt(EPSILON)', yet should still be safe to square, and you are prepared to ignore the insignificance of unity in the intermediate interval, to have a generically reproducible and recognisably constant term on the RHS of the comparison? I also wonder if there might be merit in adding if( z < sqrtl( LDBL_EPSILON ) ) return copysignl( __fast_log1pl( z ), x ); or, (if we accept the `0.5 / EPSILON' paradigm for the overflow case) if( z < (2.0L * LDBL_EPSILON) ) return copysignl( __fast_log1pl( z ), x ); to trap any potential underflow issues. ---------------------------------------------------------------------- Comment By: Keith Marshall (keithmarshall) Date: 2008-09-23 13:26 Message: > Compiler is build from package > gcc-core-3.4.5-20060117-1-src.tar.gz As is mine, using the x86-mingw32-build.sh script package, from our SF download site, with the installation directory set as $HOME/mingw32, and the cross-compiler identifying name set simply as `mingw32'. > To compile minwgex I modify configure.in Why do you do this? As a user, you should have no need to *ever* modify configure.in. > to substitute CFLAGS (also i see this reported in the tracker). Which tracker would that be, then? > Also I remove flag -nostdinc from Makefile.in, > otherwise build fail - missing include header stddef.h. You *definitely* should not need to do this. It may not be wrong to do so, but if you are referring to the tracker item I think you are, then I do not think that the solution proposed there is necessarily the correct way to address the problem. > Dunno why after modification my library return NaN for big double. Nor do I; I certainly cannot reproduce this behaviour. FWIW, here's how I build, from CVS, with the following sandbox structure:-- $HOME/mingw-sandbox | +-- build | | | +-- runtime | +-- runtime (the CVS sources are checked out to $HOME/mingw-sandbox/runtime) $ cd $HOME/mingw-sandbox/build/runtime $ ../../runtime/configure --prefix=$HOME/mingw32 \ --build=i686-pc-linux --host=mingw32 $ make CFLAGS='-s -O2 -mtune=i686 -mms-bitfields' $ make install With this, I keep my cross-compiler up to date with the current state of the MinGW CVS. Note that I do not use this build tree for creating distributable packages -- I leave that to Chris. The Makefiles are broken in this regard; they require entirely inappropriate misuse of configure's --target spec, to achieve a correctly structured MinGW package directory hierarchy. > Also for long double I would like to propose an additional > modification that deal with big long double. Ok, thanks. That's basically a useful refinement of the technique I proposed; I've no issue with adopting it, but I think the explanatory comments could be improved, to make them more acceptable in a mathematically correct sense. You may leave it with me, to come up with suitable wording. I also think it may be advantageous to handle the potential overflow similarly, within asinhf() and asinh(), rather than relying on the obfuscated method of the present implementation. ---------------------------------------------------------------------- Comment By: Roumen Petrov (rumen) Date: 2008-09-20 14:07 Message: Also I attach my new tests case build with (part from Makefile): test-asinhex-my.exe: test-asinh.c $(MINGWCC) -posix -o $@ $^ -lmingwex-my The output is (mysterious nan for asinh is still here, but output for asinhl,asinhf look good): === sizeof(double)=8 asinh(1.0000000000000002e+299)=nan ..... z=688.47294280521965 asinh(-1.0000000000000002e+299)=nan ..... z=-688.47294280521965 asinh(1.7976931348623157e+308)=nan ..... z=709.78271289338397 asinh(-1.7976931348623157e+308)=nan ..... z=-709.78271289338397 asinh(3.4028234663852886e+038)=89.415986232628299 ..... z=89.415986232628299 asinh(-3.4028234663852886e+038)=-89.415986232628299 ..... z=-89.415986232628299 asinh(0)=0 ..... z=0 asinh(-0)=-0 ..... z=-0 === sizeof(long double)=12 asinhl(1.189731495357231765021264e+4932)=11357.21655347470389507691 asinhl(-1.189731495357231765021264e+4932)=-11357.21655347470389507691 asinhl(1.797693134862315708145274e+308)=710.4758600739439420301835 asinhl(-1.797693134862315708145274e+308)=-710.4758600739439420301835 asinhl(3.402823466385288598117042e+038)=89.41598623262829834135168 asinhl(-3.402823466385288598117042e+038)=-89.41598623262829834135168 asinhl(0)=0 asinhl(-0)=-0 === sizeof(float)=4 asinhf(3.40282347e+038)=89.4159851 asinhf(-3.40282347e+038)=-89.4159851 asinhf(0)=0 asinhf(-0)=-0 File Added: test-asinh.c ---------------------------------------------------------------------- Comment By: Roumen Petrov (rumen) Date: 2008-09-20 13:57 Message: May I ask you to apply similar modification to math/asinh{l|f}.c with call to copysign{l|f}. I'm tired. I couldn't find why my self-build library return NaN for big numbers. If I add code from asinh (with call to copysign) into test program along with inline function from fastmath.h - it work as expected, i.e. without NaN for big numbers. Compiler is build from package gcc-core-3.4.5-20060117-1-src.tar.gz To compile minwgex I modify configure.in to substitute CFLAGS (also i see this reported in the tracker). Also I remove flag -nostdinc from Makefile.in, otherwise build fail - missing include header stddef.h. Dunno why after modification my library return NaN for big double. Also for long double I would like to propose an additional modification that deal with big long double. File Added: mingwex-asinhl.diff ---------------------------------------------------------------------- Comment By: Keith Marshall (keithmarshall) Date: 2008-09-19 17:20 Message: I applied this: Index: mingwex/math/asinh.c =================================================================== RCS file: /cvs/src/src/winsup/mingw/mingwex/math/asinh.c,v retrieving revision 1.1 diff -u -r1.1 asinh.c --- mingwex/math/asinh.c 6 Oct 2004 20:31:32 -0000 1.1 +++ mingwex/math/asinh.c 19 Sep 2008 17:14:04 -0000 @@ -23,6 +23,6 @@ z = __fast_log1p (z + z * z / (__fast_sqrt (z * z + 1.0) + 1.0)); - return ( x > 0.0 ? z : -z); + return copysign (z, x); } It fixes the signed zero issue; I STILL cannot reproduce your other fault, (this time just running under wine). ---------------------------------------------------------------------- Comment By: Roumen Petrov (rumen) Date: 2008-09-19 15:25 Message: May NaN is because my mingwex library is modified to use copysign instead xx> 0.0 ... Yes in the original there is no problem for the big numbers - only sign of zero. ---------------------------------------------------------------------- Comment By: Keith Marshall (keithmarshall) Date: 2008-09-19 13:21 Message: FTR, I've now taken your test case to my own Ubuntu 8.04 box, compiled with mingw32-gcc (GCC-3.4.5, mingwrt-3.15), and run it under wine-1.0. It behaves correctly; I *cannot* reproduce the fault you report in (3). ---------------------------------------------------------------------- Comment By: Keith Marshall (keithmarshall) Date: 2008-09-19 11:10 Message: Case (1) clearly isn't an issue, for the result is correct. Case (2) is purely an issue of discriminating between positive and negative zero[1]. I agree with your assessment, that `return (z > 0.0 ? z : -z);' is incorrect, for it must *always* return a zero as negative zero. Similarly, `return (z >= 0.0 ? z : -z);' is incorrect, for it must always return zero as *positive* zero. The correct solution is to either check for `if( x == 0.0 ) return x;' on entry, or to always `return copysign( z, x );' Case (3) is *not* a MinGW bug, (or at least I cannot reproduce it[2]). I am confused by your references to native vs. emulated environment; do you mean, perhaps, that you get incorrect results running on GNU/Linux under Wine? If so, then surely this problem lies in Wine, and you should follow up on their ML. Also note that uncommenting the code to avoid setting the FPU underflow exception flag would not be expected to have any effect in case (3), for it affects only very *small* values of `x', and case (3) has a fairly *large* (negative) value of `x'. [1] Using the MSVCRT implementation of printf(), I appear to see no distinction between positive and negative zero anyway; it only becomes apparent when using the mingwex printf() functions. [2] Although the current mingwex implementation appears to produce correct results for me, in each of the test cases given, it does use what seems a rather suspect method of guarding against overflow, in its internal computation of `x * x', (or `z * z' as it is actually coded); I believe that it may be worthwhile replacing this with a more mathematically robust formulation: long double my_asinhl( long double x ) { if( isfinite( x ) && (x != 0.0L) ) { long double z = fabsl( x ), zz = z * z; if( isfinite( zz = z * z ) ) /* * `z' has been squared without overflow... * Compute the inverse sinh, using the analytically correct formula, * (using log1p() for greater precision at very small `z', than * can be achieved using log()). */ return copysign( log1pl( z + sqrtl( zz + 1.0L ) - 1.0L ), x ); /* Computation of `z squared' results in overflow... * Here we may employ the approximation that, for very large `z', * `sqrt( z * z + 1 ) ~= z', yielding the approximate result */ return copysign( logl( z ) + log1pl( 1.0L ), x ); } return x; } double my_asinh( double x ) { return (double)(my_asinhl( (long double)(x) )); } ---------------------------------------------------------------------- Comment By: Roumen Petrov (rumen) Date: 2008-09-18 13:43 Message: more tests: If we change current code from "return ( x > 0.0 ? z : -z);" to "return ( x >= 0.0 ? z : -z);" the function result is not changed. If we replace return statement with "return copysign(z,x)" the result for {+/-}0.0 argument is as expected in native and emulated environment. But if we try to use function _copysign the result is ok in native environment but fail in emulated. Also with similar changes for functions asinhf/asinhl with respective copysign, i.e. copysignf/copysignl, I see expected result. Did someone know how to resolve problem with big exponent ? ---------------------------------------------------------------------- Comment By: Roumen Petrov (rumen) Date: 2008-09-18 12:13 Message: more test vectors: [-]1.0000000000000002e+99 ->[-]228.64907138697046 //ok [-]1.0000000000000001e+199->[-]458.90758068637501 //but mingwex-nan, native - ok [-]1.0000000000000002e+299->[-]689.16608998577965 //but mingwex-nan, native - ok [-]1.6025136110019349e+308->[-]710.3609292261076 //but mingwex-nan, native - ok ---------------------------------------------------------------------- Comment By: Roumen Petrov (rumen) Date: 2008-09-18 10:53 Message: also asinhl/asinhf return -0.0 too for 0.0 argument. ---------------------------------------------------------------------- You can respond by visiting: https://sourceforge.net/tracker/?func=detail&atid=102435&aid=2117590&group_id=2435 |