5
$\begingroup$

A MATLAB library seems to overcomplicate a computation:

exp( (log(a) - log(b))/b )

which is mathematically equivalent (assuming real & positive a, b) to

(a/b).^(1/b)

Another computation in log-exp cancellation makes sense for greater numeric stability, but I can't see any advantages here. What could be the motivation?


I work with Python, that implements ** as described here. I wrote a test script to compare direct ** vs exp; results are, they're very close for a**b vs a**(1/b), thus it seems best to simply choose former for readability.

enter image description here

import numpy as np
import matplotlib.pyplot as plt

def viz(h, a, out1, out2, N, reciprocal=False):
    tkw = dict(weight='bold', fontsize=16)
    if reciprocal:
        plt.title("a**(1/b) / exp((1/b)*log(a))", **tkw)
    else:
        plt.title("a**b / exp(b*log(a))", **tkw)

    plt.hist(h.ravel(), bins=500)
    annot = ("min={}\nmax={}\n** (nans, infs) = ({}, {})\nexp (nans, infs) = "
             "({}, {})").format(h.min(), h.max(),
                                np.isnan(out1).sum(), np.isinf(out1).sum(),
                                np.isnan(out2).sum(), np.isinf(out2).sum())
    plt.annotate(annot, xycoords='axes fraction', xy=(.65, .8),
                 weight='bold', fontsize=14)
    plt.show()


def test_exp(a, b, reciprocal=False):
    a, b = a.reshape(-1, 1), b.reshape(-1, 1)

    out1 = a ** b.T
    out2 = np.exp(np.log(a) * b.T)
    h = out2 / out1
    h[np.isnan(h) | np.isinf(h)] = 1

    viz(h, a, out1, out2, len(a), reciprocal)

a = np.linspace(1e-6, 40, 2000)
b = np.linspace(1e-6, 40, 2000)
test_exp(a, b)
test_exp(a, 1/b, reciprocal=True)
$\endgroup$
2
  • $\begingroup$ I think that ‘(a/b)^(1/b)‘ is implemented exactly using that logarithmic trick. See also stackoverflow.com/questions/1375953/… and related questions @OverLordGoldDragon $\endgroup$ – VoB Jan 19 at 9:24
  • 1
    $\begingroup$ Maybe this was ported from a language that did not have a working pow(float, float) without thinking too much? $\endgroup$ – Federico Poloni Jan 19 at 15:52
4
$\begingroup$

While $a^{b}$ is usually computed via $\exp(b \log (a))$, in real-life high-quality implementations of functions like pow() the logarithm is computed with extra precision to counter the error magnification effect of $\exp$, as demonstrated in this answer for example. This suggests that we should favor computation via a built-in general exponentiation operator or function like pow(). In many cases this will give more accurate results, and even for implementations of general exponentiation that are not state of the art it is unlikely to give results that are worse. For a quick experiment along those lines, see the ISO-C99 code below. On my platform, it prints:

max ulp error using expf(): 3.72998984e+005
max ulp error using powf(): 2.16050157e+004

Even so, the computation of reciprocal powers via general exponentiation likewise suffers from error magnification of the error incurred in the computation of $\frac{1}{b}$, so where specific functions are available to compute, say, the cube root (i.e., $b=3$), it is advisable to use those, as discussed in this answer. Some programming environments offer a rootn() function for the computation of the $n$-th root, for example.

#include <stdio.h>
#include <stdlib.h>
#include <stdint.h>
#include <string.h>
#include <math.h>

#define N  (10000000)

float func (float a, float b)
{
    return expf ((logf (a) - logf (b)) / b);
}

float func_alt (float a, float b)
{
    return powf (a / b, 1.0f / b);
}

double func_ref (double a, double b)
{
    return pow (a / b, 1.0 / b);
}

uint32_t float_as_uint32 (float a)
{
    uint32_t r;
    memcpy (&r, &a, sizeof r);
    return r;
}

float uint32_as_float (uint32_t a)
{
    float r;
    memcpy (&r, &a, sizeof r);
    return r;
}

uint64_t double_as_uint64 (double a)
{
    uint64_t r;
    memcpy (&r, &a, sizeof r);
    return r;
}

