2e:null 2f:["$","title",null,{"children":"Computing π with the Bailey-Borwein-Plouffe Formula | Ricky Reusser | Observable"}] 34:T45a,md` The [Baily-Borwein-Plouffe](https://en.wikipedia.org/wiki/Bailey%E2%80%93Borwein%E2%80%93Plouffe_formula) (BBP) formula is a remarkable formula for computing the hexadecimal digits of ${tex`\pi`}, starting at the ${tex`n^{th}`} digit, *without first computing preceeding digits*! Such algorithms are called *spigot algorithms*. The formula was discovered experimentally in 1995 via the [PSLQ algorithm](https://www.davidhbailey.com/dhbpapers/pslq-comp-alg.pdf), which itself was named one of the [Top Ten Algorithms of the Century](https://archive.siam.org/pdf/news/637.pdf). That such a formula could exist is a bit shocking—though the initial shock is tempered by the realization that it does, of course, take longer as you start at later digits (${tex`O(n \log n)`} as a function of the starting digit, to be precise). Additionally, the hexadecimal digits cannot be converted to decimal without first computing all of the leading digits. What it does enable though, is to compute digits deep into ${tex`\pi`} with hardly any memory consumption and nothing more than a basic double precision data type. `35:T96a,table = { const nDigits = Math.floor((width / 954) * 128); const ref = reference.digits.slice(startingDigit, startingDigit + nDigits); function compareStrings(a, b, offset) { const out = DOM.element("span"); const leading = DOM.element("span"); leading.textContent = startingDigit === 0 && offset > 0 ? ref[0] + "." + ref.slice(1, offset) : ref.slice(0, offset); out.appendChild(leading); out.style.color = "#ccc"; for (let i = 0; i < ref.length - offset; i++) { const s = DOM.element("span"); s.textContent = b[i]; out.appendChild(s); const color = a[i] === b[i] ? "rgb(32, 165, 186)" : "#c12"; s.style.color = color; s.style.fontWeight = 900; if (startingDigit === 0 && offset === 0 && i === 0) { const dot = DOM.element("span"); dot.textContent = "."; dot.style.color = color; out.appendChild(dot); } } const trailing = DOM.element("span"); trailing.textContent = ref.slice(offset + digits, ref.length); out.appendChild(trailing); out.style.color = "#bbb"; return out; } const lines = []; const refStr = startingDigit === 0 ? ref[0] + "." + ref.slice(1) : ref; for (let i = startingDigit; i < startingDigit + ref.length; i += 4) { const line = DOM.element("code"); line.style.fontSize = "inherit"; line.style.color = "#ccc"; line.textContent = refStr; lines.push(line); } const children = lines.map( (l) => html`${l}
` ); const fontSize = 11; // Math.min(14, 1.5 * (width / 128)) const output = html`

${children}

`; const isFloat = method === "Floating point"; const piDigits = isFloat ? computePi : BBPbigint; for ( let j = 0, i = startingDigit; i < startingDigit + ref.length; j++, i += 4 ) { const t1 = performance.now(); const pi = piDigits(i, digits); const t2 = performance.now(); const str = pi.toString(16).padStart(digits, "0"); const line = compareStrings( reference.digits.slice(i, i + digits), str, i - startingDigit ); const time = DOM.element("code"); time.style.fontSize = "inherit"; time.textContent = ` (${(t2 - t1).toFixed(1)}ms)`; line.appendChild(time); lines[j].replaceWith(line); yield output; } }36:T11c6,md`## Using the formula I'll be honest. I've seen this formula probably ten or twenty times, and apart from trivially plugging in some numbers, I'd never stopped to work through how to properly use it. So let's walk through using the identity to compute ${tex`\pi`} starting at an arbitrary location. I'll largely repeat the eponymous David Bailey's explanation from *[The BBP Algorithm for Pi](https://www.davidhbailey.com/dhbpapers/bbp-alg.pdf)*, filling in some details I learned along the way. The implemenation is also heavily influenced by the version from [literateprograms.org](https://literateprograms.org/pi_with_the_bbp_formula__python_.html). Our starting point is the BBP formula, ${tex.block` \pi =\sum_{k=0}^{\infty }\left[{\frac {1}{16^{k}}}\left({\frac {4}{8k+1}}-{\frac {2}{8k+4}}-{\frac {1}{8k+5}}-{\frac {1}{8k+6}}\right)\right]. `} We first multiply through by ${tex`16^d`}, then take the fractional part (modulo 1, if you prefer) which we denote by ${tex`\{\cdot\}.`} ${tex.block` \{16^d\pi\} = \left\{\sum _{k=0}^{\infty }\left( {\frac {16^{d-k} 4}{8k+1}}- {\frac {16^{d-k} 2}{8k+4}}- {\frac {16^{d-k}}{8k+5}}- {\frac {16^{d-k}}{8k+6}} \right)\right\} `} Consider ${tex`d = 4`} (starting at the fourth digit): ${tex.block` \begin{aligned} \pi &\approx \texttt{3.243f6a8885a368d31\ldots}_{16} \\ 16^4 \pi &\approx \texttt{3243f.6a8885a368d31\ldots}_{16} \\ \left\{16^d \pi\right\} = \left\{16^4 \pi\right\} &\approx \phantom{0000\,}\texttt{0.6a8885a368d31\ldots}_{16} \end{aligned} `} In hexadecimal, ${tex`16^d`} shifts the value ${tex`d`} places to the left (analogous to multiplication by ten in base ten), while the fractional part operation chops off any leading digits to the left of the [radix point](https://en.wikipedia.org/wiki/Radix_point). Thus, it pretty directly gives us the hexadecimal representation of ${tex`\pi`} starting at the ${tex`d^{th}`} digit. If we then multiply by ${tex`16^n`} and take the integer part, we obtain ${tex`n`} digits of ${tex`\pi`}, starting at digit ${tex`d`}. Consider ${tex`n = 6`}: ${tex.block` \begin{aligned} \left\{16^4 \pi\right\} &\approx \texttt{0.6a8885a368d31\ldots}_{16} \\ 16^6 \left\{16^4 \pi\right\} &\approx \phantom{0\;\;} \texttt{6a8885.a368d31\ldots}_{16} \\ \lfloor 16^n \left\{16^d \pi\right\} \rfloor = \lfloor 16^6 \left\{16^4 \pi\right\} \rfloor &\approx \phantom{0\;\;} \texttt{6a8885}_{16} \end{aligned} `} We could evaluate the formula and compute this result directly, but we'd pretty quickly find ourselves dealing with a very large number of digits, only to throw almost all of them away with the final modulo 1. Our key observation is that, in order to keep numbers small, the fractional part operation can be distributed into each term of the summation. This works since, for example, ${tex`\left\{3.4 + 5.1\right\} = \left\{\left\{3.4\right\} + \left\{5.1\right\}\right\} = \left\{0.4 + 0.1\right\} = 0.5`}, and yields ${tex.block` \{16^d\pi\} = \left\{\sum _{k=0}^{\infty }\left( \left\{\frac {16^{d-k} 4}{8k+1}\right\}- \left\{\frac {16^{d-k} 2}{8k+4}\right\}- \left\{\frac {16^{d-k}}{8k+5}\right\}- \left\{\frac {16^{d-k}}{8k+6}\right\} \right)\right\}. `} To simplify—and realizing that we can move factors like ${tex`2`} and ${tex`4`} around as well—we break this up and write ${tex.block` \left\{16^d \pi\right\} = \left\{ 4\left\{ 16^d S_1 \right\} - 2\left\{ 16^d S_4 \right\} - \left\{ 16^d S_5 \right\} - \left\{ 16^d S_6 \right\} \right\}, `} where ${tex.block` S_j = \sum_{k=0}^\infty \frac{1}{16^k(8k+j)}. `} To separate the positive (growing) and negative (quickly shrinking) powers of ${tex`16`}, we can split each summation into two parts, from ${tex`0`} to ${tex`d`} and from ${tex`d + 1`} to ${tex`\infty`}. ${tex.block` \begin{aligned} \left\{16^d S_j \right\} &= \left\{\left\{ \sum_{k=0}^d \frac{16^{d-k}}{8k + j}\right\} + \sum_{k=d+1}^\infty \frac{16^{d-k}}{8k+j}\right\} \end{aligned} `} Finally, to reduce the size of the numerator in the first term, we insert a modulo, which we're able to do since it leaves the fractional part we care about unchanged. ${tex.block` \begin{aligned} \left\{16^d S_j\right\} &= \left\{\left\{ \sum_{k=0}^d \frac{16^{d-k} \operatorname{mod} 8k+j}{8k + j}\right\} + \sum_{k=d+1}^\infty \frac{16^{d-k}}{8k+j}\right\} \end{aligned} `} The second term gets small quickly as ${tex`k`} increases. As for the first term, expressions of the form ${tex`b^e \operatorname{mod} m`} are called [*modular exponentiation*](https://en.wikipedia.org/wiki/Modular_exponentiation). `37:T2174,reference = ({ // pi, in hexadecimal, as measured with a ruler ;) digits: ` 3243f6a8885a308d313198a2e03707344a4093822299f31d0082efa98ec4e6c8 9452821e638d01377be5466cf34e90c6cc0ac29b7c97c50dd3f84d5b5b547091 79216d5d98979fb1bd1310ba698dfb5ac2ffd72dbd01adfb7b8e1afed6a267e9 6ba7c9045f12c7f9924a19947b3916cf70801f2e2858efc16636920d871574e6 9a458fea3f4933d7e0d95748f728eb658718bcd5882154aee7b54a41dc25a59b 59c30d5392af26013c5d1b023286085f0ca417918b8db38ef8e79dcb0603a180 e6c9e0e8bb01e8a3ed71577c1bd314b2778af2fda55605c60e65525f3aa55ab9 45748986263e8144055ca396a2aab10b6b4cc5c341141e8cea15486af7c72e99 3b3ee1411636fbc2a2ba9c55d741831f6ce5c3e169b87931eafd6ba336c24cf5 c7a325381289586773b8f48986b4bb9afc4bfe81b6628219361d809ccfb21a99 1487cac605dec8032ef845d5de98575b1dc262302eb651b8823893e81d396acc 50f6d6ff383f442392e0b4482a484200469c8f04a9e1f9b5e21c66842f6e96c9 a670c9c61abd388f06a51a0d2d8542f68960fa728ab5133a36eef0b6c137a3be 4ba3bf0507efb2a98a1f1651d39af017666ca593e82430e888cee8619456f9fb 47d84a5c33b8b5ebee06f75d885c12073401a449f56c16aa64ed3aa62363f770 61bfedf72429b023d37d0d724d00a1248db0fead349f1c09b075372c980991b7 b25d479d8f6e8def7e3fe501ab6794c3b976ce0bd04c006bac1a94fb6409f60c 45e5c9ec2196a246368fb6faf3e6c53b51339b2eb3b52ec6f6dfc511f9b30952 ccc814544af5ebd09bee3d004de334afd660f2807192e4bb3c0cba85745c8740 fd20b5f39b9d3fbdb5579c0bd1a60320ad6a100c6402c7279679f25fefb1fa3c c8ea5e9f8db3222f83c7516dffd616b152f501ec8ad0552ab323db5fafd23876 053317b483e00df829e5c57bbca6f8ca01a87562edf1769dbd542a8f6287effc 3ac6732c68c4f5573695b27b0bbca58c8e1ffa35db8f011a010fa3d98fd2183b 84afcb56c2dd1d35b9a53e479b6f84565d28e49bc4bfb9790e1ddf2daa4cb7e3 362fb1341cee4c6e8ef20cada36774c01d07e9efe2bf11fb495dbda4dae90919 8eaad8e716b93d5a0d08ed1d0afc725e08e3c5b2f8e7594b78ff6e2fbf2122b6 48888b812900df01c4fad5ea0688fc31cd1cff191b3a8c1ad2f2f2218be0e177 7ea752dfe8b021fa1e5a0cc0fb56f74e818acf3d6ce89e299b4a84fe0fd13e0b 77cc43b81d2ada8d9165fa2668095770593cc7314211a1477e6ad206577b5fa8 6c75442f5fb9d35cfebcdaf0c7b3e89a0d6411bd3ae1e7e4900250e2d2071b35 e226800bb57b8e0af2464369bf009b91e5563911d59dfa6aa78c14389d95a537 f207d5ba202e5b9c5832603766295cfa911c819684e734a41b3472dca7b14a94 a1b5100529a532915d60f573fbc9bc6e42b60a47681e6740008ba6fb5571be91 ff296ec6b2a0dd915b6636521e7b9f9b6ff34052ec585566453b02d5da99f8fa 108ba47996e85076a4b7a70e9b5b32944db75092ec4192623ad6ea6b049a7df7 d9cee60b88fedb266ecaa8c71699a17ff5664526cc2b19ee1193602a575094c2 9a0591340e4183a3e3f54989a5b429d656b8fe4d699f73fd6a1d29c07efe830f 54d2d38e6f0255dc14cdd20868470eb266382e9c6021ecc5e09686b3f3ebaefc 93c9718146b6a70a1687f358452a0e286b79c5305aa5007373e07841c7fdeae5 c8e7d44ec5716f2b8b03ada37f0500c0df01c1f040200b3ffae0cf51a3cb574b 225837a58dc0921bdd19113f97ca92ff69432477322f547013ae5e58137c2dad cc8b576349af3dda7a94461460fd0030eecc8c73ea4751e41e238cd993bea0e2 f3280bba1183eb3314e548b384f6db9086f420d03f60a04bf2cb8129024977c7 95679b072bcaf89afde9a771fd9930810b38bae12dccf3f2e5512721f2e6b712 4501adde69f84cd877a5847187408da17bc9f9abce94b7d8cec7aec3adb851df a63094366c464c3d2ef1c18473215d908dd433b3724c2ba1612a14d432a65c45 150940002133ae4dd71dff89e10314e5581ac77d65f11199b043556f1d7a3c76 b3c11183b5924a509f28fe6ed97f1fbfa9ebabf2c1e153c6e86e34570eae96fb 1860e5e0a5a3e2ab3771fe71c4e3d06fa2965dcb999e71d0f803e89d65266c82 52e4cc9789c10b36ac6150eba94e2ea78a5fc3c531e0a2df4f2f74ea7361d2b3 d1939260f19c279605223a708f71312b6ebadfe6eeac31f66e3bc4595a67bc88 3b17f37d1018cff28c332ddefbe6c5aa56558218568ab9802eecea50fdb2f953 b2aef7dad5b6e2f841521b62829076170ecdd4775619f151013cca830eb61bd9 60334fe1eaa0363cfb5735c904c70a239d59e9e0bcbaade14eecc86bc60622ca 79cab5cabb2f3846e648b1eaf19bdf0caa02369b9655abb5040685a323c2ab4b 3319ee9d5c021b8f79b540b19875fa09995f7997e623d7da8f837889a97e32d7 711ed935f166812810e358829c7e61fd696dedfa17858ba9957f584a51b22726 39b83c3ff1ac24696cdb30aeb532e30548fd948e46dbc312858ebf2ef34c6ffe afe28ed61ee7c3c735d4a14d9e864b7e342105d14203e13e045eee2b6a3aaabe adb6c4f15facb4fd0c742f442ef6abbb5654f3b1d41cd2105d81e799e86854dc 7e44b476a3d816250cf62a1f25b8d2646fc8883a0c1c7b6a37f1524c369cb749 247848a0b5692b285095bbf00ad19489d1462b17423820e0058428d2a0c55f5e a1dadf43e233f70613372f0928d937e41d65fecf16c223bdb7cde3759cbee746 04085f2a7ce77326ea607808419f8509ee8efd85561d99735a969a7aac50c06c 25a04abfc800bcadc9e447a2ec3453484fdd567050e1e9ec9db73dbd3105588c d675fda79e3674340c5c43465713e38d83d28f89ef16dff20153e21e78fb03d4 ae6e39f2bdb83adf7e93d5a68948140f7f64c261c94692934411520f77602d4f 7bcf46b2ed4a20068d40824713320f46a43b7d4b7500061af1e39f62e9724454 614214f74bf8b88404d95fc1d96b591af70f4ddd366a02f45bfbc09ec03bd978 57fac6dd031cb850496eb27b355fd3941da2547e6abca0a9a28507825530429f 40a2c86dae9b66dfb68dc1462d7486900680ec0a427a18dee4f3ffea2e887ad8 cb58ce0067af4d6b6aace1e7cd3375fecce78a399406b2a4220fe9e35d9f385b 9ee39d7ab3b124e8b1dc9faf74b6d185626a36631eae397b23a6efa74dd5b433 26841e7f7ca7820fbfb0af54ed8feb397454056acba48952755533a3a20838d8 7fe6ba9b7d096954b55a867bca1159a58cca9296399e1db33a62a4a563f3125f 95ef47e1c9029317cfdf8e80204272f7080bb155c05282ce395c11548e4c66d2 248c1133fc70f86dc07f9c9ee41041f0f404779a45d886e17325f51ebd59bc0d 1f2bcc18f41113564257b7834602a9c60dff8e8a31f636c1b0e12b4c202e1329 eaf664fd1cad181156b2395e0333e92e13b240b62eebeb92285b2a20ee6ba0d9 9de720c8c2da2f728d012784595b794fd647d0862e7ccf5f05449a36f877d48f ac39dfd27f33e8d1e0a476341992eff743a6f6eabf4f8fd37a812dc60a1ebddf 8991be14cdb6e6b0dc67b55106d672c372765d43bdcd0e804f1290dc7cc00ffa 3b5390f92690fed0b667b9ffbcedb7d9ca091cf0bd9155ea3bb132f88515bad2 47b9479bf763bd6eb37392eb3cc1159798026e297f42e312d6842ada7c66a2b3 b12754ccc782ef11c6a124237b79251e706a1bbe64bfb63501a6b101811caedf a3d25bdd8e2e1c3c9444216590a121386d90cec6ed5abea2a64af674eda86a85 fbebfe98864e4c3fe9dbc8057f0f7c08660787bf86003604dd1fd8346f6381fb 07745ae04d736fccc83426b33f01eab71b08041873c005e5f77a057bebde8ae2 455464299bf582e614e58f48ff2ddfda2f474ef388789bdc25366f9c3c8b38e7 4b475f25546fcd9b97aeb26618b1ddf84846a0e79915f95e2466e598e20b4577 08cd55591c902de4cb90bace1bb8205d011a862487574a99eb77f19b6e0a9dc0 9662d09a1c4324633e85a1f0209f0be8c4a99a0251d6efe101ab93d1d0ba5a4d fa186f20f2868f169dcb7da83573906fea1e2ce9b4fcd7f5250115e01a70683f aa002b5c40de6d0279af88c27773f8641c3604c0661a806b5f0177a28c0f586e 0006058aa30dc7d6211e69ed72338ea6353c2dd94c2c21634bbcbee5690bcb6d eebfc7da1ce591d766f05e4094b7c018839720a3d7c927c2486e3725f724d9db 91ac15bb4d39eb8fced54557808fca5b5d83d7cd34dad0fc41e50ef5eb161e6f 8a28514d96c51133c6fd5c7e756e14ec4362abfceddc6c837d79a32349263821 2670efa8e406000e03a39ce37d3faf5cfabc277375ac52d1b5cb0679e4fa3374 2d382274099bc9bbed5118e9dbf0f7315d62d1c7ec700c47bb78c1b6b21a1904 5b26eb1be6a366eb45748ab2fbc946e79c6a376d26549c2c8530ff8ee468dde7 dd5730a1d4cd04dc62939bbdba9ba4650ac9526e8be5ee304a1fad5f06a2d519 a63ef8ce29a86ee22c089c2b843242ef6a51e03aa9cf2d0a483c061ba9be96a4 d8fe51550ba645bd62826a2f9a73a3ae14ba99586ef5562e9c72fefd3f752f7d a3f046f6977fa0a5980e4a91587b086019b09e6ad3b3ee593e990fd5a9e34d79 72cf0b7d9022b8b5196d5ac3a017da67dd1cf3ed67c7d2d281f9f25cfadf2b89 b5ad6b4725a88f54ce029ac71e019a5e647b0acfded93fa9be8d3c48d283b57c cf8d5662979132e28785f0191ed756055f7960e44e3d35e8c15056dd488f46db a03a161250564f0bdc3eb9e153c9057a297271aeca93a072a1b3f6d9b1e6321f 5f59c66fb26dcf3197533d928b155fdf5035634828aba3cbb28517711c20ad9f 8abcc5167ccad925f4de817513830dc8e379d58629320f991ea7a90c2fb3e7bc e5121ce64774fbe32a8b6e37ec3293d4648de53696413e680a2ae0810dd6db22 469852dfd09072166b39a460a6445c0dd586cdecf1c20c8ae5bbef7dd1b588d4 0ccd2017f6bb4e3bbdda26a7e3a59ff453e350a44bcb4cdd572eacea8fa6484b b8d6612aebf3c6f47d29be463542f5d9eaec2771bf64e6370740e0d8de75b135 7f8721671af537d5d4040cb084eb4e2cc34d2466a0115af84e1b0042895983a1 d06b89fb4ce6ea0486f3f3b823520ab82011a1d4b277227f8611560b1e7933fd cbb3a792b344525bda08839e151ce794b2f32c9b7a01fbac9e01cc87ebcc7d1f 6cf0111c3a1e8aac71a908749d44fbd9ad0dadecbd50ada380339c32ac691366 78df9317ce0b12b4ff79e59b743f5bb3af2d519ff27d9459cbf97222c15e6fc2 a0f91fc719b941525fae59361ceb69cebc2a8645912baa8d1b6c1075ee3056a0 c10d25065cb03a442e0ec6e0e1698db3b4c98a0be3278e9649f1f9532e0d392d fd3a0342b8971f21e1b0a74414ba3348cc5be7120c37632d8df359f8d9b992f2 ee60b6f470fe3f11de54cda541edad891ce6279cfcd3e7e6f1618b166fd2c1d0 5848fd2c5f6fb2299f523f357a632762393a8353156cccd02acf081625a75ebb 56e16369788d273ccde96629281b949d04c50901b71c65614e6c6c7bd327a140 a45e1d006c3f27b9ac9aa53fd62a80f00bb25bfe235bdd2f671126905b204022 2b6cbcf7ccd769c2b53113ec01640e3d338abbd602547adf0ba38209cf746ce7 677afa1c52075606085cbfe4e8ae88dd87aaaf9b04cf9aa7e1948c25c02fb8a8 c01c36ae4d6ebe1f990d4f869a65cdea03f09252dc208e69fb74e6132ce77e25` .split('\n') .join('') })38:[] 32:{"id":"4d5f695dcfb9d8a4","slug":"computing-with-the-bailey-borwein-plouffe-formula","title":"Computing π with the Bailey-Borwein-Plouffe Formula","thumbnail":"e86f7408a08609c4ed1bdfacc2d0645e7c57b939827f65c847ac2a571dafa781","access_level":"public","create_time":"$D2021-05-04T03:18:13.848Z","update_time":"$D2026-01-13T17:25:53.600Z","view_time":null,"first_public_time":"$D2021-05-22T05:13:19.347Z","has_unpublished_edits":null,"fork":false,"stars":20,"starred":false,"owner":{"type":"team","id":"f1b8352b963a77b4","login":"rreusser","name":"Ricky Reusser","avatar_url":"https://avatars-next.observableusercontent.com/49c18a2ac0092176d4ea8348a1dd538e0b460b086b89d602b8d21ffe0821d92e","tier":"starter_2024"},"version":1832,"publish_version":null,"role":null,"authored":false,"authors":[{"type":"individual","id":"908c4be2279a714e","login":"rreusser","name":"Ricky Reusser","avatar_url":"https://avatars-next.observableusercontent.com/fa41989f47bb986723af5b9c22176cb503748426f56e565000ea93acd9eab86c","tier":null}],"members":[{"type":"individual","id":"908c4be2279a714e","login":"rreusser","name":"Ricky Reusser","avatar_url":"https://avatars-next.observableusercontent.com/fa41989f47bb986723af5b9c22176cb503748426f56e565000ea93acd9eab86c","tier":null,"role":"editor"}],"comments":[{"id":"b3117760a84493d7","reply_to_id":null,"content":"this visualisation is beautiful !","create_time":"$D2021-05-22T13:32:47.198Z","update_time":null,"resolved":false,"user":{"type":"individual","id":"7affeb11879b11eb","login":"stardisblue","name":"Fati CHEN","avatar_url":"https://avatars-next.observableusercontent.com/c3f14cf2d58aebca19459fc111e35113235be85745ae353f37b7731ea9a84d79","tier":null}}],"custom_thumbnail":true,"trashed":false,"lapsed":null,"subscribed":null,"forks":3,"collection":{"id":"dcd095ad05c2b722","slug":"writeups","title":"Fun Notebooks","owner":{"id":"f1b8352b963a77b4","login":"rreusser","name":"Ricky Reusser"}},"type":"notebook","body":{"title":"Computing π with the Bailey-Borwein-Plouffe Formula","cells":[{"id":301,"value":"md`# Computing π with the Bailey-Borwein-Plouffe Formula`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":1830,"value":"*[Read this notebook ➡️](https://rreusser.github.io/notebooks/computing-pi-with-the-bailey-borwein-plouffe-formula/) on my personal website via [Observable Notebook Kit](https://observablehq.com/notebook-kit/kit).*","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined","output":"","database":"$undefined"},{"id":1220,"value":"tex.block` \\pi =\\sum _{k=0}^{\\infty }\\left[{\\frac {1}{16^{k}}}\\left({\\frac {4}{8k+1}}-{\\frac {2}{8k+4}}-{\\frac {1}{8k+5}}-{\\frac {1}{8k+6}}\\right)\\right].`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":1224,"value":"$34","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":1193,"value":"viewof dummy = {\n if (!this) {\n yield html`Click to compute π starting at the millionth hexadecimal digit.`;\n return;\n }\n recompute;\n yield html`Computing…`;\n const t1 = performance.now();\n const digits = 10;\n const ans = computePi(1e6, digits);\n const t2 = performance.now();\n yield html`

\n Hexadecimal digits of π starting at the millionth (computed in ${(\n (t2 - t1) /\n 1000\n ).toFixed(2)} seconds):
${ans\n .toString(16)\n .padStart(digits, \"0\")}\n

`;\n}","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":1242,"value":"viewof recompute = Button('Compute')","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":1586,"value":"md`The table below shows hexadecimal digits of ${tex`\\pi`}. Each line represents an independent run of the algorithm starting at a slightly different digit. Green digits are correctly computed digits, while red represents incorrect digits. Due to our termination criteria, we usually get a couple fewer digits than requested.\n\n${\n supportsBigInt\n ? `\nChange from double precision floating point to \\`BigInt\\` and observe that, although starting at the same digit is at least a factor of ten slower (*update: BigInts have gotten much, much faster!*), it doesn't cost too much extra to compute many more digits.\n`\n : `\nYour browser does not support \\`BigInt\\`, but on a device which supports \\`BigInt\\`, you may compare floating point performance to an equivalent computation using \\`BigInt\\`. Although starting at the same digit is at least a factor of ten slower (*update: BigInts have gotten much, much faster!*), it doesn't cost too much extra to compute many more digits.\n`\n}`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":805,"value":"viewof method = {\n if (supportsBigInt) {\n return Radio(['Floating point', 'BigInt'], {\n label: 'Data type',\n value: 'Floating point'\n });\n } else {\n const ret = html``;\n ret.value = 'Floating point';\n return ret;\n }\n}","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":333,"value":"viewof startingDigit = Range([0, reference.digits.length - 128], {\n value: 0,\n label: 'Starting digit',\n step: 1\n})","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":437,"value":"viewof digits = Range([2, method === 'Floating point' ? 14 : 128], {\n label: 'Digits to compute',\n value: method === 'Floating point' ? 14 : 128,\n step: 1\n})","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":478,"value":"$35","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":1460,"value":"$36","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":1480,"value":"md`## Modular exponentiation\n\nWorking out an implementation requires a small digression into modular exponentiation. For a start, let's consider *binary exponentiation*. Consider evaluating\n\n${tex.block`\n3^{17} = 3 \\cdot 3 \\cdot 3 \\cdot 3 \\cdot 3 \\cdot 3 \\cdot 3 \\cdot 3 \\cdot 3 \\cdot 3 \\cdot 3 \\cdot 3 \\cdot 3 \\cdot 3 \\cdot 3 \\cdot 3 \\cdot 3.\n`}\n\nWe can multiply all the ${tex`3`}'s, but it's rather wasteful. Instead, observe that\n\n${tex.block`\n\\begin{aligned}\n3^{17} &= 3^{16} \\cdot 3 \\\\\n&= (3^8)^2 \\cdot 3 \\\\\n&= ((3^4)^2)^2 \\cdot 3 \\\\\n&= (((3^2)^2)^2)^2 \\cdot 3 \\\\\n\\end{aligned}\n`}\n\nBy repeatedly squaring, we can immediately reduce the number of multiplications from sixteen to five. As an algorithm, we repeatedly square the base as we step through the bits of the exponent. Here, the binary representation of ${tex`17`} is ${tex`10001_2`}. If the ${tex`n^{th}`} bit is a 1, we multiply by the ${tex`n`}-times squared base.`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":1498,"value":"function binaryPow(b, e) {\n let y = 1;\n while (e > 0) {\n if (e & 1) y *= b;\n b *= b;\n e >>= 1;\n }\n return y;\n}","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":1685,"value":"binaryPow(3, 17)","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":1017,"value":"md`\n*[Modular exponentiation](https://en.wikipedia.org/wiki/Modular_exponentiation)* computes the exponentiation ${tex`\\operatorname{mod} m`},\n\n${tex.block`\nb^e \\operatorname{mod} m.\n`}\n\nWhile ${tex`b^e`} may grow very large, ${tex`b^e \\operatorname{mod} m`} will always be less than ${tex`m`}. An efficient algorithm, based on binary exponentiation, makes use of [the fact that](https://en.wikipedia.org/wiki/Modular_exponentiation#Memory-efficient_method)\n\n${tex.block`\n(a \\cdot b) \\operatorname{mod} m = [(a \\operatorname{mod} m) \\cdot (b \\operatorname{mod} m)] \\operatorname{mod} m.\n`}\n\nPut simply, we just load the \\`binaryPow\\` function up with ${tex`\\operatorname{mod} m`} and keep our result between ${tex`0`} and ${tex`m`}.`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":832,"value":"function modPow(b, e, m) {\n b = b % m;\n let y = 1;\n while (e > 0) {\n if (e & 1) y = (y * b) % m;\n b = (b * b) % m;\n e >>= 1;\n }\n return y;\n}","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":1700,"value":"md`Consider ${tex`3^{17} \\operatorname{mod} 50`}:`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":1695,"value":"modPow(3, 17, 50)","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":1706,"value":"md`Compare with`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":1697,"value":"3 ** 17 % 50","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":1487,"value":"md`\nAs one final prelude, we implement an alternative to the JavaScript \\`%\\` operator since its behavior with negative numbers is somewhat undesirable. \\`-1.6 % 1\\` is equal to ${tex`-0.6`}, but we prefer the range ${tex`[0, 1)`} instead.\n`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":829,"value":"function mod(x, n) {\n x %= n;\n return x < 0 ? x + n : x;\n}","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":1501,"value":"md`## Implementation\n\n\nThe function below directly implements the above summation ${tex`S_j`} and the function \\`piBBP(d, n)\\` which starts at digit ${tex`d`} and computes ${tex`n`} digits of ${tex`\\pi`}. It returns the result as an integer JavaScript \\`Number\\` so that we'll need to run \\`.toString(16)\\` to convert the output to hexadecimal.\n`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":792,"value":"computePi = {\n function S(j, n) {\n // Lefthand summation. Add a number of terms equal\n // to the starting digit. Most of the time is spent\n // in this loop.\n let left = 0;\n for (let k = 0; k <= n; k++) {\n let r = 8 * k + j;\n left = mod(left + modPow(16, n - k, r) / r, 1);\n }\n\n // Righthand summation. Add a small number of terms\n // until the result stops changing.\n let right = 0;\n for (let k = n + 1; ; k++) {\n let rnew = right + 16 ** (n - k) / (8 * k + j);\n if (right === rnew) break;\n right = rnew;\n }\n\n return left + right;\n }\n\n // Start at digit d and compute n digits\n return function piBBP(d, n) {\n // Seems to be the convention for including the leading 3\n d -= 1;\n // Shift n digits to the left of the radix point to obtain our final\n // result as a integer\n return Math.floor(\n 16 ** n * mod(4 * S(1, d) - 2 * S(4, d) - S(5, d) - S(6, d), 1)\n );\n };\n}","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":1768,"value":"computePi(0, 10).toString(16)","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":1512,"value":"md`## BigInt implementation`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":1565,"value":"md`With minor modifications, we can repeat the above computation with integer arithmetic. The code below repeats the above implementation, except that everything is written with [BigInt](https://developer.mozilla.org/en-US/docs/Web/JavaScript/Reference/Global_Objects/BigInt)s rather than double precision numbers. This means we can compute an arbitrary number of digits at once.\n\nA small downside is that we need to generate the code dynamically so that we don't pollute all downstream Observable cells with syntax errors in the event that \\`BigInt\\` is not supported.`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":1347,"value":"supportsBigInt = {\n try {\n new Function('return 0n;');\n return true;\n } catch (e) {\n return false;\n }\n}","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":0,"value":"BBPbigint = {\n if (!supportsBigInt) return () => () => 0;\n\n return new Function(`\nfunction modPow(b, e, m) {\n if (!e) return 1n;\n b = b % m;\n let y = 1n;\n while (true) {\n if (e & 1n) y = (y * b) % m;\n e >>= 1n;\n if (!e) return y;\n b = (b * b) % m;\n }\n}\n\nfunction S(j, n, d, mask) {\n const shift = d << 2n;\n let left = 0n;\n for (let k = 0n; k <= n; k++) {\n let r = k * 8n + j;\n left = (left + (modPow(16n, n - k, r) << shift) / r) & mask;\n }\n let right = 0n;\n for (let k = n + 1n; ; k++) {\n let rnew = right + 16n ** (d + n - k) / (k * 8n + j);\n if (right === rnew) break;\n right = rnew;\n }\n return left + right;\n}\n\nreturn function computePi(dd, nn) {\n const n = BigInt(nn);\n const d = BigInt(dd) - 1n;\n const mask = 16n ** n - 1n;\n return (\n 4n * S(1n, d, n, mask) -\n 2n * S(4n, d, n, mask) -\n S(5n, d, n, mask) -\n S(6n, d, n, mask)\n ) & mask;\n};`)();\n}","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":1739,"value":"md`The article from literateprograms.org, [Pi with the BBP formula](https://literateprograms.org/pi_with_the_bbp_formula__python_.html), notes that there's an additional, more efficient variation that generates digits sequentially, but since it requires a rational data type, I've not implemented it.`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":1521,"value":"md`## References\n\n

${md`\n0. Bailey, D. H., *[The BBP Algorithm for Pi](https://www.davidhbailey.com/dhbpapers/bbp-alg.pdf).*\n0. Bailey, D. H., [PSLQ: An Algorithm to Discover Integer Relations](https://www.davidhbailey.com/dhbpapers/pslq-comp-alg.pdf).\n0. *[Pi with the BBP formula](https://literateprograms.org/pi_with_the_bbp_formula__python_.html)*, literateprograms.org.\n0. *[Modular exponentiation](https://en.wikipedia.org/wiki/Modular_exponentiation)*, wikipedia.org.\n`}`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":1515,"value":"md`## Imports`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":150,"value":"$37","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":285,"value":"import { Button, Range, Radio } from '@observablehq/inputs'","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":1553,"value":"html``","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"}],"files":"$W38","stdlib":"1","titling":"auto"}} 33:"$32:body"