Groups | Search | Server Info | Keyboard shortcuts | Login | Register [http] [https] [nntp] [nntps]

Groups > comp.lang.python > #111467 > unrolled thread

math.frexp

Back to article view | Back to comp.lang.python

Contents

  math.frexp Steven D'Aprano <steve@pearwood.info> - 2016-07-15 21:39 +1000
    Re: math.frexp Chris Angelico <rosuav@gmail.com> - 2016-07-15 21:48 +1000
      Re: math.frexp Steven D'Aprano <steve@pearwood.info> - 2016-07-16 02:32 +1000
        Re: math.frexp Chris Angelico <rosuav@gmail.com> - 2016-07-16 02:44 +1000
        Re: math.frexp Random832 <random832@fastmail.com> - 2016-07-15 14:28 -0400
        Re: math.frexp Paul Rubin <no.email@nospam.invalid> - 2016-07-15 13:24 -0700
          Re: math.frexp Steven D'Aprano <steve@pearwood.info> - 2016-07-16 20:21 +1000
            Re: math.frexp Paul Rubin <no.email@nospam.invalid> - 2016-07-16 12:49 -0700
    Re: math.frexp Nobody <nobody@nowhere.invalid> - 2016-07-15 17:39 +0100
    Re: math.frexp Vlastimil Brom <vlastimil.brom@gmail.com> - 2016-07-17 01:18 +0200

#111467 — math.frexp

I'm experimenting with various implementations of product(values). Here's a
naive version that easily suffers from overflow:

from functools import reduce
from operator import mul

def product_naive(values):
    return reduce(mul, values)


py> product_naive([2, 4, 5])
40
py> product_naive([1e200, 1e200, 1e200])
inf


So I started experimenting with math.frexp, but I'm having trouble
understanding what I'm doing wrong here.


Help on built-in function frexp in module math:

frexp(...)
    frexp(x)

    Return the mantissa and exponent of x, as pair (m, e).
    m is a float and e is an int, such that x = m * 2.**e.
    If x is 0, m and e are both 0.  Else 0.5 <= abs(m) < 1.0.



If x and y are both floats, then given:

m1, e1 = math.frexp(x)
m2, e2 = math.frexp(y)

Then x*y = m1*m2 * 2.0**(e1 + e2). We can test that:


py> from math import frexp
py> x = 2.5
py> y = 3.5
py> x*y
8.75
py> m1, e1 = math.frexp(x)
py> m2, e2 = math.frexp(y)
py> m1*m2 * 2.0**(e1 + e2)
8.75


Looks good to me. So let's try a less naive version of product():


def product_scaled(values):
    scale = 0
    prod = 1.0
    for a in values:
        m1, e1 = math.frexp(a)
        m2, e2 = math.frexp(prod)
        scale += (e1 + e2)
        prod *= (m1*m2)
    return (prod * 2.0**scale)


py> product_scaled([2.5, 3.5])  # expected 8.75
2.734375



What am I doing wrong?




-- 
Steven
“Cheer up,” they said, “things could be worse.” So I cheered up, and sure
enough, things got worse.

[toc] | [next] | [standalone]

#111469

On Fri, Jul 15, 2016 at 9:39 PM, Steven D'Aprano <steve@pearwood.info> wrote:
> py> from math import frexp
> py> x = 2.5
> py> y = 3.5
> py> x*y
> 8.75
> py> m1, e1 = math.frexp(x)
> py> m2, e2 = math.frexp(y)
> py> m1*m2 * 2.0**(e1 + e2)
> 8.75
>
>
> Looks good to me. So let's try a less naive version of product():
>
>
> def product_scaled(values):
>     scale = 0
>     prod = 1.0
>     for a in values:
>         m1, e1 = math.frexp(a)
>         m2, e2 = math.frexp(prod)
>         scale += (e1 + e2)
>         prod *= (m1*m2)
>     return (prod * 2.0**scale)
>
>
> py> product_scaled([2.5, 3.5])  # expected 8.75
> 2.734375
>

You're chaining your product twice. (Also your scale, although that
appears to be correct.) Changing it to "prod = m1 * m2" gives 8.75.

But what do you gain by this? You're still stuffing the result back
into a float at the end, so all you do is change from getting
float("inf") to getting OverflowError. How can you make it not
overflow?

ChrisA

[toc] | [prev] | [next] | [standalone]

#111481

On Fri, 15 Jul 2016 09:48 pm, Chris Angelico wrote:

