Математические функции - C (СИ)

#include <stdio.h>
#include <stdlib.h>
#define M_PI ((double)3.141592653589793)
#define M_PI12 (M_PI / 12.)
#define M_PI6 (M_PI / 6.)
#define M_PI2 (M_PI / 2.)
/* квадратный корень из 3 */
#define SQRT3 ((float)1.732050807569)
double power(double, double);
 
double mfmod(double a, double b)
{
    int i = 0;
    while(b * i <= a)
        i++;
    return (a - b * (i - 1));
}
 
double mexp(double x)
{
    int i = 1;
    double res = 1, p = 1, eps = 1e-5;
    while(p > eps)
    {
        p *= x/i;
        res += p;
        i++;
    }
    return res;
}
 
double mln(double x)
{
    double res, eps = 1e-7, tmp = 1;
    int i = 2;
    x = x / (x - 1.);
    res = 1 / x;
 
    while(tmp > eps)
    {
        tmp = 1 / (i * power(x, (double)i));
        res += tmp;
        i++;
    }
    return res;
}
 
double power(double a, double b)
{
    double res = 1, fract = 0;
    fract = b -(int)b;
    if (fract) b = (int)b;
    while(b-- > 0)
        res *= a;
    if (fract) 
        res = res * mexp(fract * mln(a));
    return res;
}
 
double factorial(double x)
{
    if (x > 1)
        return x * factorial(x - 1);
    else
        return 1; 
}
 
double mcos(double x)
{
    int i, z = 1;
    double res = 1, x0, e = 0.00001, tmp = 1;
    x = mfmod(x, M_PI * 2);
    x0 = x;
    x *= x;
    i = 1;
 
    while(tmp > e)
    {
        z = -z;         
        tmp = x / factorial((double)(i * 2));
        res += z *tmp;
        x *= x0 * x0;
        i++; 
    }
    return res;
}
 
double msin(double x)
{
    int i = 1, z = 1;
    double res, x0, tmp = 1, e = 0.0001;
    x = mfmod(x, M_PI * 2);
    res = x0 = x;
    x *= x;
    while(tmp > e)
    {
        z = -z;
        tmp = x / factorial(i * 2 + 1);
        res += z * tmp;
        x *= x0 * x0; 
        i++;
    }
    return res;
}
 
double mtan(double x)
{
    if (mcos(x))
        return (msin(x) / mcos(x));
    else 
        exit(-1);
}
 
double matan(double x) 
{
    int sta = 0,sp = 0;
    double x2, a;
    /* проверка смены знака */
    if(x < 0.)
    {
        x = -x;
        sta |= 1;
    }
    /* проверка инверсии */
    if(x>1.) 
    {
        x = 1. / x;
        sta |= 2;
    }
    /* сжатие области, пока не x<PI/12 */
    while(x > M_PI12)
    {
        sp++;
        a = x + SQRT3; 
        a = 1. / a;
        x *= SQRT3; 
        x -= 1.; 
        x *= a;
    }
    /* вычисления */
    x2 = x * x;  
    a = x2 + 1.4087812;
    a = 0.55913709 / a;
    a += 0.60310579;
    a -= 0.05160454F * x2;
    a *= x;
    /* пока не sp=0 */
    while(sp > 0)
    {
        a += M_PI6;
        sp--;
    }
    if(sta & 2)
        a = M_PI2 - a;
    if(sta & 1) 
        a = -a;
 
    return a;
}
 
double msqrt(double a)
{
    return power(a, .5);
}
double masin(double x)
{
    return matan(x / msqrt(1 - x * x));
}
 
double macos(double x)
{
    return (2 * matan(msqrt((1 - x) / (1 + x))));
}
 
int mabs(int a)
{
    return (a > 0) ? a : -a;
}
 
double mfabs(double a)
{
    return (a > 0.) ? a : -a;
}
 
int main(void)
{
    double a = 2., b = 2.2;
    printf("%lg^%lg = %lf", a, b, power(a, b));
    getchar();
    return 0;
}