Series_L_Upper_Lower_Bounds.py

# series_L_Upper_Lower_Bounds22b.py --> Series_L_Upper_Lower_Bounds.py

# try to set up an option to set that shows a new result only when there is a new idx
# Did that.
# Also set up an easy way to configure to compute for only a single n
#
#===========================================

# "The Euler Series" on p.31 in "The Pleasures of Pi,e" is

# 1/1**2 + 1/2**2 + 2/3**2 + ... + 1/n**2 + ... --> pi**2/6

# The partial sums converge very slowly to pi**2/6, so I use Hans Lundmark's post on http://math.stackexchange.com, at
# http://math.stackexchange.com/questions/8337/different-methods-to-compute-sum-n-1-infty-frac1n2

# There's an image of the last inequality in his proof at http://www.rcblue.com/images/PinchingForEuler.jpg

# I use this inequality after multiplying the 3 expressions by (math.pi**2/(4**(n+1)))
# Thus upper_bound_partial_sum = (2*(4**n - 1)/3) * (math.pi**2/(4**(n+1)))
# and lower_bound_partial_sum = (2*(4**n - 1)/3 - (2**n - 1)) * (math.pi**2/(4**(n+1)))

# The middle expression becomes lines 49-56 below,
# which compute the euler_series_partial_sum for a given n

# NOTE: the above 2 lines are no longer true. The middle term takes so long to compute that
# the demo is not very intersting. So I decided to just show that the upper and lower bounds on the Euler
# Series partial sums get VERY close to one another, and in fact are identical
# for a large number of their initial digits.

#For singleton_first_n = 10,000,000 the output is
# (for singleton_first_n, see main() )

#=================================================
#1.6449340668482264364724151666...351904003992141795698445313293*7211734124
#For n = 10,000,000  idx = 3,010,301
#the first 3,010,300 digits of upper and lower match up to the asterisk
#That's 4.8e+00 miles at 10 digits/inch.
#1.6449340668482264364724151666...351904003992141795698445313293*4485269029

#time was 23.4294 seconds  or  00:00:23.43
#=====================================================

# Thus this script demos the "pinching down" (for increasing values of n)
# on the middle expression's value by the upper_bound_partial_sum from above,
# and the lower_bound_partial_sum from below.
# Because the limits of both the upper_bound_partial_sum's and the lower_bound_partial_sum's are pi^2/6,
# the limit of the euler_series_partial_sum's must also be pi^2/6. Q.E.D.

from mpmath import *; mp.pretty = True  # mp.dps is set in main() for each n

from time import clock
import sys

def tim(t, digits, p=1):
    time_  = clock() - t
    if digits:
        seconds = round(time_, digits)
        minutes = round(time_/60, digits)
        hours = round(time_/3600, digits)
    else:
        seconds = int(round(time_, digits))
        minutes = int(round(time_/60, digits))
        hours = int(round(time_/3600, digits))
    if p and seconds > 3600:
        print("time was", hours, "hours  or ", secsToHMS(time_))
    elif p and seconds > 60:
        print("time was", minutes, "minutes  or ", secsToHMS(time_))
    elif p:
        print("time was", seconds, "seconds  or ", secsToHMS(time_))
    else:
        return secsToHMS(time_)

def secsToHMS(seconds):
    """
    Convert seconds to hours:minutes:seconds, with seconds rounded to hundredths of a second.
    """
    minutes, seconds = divmod(seconds, 60)
    hours, minutes = divmod(minutes, 60)
    return "%02d:%02d:%05.2f" % (hours, minutes, seconds)


def int_commas(n):
    """
    inserts commas into integers and returns the string

    E.g. -12345678 -> -12,345,789
    """
    return format(n, ',d')

def compareSequences(seq1, seq2):
    """
    Returns first index at which two sequences differ. If one is longer than the other,
    ignores the elements of the longer for which there are no corresponding elements in the shorter.
    """
    if len(seq2) > len(seq1):
        seq1, seq2 = seq2, seq1
    for index in range(0, len(seq2)):
        if seq1[index] != seq2[index]:
            return index

def num_digits2int_length(num_digits, dpi = 10):
    """
    Given number of digits of (large) integer n,
    returns length of n in miles

    For my laptop monitor, there are
    about 10 digits per inch (dpi = 10)
    """
    #digits = int_digits(n)
    miles = num_digits/dpi/12/5280
    return miles

