Arithmetic operator of power in C

The one operator C doesn't have is an exponent operator; you'll either have to use one of the pow functions:

#include <math.h>

double pow(double x, double y);
float powf(float x, float y);
long double powl(long double x, long double y);

or, multiply it out by hand:

y = x * x * x * ...;

or, use a loop:

result = 1;
for ( size_t i = 0; i < power; i++ )
  result *= x;

You'd just have to roll your own. This seems to work (no Taylor series needed):

// Assumes IEEE 754 format for double and Intel 80-bit for long double, and
// little-endian storage.

long double ldsqrt(long double x) {
    if (x < 0)
        return 0/0.;
    if (x == x + x || x != x)
        return x;
    long double r = x, s;
    short *b = (short *)&r;
    if (sizeof(long double) > 8)
        b[4] = (b[4] >> 1) + 8191;
    else
        b[3] = (((b[3] >> 5) + 511) << 4) + (b[3] & 15);
    do {
        s = r;
        r = ((x / r) + r) / 2;
    } while (r != s);
    return r;
}

double dpow(double x, double y) {
    if (y == 0)
        return 1;
    if (x == 0)
        return 0;
    if (x == 1)
        return 1;
    int s = y < 0;
    if (s)
        y = -y;
    if (y > 3200000000000000000)
        return x < 0 ? 0/0. : x < 1 ? 0 : 1/0.;
    long long i = y;
    y -= i;
    if ((y || i > 9007199254740991) && x < 0)
        return 0/0.;
    long double p = 1, q = x;
    for (;;) {
        if (i & 1)
            p *= q;
        i >>= 1;
        if (i == 0)
            break;
        q *= q;
    }
    q = x;
    while (y) {
        q = ldsqrt(q);
        y *= 2;
        if (y >= 1) {
            p *= q;
            y -= 1;
        }
    }
    return s ? p ? 1 / p : 1/0. : p;
}

In case you have an unsigned integer exponent, you can use this fast exponential:

#include <stdio.h>

double fast_pow(double x, unsigned n)
{
    double result = 1.0;
    for (   unsigned msk = ~((~0U) >> 1); 
            msk; 
            msk >>= 1)
    {
        result *= result; /* square */
        if (msk & n) result *= x; /* multiply by x */
    }
    return result;
} /* fast_pow */

#define P(expr) printf(#expr" = %.10f\n", expr)

int main()
{
    P(fast_pow(1.0000118, 6430));
} /* main */

which produces:

$ a.out
fast_pow(1.0000118, 6430) = 1.0788261505
$ _