double floatUlpErr (float res, double ref)
{
    uint64_t i, j, err, refi;
    int expoRef;
    
    /* ulp error cannot be computed if either operand is NaN, infinity, zero */
    if (isnan (res) || isnan (ref) || isinf (res) || isinf (ref) ||
        (res == 0.0f) || (ref == 0.0f)) {
        return 0.0;
    }
    /* Convert the float result to an "extended float". This is like a float
       with 56 instead of 24 effective mantissa bits.
    */
    i = ((uint64_t)float_as_uint32(res)) << 32;
    /* Convert the double reference to an "extended float". If the reference is
       >= 2^129, we need to clamp to the maximum "extended float". If reference
       is < 2^-126, we need to denormalize because of the float types's limited
       exponent range.
    */
    refi = double_as_uint64(ref);
    expoRef = (int)(((refi >> 52) & 0x7ff) - 1023);
    if (expoRef >= 129) {
        j = 0x7fffffffffffffffULL;
    } else if (expoRef < -126) {
        j = ((refi << 11) | 0x8000000000000000ULL) >> 8;
        j = j >> (-(expoRef + 126));
    } else {
        j = ((refi << 11) & 0x7fffffffffffffffULL) >> 8;
        j = j | ((uint64_t)(expoRef + 127) << 55);
    }
    j = j | (refi & 0x8000000000000000ULL);
    err = (i < j) ? (j - i) : (i - j);
    return err / 4294967296.0;
}

// Fixes via: Greg Rose, KISS: A Bit Too Simple. http://eprint.iacr.org/2011/007
static unsigned int z=362436069,w=521288629,jsr=362436069,jcong=123456789;
#define znew (z=36969*(z&0xffff)+(z>>16))
#define wnew (w=18000*(w&0xffff)+(w>>16))
#define MWC  ((znew<<16)+wnew)
#define SHR3 (jsr^=(jsr<<13),jsr^=(jsr>>17),jsr^=(jsr<<5)) /* 2^32-1 */
#define CONG (jcong=69069*jcong+13579)                     /* 2^32 */
#define KISS ((MWC^CONG)+SHR3)

int main (void)
{
    float a, b, res1, res2;
    double ref, ulp1, ulp2, maxulp1 = 0, maxulp2 = 0;

    for (int i = 0; i < N; i++) {
        do { 
            a = uint32_as_float (KISS & 0x7fffffff); // ensure positive
        } while (!isnormal (a));
        do {
            b = uint32_as_float (KISS & 0x7fffffff); // ensure positive
        } while (!isnormal (a));
        res1 = func (a, b);
        res2 = func_alt (a, b);
        ref = func_ref ((double)a, (double)b);
        ulp1 = floatUlpErr (res1, ref);
        ulp2 = floatUlpErr (res2, ref);
        maxulp1 =  fmax (ulp1, maxulp1);
        maxulp2 =  fmax (ulp2, maxulp2);
    }
    printf ("max ulp error using expf(): %15.8e\n", maxulp1);
    printf ("max ulp error using powf(): %15.8e\n", maxulp2);
    return EXIT_SUCCESS;
}
```
$\endgroup$
4
$\begingroup$

It's not overcomplicated: That's how $x^y$ is actually defined -- via $$ x^y = \exp(y \log x). $$ The point is that for non-integer $y$, it's not at all obvious what $x^y$ actually means: It is not just a product of $x$ with itself, repeated $y$ times. The definition of $x^y$ then is done by pulling the expression back to functions we have previously defined, namely the exponential and logarithmic functions, which are in terms defined via Taylor series.

$\endgroup$
4
  • $\begingroup$ Interesting, though it does considerably degrade readability; is there any numeric stability or precision advantage from this? I code it in Python, which computes via Numpy and C, so the answer lies there, but is there even an 'answer' to search $\endgroup$ – OverLordGoldDragon Jan 19 at 17:21
  • 1
    $\begingroup$ What I'm saying is that whether you write $x^y$ or $\exp(y \ln x)$ doesn't matter if $y$ is a floating point number. If you write the former, then any compiler will have to translate it to the latter internally because that's how this is defined. So there is no difference in stability or precision. $\endgroup$ – Wolfgang Bangerth Jan 19 at 18:59
  • $\begingroup$ That's maybe not quite correct. One could interpret $x^y$ as $x^{\lfloor y \rfloor} x^{(y - \lfloor y \rfloor})$ and only convert the latter part. I don't know whether compilers do this transformation. $\endgroup$ – Wolfgang Bangerth Jan 19 at 19:01
  • $\begingroup$ Opened one here to be sure. $\endgroup$ – OverLordGoldDragon Jan 19 at 21:12

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.