>> py> product_scaled([2.5, 3.5])  # expected 8.75
>> 2.734375
>>
> 
> You're chaining your product twice. (Also your scale, although that
> appears to be correct.) Changing it to "prod = m1 * m2" gives 8.75.

D'oh!

Thanks! I needed some fresh eyes on that, because I have been staring at it
for two days.


> But what do you gain by this? You're still stuffing the result back
> into a float at the end, so all you do is change from getting
> float("inf") to getting OverflowError. How can you make it not
> overflow?

If the result is too big to be represented as a float at the end of the
product, then of course it will overflow. But this can give some protection
against overflow of intermediate values. Consider multiplying:

2.0, 1e200, 1e200, 1e-200, 1e-200, 3.0


Mathematically, the answer should be 6. In principle, by rescaling when
needed to prevent overflow (or underflow), product() should be able to get
something very close to 6, if not exactly 6.

But I'm not actually writing a general product() function. I'm doing this
for geometric mean, so I return the scaling exponent and the mantissa[1]
separately, and then take the nth-root of them individually, before
combining them into the final result. Here's an example from my test suite.
I can take the geometric mean of: 

[2.0**900, 2.0**920, 2.0**960, 2.0**980, 2.0**990]

and get to within a relative error of 1e-14 of the correct answer, 2**950:

py> geometric_mean([2.0**900, 2.0**920, 2.0**960, 2.0**980, 2.0**990])
9.516908214257811e+285
py> 2.0**950
9.516908214257812e+285


Here's one that's exact:

py> geometric_mean([2e300, 4e300, 8e300, 16e300, 32e300])
8e+300






[1] I don't care what it's actually called :-)


-- 
Steven
“Cheer up,” they said, “things could be worse.” So I cheered up, and sure
enough, things got worse.

[toc] | [prev] | [next] | [standalone]

#111483

On Sat, Jul 16, 2016 at 2:32 AM, Steven D'Aprano <steve@pearwood.info> wrote:
> If the result is too big to be represented as a float at the end of the
> product, then of course it will overflow. But this can give some protection
> against overflow of intermediate values. Consider multiplying:
>
> 2.0, 1e200, 1e200, 1e-200, 1e-200, 3.0
>
>
> Mathematically, the answer should be 6. In principle, by rescaling when
> needed to prevent overflow (or underflow), product() should be able to get
> something very close to 6, if not exactly 6.

I was under the impression that appropriate reordering of the elements
could prevent over/underflow and maximize accuracy, but that may be a
mistaken memory. However...

> But I'm not actually writing a general product() function. I'm doing this
> for geometric mean, so I return the scaling exponent and the mantissa[1]
> separately, and then take the nth-root of them individually, before
> combining them into the final result.

... this makes a lot of sense. In effect, I *think*, you're basically
doing the multiplication on something rather larger than a 64-bit IEEE
float, then taking the nth-root, and then combine them and convert
back to float. Is this about right?

ChrisA

[toc] | [prev] | [next] | [standalone]

#111484

On Fri, Jul 15, 2016, at 12:32, Steven D'Aprano wrote:
> I can take the geometric mean of: 
> 
> [2.0**900, 2.0**920, 2.0**960, 2.0**980, 2.0**990]
> 
> and get to within a relative error of 1e-14 of the correct answer,
> 2**950:
> 
> py> geometric_mean([2.0**900, 2.0**920, 2.0**960, 2.0**980, 2.0**990])
> 9.516908214257811e+285
> py> 2.0**950
> 9.516908214257812e+285
>
> Here's one that's exact:
> 
> py> geometric_mean([2e300, 4e300, 8e300, 16e300, 32e300])
> 8e+300

My own attempt gave an exact result for both test cases.

from math import frexp as _frexp, ldexp

def frexp(x):
    if type(x) is tuple:
        m, e1 = x
        m, e2 = _frexp(m)
        return m, e1 + e2
    else:
        return _frexp(x)

def _mul(a, b):
    ma, ea = frexp(a)
    mb, eb = frexp(b)
    return frexp((ma * mb, ea + eb))

def _product(values):
    return reduce(_mul, values)

def geometric_mean(values):
    n = len(values)
    pm, pe = _product(values)
    me, re = divmod(pe, n)
    pm = ldexp(pm, re)
    return ldexp(pm**(1./n), me)

[toc] | [prev] | [next] | [standalone]

#111487

