Probability sequential terms of a linear congruential generator are coprime

450 Views Asked by At

I was following section 3.1.2 of the Structure and Interpretation of Computer Programming where they calculate $\pi$ using Cesàro's theorem that the probability two randomly chosen numbers are coprime is $6 \over \pi^2$.

When I implemented this using a linear congruential generator, I found the probability looked close to $6 \over e^2$ for a variety of seeds and a couple of generators. (With a Blum Blum Shub generator I got $6 \over \pi^2$ as expected).

Is the probability that two sequential terms of a (suitable) linear congruential generator coprime $6 \over e^2$? Is there a proof/disproof?

Here is the Scheme code for numerical validation:

(define (lcg modulus multiplier increment)
  (lambda (x)
    (remainder (+ (* multiplier x) increment) modulus)))
;From Knuth's MMIX
(define rand-update (lcg (expt 2 64) 6364136223846793005 1442695040888963407))
(define random-init 1)
;Below is copied from SICP
(define rand
  (let ((x random-init))
    (lambda ()
      (set! x (rand-update x))
      x)))
(define (estimate-pi trials)
  (sqrt (/ 6 (monte-carlo trials cesaro-test))))
(define (cesaro-test)
   (= (gcd (rand) (rand)) 1))
(define (monte-carlo trials experiment)
  (define (iter trials-remaining trials-passed)
    (cond ((= trials-remaining 0)
           (/ trials-passed trials))
          ((experiment)
           (iter (- trials-remaining 1) (+ trials-passed 1)))
          (else
           (iter (- trials-remaining 1) trials-passed))))
  (iter trials 0))
;(estimate-pi 100000) = 2.72268....
1

There are 1 best solutions below

1
On BEST ANSWER

I get that the probability that two numbers of different parities are relatively primes is $p=\frac{8}{\pi^2}$.(*) If you then did your above computation:

$$\sqrt{\frac{6}{p}} = \frac{\sqrt{3}}2 \pi\approx 2.7207$$

which is close to what you got. The fact that this is somewhat close to $e$ is likely a coincidence.

While it is easier in an LCG to use a power of $2$ as your modulus, it almost always gives parity problems like this if you do. In this particular example, it is better to use a large prime, I think, for the modulus.

(*) Technically not a probability, but a limit probability or a density, but we'll call it a probability here.