Возведение числа в степень Си

Question

Нужно сделать свою pow, вот код:

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <stdbool.h>
#define s21_PI 3.141592654
#define eps 0.0000001

long double s21_pow(double base, double exp);

int main() {
    printf("%Lf\n", s21_pow(2, 2.999));
    printf("%f\n", pow(2, 2.999));

    return 0;
}


long double s21_pow(double base, double exp) {
    long double celoe = floor(fabs(exp));
    long double drob = fabs(exp) - celoe;
    long double znamenatel = 1;

    long double returnValue = 1.0;
    long double ValueDrob = 1.0;
    long double startX = 1;
    for (double i = 1; i <= celoe; i++) {
        returnValue = returnValue * base;
    }
    if (drob > eps) {
        long double znamenatel = 1;
        while ((drob) - floor(drob) > eps && (drob) - floor(drob) < 1 - eps){
            drob = drob * 10;
            znamenatel = znamenatel * 10;
        }
        while (fmod(drob, 2) < eps && fmod(znamenatel, 2) < eps) {
            drob = drob / 2;
            znamenatel = znamenatel / 2;
        }
        while (fmod(drob, 5) < eps && fmod(znamenatel, 5) < eps) {
            drob = drob / 5;
            znamenatel = znamenatel / 5;
        }
        long double forBase = s21_pow(base, drob);
        long double startX = forBase/znamenatel;
        long double oldX;
        bool end = false;
        while (!end) {
            oldX = startX;
            startX = 1/znamenatel * ((znamenatel - 1) * startX + forBase/s21_pow(startX, znamenatel - 1));
            if (!(oldX - startX > eps)) {
                end = true;
            }
        }
        returnValue = returnValue * startX;
    }
    if (exp < 0) {
        returnValue = 1/returnValue;
    } 
    if (fabs(exp) < eps) {
        returnValue = 1;
    }
    return returnValue;
}

Не понимаю, как правильно реализовать, если у степени есть десятичная часть. Точнее я сделал это, но при 3х знаках после запятой у степени, там уже получается возведение в 1000 степень и на это тратится секунд 5. При 4х знаках после запятой уже ловлю seg fault.

Answer 1

a^b = e^b·log(a):

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

double s21_pow(double a, double b) {
    return exp(b * log(a));
}

int main() {
    double a, b;
    if (scanf("%lf%lf", &a, &b) == 2) {
        printf("%lf\n", s21_pow(a, b));
    }
}

$ gcc -std=c11 -pedantic -Wall -Wextra -Werror temp.c -lm 

$ echo 100 0.5 | ./a.out 
10.000000

$ echo 100 1 | ./a.out 
100.000000

$ echo 100 2 | ./a.out 
10000.000000

Это решение кажется вам слишком простым? Согласен. Решим задачу без использования <math.h>.

Идею подсказал eri. Я ему благодарен.

Например 1.5^10.625. Переведём показатель в двоичную систему счисления: 10.625₁₀ = 1010.101₂:

 i - номер   значение  1.5^(2^i)
 бита        бита
 3           1          25.62890625   (= 5.0625 * 5.0625)
 2           0           5.0625       (= 2.25 * 2.25)
 1           1           2.25         (= 1.5 * 1.5)
 0           0           1.5
             .
-1           1           1.22474487   (= sqrt(1.5))
-2           0           1.10668192   (= sqrt(1.22474487))   
-3           1           1.05198951   (= sqrt(1.10668192))

Множители в правом столбце получаются последовательным возведением основания в квадрат для положительных бит. Для отрицательных бит надо последовательно извлекать квадратные корни.

Когда таблица составлена, надо перемножить числа соответствующие единичным битам:

1.5^10.625 = 25.62890625 · 2.25 · 1.22474487 · 1.05198951 = 74.29671764.

Осталось сосчитать корень квадратный. Это можно сделать двоичным поиском, если учесть что корень - функция обратная к возведению в квадрат. В итоге получается программа в которой используется сложение, вычитание, умножение и деление. Но деление только на два (что для двоичной вещественной арифметики - особый случай).

Встречайте:

#include <stdio.h>

double sqrt(double a) {
    double low = 0;
    double high = (a < 1) ? 1 : a;
    for (; ; ) {
        double middle = (low + high) / 2;
        if (middle * middle <= a) {
            if (middle == low) {
                return middle;
            }
            low = middle;
        } else {
            if (middle == high) {
                return middle;
            }
            high = middle;
        }
    }
}

typedef struct {
    double b;
    double pow;
} pow_t;

void pow_integer(pow_t *p, double a, double b) {
    if (b <= p->b) {
        pow_integer(p, a * a, 2 * b);
    }
    if (b <= p->b) {
        p->b -= b;
        p->pow *= a;
    }
}

void pow_fractional(pow_t *p, double a, double b) {
    while (p->b > 0) {
        if (b <= p->b) {
            p->b -= b;
            p->pow *= a;
        }
        a = sqrt(a);
        b /= 2;
    }
}

double pow(double a, double b) {
    pow_t p = {b, 1};
    pow_integer(&p, a, 1);
    pow_fractional(&p, a, 1);
    return p.pow;
}

int main() {
    double b, e;
    if (scanf("%lf%lf", &b, &e) == 2) {
        printf("%lf\n", pow(b, e));
    }
}

БЛОГ НА HUSL

Возведение числа в степень Си

Ответы (1 шт):