60 lines
1.4 KiB
Matlab
60 lines
1.4 KiB
Matlab
function pi_str = piSpigot(N)
|
|
% Return N digits of pi using Gibbons's first spigot algorithm.
|
|
% If N is omitted, the digits are printed ad infinitum.
|
|
% Uses the expansion
|
|
% pi = sum_{i=0} (i!)^2 2^{i+1} /(2i+1)!
|
|
% = 2 + 1/3 * ( 2 + 2/5 * (2 + 3/7 * ( 2 + 4/9 * ( ..... )))))
|
|
% = (2 + 1/3 *)(2 + 2/5 *)(2 + 3/7 *)...
|
|
% where the terms in the last expression represent Linear Fractional
|
|
% Transforms (LFTs).
|
|
%
|
|
% Requires the Variable Precision Integer (vpi) Toolbox
|
|
%
|
|
% Reference:
|
|
% "Unbounded Spigot Algorithms for the Digits of Pi" by J. Gibbons, 2004
|
|
% American Mathematical Monthly, vol. 113.
|
|
if nargin < 1
|
|
N = Inf;
|
|
lineLength = 50;
|
|
else
|
|
pi_str = repmat(' ',1,N);
|
|
end
|
|
|
|
q = vpi(1);
|
|
r = vpi(0);
|
|
t = vpi(1);
|
|
k = 1; % If printing more than 3E15 digits, use k = vpi(1);
|
|
|
|
i = 1;
|
|
first_digit = true;
|
|
while i <= N
|
|
threeQplusR = 3*q + r;
|
|
n = double(threeQplusR / t);
|
|
if q+threeQplusR < (n+1)*t
|
|
d = num2str(n);
|
|
if isinf(N)
|
|
fprintf(1,'%s', d);
|
|
if first_digit
|
|
fprintf(1,'.');
|
|
first_digit = false;
|
|
i = i+1;
|
|
end
|
|
if i == lineLength
|
|
fprintf(1,'\n');
|
|
i = 0;
|
|
end
|
|
else
|
|
pi_str(i) = d;
|
|
end
|
|
q = 10*q;
|
|
r = 10*(r-n*t);
|
|
i = i + 1;
|
|
else
|
|
t = (2*k+1)*t;
|
|
r = (4*k+2)*q + (2*k+1)*r;
|
|
q = k*q;
|
|
k = k + 1;
|
|
end
|
|
end
|
|
end
|