Aaron S. Hawley (aaronhawley) wrote,
Aaron S. Hawley
aaronhawley

Frobbing primes with Emacs Lisp

Yoni Rabkin wrote a short post called waaaaa? pi? again?. He read that the probability of integer pairs being coprime is equal to 6 over pi squared (6 / pi^2).

I wrote some Emacs Lisp to test this out in GNU Emacs 22, and it seems to get the same result.

I wrote my own greatest common divisor (GCD) function, but just used the one in Emacs.

  (require 'calc)
  (require 'calc-ext)

  ;; greatest common divisor
  (defun my-gcd (x y)
    (if (eq y 0)
        x
      (my-gcd y (% x y))))

  ;; x and y are coprime if gcd(x, y) = 1
  (defun coprimep (x y)
    (eq 1 (math-gcd x y)))

  (defun coprime-probability (n maxint)
    "Probability of coprime for N random numbers of MAXINT."
    (random t)
    (/
     (apply '+
            (mapcar (function
                     (lambda (p)
                       (if (coprimep (car p) (cdr p))
                           1.0 0.0)))
                    (mapcar* 'cons
                             (mapcar 'random
                                     (make-vector n maxint))
                             (mapcar 'random
                                     (make-vector n maxint)))))
     n))

Here's an example of running the coprime generation using a million pairs and using as much precision Emacs can handle for integers, and then using the result to calculate pi.

(let* ((n #xfffff) ;; 1048575
       (maxint #xfffffff) ;; 268435455
       (p (coprime-probability n maxint))
       (pi-approx (if (eq 0 p) 0.0 (sqrt (/ 6 p)))))
  (princ (format
          "Probability that %d pairs of integers less than %d are prime: %f\n"
          n maxint p))
  (princ (format "pi (%s) is approximately "
                 (math-format-number (math-pi))))
  pi-approx)

Here's a result of running this once:

Probability that 1048575 pairs of integers less than 268435455 are prime: 0.607957

pi (3.14159265359) is approximately 3.141516663694778

It only takes a few seconds to run. Though, if you increase the number of pairs, Emacs quickly started choking my computer.

What was most intriguing to me about this factoid, was not the relationship to pi. I was surprised that over 60 percent of all points in two-dimensional space are coprime. I didn't expect that to be that high.

  (require 'calc)
  (require 'calc-ext)
  (math-format-number (math-div 6 (math-pow (math-pi) 2)))
  >> "0.607927101854"

The probability of selecting a generally prime number is zero. "Apples to oranges", I know.

Tags: emacs, free software, programming
Subscribe

  • User liberation: New video from the FSF

    from fsf.org community blog The last 45 seconds is pretty cool. There's a build of Gstreamer, interspersed with screenshots of Gnome,…

  • Big Emacs reference card updated

    With the release of Emacs 24.3 last month and the big changes at EmacsWiki, I've posted an updated version of the giant Emacs reference card. It…

  • M-x in Emacs 24.3 is now in Lisp

    It didn't make the NEWS file for Emacs 24.3, but Emacs now ships with an ` M-x' (` execute-extended-command') that is written in Lisp. It is no…

  • Post a new comment

    Error

    default userpic

    Your reply will be screened

    Your IP address will be recorded 

    When you submit the form an invisible reCAPTCHA check will be performed.
    You must follow the Privacy Policy and Google Terms of use.
  • 4 comments