NORX constants

"""Calculate initialisation constants of NORX

Author: Christian Heimes <[email protected]>
License: CC0
"""
from __future__ import print_function
from decimal import Decimal, localcontext, getcontext

# http://www.geom.uiuc.edu/~huberty/math5337/groupe/digits.html
PI = "3.141592653589793238462643383279502884197169399375105820974944592307816406286208998"

NORX_64 = (0x243F6A8885A308D3, 0x13198A2E03707344, 0xA4093822299F31D0,
           0x082EFA98EC4E6C89, 0xAE8858DC339325A1, 0x670A134EE52D7FA6,
           0xC4316D80CD967541, 0xD21DFBF8B630B762, 0x375A18D261E7F892,
           0x343D1F187D92285B)

NORX_32 = (0x243F6A88, 0x85A308D3, 0x13198A2E, 0x03707344, 0x254F537A,
           0x38531D48, 0x839C6E83, 0xF97A3AE5, 0x8C91D88C, 0x11EAFB59)


def approx_pi():
    """Compute Pi to the current precision.

    >>> print(pi())
    3.141592653589793238462643383

    https://docs.python.org/3/library/decimal.html#recipes
    """
    getcontext().prec += 2  # extra digits for intermediate steps
    three = Decimal(3)      # substitute "three=3.0" for regular floats
    lasts, t, s, n, na, d, da = 0, three, 3, 1, 0, 0, 24
    while s != lasts:
        lasts = s
        n, na = n+na, na+8
        d, da = d+da, da+32
        t = (t * n) / d
        s += t
    getcontext().prec -= 2
    return +s               # unary plus applies the new precision


def norx_pi(w=64):
    with localcontext() as ctx:
        ctx.prec = 80
        # pi = Decimal(PI)
        pi = approx_pi()
        frac = pi - Decimal(3)
        ft = int(frac * 2**(4*w))
    return tuple((ft >> i*w) & (2**w-1) for i in (3, 2, 1, 0))


def g(a, b, c, d, w=64):
    r1, r2, r3, r4 = (8, 19, 40, 63) if w == 64 else (8, 11, 16, 31)
    mask = 2**w - 1

    def rotr(x, r):
        return ((x >> r) | (x << (w - r))) & mask

    def h(x, y):
        return ((x ^ y) ^ ((x & y) << 1)) & mask

    a = h(a, b)
    d = rotr(a ^ d, r1)
    c = h(c, d)
    b = rotr(b ^ c, r2)
    a = h(a, b)
    d = rotr(a ^ d, r3)
    c = h(c, d)
    b = rotr(b ^ c, r4)
    return a, b, c, d


def norx32():
    norx32_0 = norx_pi(32)
    norx32_1 = g(*norx32_0, w=32)
    norx32_2 = g(*norx32_1, w=32)
    return norx32_0 + norx32_1 + norx32_2[:2]


def norx64():
    norx64_0 = norx_pi(64)
    norx64_1 = g(*norx64_0, w=64)
    norx64_2 = g(*norx64_1, w=64)
    return norx64_0 + norx64_1 + norx64_2[:2]


if __name__ == "__main__":
    for number in (NORX_64, norx64()):
        print(" ".join("%0.16x" % n for n in number))
    for number in (NORX_32, norx32()):
        print(" ".join("%0.8x" % n for n in number))