A short pow(float x, float y) = x^y approximation that doesn't use any libraries (not even math.h.)

powff.c

#include <stdio.h>

// Iterations for exp(x) approximation. Higher = more accurate, but higher likelihood of running into inf.
#define EXP_ITERATIONS 100

// Iterations for ln(x) approximation. Higher = more accurate, but slower.
#define LN_ITERATIONS 10000

// Returned if invalid input entered.
#define ERROR_RESULT -999999999

// Super fast x^p (integer power) using Exponentiation by Squaring.
double powfi(double x, unsigned p) {
    double result = 1;

    while (p) {
        if (p & 0x1) {
            result *= x;
        }
        x *= x;
        p >>= 1;
    }

    return result;
}

// Factorial
double fact2(unsigned n) {
    if (n < 0 || n >= 172) return ERROR_RESULT; // error out (fact(172) = inf)

    double result = 1;
    for (int i = 2; i <= n; i++) {
        result *= i;
    }
    return result;
}

// Approximate e^x calculation using Taylor series.
// e^x = sum_{n=0}^{inf} (x^n / n!)  for x >= 0(?)
double exp2(double x) {
    if (x < 0) return ERROR_RESULT; // error out (note: it might work for negative numbers too so this line might be useless)

    double result = 0;
    for (int n = 0; n < EXP_ITERATIONS; n++) {
        result += powfi(x, n) / fact2(n);
    }
    return result;
}

// Approximate ln(x) calculation using Taylor series
// ln(x) = 2 * sum_{n=0}^{inf} (((x-1)/(x+1))^(2n-1))/(2n-1)  for x > 0
double ln2(double x) {
    if (x <= 0) return ERROR_RESULT; // error out

    double result = 0;
    for (int n = 1; n < LN_ITERATIONS; n++) {
        result += powfi((x-1)/(x+1), 2*n-1)/(2*n-1);
    }
    return 2*result;
}

// Finally, the custom pow(float x, float y) function.
// Uses the fact that pow(x,y) = exp(y*log(x))
double powff(double x, double p) {
    if (x <= 0) return ERROR_RESULT; // error out

    // split the calculation into two to discourage inf
    // y = y - (int)y + (int)y
    // so
    // x^y = x^((int)y) * x^(y-((int)y))

    // first calculate x^((int)y)
    int int_power = (unsigned)p;
    double int_pow_result = powfi(x, (unsigned)int_power);

    // then calculate x^(y-((int)y))
    double remaining_power = p - int_power;
    double remaining_result = exp2(remaining_power*ln2(x));

    // Then multiply them together!
    return int_pow_result * remaining_result;
}

// Example
int main() {
    float a = 1305.3968;
    float b = 12.530006;
    printf("%f ^ %f = %f\n", a, b, powff(a, b));

    return 0;
}