Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Add C99's log2() function to the math library #56097

Closed
rhettinger opened this issue Apr 20, 2011 · 35 comments
Closed

Add C99's log2() function to the math library #56097

rhettinger opened this issue Apr 20, 2011 · 35 comments
Assignees
Labels
type-feature A feature request or enhancement

Comments

@rhettinger
Copy link
Contributor

BPO 11888
Nosy @rhettinger, @jcea, @mdickinson, @vstinner
Files
  • issue11888.patch
  • issue11888-part2.patch
  • Note: these values reflect the state of the issue at the time it was migrated and might not reflect the current state.

    Show more details

    GitHub fields:

    assignee = 'https://github.com/mdickinson'
    closed_at = <Date 2011-05-11.20:08:22.426>
    created_at = <Date 2011-04-20.15:58:48.806>
    labels = ['type-feature']
    title = "Add C99's log2() function to the math library"
    updated_at = <Date 2011-05-12.07:29:45.446>
    user = 'https://github.com/rhettinger'

    bugs.python.org fields:

    activity = <Date 2011-05-12.07:29:45.446>
    actor = 'mark.dickinson'
    assignee = 'mark.dickinson'
    closed = True
    closed_date = <Date 2011-05-11.20:08:22.426>
    closer = 'vstinner'
    components = []
    creation = <Date 2011-04-20.15:58:48.806>
    creator = 'rhettinger'
    dependencies = []
    files = ['21861', '21918']
    hgrepos = []
    issue_num = 11888
    keywords = ['patch']
    message_count = 35.0
    messages = ['134159', '134197', '134198', '134200', '134202', '134203', '134427', '134428', '134946', '134998', '135001', '135201', '135352', '135426', '135427', '135428', '135549', '135569', '135572', '135575', '135576', '135577', '135578', '135579', '135583', '135585', '135586', '135733', '135734', '135735', '135737', '135742', '135748', '135800', '135822']
    nosy_count = 6.0
    nosy_names = ['rhettinger', 'jcea', 'mark.dickinson', 'vstinner', 'rpetrov', 'python-dev']
    pr_nums = []
    priority = 'normal'
    resolution = 'fixed'
    stage = None
    status = 'closed'
    superseder = None
    type = 'enhancement'
    url = 'https://bugs.python.org/issue11888'
    versions = ['Python 3.3']

    @rhettinger
    Copy link
    Contributor Author

    The three most popular logarithm bases are 10, e, and 2. The math library has direct function calls for the first two but not the latter which is important in informatics.

    Since a direct call can use a custom algorithm or native hardware support (such as the FLDLN2 fpu instruction), it provides better speed and accuracy than our existing math.log(x, 2) option.

    @rhettinger rhettinger added the type-feature A feature request or enhancement label Apr 20, 2011
    @mdickinson
    Copy link
    Member

    See also bpo-3724.

    I'm -0 on this: between log(x, 2) and int.bit_length, there's not much need for log2. log(x, 2) should be plenty accurate enough for most numerical needs; the exception is when you're taking log base 2 of an integer and need a guarantee of exact results for powers of 2, and int.bit_length generally solves that problem.

    The main issue is that we'd have to provide (and maintain) our own implementation of log2 for Windows (and other OSs that don't have all the C99 support. Solaris?) That implementation should, ideally:

    • provide exact values for powers of 2, and
    • be monotonic.

    and that's not trivial. As Raymond points out, on x86 / x64 we might be able to use inline assembly directly; that would probably cover us for Windows.

    @vstinner
    Copy link
    Member

    The main issue is that we'd have to provide (and maintain) our own
    implementation of log2 for Windows (and other OSs that don't have all
    the C99 support. Solaris?)

    No, we don't have to. Python has already a lot of optional functions, see for example the os module. We can provide log2() only if the C library has this function.

    @mdickinson
    Copy link
    Member

    We can provide log2() only if the C library has this function.

    Big -1 from me: I'd hate to see working Python scripts written on Unix fail on Windows because of a missing log2.

    @mdickinson
    Copy link
    Member

    Rather than reinventing the wheel, it may be worth looking at what numpy does here.

    @mdickinson
    Copy link
    Member

    it may be worth looking at what numpy does here.

    ... or it may not. NumPy just uses (approximation to 1/log(2)) * log(x) when log2 doesn't already exist. And indeed, on Windows:

    Python 2.7.1 |EPD 7.0-2 (64-bit)| (r271:86832, Dec  2 2010, 10:23:25) [MSC v.1500 64 bit (AMD64)] on win32
    Type "help", "copyright", "credits" or "license" for more information.
    >>> import numpy
    >>> numpy.log2(8.0)
    2.9999999999999996

    I think we should be able to do better than this. :-)

    @vstinner
    Copy link
    Member

    The main issue is that we'd have to provide (and maintain) our own
    implementation of log2 for Windows (and other OSs that don't have all
    the C99 support. Solaris?)

    Can't we simply use (approximation to 1/log(2)) * log(x)? Is it worse than reimplementing it using log(x)/log(2) in Python?

    That implementation should, ideally: ...

    Because some platforms are less accurate, you prefer to not provide a more accurate function when log2() is available? Can't we start with something simple and improve it later?

    It can be documented than the Python log2 function may be the same than log(x)/log(2) if the platform/CPU doesn't provide a C log2 function.

    --

    "... And indeed, on Windows ...

    >>> numpy.log2(8.0)
    2.9999999999999996"

    Oh. Python is better on Linux:

    Python 2.6.6 (r266:84292, Dec 26 2010, 22:31:48) 
    [GCC 4.4.5] on linux2
    Type "help", "copyright", "credits" or "license" for more information.
    >>> import math
    >>> math.log(8) / math.log(2)
    3.0

    @vstinner
    Copy link
    Member

    Can't we simply use (approximation to 1/log(2)) * log(x)?
    Is it worse than reimplementing it using log(x)/log(2) in Python?

    Hum. With a x86 and the right compiler optimization level, log(x)/log(2) in C can be more accurate than log(x)/log(2) in Python, because the FPU works with 80 bits float internally, and the result is only "truncated" to 64 bits float at the end. In Python, the result is truncated to 64 bits on each Python instruction.

    I don't know if it should be called a feature or a bug. In PHP world, it would be called a bug :-D http://bugs.php.net/bug.php?id=53632

    @vstinner
    Copy link
    Member

    vstinner commented May 1, 2011

    Oh... math.log() has an optional second argument: base. math.log(x, 2). But it is equivalent as math.log(x) / math.log(2) in Python. math.log(x, 2) is implemented as:
    num=math.log(x)
    den=math.log(2)
    return num / den
    where num and den are Python floats (64 bits).

    So we don't benefit from 80 bits float used internally in x87.

    @mdickinson
    Copy link
    Member

    Here's a patch implementing log2. Still to do: use the system log2 where available.

    @rhettinger
    Copy link
    Contributor Author

    Wow Mark, that is really nice work. Thanks.

    @vstinner
    Copy link
    Member

    vstinner commented May 5, 2011

    Updated patch to use the system log2() if it is available. The test pass with the system log2() on Linux (Debian Sid, eglibc 2.11.2).

    @mdickinson
    Copy link
    Member

    Thanks, Victor. I suspect we're going to need to be a bit more careful, though: when the extra tests were added for math.log, it turned out that it had all sorts of strange special-case behaviour on various platforms.

    So I suspect that even on platforms that have log2 natively, it'll be necessary to factor out special cases and deal with those first, only passing positive finite floats onto the system log2. Take a look at m_log and the comment directly above it to see how that works.

    I'd also like to check in the non-system version first, just to give it a thorough test on the buildbots, before adding in the version that uses the system log2 when available.

    @vstinner
    Copy link
    Member

    vstinner commented May 7, 2011

    Thanks, Victor. I suspect we're going to need to be a bit more
    careful, though: when the extra tests were added for math.log, it
    turned out that it had all sorts of strange special-case behaviour on
    various platforms.

    So I suspect that even on platforms that have log2 natively, it'll be
    necessary to factor out special cases and deal with those first, only
    passing positive finite floats onto the system log2. Take a look at
    m_log and the comment directly above it to see how that works.

    Oh, I see:

    /*
    Various platforms (Solaris, OpenBSD) do nonstandard things for
    log(0),
    log(-ve), log(NaN). Here are wrappers for log and log10 that deal
    with
    special values directly, passing positive non-special values through
    to
    the system
    log/log10.
    */

    I'd also like to check in the non-system version first, just to give
    it a thorough test on the buildbots, before adding in the version that
    uses the system log2 when available.

    Yes, we can use log2() only for the "x > 0.0" case. My secret plan was
    to check system log2() using the buildbots. But if we know that system
    log() only works correctly with strictly positive numbers, it's faster
    to directly only use system log2() for x > 0.0.

    Updated patch (version 3) implements that.

    @vstinner
    Copy link
    Member

    vstinner commented May 7, 2011

    (Oh, I hit the wrong keyboard shortcut and sent my email too fast)

    You can commit bpo-11888.patch first if you would like to test it.

    In this case, here is a patch to use system log2(), patch to apply *after* bpo-11888.patch. It only uses log2() if x > 0.0.

    @vstinner
    Copy link
    Member

    vstinner commented May 7, 2011

    By the way, bpo-11888.patch is just fine: you can commit it. I like your frexp "trick" to improve the accuracy.

    @python-dev
    Copy link
    Mannequin

    python-dev mannequin commented May 8, 2011

    New changeset 6d1cbfcee45a by Victor Stinner in branch 'default':
    Issue bpo-11888: Add log2 function to math module. Patch written by Mark
    http://hg.python.org/cpython/rev/6d1cbfcee45a

    @mdickinson
    Copy link
    Member

    Thanks, Victor. You caught me by surprise a bit: I had some more minor changes to that patch pending, so I've committed those separately.

    Any news from the buildbots?

    @vstinner
    Copy link
    Member

    vstinner commented May 9, 2011

    Thanks, Victor. You caught me by surprise a bit

    Oh, I thought that the patch was ready to be commited.

    I had some more minor changes to that patch pending,
    so I've committed those separately.

    You should add "Issue bpo-11888: " prefix to your commit messages so a bot automatically add comments for the new commits to this issue.

    "Fix cut-and-paste typo in comment: log10 -> log2."

    Oh, I didn't notice it, but Terry Reedy

    "Fix nonunique test ids in math_testcases.txt."

    Oh, I didn't notice it: why do we need explicit identifiers? Can't we use the line number or something like that? What happens if two tests have the same identifier? Is the first test skipped?

    Any news from the buildbots?

    Yes, most of them are happy. Some of them are busy (don't have test log2() yet), others are unhappy but not because of log2(). Said differently: log2 tests pass on most buildbots, but we have to wait 12 hours or maybe one day to wait for all buildbots.

    If the test pass on all buildbots, I will commit bpo-11888-part2.patch to use the system log2 function.

    @vstinner
    Copy link
    Member

    vstinner commented May 9, 2011

    we have to wait 12 hours or maybe one day to wait for all buildbots

    Oh, it's faster than expected: test_math passed on FreeBSD 6.4 3.x buildbot. I was waiting for this one because it's an old OS and many tests fail on this buildbot (because it's old but also slow). So test_math passed on all buildbots. Let's try the system log2 :-)

    @python-dev
    Copy link
    Mannequin

    python-dev mannequin commented May 9, 2011

    New changeset 565f43f6bed4 by Victor Stinner in branch 'default':
    Issue bpo-11888: Use system log2() when available
    http://hg.python.org/cpython/rev/565f43f6bed4

    @vstinner
    Copy link
    Member

    vstinner commented May 9, 2011

    Issue bpo-11888: Use system log2() when available
    http://hg.python.org/cpython/rev/565f43f6bed4

    "I expect the system libc to use more accurate functions than Python."

    You know what? Mac OS X log2 is less accurate than Python log2! A log2 test failed on "x86 Tiger 3.x":

    http://www.python.org/dev/buildbot/all/builders/x86%20Tiger%203.x/builds/2488
    ======================================================================
    FAIL: testLog2 (test.test_math.MathTests)
    ----------------------------------------------------------------------

    Traceback (most recent call last):
      File "/Users/db3l/buildarea/3.x.bolen-tiger/build/Lib/test/test_math.py", line 658, in testLog2
        self.assertEqual(actual, expected)
    AssertionError: Lists differ: [-324.0, -323.0, -322.0, -321.... != [-324.0, -323.0, -322.0, -321....

    First differing element 69:
    -254.99999999999997
    -255.0
    ----------------------------------------------------------------------

    Should I revert my patch or should we test the system log2 in configure to check if it is as accurate or more accurate than Mark's algorithm?

    @mdickinson
    Copy link
    Member

    You know what? Mac OS X log2 is less accurate than Python log2!

    That doesn't surprise me much. Though it's probably still true that log2 from OS X is more accurate than our log2 for some other values. It's just that getting the answer wrong for a power of 2 is a 'high-visibility' error.

    Testing log2 sounds long-winded and error-prone. I'd suggest just marking that test as an expected failure on Tiger. BTW, I don't see any such failure on Snow Leopard.

    @mdickinson
    Copy link
    Member

    One other thought: we should check that it's not pow that's at fault here, rather than log2. The test uses math.log2(2.0**n). It would probably be better off using math.log2(ldexp(1.0, n)), or similar: the libm pow operation is also notorious for inaccuracies (due to poor implementations or otherwise) on various platforms.

    @vstinner
    Copy link
    Member

    vstinner commented May 9, 2011

    we should check that it's not pow that's at fault here

    Some tests on Mac OS X Tiger:

    >>> (2.0 ** -255).hex()
    '0x1.0000000000000p-255'

    => pow is correct

    >>> import ctypes; import ctypes.util, math
    >>> libc = ctypes.cdll.LoadLibrary(ctypes.util.find_library('c'))
    >>> clog2=libc.log2
    >>> clog2.restype=ctypes.c_double
    >>> clog2.argtypes=(ctypes.c_double,)
    >>> clog2(2.0**-255)
    -254.99999999999997
    >>> math.log(2.0**-255) / math.log(2.0)
    -255.0
    
    >>> math.log(2.0**-255)
    -176.75253104278605
    >>> math.log(2.0**-255).hex()
    '-0x1.61814bbfb3fb5p+7'
    >>> math.log(2.0)
    0.6931471805599453
    >>> math.log(2.0).hex()
    '0x1.62e42fefa39efp-1'
    
    >>> clog2(2.0**-255).hex()
    '-0x1.fdfffffffffffp+7'
    >>> (math.log(2.0**-255) / math.log(2.0)).hex()
    '-0x1.fe00000000000p+7'

    clog2() is wrong for 2^-255.

    @mdickinson
    Copy link
    Member

    Okay, thanks. We should still be using ldexp rather than 2.0**... in the tests, though; I've fixed this, and also fixed the incorrect (too small) range for those tests, so that all representable powers of 2 are now covered.

    @mdickinson
    Copy link
    Member

    Grr. Got the issue number wrong in the commit message; see msg135584.

    New changeset 1f23d63b578c by Mark Dickinson in branch 'default':
    Issue bpo-11188: In log2 tests, create powers of 2 using ldexp(1, n) instead of the less reliable 2.0**n.
    http://hg.python.org/cpython/rev/1f23d63b578c

    @rpetrov
    Copy link
    Mannequin

    rpetrov mannequin commented May 10, 2011

    Why configure script check two times for log2 function ?

    @python-dev
    Copy link
    Mannequin

    python-dev mannequin commented May 10, 2011

    New changeset d3f9895e2e19 by Mark Dickinson in branch 'default':
    Issue bpo-11888: remove duplicate check for log2 in configure.in.
    http://hg.python.org/cpython/rev/d3f9895e2e19

    @mdickinson
    Copy link
    Member

    Thanks, Roumen. Fixed.

    @mdickinson
    Copy link
    Member

    Victor, what do you think about simply #undefining HAVE_LOG2 on Tiger (e.g. in pyport.h), so that the fallback log2 version is used there instead of the system version?

    Does anyone know the appropriate preprocessor check for OS X <= 10.4? I can get as far as "#ifdef __APPLE__", but don't know how to check for specific versions of OS X.

    @python-dev
    Copy link
    Mannequin

    python-dev mannequin commented May 10, 2011

    New changeset 34871c3072c9 by Victor Stinner in branch 'default':
    Issue bpo-11888: skip some log2 tests on Mac OS X Tiger
    http://hg.python.org/cpython/rev/34871c3072c9

    @vstinner
    Copy link
    Member

    New changeset 34871c3072c9 by Victor Stinner in branch 'default':
    Issue bpo-11888: skip some log2 tests on Mac OS X Tiger

    Oh... I realized that the test doesn't fail on Mac OS X Tiger PPC, only on Mac OS X Tiger x86. But I am too lazy to patch the test. Or should I do it?

    I wait for the following build to close this issue.
    http://www.python.org/dev/buildbot/all/builders/x86%20Tiger%203.x/builds/2507

    @vstinner
    Copy link
    Member

    I wait for the following build to close this issue.
    http://www.python.org/dev/buildbot/all/builders/x86%20Tiger%203.x/builds/2507

    Oh, it's the wrong build. The correct build is:
    http://www.python.org/dev/buildbot/all/builders/x86%20Tiger%203.x/builds/2508

    And it passed so I close this issue.

    @mdickinson
    Copy link
    Member

    Thanks, Victor.

    @ezio-melotti ezio-melotti transferred this issue from another repository Apr 10, 2022
    Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
    Labels
    type-feature A feature request or enhancement
    Projects
    None yet
    Development

    No branches or pull requests

    3 participants