``` 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37``` ```function [out, terms] = mp_BBP(order, precision) %# function [out, last] = mp_BBP(order, precision) This is a %# straight implementation of the Bailey-Borwein-Plouffe formula, %# aiming at computing pi approximation. The inputs are the number %# of terms to use, and the precision. Output is the sum of terms, %# converging to pi; and the smallest of the series, to get an idea %# of the order of the residuals. %# See http://en.wikipedia.org/wiki/Bailey%E2%80%93Borwein%E2%80%93Plouffe_formula if (0 == nargin) precision = 0; end if (precision < 1) %# better not to try mp(0, -something) mp_defaults precision = default_precision; end %# compute the last term in the sequence ii = order; i8 = 8 * mp (ii, precision); x = 4/(i8 + 1)-2/(i8 + 4)-1/(i8 + 5)-1/(i8 + 6); x = x * power (mp (16, precision), -ii); [out, terms(1, 1+order)] = deal (x); %# compute the rest for ii = (order-1:-1:0) i8 = 8 * mp (ii, precision); x = 4/(i8 + 1)-2/(i8 + 4)-1/(i8 + 5)-1/(i8 + 6); x = x * power (mp (16, precision), -ii); out = out + x; terms(1+ii) = x; end end ```