Self-Locating Strings In Pi

January 4, 2019

We have written several programs that compute the digits of π, including this one just a week after the inception of Programming Praxis, based on a spigot algorithm of Jeremy Gibbons:

(define (pi-spigot z)
  (let loop ((z z) (ds '()) (q 1) (r 0) (t 1) (k 1) (n 3) (m 3))
    (cond ((zero? z) (reverse ds))
          ((< (+ q q q q r (- t)) (* n t))
            (loop (- z 1) (cons n ds) (* 10 q) (* 10 (- r (* n t)))
                  t k (- (quotient (* 10 (+ q q q r)) t) (* 10 n)) m))
          (else (loop z ds (* q k) (* (+ q q r) m) (* t m) (+ k 1)
                      (quotient (+ (* q (+ k k k k k k k 2)) (* r m))
                      (* t m)) (+ m 2))))))

It takes about ten minutes to compute the first 50,000 digits of π:

> (define pi (pi-spigot 50000))

Then it’s just a matter of indexing through the digits and checking for self-locating indices:

> (do ((n 0 (+ n 1)) (ps (cdr pi) (cdr ps))) ((null? ps))
    (let* ((ds (digits n)) (len (length ds)))
      (when (equal? ds (take len ps)) (display n) (newline))))
0
6
27
13598
43611

The initial 0 surprised me. When n = 0, digits returns (). But then the length of the digit string is 0, so take returns a null list, which is equal? to ds, so 0 is a proper result. Which it isn’t. I guess that’s a design bug in digits.

You can run the program at https://ideone.com/TnoBhB.

Advertisements

Pages: 1 2

7 Responses to “Self-Locating Strings In Pi”

  1. I’m glad your getting better, back at work, and had family time at Christmas.

    I came across your website because mine’s called https://www.praxismachinecodinglessons.com

    Adapting your coding challenges to developing curriculum for high school computer science teachers and students will be a main theme of mine in 2019.

    I’ll keep you posted,
    Peter L.
    ________________________________

  2. Paul said

    Hi Phil, welcome back.
    I cheated a little with this exercise and downloaded 1 million digits from the internet. Results were the same as from @programmingpraxis. Running time is about 4 seconds.

    from collections import deque
    
    REVERSE_SLICE = slice(None, None, -1)
    DDIG = dict(zip("0123456789", range(10)))
    
    def digits(number, reverse=False):
        s = str(abs(number))
        if reverse:
            s = s[REVERSE_SLICE]
        return (DDIG[c] for c in s)
    
    def pidigs(fname):
        with open(fname) as f:
            yield from (int(char) for char in f.readline() if char != ".")
            for line in f:
                yield from (int(char) for char in line)
    
    def self_locating_strings(fname, NMAX):
        pidigits = pidigs(fname)
        nd = len(str(NMAX))  # make the deque long enough
        D = deque(it.islice(pidigits, nd), nd)
        for i in range(NMAX-10):
            D.append(next(pidigits))
            if all(d1 == d2 for d1, d2 in zip(digits(i), D)):
                print(i)
    
    t0 = clock()
    self_locating_strings(r"D:\Data\Development\python\src\pywork\prpr2\pi1000000.txt", 1000000)
    print(clock()-t0)
    
  3. Steve said

    $ date && awk ‘END { for (i = 0; i <= 999999; i++) { len = log(i)/log(10); len = int(len) + 1; str = substr($0,i+1,len); if (str ==i) print str } }’ value_of_pi.txt && date

    Fri, Jan 4, 2019 10:32:49 AM

    6
    27
    13598
    43611

    Fri, Jan 4, 2019 10:32:50 AM

  4. Steve said

    AWK version

    $ date && awk 'END { for (i = 0; i <= 999999; i++) { len = log(i)/log(10); len = int(len) + 1; str = substr($0,i+1,len); if (str ==i) print str } }' value_of_pi.txt && date
    
    Fri, Jan  4, 2019 10:32:49 AM
    
    6
    27
    13598
    43611
    
    Fri, Jan  4, 2019 10:32:50 AM
    
    ;; Less than 2 seconds
    ;; Amazes me that AWK supports variables that can hold 1,000,000 characters!
    
    
  5. chaw said

    Welcome back @programmingpraxis, and Hacky New Year to all.

    For the digits of pi, I cheated and used Simon Tatham’s excellent
    ‘spigot’ program to do the real work, so my main code isn’t much worth
    posting. However, apropos the digits bug you mention, here’s the
    simple implementation in R7RS Scheme and SRFI 8 (for receive) from my
    homegrown utility library.

    (define digits
      (case-lambda
        ((n) (digits n 10)) 
        ((n radix)
         (unless (> radix 1)
           (error "radix must be greater than 1" radix))
         (if (zero? n)
             '(0)
             (let f ((n n)
                     (digits '()))
               (if (zero? n)
                   digits
                   (receive (quot rem) (floor/ n radix)
                     (f quot
                        (cons rem digits)))))))))
    

  6. matthew said

    Great to have you back again! Hope recovery continues to go well.

    Interesting problem, never really looked at generation of digits of pi before. The Gibbons spigot is neat but seems to slow down drastically as we get further into the number – lots of very large numbers in the internal state that presumably are needed for later digits. The Machin formula is neat as we can just sum lot of reciprocals to get our approximation, using some sort of fixed-point representation. I tried some of the “optimal” Machin formulas, but there wasn’t an enormous speedup that I saw (maybe 20-30%) so stuck with the original 1706 method:

    """
    Taylor series for arctan:
    
    arctan(z) = z - z^3/3 + z^5/5 - z^7/7 + ...
    
    And Machin's 1706 formula:
    
    pi/4 = 4*arctan(1/5) - arctan(1/239)
    
    So we can compute pi just by summing (integer multiples of) reciprocals.
    
    Use fixed point arithmetic with K = 10**N places, then X/n is approximately
    m/K where m = (XK + n//2)//n (we add n//2 to round to closest).
    
    So we can just generate the series above until the term absolute value is
    less than K (the steps are decreasing and alternating sign, so each term
    bounds the sum of the remaining sequence).
    
    The last few digits of the output are likely to be wrong due to rounding error
    (eg. for 50000 digits, the last 3 are wrong).
    """
    
    import sys
    
    def sumarctan(z,A):
        total = 0
        k,t,s = 1,z,1
        while True:
            x = k*t
            if x > A: break
            total += (A+x//2)//x*s
            k,t,s = k+2,t*z*z,-s
        return total
        
    def pi(K): return (sumarctan(5,16*K) - sumarctan(239,4*K))
    
    if __name__ == '__main__':
        N = 1000 if len(sys.argv) == 1 else int(sys.argv[1])
        K = 10**N
        # mod to strip off integer part
        pidigits = str(pi(K)%K)
        #print(pidigits)
        for i in range(1,N):
            s = str(i)
            if s == pidigits[i:i+len(s)]: print(s)
    

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out /  Change )

Google photo

You are commenting using your Google account. Log Out /  Change )

Twitter picture

You are commenting using your Twitter account. Log Out /  Change )

Facebook photo

You are commenting using your Facebook account. Log Out /  Change )

Connecting to %s

%d bloggers like this: