Возведение вещественного числа в вещественную степень

#pragma once
#include <stdio.h>
#include <time.h>
 
const double EPS = 1.e-15;
 
double Simpson(double f_sh(double), double a, double b, int n = 64);//формула Симпсона (интегрирование)
 
double FAbs(double x);//модуль числа
 
double power(double a, double b);//Возведение вещественного числа в вещественную степень (кроме отрицательных, их только в целые)
double ln(double x);//логарифм
double powerd(double a, long b);//возведение в целые степени

#include "funcs.h"
 
double FAbs(double x)//модуль
{
    if(x > 0)
        return x;
    return -x;
}
 
double Simpson(double f_sh(double), double a, double b, int n)//h*(f0+f2n+2*(f2k)+4*(f2k-1))/3
{//интеграл от функции
    double integ = 0, h = (b - a) / (2 * n), x0 = a;
    integ = f_sh(b) + f_sh(a);
    for (int i = 1; i < 2 * n; i++)
    {
        x0 += h;
        integ += 2 * f_sh(x0) * ((i + 1) % 2) + 4 * f_sh(x0) * (i % 2);
    }
    return integ * h / 3;
}
 
double powerd(double a, long b)
{   
    double S = 1.0;
    while(b != 0)
    {
        if((b & 1) != 0)
            S *= a;
        a *= a;
        b >>= 1;
    }
    return S;
}
 
double ln(double x)
{
    x = (x - 1)/(1 + x);
    double el = x, S = 0;
    x *= x;
    for(int i = 1; FAbs(el) > EPS; i += 2)
    {
        S += el / i;
        el *= x;
    }
    return S * 2;
}
 
double power(double a, double b)
{
    
    if(a == 0 || a == 1 || a == -1)
        return (int)b % 2 ? a : a * a;
    
    double tmp = FAbs(b), sum2 = 0, sum1 = powerd(a, (long)tmp);
    tmp -= (long) tmp;
    
    if(tmp > EPS)
    {
        double LN = (FAbs(a) > 1 ? ln(a) : -ln(1.0 / a)) * tmp;
        double el = 1.0;
        for(int i = 1; FAbs(el) > EPS; i++)
        {
            sum2 += el;
            el *= LN / i;
        }
    }   
    else
        sum2 = 1;
    
    sum2 *= sum1;
    if(b > 0)
        return sum2;
    return 1.0 / sum2;
}

#include "funcs.h"
#include <math.h>
 
int main()
{   
    for(double a = -2.4; a < 2.6; a += 1)
    {
        for(double b = -1.0; b < 3.1; b += 3)
            printf("%g^%g = %3.17lf ", a, b, power(a, b));
        puts("\n");
    }
    
    for(double a = -2.4; a < 2.6; a += 1)
    {
        for(double b = -1.0; b < 3.1; b += 3)
            printf("%g^%g = %3.17lf ", a, b, pow(a, b));
        puts("\n");
    }
    
    time_t st, fin;
    
    st = clock();
    for(int i = 0; i < 100000; i++)
        power(2, 10);
    fin = clock();
    printf("\n%d\n", fin - st);
    
    for(int i = 0; i < 100000; i++)
        pow(2, 10);
    fin = clock();
    printf("\n%d\n", fin - st);
    
    st = clock();
    for(int i = 0; i < 100000; i++)
        power(2, 10.2);
    fin = clock();
    printf("\n%d\n", fin - st);
    
    for(int i = 0; i < 100000; i++)
        pow(2, 10.2);
    fin = clock();
    printf("\n%d\n", fin - st);
    
            
    getchar();
    return 0;
}