why is C's pow() more precise than cpow()?

When I compile the following code

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


double complex testVar;

int main(void) {
    testVar = 10;
    double complex a = 9;
    double complex b = 6;
    double aDouble = 9;
    double bDouble = 6;
    int aInt = 9;
    int bInt = 6;
    printf("Result of cpow() with double complex arguments: %.15f\n", creal(cpow(a, b)));
    printf("Result of cpow() with double arguments: %.15f\n", creal(cpow(aDouble, bDouble)));
    printf("Result of cpow() with int arguments: %.15f\n", creal(cpow(aInt, bInt)));
    printf("Result of cpow() no variables: %.15f\n", creal(cpow(9, 6)));
    printf("Result of pow(): %.15f\n", pow(aDouble, bDouble));
    }

with gcc test_cpow.c -lm -Wall (version 11.4.0 on ubuntu) and I run it I get:

Result of cpow() with double complex arguments: 531441.000000000116415
Result of cpow() with double arguments: 531441.000000000116415
Result of cpow() with int arguments: 531441.000000000116415
Result of cpow() no variables: 531441.000000000000000
Result of pow(): 531441.000000000000000

Is there a bug in cpow()?

I would expect its behavior to have the same precision as pow() if the numbers are real, although they are represented by a double complex type.

Note: this question is not answered by this, because pow() does a better job than cpow()

Solution

This is actually related to glibc, not gcc.

If we look at the glibc code for cpow, specifically at math/s_cpow_template.c, we find the following:

CFLOAT
M_DECL_FUNC (__cpow) (CFLOAT x, CFLOAT c)
{
  return M_SUF (__cexp) (c * M_SUF (__clog) (x));
}

So it looks like cpow is using the "x^y = e^{y ln x}" equivalency.

If we do the equivalent operation with real types:

printf("Result of e^(y ln x)(): %.15f\n", exp(bDouble * log(aDouble)));

The result matches what the complex types comes up with:

Result of e^(y ln x)(): 531441.000000000116415

On the other end, if we look at math/w_pow_template.c we see that pow contains the following code:

M_DECL_FUNC (__pow) (FLOAT x, FLOAT y)
{
  FLOAT z = M_SUF (__ieee754_pow) (x, y);
  /* checks for infinity / nan, removed for brevity */
  return z;
}

Which looks to be calling an intrinsic function.

So the difference between the two is that one uses "x^y = e^{y ln x}" to get the result, while the other doesn't.