The simplest way to pow using only simple arithmetic

Question

i want to get function $f(x, a) = x^a$, for both x and a - real numbers, that uses only + - * /. So only way I found is: get taylor series for $$x^a = \sum_{n = 0}^{\infty}\frac{{a}^{n}{log}^{n}(x)}{n!}$$, use taylor series for log and found an answer, but it has complexity of $O(n^2)$ and seems be a bad solution.

Answer 1

So I found my own solution for any a and b. But your post was very helpful. It gave me this idea.

Everybody knows that $${a}^{b} = {e}^{b*ln(a)}$$ We also know that $${a}^{b} = {a}^{n}{a}^{f}, n \in Z, f \in (0..1)$$ Idea is: we can find easily ${a}^{n}$, after we can find $${a}^{f} = {e}^{f*ln(a)}$$ $ln(a)$ can be also easily find by integrating $$\int_{1}^{a}\frac{dt}{t}$$ by Simpson's rule. And exponent can be easily find by taylor series $${e}^{x} = \sum \frac{{x}^{n}}{n!}$$

so here is my C++ realisation:

#include <iostream>
#include <cmath>

using namespace std;


double simpson(double(*f)(double), double a, double b)
{
    double result = (b - a) / 6;
    result *= f(a) + 4 * f((a + b) / 2) + f(b);
    return result;
}

double f(double x)
{
    return 1 / x;
}

double my_ln(double x)
{
    if (x == 1)
        return 0;
    if (x <= 0)
        return NAN;
    double result = 0;
    double step = (x - 1) / 100;
    for (double a = 1; a < x; a += step)
        result += simpson(f, a, a + step);
    return result;
}

double my_exp(double x, double epsilon)
{
    double result = 1;
    double m = x;
    for (int i = 2; abs(m) >= epsilon; i++)
    {
        result += m;
        m *= x / i;
    }
    return result;
}

double my_pow(double a, double b, double epsilon)
{
    if (b < 0)
        return my_pow(1 / a, -b, epsilon);
    int div = (int)b;
    double integerPart = 1;
    for (int i = 0; i < div; i++)
        integerPart *= a;
    return integerPart * my_exp((b - div) * my_ln(a), epsilon);
}

int main()
{   
    cout << my_ln(15) << endl;
    cout << my_exp(3, 0.001) << endl;
    cout << my_pow(2, 3, 0.001) << endl;
}

---

easier way is a bit more compact:

#include <iostream>
#include <cmath>

using namespace std;


double simpson(double(*f)(double), double a, double b)
{
    double result = (b - a) / 6;
    result *= f(a) + 4 * f((a + b) / 2) + f(b);
    return result;
}

double f(double x)
{
    return 1 / x;
}

double my_ln(double x)
{
    if (x == 1)
        return 0;
    if (x <= 0)
        return NAN;
    double result = 0;
    double step = (x - 1) / 100;
    for (double a = 1; a < x; a += step)
        result += simpson(f, a, a + step);
    return result;
}

double my_pow(double a, double b, double epsilon)
{
    double result = 1, tmp = 1;
    double m = b * my_ln(a);
    for (int i = 1; abs(tmp) >= epsilon; i++)
    {
        tmp *= m / i;
        result += tmp;
    }
    return result;
}

int main()
{
    cout << my_ln(15) << endl;
    cout << my_pow(2, 3, 0.001) << endl;
}

from formulas $$x^a = \sum_{n = 0}^{\infty}\frac{{a}^{n}{log}^{n}(x)}{n!}$$ I dunno why I didn't use it first time.

Answer 2

Answering for $a\in\mathbb{Q}$ (the set of rational numbers):

• Let $\displaystyle a=\frac{N}{D}$

• $\displaystyle x^a=\sqrt[D]{x^N}$

• $\displaystyle x^N$ can be calculated easily

• $\displaystyle \sqrt[D]{x^N}$ can be calculated with the following (converging) sequence:

  • $\displaystyle S_0=x^N$

  • $\displaystyle S_{n+1}=S_n-\frac{{S_n}^D-x^N}{D\cdot{S_n}^{D-1}}$

---

Here is a method for calculating the numerator and the denominator of a rational value:

Function (Input value, Output numerator, Output denominator):

    Set numerator   = 0
    Set denominator = 1

    If value > 0:
        Set sign = 1
    Else:
        Set sign  = -1
        Set value = -value

    While value > 0:
        Set intValue    = Int(value)
        Set value       = value - intValue
        Set numerator   = numerator + intValue
        Set value       = value * 2
        Set numerator   = numerator * 2
        Set denominator = denominator * 2

    If numerator > denominator:
        Set hi = numerator
        Set lo = denominator
    Else:
        Set hi = denominator
        Set lo = numerator

    While True:
        Set rem = hi mod lo
        If rem == 0:
            Break
        Set hi = lo
        Set lo = rem

    Set numerator   = numerator   / lo
    Set denominator = denominator / lo
    Set denominator = denominator * sign

---

In addition, you might find this C++ implementation useful.