[Bug math/14473] New: Inaccurate CPOWF function on x86_32 and x86_64

[Bug math/14473] New: Inaccurate CPOWF function on x86_32 and x86_64

[liubov.dmitrieva at gmail dot com]

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

             Bug #: 14473
           Summary: Inaccurate CPOWF  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


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

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

#define ireal(a)    (((int*)&(a))[0])
#define iimag(a)    (((int*)&(a))[1])
#define  real(a)    (((float*)&(a))[0])
#define  imag(a)    (((float*)&(a))[1])

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

    ireal(a)=0xbf0e38e4;    iimag(a)=0xbde38e3e;    //argument
    ireal(b)=0x3f0e38e3;    iimag(b)=0xbde38e3e;    //argument

    fesetround(FE_TONEAREST);

    r1 = cpowf(a,b);  //actual result
    r2 = cpow(a,b);   //actual result 2

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

    return 0;
}

Result x86_32 (core2):

inputs: a.real = (0xbf0e38e4)   -5.555556e-01
        a.imag = (0xbde38e3e)   -1.111111e-01

inputs: b.real = (0x3f0e38e3)   5.555555e-01
        b.imag = (0xbde38e3e)   -1.111111e-01

actual.real    = (0xba6f8923)   -9.1375614283606e-04
expected.real  = (0xba6f8d77)   -9.1382063692436e-04
error.real     = 1108 ulp

actual.imag    = (0xbf069c6d)   -5.2582436800003e-01
expected.imag  = (0xbf069c6d)   -5.2582436800003e-01
error.imag     = 0 ulp


Result x86_64 (core2):

inputs: a.real = (0xbf0e38e4)   -5.555556e-01
        a.imag = (0xbde38e3e)   -1.111111e-01

inputs: b.real = (0x3f0e38e3)   5.555555e-01
        b.imag = (0xbde38e3e)   -1.111111e-01

actual.real    = (0xba6f8923)   -9.1375614283606e-04
expected.real  = (0xba6f8d77)   -9.1382063692436e-04
error.real     = 1108 ulp

actual.imag    = (0xbf069c6d)   -5.2582436800003e-01
expected.imag  = (0xbf069c6d)   -5.2582436800003e-01
error.imag     = 0 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/14473] Inaccurate CPOWF function on x86_32 and x86_64

[joseph at codesourcery dot com]

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

--- Comment #1 from joseph at codesourcery dot com <joseph at codesourcery dot com> 2012-08-15 15:22:55 UTC ---
Again, please give testcases that use hex floats rather than using type 
punning, so that the test inputs can more readily be extracted for use 
elsewhere.

cpowf, cpow, cpowl are generically inaccurate for all architectures and 
floating-point formats - there's nothing x86-specific.  They are also by 
far the hardest complex functions to make accurate; the errors in this 
example are tiny compared to those you can get with properly chosen 
examples.  In general for long double you need to compute log and atan2 
results to over 16000 places to get good results for cpowl.  Consider 
cpowl (LDBL_MAX, I * LDBL_MAX) for example - you need to compute LDBL_MAX 
* log (LDBL_MAX), reduced mod 2pi.  Many cases will overflow or underflow 
but you still need to determine the signs of zero or infinity in those 
cases.

-- 
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/14473] Inaccurate CPOWF function on x86_32 and x86_64

[liubov.dmitrieva at gmail dot com]

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

--- Comment #2 from Liubov Dmitrieva <liubov.dmitrieva at gmail dot com> 2012-08-16 07:51:45 UTC ---
Reformatted test case:

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

#define  real(a)    (((float*)&(a))[0])
#define  imag(a)    (((float*)&(a))[1])

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

    a=-0x1.1c71c8p-1 + I * -0x1.c71c7cp-4; //argument
    b= 0x1.1c71c6p-1 + I * -0x1.c71c7cp-4; //argument

    fesetround(FE_TONEAREST);

    r1 = cpowf(a,b);  //actual result
    r2 = cpow(a,b);   //actual result 2

    printf("inputs: a = %15.6a + I * %15.6a\n", creal(a), cimag(a));
    printf("        b = %15.6a + I * %15.6a\n", creal(b), cimag(b));
    printf("actual    = %15.6a + I * %15.6a\n", creal(r1),cimag(r1));
    printf("expected  = %15.6a + I * %15.6a\n", creal(r2),cimag(r2));
    printf("error     = %15f + I * %15f ulp\n", fabs(creal(r1) - creal(r2)) *
exp2(23.0-logbf(creal(r1))),
                                                fabs(cimag(r1) - cimag(r2)) *
exp2(23.0-logbf(cimag(r1))) );
    return 0;
}


Results (same x86 and x86_64)

inputs: a =  -0x1.1c71c8p-1 + I *  -0x1.c71c7cp-4
        b =   0x1.1c71c6p-1 + I *  -0x1.c71c7cp-4
