Is there a way to calculate decimal powers using only addition, subtraction, multiplication and division?

Question

I am a programmer who is trying to build an arbitrary precision decimal class library for C# and was successful in implementing all the basic operations like addition, subtraction, multiplication and division but was stopped when trying to calculate powers to decimal (and negative decimal) exponents. Calculating Log and Ln was also not possible since the depend on calculating powers.

So here is what I CAN do (with arbitrary accuracy):

• Calculate the result of ALL operators if both numbers are integers.

• Calculate the result of addition, subtraction, multiplication and division if the numbers are mixed (decimals and integers).

• Calculate the result of the power of a decimal if the exponent is integer.

Here is what I CAN'T do (and I hope to find the answer here):

• Calculate the power of decimals or integers to fractional, negative or decimal exponents.

Some (programmers) might point out that I can use functions already in C# to calculate these but this would not work as I need to build the library to be of arbitrary precision so I have to do the calculation manually.

So, is there a way to do it with the current "tools" I have? Iterative methods are not a problem since I will not be calculating the powers by hand, the CPU will be doing the counting.

Answer 1

The simplest way would be to implement log and exp via Taylor series (using the Taylor remainder theorem in order to bound the error to the desired precision) and then just rewrite

$$a^b = \exp(b \log a)$$

However, this is probably not the ideal way of doing things. For instance, you could use the Sasaki-Kanada algorithm to obtain the logarithm.

A lot of this comes down to the kind of tradeoffs you're trying to make. Do you just want something that kind of functions, or do you want to be able to make guarantees about the accuracy of your calculations? Do you need speed, or can it take all the time it wants? Is there a solid reason not to use existing libraries for arbitrary-precision arithmetic, and if so what is it?

There's an entire field of mathematics studying these questions, namely numerical analysis, and many of these problems have been solved (or at least have best-practice solutions). But there's often no simple, one-size-fits-all answer in a lot of cases.

To maybe illustrate a bit of the complexity here: given numbers $a$ and $b$, we know that $a^b = e^{b \log a}$. But if you calculate $\log a$ to a particular precision, and then calculate the product with $b$ to the same precision, and finally exponentiate that to the same precision, then you will not have calculated $a^b$ to that precision -- every time you do an arithmetic calculation you lose a bit of accuracy, and this can compound dramatically.

Answer 2

Here is what you can do, roughly:

First, we are trying to implement the exponential function $\exp$, which is given by the power series:

$$\exp x =\sum_{k=0}^\infty \frac{x^n}{n!}$$

We are going to use a few of the terms of the series to approximate $\exp x$, say the first four ones:

$$\exp x \approx 1 + x + \frac{x^2}{2} + \frac{x^3}{6}$$

This formula works well for small $|x|$ and increasingly badly for large $|x|$. That is why we reduce all the calculation of $\exp x$ to those with small $|x|$: We have $\exp x = \exp(x\frac{k}{k}) = \exp(\frac{x}{k})^k$. You choose (somehow) a large enough $k$ and calculate $\exp x$ that way, where $\exp(\frac{x}{k})$ is calculated by the aforementoined formula.

Next, observe, that $a^b = \exp( \log a \cdot b)$, if it is defined. Calculating $\log a$ is a little bit more complicated (afaik), you can read up on that somewhere else too. It get's very difficult for $a$ close to $0$ (because $\log a$ goes to $-\infty$). Here is an idea on how to do it: We have the power series:

$$\log(x + 1) = -\sum_{k=1}^\infty \frac{(-1)^k\cdot x^k}{k}$$

which works for $|x|<1$. We use a few terms:

$$\log(x + 1) \approx x - \frac{x^2}{2} + \frac{x^3}{3}$$

(for $|x|<1$)

If $a>1$, then we want to calculate $\log(x+1)$ somehow for an $x$ smaller then $1$, but better close to $0$. Pick the largest $k$ with $\frac{a}{2^k} > 1$. Consider: $$\log a = \log(\frac{a}{2^k}2^k) = \log(x+1)+\log(2^k) = \log(x+1) + k\log 2$$ for $x = \frac{a}{2^k} - 1$. Now you can calculate $\log(x+1)$ with the above formula and remember the value of $\log 2$.

For $a<1$, I really don't know what to do.

Answer 3

Suppose you have a class RationalNum which implements all the other operations.

Suppose that every RationalNum instance has the following attributes:

• bool m_bSign
• NaturalNum m_cEnumerator
• NaturalNum m_cDenominator

Where class NaturalNum implements all the relevant arithmetic operations.

You can use Newton-Raphson method in order to calculate the $n$th root as follows:

RationalNum RationalNum::Root(unsigned int iExponent,unsigned int iAccuracy) const
{
    if (iExponent == 0)
        throw "Infinite Result";
    if (m_bSign && (iExponent%2) == 0)
        throw "Undefined Result";

    if (m_cEnumerator == 0)
        return 0;

    RationalNum cRes;

    cRes.m_bSign = m_bSign;
    cRes.m_cEnumerator  = 1;
    cRes.m_cDenominator = 1;

    cRes.m_cEnumerator  <<= m_cEnumerator.BitCount()/m_cDenominator.BitCount()/iExponent;
    cRes.m_cDenominator <<= m_cDenominator.BitCount()/m_cEnumerator.BitCount()/iExponent;

    for (unsigned int i=0; i<iAccuracy; i++)
    {
        RationalNum cTmp = cRes.Power(iExponent-1);
        cRes -= (cRes*cTmp-*this)/(iExponent*cTmp);
    }

    return cRes;
}