testing non-uniformity with trunc(R2R1R0 . 0,X1Y1Y0X0) computed with neglect of R0.X0

maxbad.py

#!/bin/env/python

"""
This is occupied with finding impact on non-uniformity of evaluating

R times 0,X1Y1Y0X0 in an approximate way which neglects the lowest
order term R0 * X0

Compared to absolute bound 2**(-24) of exact formula near 2**32,
we can not make it much worse because the worst we found is

R=4287201279 N=1388493662 the non-uniformity raises to
1.042290287092328 times 2**(-24).

But there was a bug in the handling of N=0 case, so a worse example
is

R=4282580991 N=0 giving 1.0542947321664542 times 2**(-24)

Compared to R/2**56 however we obtained in the 2**28, 2**32 range:

R=271466495 N=123575687 which gives a non-uniformity
equal to 1.8409720359781416 times R/2**56

(theoretical analysis says it is less than < 2 * R/2**56)

After fixing the N=0 bug I obtained

R=272662527 N=0 which gives 1.9687537664462413 times R/2**56

For R < 2**28, and again comparing to R/2**56 we can obtain
even more dramatic increase as in this example

R=147455, N=120403 non-uniformity is 1944.5611949408294 * R/2**56

But compared to absolute value 2**(-24) we did not find worse than

R=266371071, N=13945551 which gives 0.120624631177634 * 2**(-24)
and this is about 1.93 times 2**(-28). Theoretical guaranteed bound
is 3 times 2**(-28) in this R<2**28 range.

Jeudi 10 mai 2018 à 16:11:07

I have edited Jeudi 10 mai 2018 à 17:16:47 because 
the formula in Fless(N,R) should have been restricted to N>0
and return 0 for N=0 always. Hence maximizebad3() and
maximizeotherbad3() procedures were buggy in former version.

J.-F. B.
"""
def ifloor(a,b):
    """b is positive integer and a is a non-negative integer
    ifloor(a, b) is the maximal integer n with n*b <= a

    """
    return a//b

