The distributions of the integers of a continued fraction given by Gauss-Kuzmin is different from I had worked out some years ago. They give P(an = k) = −ln(1 − (k + 1)−2). I got the simpler result P(an = k) = 6/(πk)2
Guess what! Gauss was right. This program contrasts my distribution with Gauss’:
#include <stdio.h> #include <math.h> typedef double f; static f sq(f x){return x*x;} static f lg(f x){return log(x)/log(2);} int main(){f ms=0, gs=0; f pi = 3.141592653589793238; f nrm = pi*pi/6; int j; for(j = 1; j<20000; ++j) {f m = 1/(sq(j)*nrm), g = -lg(1-1/sq(j+1)); ms += m; gs += g; if(j<21) printf("%d %19.16f %19.16f\n", j, m, g);} printf("%19.16f %19.16f\n", ms, gs); return 0;}which yields:
1 0.6079271018540267 0.4150374992788438 2 0.1519817754635067 0.1699250014423125 3 0.0675474557615585 0.0931094043914815 4 0.0379954438658767 0.0588936890535686 5 0.0243170840741611 0.0406419844973459 6 0.0168868639403896 0.0297473433940521 7 0.0124066755480414 0.0227200765000835 8 0.0094988609664692 0.0179219079972625 9 0.0075052728623954 0.0144995696951151 10 0.0060792710185403 0.0119726416660759 11 0.0050241909244134 0.0100536646639229 12 0.0042217159850974 0.0085620135034241 13 0.0035972017861185 0.0073795303655975 14 0.0031016688870103 0.0064262691594330 15 0.0027018982304623 0.0056465631411421 16 0.0023747152416173 0.0050006810583665 17 0.0021035539856541 0.0044596481906998 18 0.0018763182155988 0.0040019305574962 19 0.0016840085923934 0.0036112535523788 20 0.0015198177546351 0.0032751320328610 0.9999696028849843 0.9999278670512570and this Algol68 program, using 4000 digits of precision:
INT lm = 100; INT big := 0; MODE B = LONG LONG REAL; [lm] INT h; FOR j TO lm DO h[j] := 0 OD; B x := long long pi; TO 10000 DO INT q = SHORTEN SHORTEN ENTIER x; IF q <= lm THEN h[q] +:= 1 ELSE big +:= 1 FI; x := 1/(x-q) OD; FOR j TO lm DO print((j, h[j], newline)) OD; print(big)considers the first 10000 terms for π and reports
+1 +4182 +2 +1675 +3 +939 +4 +595 +5 +404 +6 +286 +7 +254 +8 +181 +9 +123 +10 +129 +11 +77 +12 +88 +13 +76 +14 +71 +15 +58 +16 +46 +17 +43 +18 +43 +19 +30 +20 +31 +21 +26 +22 +22 +23 +30 +24 +28 +25 +17 +26 +18 +27 +14 +28 +11 +29 +17 +30 +17 +31 +12 +32 +20 +33 +14 +34 +10 +35 +8 +36 +14 +37 +12 +38 +4 +39 +8 +40 +9 +41 +10 +42 +9 +43 +9 +44 +5 +45 +6 +46 +5 +47 +11 +48 +7 +49 +3 +50 +11 +51 +3 +52 +6 +53 +11 +54 +5 +55 +10 +56 +3 +57 +6 +58 +2 +59 +6 +60 +5 +61 +4 +62 +4 +63 +6 +64 +1 +65 +4 +66 +4 +67 +4 +68 +2 +69 +1 +70 +3 +71 +1 +72 +3 +73 +7 +74 +1 +75 +0 +76 +0 +77 +2 +78 +2 +79 +2 +80 +5 +81 +0 +82 +0 +83 +0 +84 +6 +85 +1 +86 +2 +87 +1 +88 +4 +89 +0 +90 +0 +91 +1 +92 +2 +93 +2 +94 +4 +95 +3 +96 +0 +97 +1 +98 +1 +99 +2 +100 +0 +139