[Bug math/14469] New: Inaccurate j0f function on x86_32 and x86_64

[Bug math/14469] New: Inaccurate j0f function on x86_32 and x86_64

[liubov.dmitrieva at gmail dot com]

http://sourceware.org/bugzilla/show_bug.cgi?id=14469

             Bug #: 14469
           Summary: Inaccurate j0f function on x86_32 and x86_64
           Product: glibc
           Version: 2.17
            Status: NEW
          Severity: normal
          Priority: P2
         Component: math
        AssignedTo: unassigned@sourceware.org
        ReportedBy: liubov.dmitrieva@gmail.com
    Classification: Unclassified


Test case example:

// gcc main.c -std=c99 -lm <-m32>

#include <stdio.h>
#include <fenv.h>
#include <stdlib.h>  // abs function for int argument

extern float j0f(float);
extern double j0(double);

#define i32(a)    (((int*)&(a))[0])

int main() {
    float a,r1,r2;

    i32(a)=0x413caa1f;    //argument

    fesetround(FE_TONEAREST);

    r1 = j0f(a);  //actual result
    r2 = j0(a);   //actual result 2

    printf("inputs: a = (0x%08x)\t%e\n", i32(a), a);
    printf("actual    = (0x%08x)\t%.13e\n", i32(r1), r1);
    printf("expected  = (0x%08x)\t%.13e\n", i32(r2), r2);
    printf("error     = %d ulp\n", abs(i32(r1) - i32(r2)));

    return 0;
}

Results x86_32:

inputs: a = (0x413caa1f)        1.179153e+01
actual    = (0xb471e828)        -2.2529331999976e-07
expected  = (0xb471d452)        -2.2522115727952e-07
error     = 5078 ulp

Results x86_64:

inputs: a = (0x413caa1f)        1.179153e+01
actual    = (0xb471d9d6)        -2.2524122300638e-07
expected  = (0xb471d452)        -2.2522115727952e-07
error     = 1412 ulp

-- 
Configure bugmail: http://sourceware.org/bugzilla/userprefs.cgi?tab=email
------- You are receiving this mail because: -------
You are on the CC list for the bug.

[Bug math/14469] Inaccurate j0f function on x86_32 and x86_64

[liubov.dmitrieva at gmail dot com]

http://sourceware.org/bugzilla/show_bug.cgi?id=14469

--- Comment #1 from Liubov Dmitrieva <liubov.dmitrieva at gmail dot com> 2012-08-15 14:24:41 UTC ---
*** Bug 14468 has been marked as a duplicate of this bug. ***

-- 
Configure bugmail: http://sourceware.org/bugzilla/userprefs.cgi?tab=email
------- You are receiving this mail because: -------
You are on the CC list for the bug.

[Bug math/14469] Inaccurate j0f function on x86_32 and x86_64

[joseph at codesourcery dot com]

http://sourceware.org/bugzilla/show_bug.cgi?id=14469

--- Comment #2 from joseph at codesourcery dot com <joseph at codesourcery dot com> 2012-08-15 15:11:00 UTC ---
For reference for anyone working on improving accuracy for Bessel 
functions (an issue for all the versions for all floating-point formats, 
not just x86 / x86_64): http://www.cl.cam.ac.uk/~jrh13/papers/bessel.ps.gz 
discusses how to achieve correct rounding for Bessel functions of any 
fixed constant order (jn/yn will be harder).

-- 
Configure bugmail: http://sourceware.org/bugzilla/userprefs.cgi?tab=email
------- You are receiving this mail because: -------
You are on the CC list for the bug.

[Bug math/14469] Inaccurate j0f function on x86_32 and x86_64

[liubov.dmitrieva at gmail dot com]

http://sourceware.org/bugzilla/show_bug.cgi?id=14469

--- Comment #3 from Liubov Dmitrieva <liubov.dmitrieva at gmail dot com> 2012-08-16 08:06:33 UTC ---
>> please give testcases that use hex floats

New test case for J0f Bessel function. 

#include <stdio.h>
#include <math.h>
#include <fenv.h>

extern float j0f(float);
extern double j0(double);

int main() {
    float a,r1;
    double r2;

    a=0x1.79543ep3;    //argument

    fesetround(FE_TONEAREST);

    r1 = j0f(a);  //actual result
    r2 = j0(a);   //actual result 2

    printf("inputs: a = %.6a\n", a);
    printf("actual    = %.6a\n", r1);
    printf("expected  = %.6a\n", r2);
    printf("error     = %f ulp\n", fabs(r1 - r2) * exp2(23.0-logbf(r1)));
    return 0;
}

64 bit result:
inputs: a = 0x1.79543ep+3
actual    = -0x1.e3b3acp-23
expected  = -0x1.e3a8a3p-23
error     = 1412.411834 ulp

32 bit result:
inputs: a = 0x1.79543ep+3
actual    = -0x1.e3d050p-23
expected  = -0x1.e3a8a3p-23
error     = 5078.411826 ulp

-- 
Configure bugmail: http://sourceware.org/bugzilla/userprefs.cgi?tab=email
------- You are receiving this mail because: -------
You are on the CC list for the bug.

[Bug math/14469] Inaccurate j0f function

[jsm28 at gcc dot gnu.org]

https://sourceware.org/bugzilla/show_bug.cgi?id=14469

Joseph Myers <jsm28 at gcc dot gnu.org> changed:

           What    |Removed                     |Added
----------------------------------------------------------------------------
            Summary|Inaccurate j0f function on  |Inaccurate j0f function
                   |x86_32 and x86_64           |

-- 
You are receiving this mail because:
You are on the CC list for the bug.

[Bug math/14469] Inaccurate j0f function

[fweimer at redhat dot com]

https://sourceware.org/bugzilla/show_bug.cgi?id=14469

Florian Weimer <fweimer at redhat dot com> changed:

           What    |Removed                     |Added
----------------------------------------------------------------------------
              Flags|                            |security-

-- 
You are receiving this mail because:
You are on the CC list for the bug.

[Bug math/14469] Inaccurate j0f function

[zimmerma+gcc at loria dot fr]

https://sourceware.org/bugzilla/show_bug.cgi?id=14469

--- Comment #11 from Paul Zimmermann <zimmerma+gcc at loria dot fr> ---
with glibc-2.32, the maximum error is now 6177902 ulps on x86_64 (for
x=3.153646966e+38)

-- 
You are receiving this mail because:
You are on the CC list for the bug.

[Bug math/14469] Inaccurate j0f function

[zimmerma+gcc at loria dot fr]

https://sourceware.org/bugzilla/show_bug.cgi?id=14469

Paul Zimmermann <zimmerma+gcc at loria dot fr> changed:

           What    |Removed                     |Added
----------------------------------------------------------------------------
   Target Milestone|---                         |2.34
             Status|NEW                         |RESOLVED
         Resolution|---                         |FIXED

--- Comment #12 from Paul Zimmermann <zimmerma+gcc at loria dot fr> ---
fixed by commit 9acda61

-- 
You are receiving this mail because:
You are on the CC list for the bug.