def iceil(a,b):
    """b is positive integer and a is a non-negative integer
    iceil(a,b) is the number of non-negative integers n with n*b < a

    """
    return -((-a)//b)

twotothefiftysix = 2**56
twotothefourtytwo = 2**42
twotothefourteen = 2**14

def Fless(N, R):
    """With N <= R+1 integers, N non-negative,
    this counts how many integers n = 2**14*x + y there are verifying

    0<= x < 2**42, 
    0<= y < 2**14,

    and

    x*R/2**42 + y*(R-R0)/2**56 < N

    Here R0 = R%2**14, so R-R0 is a multiple of 2**14.

    """
    # attention to special case N=R+1, we want to return 2**56 then
    if N > R:
        return twotothefiftysix

    # attention to special case N=0, we want to return 0 then
    if N == 0:
        return 0

    # if N is at most R, this A will be at most 2**42,
    # And A-1 < 2**42 will hold.
    A = iceil(N*twotothefourtytwo, R)

    return (A - 1) * twotothefourteen + \
        iceil(twotothefourtytwo * N - (A - 1) * R, ifloor(R, twotothefourteen))


def Fequal(N, R):
    """With N <= R integers, N non-negative,
    this counts how many integers n = 2**14*x + y there are verifying

    0<= x < 2**42, 
    0<= y < 2**14,
    
    and

    N = trunc(x*R/2**42 + y*(R-R0)/2**56)

    Here R0 = R%2**14, so R-R0 is a multiple of 2**14.

    """
    return Fless(N+1, R) - Fless(N, R)


def badness(N, R):
    """With N <= R integers, N non-negative, the badness(N, R) is defined
    as a float by this formula:

    Fequal(N, R) = 2**56/R * (1 + badness(N, R) * 2**(-24))

    It measures the non-unformity rescaled by 2**24 for the probability
    of obtaining value N by trunc(R * 0,X_1Y_1Y_0X_0), R=R_2R_1R_0 where the
    R_0 * X_0 term is neglected.

    We haven't incorporated a shift here.

    """
    return abs(Fequal(N, R) * R - 2**56)/2**32


def otherbadness(N, R):
    """With N <= R integers, N non-negative, the otherbadness(N, R) is defined
    as a float by this formula:

    Fequal(N, R) = 2**56/R * (1 + otherbadness(N, R) * R / 2**56)

    It measures the non-unformity for the probability
    of obtaining value N by trunc(R * 0,X_1Y_1Y_0X_0), R=R_2R_1R_0 where the
    R_0 * X_0 term is neglected, compared to R / 2**56 formula

    We haven't incorporated a shift here.

    """
    return abs(Fequal(N, R) * R - 2**56)/R


maxbadness = 0
Rmaxbad = 0
Nmaxbad = 0

maxotherbadness = 0
Rmaxotherbad = 0
Nmaxotherbad = 0

import random

def maximizebad(reps, Rmin, Rmax):
    """This tries to maximize the badness by randomtrials.

    After running about 10**7 times with Rmin=2**28, Rmax=2**32
    I obtained for example

    4252038899 3491140756 1.021719359094277
    """
    global maxbadness, Rmaxbad, Nmaxbad
    for k in range(reps):
        R = random.randrange(Rmin, Rmax)
        N = random.randrange(R)
        b = badness(N, R)
        if b > maxbadness:
            Rmaxbad = R
            Nmaxbad = N
            maxbadness = b
            print(R, N, b)


def maximizebad2(reps, Rmin, Rmax):
    """This tries to maximize the badness by randomtrials.

    But we focus on R's with R0 = 2**14 -1 only

    I obtained with Rmin=2**28, Rmax=2**32

    4287201279 1388493662 1.042290287092328

    (a worse non-uniformity was found by maximizebad3(), see its docstring)

    Then with Rmin=2**21, Rmax=2**28

    266371071 13945551 0.120624631177634

    Then with Rmin=2**17, Rmax=2**21

    147455 66949 0.06676075770519674
    """
    Rmin2 = Rmin//2**14
    Rmax2 = Rmax//2**14
    global maxbadness, Rmaxbad, Nmaxbad
    for k in range(reps):
        R = (random.randrange(Rmin2, Rmax2) + 1) * twotothefourteen - 1
        N = random.randrange(R)
        b = badness(N, R)
        if b > maxbadness:
            Rmaxbad = R
            Nmaxbad = N
            maxbadness = b
            print(R, N, b)


def maximizebad3(reps, Rmin, Rmax):
    """This tries to maximize the badness by randomtrials with only N=0.

    But we focus on R's with R0 = 2**14 -1 only and only
    test the value N = 0.

    I obtained with Rmin=2**28, Rmax=2**32 and 10**7 trials

    4282580991 0 1.0542947321664542

    """
    Rmin2 = Rmin//2**14
    Rmax2 = Rmax//2**14
    global maxbadness, Rmaxbad, Nmaxbad
    for k in range(reps):
        R = (random.randrange(Rmin2, Rmax2) + 1) * twotothefourteen - 1
        N = 0
        b = badness(N, R)
        if b > maxbadness:
            Rmaxbad = R
            Nmaxbad = N
            maxbadness = b
            print(R, N, b)


def maximizeotherbad(reps, Rmin, Rmax):
    """This tries to maximize the otherbadness by randomtrials.

    After running about 10**7 times with Rmin=2**28, Rmax=2**32
    I obtained for example

    268483913 178568437 1.7457040638408752
    """
    global maxotherbadness, Rmaxotherbad, Nmaxotherbad
    for k in range(reps):
        R = random.randrange(Rmin, Rmax)
        N = random.randrange(R)
        b = otherbadness(N, R)
        if b > maxotherbadness:
            Rmaxotherbad = R
            Nmaxotherbad = N
            maxotherbadness = b
            print(R, N, b)


def maximizeotherbad2(reps, Rmin, Rmax):
    """This tries to maximize the otherbadness by randomtrials.

    But we focus on R's with R0 = 2**14 -1 only

    I obtained with Rmin=2**28, Rmax=2**32 after about 10**7 trials

    271466495 123575687 1.8409720359781416

    But with N=0 only a larger example was found by maximizeotherbad3()

    272662527 0 1.9687537664462413

    Then with Rmin=2**17, Rmax=2**28

    147455 48744 1944.5611949408294

    Then Rmin=2**17, Rmax=2**21 found

    147455 120403 1944.5611949408294

    which is same worse R and another N

    """
    Rmin2 = Rmin//2**14
    Rmax2 = Rmax//2**14
    global maxotherbadness, Rmaxotherbad, Nmaxotherbad
    for k in range(reps):
        R = (random.randrange(Rmin2, Rmax2) + 1) * twotothefourteen - 1
        N = random.randrange(R)
        b = otherbadness(N, R)
        if b > maxotherbadness:
            Rmaxotherbad = R
            Nmaxotherbad = N
            maxotherbadness = b
            print(R, N, b)


def maximizeotherbad3(reps, Rmin, Rmax):
    """This tries to maximize the otherbadness by randomtrials with N=0 only

    But we focus on R's with R0 = 2**14 -1 only. Furthermore we only
    test the value N = 0.

    I obtained with Rmin=2**28, Rmax=2**32

    272662527 0 1.9687537664462413
    """
    Rmin2 = Rmin//2**14
    Rmax2 = Rmax//2**14
    global maxotherbadness, Rmaxotherbad, Nmaxotherbad
    for k in range(reps):
        R = (random.randrange(Rmin2, Rmax2) + 1) * twotothefourteen - 1
        N = 0
        b = otherbadness(N, R)
        if b > maxotherbadness:
            Rmaxotherbad = R
            Nmaxotherbad = N
            maxotherbadness = b
            print(R, N, b)