actual    = -0x1.df1246p-11 + I *  -0x1.0d38dap-1
expected  = -0x1.df1aefp-11 + I *  -0x1.0d38d9p-1
error     =  1108.410103 + I *  0.495932 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/14473] Inaccurate CPOWF function on x86_32 and x86_64

[glisse at gcc dot gnu.org]

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

Marc Glisse <glisse at gcc dot gnu.org> changed:

           What    |Removed                     |Added
----------------------------------------------------------------------------
                 CC|                            |glisse at gcc dot gnu.org

--- Comment #3 from Marc Glisse <glisse at gcc dot gnu.org> 2013-01-25 16:48:59 UTC ---
Hello,

as an additional data point, using cpow to compute i^2 yields a number with a
non-zero imaginary part, as was noticed in later messages of the conversation
started here:

http://gcc.gnu.org/ml/libstdc++/2013-01/msg00058.html

-- 
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/14473] Inaccurate CPOW* function on all machines.

[carlos at redhat dot com]

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

Carlos O'Donell <carlos at redhat dot com> changed:

           What    |Removed                     |Added
----------------------------------------------------------------------------
   Last reconfirmed|                            |2013-04-04
                 CC|                            |carlos at redhat dot com
            Summary|Inaccurate CPOWF  function  |Inaccurate CPOW*  function
                   |on x86_32 and x86_64        |on all machines.

--- Comment #4 from Carlos O'Donell <carlos at redhat dot com> 2013-04-04 16:12:37 UTC ---
Related LSB test suite failure "/olver/math/cexp/tests/cpow_spec"
https://lsbbugs.linuxfoundation.org/show_bug.cgi?id=2912

-- 
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/14473] Inaccurate CPOW* function on all machines.

[carlos at redhat dot com]

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

Carlos O'Donell <carlos at redhat dot com> changed:

           What    |Removed                     |Added
----------------------------------------------------------------------------
           See Also|                            |https://lsbbugs.linuxfounda
                   |                            |tion.org/show_bug.cgi?id=29
                   |                            |12

-- 
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/14473] Inaccurate CPOW* function on all machines.

[licquia at linuxfoundation dot org]

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

Jeff Licquia <licquia at linuxfoundation dot org> changed:

           What    |Removed                     |Added
----------------------------------------------------------------------------
                 CC|                            |licquia at linuxfoundation dot org

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

[Bug math/14473] Inaccurate CPOW* function on all machines.

[mwelinder at gmail dot com]

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

M Welinder <mwelinder at gmail dot com> changed:

           What    |Removed                     |Added
----------------------------------------------------------------------------
                 CC|                            |mwelinder at gmail dot com

--- Comment #5 from M Welinder <mwelinder at gmail dot com> ---
There are at least two classes of accuracy issues with cpow(z,z2):

(1) Insufficient precision in intermediate results, affecting primarily "large"
    arguments

(2) Bad handling of |z1| == 1.  Sample program below.  This is actually
    easy to solve: use atan2pi, sinpi, and cospi instead of atan2, sin, and
    cos.


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

int dont_optimize = 0;

int
main (int argc, char **argv)
{
  double complex z1, z2, z;

  /*
   * z here should be real.  I see an imaginary part of
   * 2.54571e-17, except when gcc is allowed to optimize
   * the whole thing at compile time.  Exact result should
   * be -exp(-pi/2) =~ -0.20788
   *
   * Surprise: Mathematica gets it wrong too.  I see the
   * same value as gcc for N[I^(I+2)].
   */
  z1 = dont_optimize + I;
  z2 = 2 + I;
  z = cpow (z1, z2);
  printf ("%g + %g i\n", creal (z), cimag (z));

  /*
   * Really scary version.  Exact result is exp(-pi/2)
   * =~ 0.20788.  I see "0.177058 + -0.108924 i".
   */
  z2 = ldexp (1,53) + I;
  z = cpow (z1, z2);
  printf ("%g + %g i\n", creal (z), cimag (z));

  return 0;
}

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

[Bug math/14473] Inaccurate CPOW* function on all machines.

[mwelinder at gmail dot com]

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

--- Comment #6 from M Welinder <mwelinder at gmail dot com> ---
Ok, make that three classes of inaccuracy:

(3) Bad handling of real arguments.  I see different values in this
    case:

  z1 = dont_optimize + 607/128;
  z2 = 1234567 / 8192;
  z = cpow (z1, z2);
  printf ("%.20g + %g i\n", creal (z), cimag (z));
  printf ("%.20g\n", pow (creal (z1), creal (z2)));

I get:

2.0370359763344673153e+90 + 0 i
2.0370359763344860863e+90


For real arguments cpow should defer to pow.  (Or at least use a method
that is no less accurate than pow; I think the current evidence is that
pow is better.)

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

[Bug math/14473] Inaccurate CPOW* function on all machines.

[fweimer at redhat dot com]

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

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.