def sig_digits(n, d=6):
    """
    Return any real number n to d significant digits, in scientific notation.
    """
    return '%.*e' % (d-1, n)

def demo(n, complete, first_n, idx, prev_idx):
    n = mpf(n)
    prev_idx = idx
    upper_bound_partial_sum = (2*(4**n - 1)/3) * (pi()**2/(4**(n+1)))
    lower_bound_partial_sum = (2*(4**n - 1)/3 - (2**n - 1)) * (pi()**2/(4**(n+1)))

    idx = compareSequences(str(upper_bound_partial_sum), str(lower_bound_partial_sum))

    if n > first_n and idx != prev_idx:
        if n < 384 or complete == True:
            #if complete == True:
            upper = str(upper_bound_partial_sum)[:idx] + '*' + str(upper_bound_partial_sum)[idx:idx+10]
            lower = str(lower_bound_partial_sum)[:idx] + '*' + str(lower_bound_partial_sum)[idx:idx+10]
            print(upper)
            print("for n =", int_commas(int(n)), " idx =", int_commas(idx))
            print("the first", idx-1, "digits of upper and lower match up to the asterisk")
            print("for n =", int_commas(int(n)), " idx =", int_commas(idx))
            print(lower)

        else:
            upper1 = str(upper_bound_partial_sum)[:30]
            upper2 = str(upper_bound_partial_sum)[idx-30:idx] + '*' + str(upper_bound_partial_sum)[idx:idx+10]
            upper = upper1 + '...' + upper2
            lower = str(lower_bound_partial_sum)[idx-50:idx] + '*' + str(lower_bound_partial_sum)[idx:idx+10]
            print(upper)
            print("For n =", int_commas(int(n)), " idx =", int_commas(idx))
            print("the first", int_commas(idx-1), "digits of upper and lower match up to the asterisk")
            print("That's", sig_digits(num_digits2int_length(idx-1),2), "miles at 10 digits/inch.")
            lower1 = str(lower_bound_partial_sum)[:30]
            lower2 = str(lower_bound_partial_sum)[idx-30:idx] + '*' + str(lower_bound_partial_sum)[idx:idx+10]
            lower = lower1 + '...' + lower2
            print(lower)
        print()
    return idx


def main():
    t0 = clock()
    singleton = 1
    singleton_first_n = 10000000
    if singleton_first_n > 646456000:
        print("DAME!! n must be <= 646456000")
        print("Exiting.")
        sys.exit()

    first_n = 564
    last_n = 574
    every_nth = 1
    if singleton:
        first_n = singleton_first_n
        last_n = first_n
        first_n -= 1
        every_nth = 1
    else:
        first_n -= every_nth

    complete = 0
    idx = 0
    prev_idx = 0
    for n in range(first_n, last_n+1, every_nth):
        #this should ensure that n+1 will not be skipped
        if n == first_n:
            n = -20
            prev_idx = 0
        #make mp.dps large enough (but not too large) for small n
        if n < 150:
            mp.dps = 65
        else:
            #get mp.dps just about right for larger n
            mp.dps = int(n/2.5)
        #until this point n shouldn't be an mpf (for n = -20, and assigning mp.dpf)
        n = mpf(n)
        idx = demo(n, complete, first_n, idx, prev_idx)

    if n >= 1000000:
        tim(t0, 4)

if __name__ == '__main__':
    main()

"""
OUTPUT
C:\P32Working\Series>series_L_Upper_Lower_Bounds22a.py
80209552093784631871255223825845694405451612628725*4702464901
for n = 646,456,000  idx = 194,602,648
That's 3.1e+02 miles at 10 digits/inch.
the first 194602647 digits of upper and lower match up to the asterisk (307 miles at 10 digits/inch)
80209552093784631871255223825845694405451612628725*1426908163

time was 47.9024 minutes  or  00:47:54.14
===========================================
79692796739533168682606523239647663116122558403219*7481745350
For n = 100,000,000  idx = 30,103,000
the first 30,102,999 digits of upper and lower match up to the asterisk
That's 4.8e+01 miles at 10 digits/inch.
79692796739533168682606523239647663116122558403219*6812104970

time was 5.256 minutes  or  00:05:15.36


"""