Steven D'Aprano <steve@pearwood.info> writes:
> But this can give some protection against overflow of intermediate
> values. 

Might be simplest to just add the logarithms.  Look up Kahan summation
for how to do that while minimizing loss of precision.

[toc] | [prev] | [next] | [standalone]

#111506

On Sat, 16 Jul 2016 06:24 am, Paul Rubin wrote:

> Steven D'Aprano <steve@pearwood.info> writes:
>> But this can give some protection against overflow of intermediate
>> values.
> 
> Might be simplest to just add the logarithms.  Look up Kahan summation
> for how to do that while minimizing loss of precision.

Simplest, but least accurate, even with Kahan summation or equivalent. Even
a naive implementation of product does better:

py> from operator import mul
py> reduce(mul, [1.0, 2.0, 3.0, 4.0, 5.0])
120.0
py> math.exp(math.fsum(math.log(x) for x in [1.0, 2.0, 3.0, 4.0, 5.0]))
119.99999999999997


That second answer might be good enough for getting an astronaut to the
Moon, but it's not good enough to get them back again *wink*



-- 
Steven
“Cheer up,” they said, “things could be worse.” So I cheered up, and sure
enough, things got worse.

[toc] | [prev] | [next] | [standalone]

#111521

Steven D'Aprano <steve@pearwood.info> writes:
>>> protection against overflow of intermediate values.
>> Might be simplest to just add the logarithms.
> Simplest, but least accurate, even with Kahan summation or equivalent.

Well the idea was to avoid overflow first, then hold on to whatever
precision you have after that.  Is there a way to use extended
(e.g. 80-bit) floats in numpy/scipy/etc.?

[toc] | [prev] | [next] | [standalone]

#111482

On Fri, 15 Jul 2016 21:39:31 +1000, Steven D'Aprano wrote:

>         prod *= (m1*m2)

Should be:

	prod = m1*m2

or:
	prod *= m1

(in the latter case, there's no point in decomposing prod).

Of course, if the result overflows, it's going to overflow whether you use
the naive approach or frexp(); in the latter case, it's the 2.0**scale
which will overflow.

One advantage of processing the scale separately is that you can use e.g.

    return ((int(prod * 2**53) * (2**(scale-53))) if scale >= 53
            else prod * 2**scale)

which will return a (long) integer if the exponent is such that the
fractional bits are lost.

[toc] | [prev] | [next] | [standalone]

#111529

2016-07-15 13:39 GMT+02:00 Steven D'Aprano <steve@pearwood.info>:
> I'm experimenting with various implementations of product(values). Here's a
> naive version that easily suffers from overflow:
>
> from functools import reduce
> from operator import mul
>
> def product_naive(values):
>     return reduce(mul, values)
>...
> --
> Steven
> --
> https://mail.python.org/mailman/listinfo/python-list

Hi,
just to add another approach to already proposed ones, I thought about
using fractions for computations instead of floats where possible.
The accuracy shoud be better, - it works fine for my tests with
multiplying however, the downside is - for computing fractional power
the fractions fall back to float, as irrational numbers are expected
in the output, hence the float overflow is also possible (unless an
individual implementation of the root computation would be used in
gmean...).

The hackish code in adapt_float_exp enhances the accuracy compared to
directly building Fractions from large floats in exponential notation.
(I can realise, that this approach might be insufficient for the scale
you intend to support, due to the mentioned fallback to floats.)

Regards,
    vbr

========================

#! Python
# -*- coding: utf-8 -*-

from fractions import Fraction

def adapt_float_exp(flt):
    """Return Fraction extracted from float literal in exponential notation."""
    if "e" not in str(flt).lower():
        return Fraction(flt)
    else:
        mant_str, exp_str = str(flt).lower().split("e")
        return Fraction(mant_str) * 10 ** Fraction(exp_str)

def prod_fract(*args):
    prod = 1
    for arg in args:
        if isinstance(arg, float):
            arg = adapt_float_exp(arg)
        prod *= arg
    return prod

def gmean(*args):
    prod = prod_fract(*args)
    root = Fraction(prod) ** Fraction(1, len(args)) # fractional power
** - irrational - implemented via floats !! - overflow ...
    return root


print(repr(prod_fract(2.0, 1e200, 1e200, 1e-200, 1e-200, 3.0)))
print(repr(gmean(2.0, 1e200, 1e200, 1e-200, 1e-200, 3.0)))

[toc] | [prev] | [standalone]

Back to top | Article view | comp.lang.python

csiph-web