The First n Digits of Pi

This program is based on several C programs, which were based on a Pascal program, which was based on the article "A spigot algorithm for pi" in the "American Mathematical Monthly", Volume 102(?), Number 3, March 1995, page 195, by Stanley Rabinowitz and Stan Wagon.

/* pi(N, Pi) is true if the list Pi contains the first N digits of pi.       */
/* e.g. pi(29, [3,1,4,1,5,9,2,6,5,3,5,8,9,7,9,3,2,3,8,4,6,2,6,4,3,3,8,3,2]). */
pi(N, Pi):-
  N > 0,
  Length is (N * 10) // 3 + 1,
  copies(Length, 2, As),
  pi_1(N, Length, 0, 0, As, [], Pi0),
  reverse(Pi0, Pi).

pi_1(0, _, _, _, _, Pi, Pi):-!.
pi_1(J, Length, Nines0, Previous0, As, Pi0, Pi):-
  next_x(Length, 0, As, Bs, X),
  next_digits(X, J, J1, Nines0, Nines, Previous0, Previous, Pi0, Pi1),
  pi_1(J1, Length, Nines, Previous, Bs, Pi1, Pi).

next_x(0, Q, [], [], Q):-!.
next_x(I, Q, [A|As], [B|Bs], X):-
  T is 10 * A + Q * I,
  B is T mod (2 * I - 1),
  Q1 is T // (2 * I - 1),
  I1 is I - 1,
  next_x(I1, Q1, As, Bs, X).

next_digits(X, J, J, Nines0, Nines, Previous, Previous, Pi, Pi):-
  X mod 10 =:= 9, !, 
  Nines is Nines0 + 1.
next_digits(X, J0, J, Nines0, Nines, Previous0, Previous, Pi0, Pi):-
  J1 is J0 - 1,
  Digit is Previous0 + X // 10,
  Digit1 is 9 * (1 - X // 10),
  append_digits(J1, J, Nines0, Nines, Digit1, [Digit|Pi0], Pi),
  Previous is X mod 10.

append_digits(J0, J, N0, N, Digit, Pi0, Pi):-
  J0 > 0, 
  N0 > 0, !,
  J1 is J0 - 1,
  N1 is N0 - 1,
  append_digits(J1, J, N1, N, Digit, [Digit|Pi0], Pi).
append_digits(J, J, N, N, _, Pi, Pi).

/* copies(N, X, Xs) is true if the list Xs contains N occurrences of the     */
/*   item X.                                                                 */
copies(0, _, []):-!.  
copies(I, X, [X|Xs]):-
  I > 0,
  I1 is I - 1,
  copies(I1, X, Xs).

/* reverse(Xs, Ys) is true if Ys is the result of reversing the order of     */
/*   the elements in the list Xs.                                            */
%reverse(Xs, Ys):-reverse_1(Xs, [], Ys).

%reverse_1([], As, As).
%reverse_1([X|Xs], As, Ys):-reverse_1(Xs, [X|As], Ys).

LPA Index     Home Page