同時に投げて、やはり裏が出たコインは取り除いていく。以下、この試行を全てのコインが
取り除かれるまで繰り返すことにする。
さてこのとき、終了するまでの試行回数の期待値は如何ほど？また、４枚のコインで始め
たら、この期待値は如何ほどに？
らすかるさんからのコメントです。（平成２９年１月１３日付け）
３枚の場合、22/7回、４枚の場合、368/105回かな？
ＧＡＩ さんからのコメントです。（平成２９年１月１３日付け）
流石、らすかるさん、はやくて正確。３枚で、22/7（円周率の近似分数）だったので、もしや
と思って、４枚にあのややこしい漸化式を解いていったら、違っていたので少しがっかりしま
した。
らすかるさんからのコメントです。（平成２９年１月１３日付け）
枚数が非常に多い時、1回で約半分に減りますので、ｎ枚のとき、log2ｎ 近辺になり、
(枚数)→∞のとき、(回数)→∞ですね。
ＧＡＩ さんからのコメントです。（平成２９年１月１３日付け）
ははは、そうでしたね。100枚でも、log2 100 = 6.6438･･･ なので、７回も繰り返すと終了か！
以外と速いな～。（数学的論理値とどの位ちがうのか？）
ちなみに10枚なら、log2 10 = 3.3219･･･ に対し、論理値は、4.7255･･･なのでそう違いはな
いが．．．。
らすかるさんからのコメントです。（平成２９年１月１３日付け）
理論値は、
10枚では、4.7255593236…回(log2 10 ≒ 3.322)
100枚では、7.9838015351…回(log2 100 ≒ 6.644)
1000枚では、11.2992526972…回(log2 1000 ≒ 9.966)
となるようです。
DD++さんからのコメントです。（平成２９年１月１４日付け）
limn→∞ (理論値-log2 ｎ) が正定数に収束しそうに見えますが、どうなんですかね？
らすかるさんからのコメントです。（平成２９年１月１４日付け）
もう少し先まで計算できれば予想できるかも知れませんが、今の計算方法（ｎ枚のときの値
を求めるのに、ｎ-1枚までの結果を使う）では、100000枚も厳しいですね。
# 1000枚までの結果を見た感じでは、log2 2.5ぐらい？
DD++さんからのコメントです。（平成２９年１月１４日付け）
Σk=1～n (-1)^(k-1) C(n,k)/(1-1/2^k) で計算できます。
３枚の場合、3/(1/2) - 3/(3/4) + 1/(7/8) = 22/7
４枚の場合、4/(1/2) - 6/(3/4) + 4/(7/8) - 1/(15/16) = 368/105
さて、ｎ→∞ を考えるときに、これをどう処理するか。
（追伸） 1 + Σk=1～n (-1)^(k-1) C(n,k)/(2^k-1) の方が楽ですかね。
３枚の場合、1 + 3/1 - 3/3 + 1/7 = 22/7
４枚の場合、1 + 4/1 - 6/3 + 4/7 - 1/15 = 368/105
らすかるさんからのコメントです。（平成２９年１月１５日付け）
d(n)=(理論値(n))-log2 ｎ は、
d(10) = 1.40363122874716549508…
d(100) = 1.33994534538219529257…
d(1000) = 1.33346841261711447858…
d(10000) = 1.33281901306011801337…
となっていて、
d(10)-d(100) = 0.0636858…
d(100)-d(1000) = 0.00647693…
d(1000)-d(10000) = 0.000649399…
なので、
d(10000)-d(100000) ≒ 0.000065
d(100000)-d(1000000) ≒ 0.0000065
d(1000000)-d(10000000) ≒ 0.00000065
・・・
と予想され、収束値は、
d(10000)-0.000065-0.0000065-0.00000065-… = 1.332746…
あたりになりそうですが、この値は、
γ/log２ + 1/2 = 1.3327461772… （γはオイラー定数 0.5772156…）
に非常に近いですね。
ＧＡＩ さんからのコメントです。（平成２９年１月１５日付け）
一晩かけて、コイン数ｎ枚、理論終了回数期待値-log2 ｎ を計算させてみました。
1000;1.333468412617114478586587147
2000;1.333107805555418028654596957
3000;1.332985130663444298927243621
4000;1.332927479707884870226526403
5000;1.332891146446736095929558592
6000;1.332864926656726407210249487
7000;1.332848394599154836112717789
8000;1.332837311204125504784067225
9000;1.332827886899555623592396525
10000;1.332819013060118013377739250
これ以上のｎは、とても時間がかかりすぎていきます。らすかるさんが、これらのデータから
γ/log２+1/2 を推測されていて、これってほとんど一致しそうですね。
オイラーは、このγを、
γ=1/2(1+1/2^2+1/3^2+1/4^2+･･･)-1/3(1+1/2^3+1/3^3+1/4^3+･･･)+1/4(1+1/2^4+1/3^4+1/4^4+･･･)+･････
ヴァッカは、
γ=1/2-1/3+2(1/4-1/5+1/6-1/7)+3(1/8-1/9+･･･-1/15)+4(1/16-1/17+･･･-1/31)+･････
ポリアは、
γ=1-(1/2+1/3)+3/4-(1/5+1/6+1/7+1/8)+5/9-(1/10+1/11+1/12+1/13+1/14+1/15)+7/16-+･････
メルテンスは、
γ=lim[n→∞](Σ[p=2,p≦n]log(p/(p-1))-log(log(n))) (pはnを越えない素数)
と表しているのを重ねて見てみると、不思議さと神秘さが感じられます。無限の世界は奥が
深い！
改めて、1000刻みで階差を見てみると、
d(1000)-d(2000)=.000360607061697
d(2000)-d(3000)=.000122674891974
d(3000)-d(4000)=.000057650955559
d(4000)-d(5000)=.000036333261149
d(5000)-d(6000)=.000026219790010
d(6000)-d(7000)=.000016532057571
d(7000)-d(8000)=.000011083395029
d(9000)-d(8000)=.000009424304570
d(10000)-d(9000)=.000008873839437
更に階差を調べていくとかなり変動していることが見て取れ（プラス、マイナスが混じります）
ますが、極限値として、
γ/log2 + 1/2
は非常に近いですね。らすかるさんの計算からは、
? 1.332819013060118013377739250-65/9/10^5
%42 = 1.332746790837895791155517028
となりますが、
? gamma/log(2)+1/2
%43 = 1.332746177276867150646417519
なので、小数第6位までは一致させられますが、上記の揺れから極限値としてこれを保証す
るにはなりませんね。
DD++さんからのコメントです。（平成２９年１月１５日付け）
(理論値(n))-log2 (n+1/2) くらいだと精度良く収束しそうな予感。一体この1/2が何を意味
するかと言われてもさっぱりですが。
らすかるさんからのコメントです。（平成２９年１月１５日付け）
1日計算して、n=100000 の理論値が出ましたので再度まとめます。
以下、n=10、100、1000、10000、100000の値です。
(理論値(n))
4.725559323634527842953987445277…
7.983801535156919988312681262704…
11.299252697279201522197545434958…
14.620531392609567404859016967646…
17.942392291438066655235560413171…
上記の階差は
3.258242211522392145358693817426…
3.315451162122281533884864172254…
3.321278695330365882661471532687…
3.321860898828499250376543445525…
であり、順調に、log2 10=3.32192809…に収束しているように見えます。
(理論値(n))-log2 n は、
1.403631228747165495083668015787…
1.339945345382195292572042403725…
1.333468412617114478586587146490…
1.332819013060118013377739249688…
1.332751817001254915883963265724…
で、こちらも順調に1.3327461772…(=γ/log2+1/2)に収束しているように見えます。
（しかし階差をとると、次は1.3327461772より小さくなるようにも思えます）
(理論値(n))-log2 (n+1/2) は、
1.333241900855767554058279184097…
1.332749843977991373929665309988…
1.332747245373460347729884922219…
1.332746880111382256279325975885…
1.332744603544084098965970558518…
となり、もし収束値が1.3327461772…ならば下に超えてしまっていますので、log[2](n+1/2)よ
りも log2 n の方が良い近似なのかも知れません。
ちなみに、Σ(1/(2^k-1))とか計算してみたら、分母にlog2、分子にq-digamma関数が出て
きましたので、γ/log2+1/2 という予想が当たっている可能性が少し高まった気がします。
余談ですが、n=100000の理論値を計算していて、「後ろn/5項は計算する必要がない」こと
に気付きました。
nCk/(2^k-1)は、k=n/3あたりで最大となり、n=100000、k=33333の場合は、約3.6×10^17606
という巨大な数になります。（→結果が18付近なので少なくとも17700桁以上の計算が必要）
それに対して、n=100000、k=78000では、約2.9×10^(-600)という微小な数になり、これ以降
も小さくなる一方ですから、100桁以内の有効数字を求めるには計算不要です。
らすかるさんからのコメントです。（平成２９年１月１６日付け）
どうしても先が知りたくて、とうとう専用プログラムを作ってしまいました。専用に高速化して
いますので、n=1000は一瞬、n=10000は1秒、n=100000は11秒、n=1000000は20分で求まり
ます。
n=100000までは今までの結果と一致していました。n=1000000の理論値は、
21.264316141955791815486195154617…
理論値(1000000)-理論値(100000)（階差）は、
3.321923850517725160250634741445… (log2 10=3.3219280948…)
理論値(1000000)-log2 1000000 は、
1.332747572631617728264278577680… (γ/log2+1/2=1.3327461772…)
理論値(1000000)-log2 (1000000+1/2) は、
1.332746851284277620602573747316… (100000の結果より増加)
でした。このプログラムならば、n=10000000も2日あれば求められますので、現在計算中です。
また、
1 + Σ[k=1..n] (-1)^(k-1) C(n,k)/(2^k-1)
= 1 + Σ[k=1..∞] 1-(1-1/2^k)^n
（ = 1 + {1-(1/2)^n} + {1-(3/4)^n} + {1-(7/8)^n} + {1-(15/16)^n} + …）
と変形できることに気付いたのですが、こう変形できてもこの先は計算できませんよね…。
DD++さんからのコメントです。（平成２９年１月１６日付け）
むしろ、あまりになんともならないので、二項展開した残骸としてComb(n,k)や(-1)^(k-1)が
出現しているわけなんです。機械計算の場合、展開前の方が計算しやすいんでしょうか？
らすかるさんからのコメントです。（平成２９年１月１６日付け）
作ったプログラムでは、
nCk/(2^k-1)=nCk/2^k+nCk/2^(2k)+nCk/2^(3k)+nCk/2^(4k)+…
と変形して計算していますので、この変形から気付きました。２進計算の場合は、左辺より
右辺の方が簡単かつ高速になります。
ａｔ さんからのコメントです。（平成２９年１月１６日付け）
「ｎ枚のときの値を求めるのに、ｎ－１枚までの結果を使う」というのは、具体的にはどのよ
うな式になるのでしょうか？よろしければ教えてください。
DD++さんが導かれた期待値の理論式は、次のように考えても導出できますね。
ｎ個のコインをｍ回投げたとき、少なくとも一枚のコインが取り除かれないで残っている確
率をQ(ｍ，ｎ)とすると、包除原理により、
Q(ｍ，ｎ)=∑[k=1,n]C(n,k)*(-1)^(k+1)*((1/2)^m)^k.
ｎ個のコイン全てが取り除かれるまでに要する試行回数の期待値は、
∑[m=0,∞]Q(ｍ，ｎ)
=∑[m=0,∞]∑[k=1,n]C(n,k)*(-1)^(k+1)*((1/2)^m)^k
=∑[k=1,n]∑[m=0,∞]C(n,k)*(-1)^(k+1)*((1/2)^k)^m
=∑[k=1,n]C(n,k)*(-1)^(k+1)/(1-(1/2)^k)
=1+∑[k=1,n]C(n,k)*(-1)^(k+1)/(2^k-1)
らすかるさんからのコメントです。（平成２９年１月１６日付け）
ｎ枚の時の試行回数をａ[ｎ]とし、ａ[0]=1 と定義すると、
a[n]={Σ[k=0～n-1]nCk・a[k]}/(2^n-1)+1
という式になりました。
ａｔ さんからのコメントです。（平成２９年１月１６日付け）
教えていただき、ありがとうございます。
DD++さんからのコメントです。（平成２９年１月１７日付け）
２進計算の場合は、左辺より右辺の方が簡単かつ高速になります。
なるほど、ビットをずらすだけで済むわけですか。
ところで、表が出る（取り除かれない）確率が 1/2 でなく p （ただし 0<p<1） の場合、
(期待値-log[1/p](n)) はどうも 1/2 - γ/log(p) になりそうに見えます。
つまり、(期待値 + (γ+log(n))/log(p) ) は p の値によらず、1/2に収束するのではないかと
いう予想。つまり、これをpで偏微分したら恒等的に0なのではないか、と思ったりしましたが、
これはこれで証明大変そう……。
らすかるさんからのコメントです。（平成２９年１月１７日付け）
機械計算の場合展開前の方が計算しやすいんでしょうか？
この貴重なコメントをスルーしていましたが、よく考えてみると、
1 + Σ[k=1..∞] 1-(1-1/2^k)^n
という式の方が圧倒的に速いことに気付きました。nが10^50などの巨大な値でも、数秒で求
まります。そして、収束値が得られるのを期待して10^50まで計算してみたのですが、実に予
想外の結果になりました。以下、(理論値(10^k))-log[2](10^k) （k=1～50）の結果(小数点以
下30桁のみ)です。
この結果を見る限り、「収束しない」ように思えます。プログラムに問題があるということも
考えましたが、全く異なる式・プログラムで計算した10^6までの結果と少なくとも小数点以下
60桁以上は一致していましたので、おそらくプログラムの問題ではないと思います。
これは予想ですが、元々式が2のべき乗に依存していますので、2のべき乗に近いかどう
かによって結果が変わるのかも知れません。(つづく)
DD++さんからのコメントです。（平成２９年１月１７日付け）
この結果を見る限り、「収束しない」ように思えます。
なんと、これは驚き。10^19 から 10^50 までで、ほぼちょうど10回増減しているように見え
るので、sin(2log[10](n)) みたいな項が微小係数で存在するんでしょうか？
これは予想ですが、元々式が2のべき乗に依存していますので、2のべき乗に近いかどう
かによって結果が変わるのかも知れません。
これが正しいかどうかはコインの枚数を 10^n じゃなく 2^n で増やしてみれば検証できそう
ですね。
らすかるさんからのコメントです。（平成２９年１月１７日付け）
そこらへんをいろいろ調べていました。n=10^kではなくn=2^kとしたらどうなるかをやってみ
ると、(以下はn=2^10,2^20,2^30,…,2^100)
1.333451762795894434203328699982…
1.332748070414984732499692226281…
1.332747383104756835916276544981…
1.332747382433555335956853765744…
1.332747382432899865741949001465…
1.332747382432899225634317258435…
1.332747382432899225009212149311…
1.332747382432899225008601695103…
1.332747382432899225008601098956…
1.332747382432899225008601098374…
のように見事に収束します。
(多いので間は省略しましたが、2^1～2^100まで綺麗に単調減少しています。)
この収束値 1.332747382432899225008601098… は、「A158468」の値です。
そして、
n=3・2^kの場合 1.332744625291612590542519660… に収束
n=5・2^kの場合 1.332746560607882545186885743… に収束
n=7・2^kの場合 1.332745656041297660877197235… に収束
n=2^k-1,2^k+1の場合 n=2^kのときと同じ値に収束
などのようになります。やはり「収束しない」が正解のようです。γ/log2+1/2の値は上記収束
値の間ですから、もしかしたら、γ/log2+1/2に収束する意味のある数列があるのかも知れま
せんね。
ＧＡＩ さんからのコメントです。（平成２９年１月１７日付け）
結論は出たようですが、偏微分は不可能なので、期待値（E）は、n、pの関数よりE(n,p)とし
て、
P(n,p)=E(n,p)+(γ+log(n))/log(p)
の数値計算を、私の非力なプログラム力で実行可能な n=5000 までの1000刻みでやってみ
ました。（１６時間位かかりました。） → 計算結果
5000枚までなので極限値とは遥かに離れていますが、1/2へ向かっているとはいえフラフラ
している感じです。
DD++さんからのコメントです。（平成２９年１月１８日付け）
なるほど、pが小さいほどぶれ幅は大きくなるのですね。2/3, 3/4, 4/5 みたいに変えていく
とぶれ幅は大きくなるのか小さくなるのか、どうでしょう。
ふーむ、つまり log[2]n の小数部分がぶれに影響している感じですかね。なるほどなるほど。
らすかるさんからのコメントです。（平成２９年１月１８日付け）
このぶれについて、もう少し詳しく調べてみました。nは本来自然数ですが、以下実数として
扱います。
まず、n=2^100～2^101を100分割（2^100.00,2^100.01,2^100.02,…,2^100.99）して、それぞれ
の「理論値(n)-log[2]n」を求めてグラフを描くと、（見た目は）非常に綺麗なサインカーブになり
ます。そして、その極値の位置と値を調べてみると、log[2]nの小数部が 0.1111036249791…
のとき、最大値 1.33274775043405996376381438046…、小数部が 0.6111038362867…のと
き、最小値 1.33274460411865843517321601464…をとります。
log[2]nの小数部の差が0.500000211307…ですから、サインカーブそっくりでもサインカーブ
ではなかったことになります。
最大値と最小値の平均をとると、1.33274617727635919946851519755となり、これは、
γ/log2+1/2の1.33274617727686715064641751940…と非常に近いですが、これも微妙に異
なります。
# 上記の微差は計算誤差ではありません。n≒2^100ではn→∞のときの収束値と小数点以
下29桁程度が一致することが数値的にわかっていますので、上に記載した値は全桁有効数
字です。
サインカーブではありませんので近似式しか作れませんが、とりあえず今までのデータから
作ってみると、
(理論値(n))-log[2]n≒1.332746177277+0.0000015731581*sin(2πlog[2]n+0.872711)
この近似式による計算結果を以前出した真値と比較してみると、
n=2^kの場合
1.33274738243291305148… ← 上記近似式の計算結果
1.33274738243289922500… ← 真値
n=3・2^kの場合
1.33274462529166362305… ← 上記近似式の計算結果
1.33274462529161259054… ← 真値
n=5・2^kの場合
1.33274656060791017970… ← 上記近似式の計算結果
1.33274656060788254518… ← 真値
n=7・2^kの場合
1.33274565604122989917… ← 上記近似式の計算結果
1.33274565604129766087… ← 真値
これだけ良い近似になっていることから、サインカーブにかなり近いことがわかりますね。
また一つ面白いことを発見しました。前に「最大値と最小値の平均がγ/log2+1/2に近いが
一致しない」と書きましたが、それの関連です。試しに、「最大値と最小値の平均」でなく、
「log[2]nの小数部が0のときの値(すなわち「A158468」の値)とlog[2]nの小数部が0.5のときの
値の平均」にしてみると、1.33274738243289922500… と1.33274497212168702001… の平
均である1.33274617727729312251… となり、この値はやはりγ/log2+1/2に近いけれども
一致しない値になります（小数点以下11桁が一致）。
そこで今度は、小数部が、0,0.25,0.5,0.75の4つの平均にしてみると、
1.332747382432899225008601098373… と1.332747188426833825810130724941… と
1.332744972121687020018551808032… と1.332745166126048531748386045602… の平均
である1.332746177276867150646417419237… となり、相変わらずγ/log2+1/2には一致し
ませんが、今度は小数点以下24桁が一致します。
更に進めて、0,0.125,0.25,…,0.875の8つの平均にしてみると
1.332747382432899225008601098373899704416743982259844536579718…
1.332747744441277104972326119911514968238749189221115770554949…
1.332747188426833825810130724941413368406388004542018131127982…
1.332746040095056077554177929803387298360834559021500119723502…
1.332744972121687020018551808032416109064464926033884719957685…
1.332744610111275303020082449731014441839771078589089832304373…
1.332745166126048531748386045602779874882365834949179244539924…
1.332746314459860117039083978868498516090134673593955996501410… の平均
1.332746177276867150646417519408115535162431531026323543911193…
となり小数点以下50桁が一致、同様に16個の平均にすると小数点以下99桁が一致、のよう
に、すごい速さでγ/log2+1/2に近づきます。
# 150桁計算にして32個平均が150桁全部一致するところまで確認しました。n=2^500で150
桁まで収束します。誤差のオーダーが1/nということですね。
つまりこれは、「一周期分の平均値はγ/log2+1/2に収束する」ということを示しています。
# また、この計算によってγを任意の桁数計算できるということにもなります。
今までは「γ/log2+1/2と近い」というだけで関係ない値の可能性もあったわけですが、これ
でようやくγ/log2+1/2とつながりました。
# 元の問題は単純でしたが、つきつめていくと実に奥が深いですね。このテーマで論文が書
けそうです。
DD++さんからのコメントです。（平成２９年１月１８日付け）
おおー、素晴らしい！念のため確認ですが、計算は Σ(1-1/2^k)^n（で、n=2^x） で行ってる
ものでいいんですかね？幾つか計算式があって、nが整数の場合は同じ式ですが整数に限
らない場合は同じ保証はないので……。
大した考えもなく、なんとなく言ってみた「収束しそうに見えますが、どうなんですかね？」が
まさかこんま壮大な話になるとは全く予想していませんでした。
らすかるさんからのコメントです。（平成２９年１月１９日付け）
計算は Σ(1-1/2^k)^n（で、n=2^x） で行ってるものでいいんですかね？
はい、おっしゃる通りです。他の計算式では高速に計算できませんので、その式を使い出し
てからはその式しか使っていません。式で書くと、
lim[m→∞]∫[m～m+1]1-x+Σ[k=1～∞](1-(1-1/2^k)^(2^x)) dx=γ/log2+1/2
でしょうか。
DD++さんからのコメントです。（平成２９年１月１９日付け）
d/dx { ∫[定数～2^x*log(a)] e^t/t dt}
= e^(2^x*log(a))/(2^x*log(a)) * 2^x*log(a) * log2 = a^(2^x) * log2
を利用すると、
lim[m→∞] ∫[m～m+1] {1-x+Σ[k=1～∞](1-(1-1/2^k)^(2^x))} dx
= lim[m→∞] { (1/2)-m+Σ[k=1～∞](1- ∫[m～m+1] (1-1/2^k)^(2^x)dx) }
= lim[m→∞] { (1/2)-m+Σ[k=1～∞](1- (1/log2) ∫[2^m*log(1-1/2^k)～2^(m+1)*log(1-1/2^k)] e^t/t dt) }
となりますね。先頭に 1/2 が出て、しかも 1/log2 が出ました。積分のところは指数積分を彷
彿とさせる形なのでγの出現が期待できます。
無限和が1に収束するような数列 a[k] と b[k] を使って、積分のところが
log2 - a[k]*γ - b[k]*log(2^m) + 0収束項
と書けるということでしょうかね。
ＧＡＩ さんからのコメントです。（平成２９年１月１９日付け）
まったく関係ないかもしれませんが、[0,1]区間のあらゆる有理数を構成する
E0=[0,1]
E1=[0,1/2,1]
E2=[0,1/3,1/2.2/3,1]
E3=[0,1/4,1/3,2/5,1/2,3/5,2/3,3/4,1]
E4=[0,1/5,1/4,2/7,1/3,3/8,2/5,3/7,1/2,4/7,3/5,5/8,2/3,5/7,3/4,4/5,1]
以下En:[･･･,q1/p1,q2/p2,･･･]に対し、E(n+1):[･･･,q1/p1,(q1+q2)/(p1+p2),q2/p2,･･･]
を割り込ませていくと、あらゆる有理数x=q/p (0<x<1)が重複することなく現れる集合がとれ
ますが、目的の値xを連分数表示したもの（k/2^mなる２進有理数）と対応させて[0,1]にある
有理数全体と全単射対応を与える関数として、Minkowski's question mark function (?(x)で
表記される）ものを使って（有理数なら有限和、無理数なら無限和となる）
γ/log2+1/2
=1.33274738243289922500860109837389970441674398225984453657972
の小数部分 0.133274･･･ を計算させてみると、
?(0.133274･･･)
= (5386379163185534471414773640069251753523547958090930032120046037215359507
85759609387310293523731869555656043490188556443564947941991142751393898160
9541905675481992802284708758507609390975756931455180006556095741937)
/
(2154551665274213788565909456027700701409419183236372012850495857896952690
40008430473775781325930022962257222404349059797525275944967528606233715293
00590103961416596156942109074193054752294185849943217159109760516096)
=0.2500000000000000000000000000
=1/4
とほとんど一致します。
らすかるさんからのコメントです。（平成２９年１月１９日付け）
γ/log2+1/2の値は(理論値(2^k))-kの収束値であり、γ/log2+1/2の値ではないですが、ど
ちらの値のときに1/4になるのでしょうか。γ/log2+1/2の値は、
1.33274617727686715064641751940811553516243153102632810163149…
です。
ＧＡＩ さんからのコメントです。（平成２９年１月１９日付け）
「A158468」から張られていた「A158469」での値です。しかし、
? Minkof(0.33274617727686715064641751940811553516243153102632810163149)
%233 = 167597599124282463744675312477573076593492072757404917221330541062160
0568484339603772481442045583674151826360229331984956558278330017274417
636495547432959/670390396497129854978701249910292306373968291029619668
8861780721860882015036773488400937149083451713845015929093243025426876
941405973284973216824503042048
? %+0.
%234 = 0.2500000000000000000000000000
こちらも同様に1/4
ここに、Minkof(x)(=?(x))は次のもので構成しました。(PARIでのプログラムです。)
Minkof(x)={v=vecextract(contfrac(x),"2..")};//先頭のベクトル成分0を取り除くため
sum(i=1,#v,(-1)^(i-1)/(2^(sum(k=1,i,v[k])-1)));
らすかるさんからのコメントです。（平成２９年１月１９日付け）
Minkowski's question mark function のグラフを見てみると、1/3付近ではほとんど平らに
なっていますので、単に「1/3に近い値の場合は1/4に非常に近い値になる」というだけのこ
とではないでしょうか。
ＧＡＩ さんからのコメントです。（平成２９年１月１９日付け）
1/3前後を1/1000単位で動かすと0.25前後の値はとりますね。
991/3000;0.2499999999977973175191436894
497/1500;0.2499999999999999933863667478
997/3000;0.2500000000000000000000000000
1/3;0.2500000000000000000000000000
1003/3000;0.2500000000000000000000000000
503/1500;0.2500000000000000069659989582
1009/3000;0.2500000000030837554731988348
どちらも0.33274辺りですから、1/3-33274/100000=89/150000=0.000593333
その差は、3/5000。十分0.25範囲に収まる所にたまたま居るわけですね。
らすかるさんからのコメントです。（平成２９年２月２日付け）
ちょっと間が空いてしまいましたが、追跡調査結果です。
1 + Σ[k=1～∞] 1-(1-1/2^k)^n がlog[2]nの小数部で決まるサインカーブっぽい曲線に収
束するという話の続きですが、「サインカーブとの差をとると、どういう曲線になるんだろう」と
いう疑問がありましたので調査しました。
もちろん、サインカーブのとりかたで残る曲線は変わるのですが、Excelのグラフ機能で振
幅と位相を微妙に調整していくと、なんと「周期が半分のサインカーブ（っぽい非常に綺麗な
曲線）」（ただし振幅は最初のサインカーブと比べて非常に小さい）が残りました。
…ということは、フーリエ変換？ ということで、log[2]n=500.00、500.01、500.02、…、500.99
の100個のデータ（いずれも150桁以上収束）を使ってフーリエ級数展開してみると、
({(理論値(n))-log[2]n}の収束値)
= (γ/log2+1/2)
+ 0.00000157315770076355974575250413267281507413914124…
sin(2πlog[2]n+0.87271099890815930080058867129393842162914996787948…)
+ 0.00000000000072847098432903633382448549588022717973…
sin(4πlog[2]n+2.51702341788490728464101194634453992745579581698030…)
+ 0.00000000000000000038951270145869284493612416910423…
sin(6πlog[2]n+5.70449041268300778774992450513287834008165858559157…)
+ 0.00000000000000000000000022090586734847275937917260…
sin(8πlog[2]n+5.81254926179414763386268318796456736588363327177119…)
+ 0.00000000000000000000000000000012939194463915830520…
sin(10πlog[2]n+3.63051777405709641694886245041298806565675604994913…)
+ 0.00000000000000000000000000000000000007735201159493…
sin(12πlog[2]n+5.90659684749221331225221137882780917590439830849139…)
+ …
のようになることがわかりました。
次項の比が100万分の1以下と非常に小さくなりますので、「綺麗なサインカーブ」「周期半
分の綺麗なサインカーブ」に見えるわけです。グラフでは100万分の1の誤差は見えません
ので。
あとは、上記の係数の正体がわかれば終わりですが、これがまた難しそうです。とりあえ
ず位相の係数はおいといて、振幅の係数は、ほぼ6桁ずつ小さくなっていますので、まずは
比をとってみると、
（上記は50桁だけ書きましたので6項分の係数しかありませんが、150桁計算で出した23項
の係数の比です）
n=1～22に対する(第n項の振幅の係数)÷(第n+1項の振幅の係数)は、
2159533.78323301…
1870211.11661050…
1763251.95040767…
1707261.35977416…
1672767.67560675…
1649371.17124660…
1632454.09665057…
1619650.33742475…
1609621.44633876…
1601553.11760446…
1594921.39919034…
1589373.82968918…
1584664.56354631…
1580616.88244393…
1577100.48902825…
1574017.20066562…
1571291.62771048…
1568864.92104579…
1566690.47104577…
1564730.88244766…
1562955.80430114…
1561340.35936744…
のようになります。何だか150万近辺の値に収束しそうですが、この収束値を数値的に求め
るのは厳しそうです。
グラフにプロットしてみると、直角双曲線の形に近く、近似式を適当に作ると、
762433/(n+0.23344)+1527050
あたりになります。
よって、1527050（より少し小さい値）あたりに収束しそうです。しかし、この値も思い当たる
ものがなく、不明です。今のところここで手詰まりですが、また何か思いついたら書きます。
らすかるさんからのコメントです。（平成２９年２月３日付け）
もう少し調査が進みました。
以前書いた振幅の隣項比 2159533.78323301… ～ 1561340.35936744… の隣項の比を
とり、1引いて逆数にして近い一般項の形との差をとって、さらに逆数にして、・・・などゴリゴ
リ計算していって、逆算することにより、この値の一般項が
e^(π^2/log2)√{(n+1)/n}
と求まり、ここからさらに逆算することで、振幅の値
0.00000157315770076355974575250413267281507413914124…
0.00000000000072847098432903633382448549588022717973…
0.00000000000000000038951270145869284493612416910423…
0.00000000000000000000000022090586734847275937917260…
0.00000000000000000000000000000012939194463915830520…
0.00000000000000000000000000000000000007735201159493…
・・・
の一般項が、 2/{e^(nπ^2/log2)√(nlog2)} と求まりました。
フーリエ級数展開によって計算した上記の値は若干誤差がありますが、（初項小数第32
位以下、第2項小数第63位以下、第3項小数第93位以下、…）、これはフーリエ変換時の積
分を100区間の区分求積で代用していることによる誤差と思われます。
一般項の式が（巨大な定数などのない）かなり綺麗な形になっていますので、おそらくこれ
で間違いないでしょう。この式を使って検算したいところですが、位相の値の一般式が出て
いませんので、今のところ検算ができません。
ここまでで、
({(理論値(n))-log[2]n}の収束値)
= (γ/log2+1/2)
+ 2sin(2πlog[2]n+0.8727109989081593…)/{e^(π^2/log2)√(log2)}
+ 2sin(4πlog[2]n+2.5170234178849072…)/{e^(2π^2/log2)√(2log2)}
+ 2sin(6πlog[2]n+5.7044904126830077…)/{e^(3π^2/log2)√(3log2)}
+ 2sin(8πlog[2]n+5.8125492617941476…)/{e^(4π^2/log2)√(4log2)}
+ 2sin(10πlog[2]n+3.6305177740570964…)/{e^(5π^2/log2)√(5log2)}
+ 2sin(12πlog[2]n+5.9065968474922133…)/{e^(6π^2/log2)√(6log2)}
+ …
となりましたが、残る位相の値の一般項はどうやって調べたら良いのでしょうね…。
# ところで、WolframAlphaは便利ですね。最初に求まったのが振幅の隣項比の収束値
e^(π^2/log2)=1527020.98232550…
なのですが、これは数値的に求めて出た1527020.98232550…という値をWolframAlphaに入
力すると「Levy constantの12乗」と出て、「Levy constant」を調べて、e^(π^2/log2) とわかり
ました。
また、収束値と隣項比から数値計算される第0項の 2.40224481…の正体も、WolframAlpha
により、2/√(log2) と判明して一般項にたどり着いたという次第です。
DD++さんからのコメントです。（平成２９年２月４日付け）
1/√π * 2^(1-n*x^2) * cos(2nπ*x) を -∞ から ∞ まで積分するとこれになるようです。
n=1,2,3,4,5 でwolfram先生に尋ねるとその通りに返ってくるので。
いかにもフーリエ変換という感じの積分なので、何か関係ありそう。とはいえ 2^(1-n*x^2)
の中にも n がいるので、単純に考えていいわけでもなさそうですが。
ペントミノさんからのコメントです。（平成２９年２月９日付け）
久々の書き込みです。オイラー定数が登場して、とても奥が深そうな問題ですね。何とか
一般式が出せないものかと、考察してみました。
より一般性を持たせるために、表が出る確率を p (0<p<1) とする。コインが n 枚あるとき
の期待値（理論値）を E(n) で表すことにする。
E(n) = 1 + Σ[m=1～∞] {1 - (1-p^m)^n}
Σ[m=1～∞] の部分はあえて残したまま、下記の A(n) を定義する。
A(n) = 1 + Σ[m=1～∞] {1 - e^(-p^m n)}
まず、n が十分大きいとき、E(n) は A(n) で近似できることを示す。
0<a<1 とするとき、1-a < e^(-a) < (1-a)^(1-a) が成り立つので、
0 < e^(-p^m n) - (1-p^m)^n
< e^(-p^m n) - e^{-p^m n / (1-p^m)}
= e^(-p^m n) [1 - e^{-p^(2m)n / (1-p^m)}]
< e^(-p^m n) p^(2m)n / (1-p^m)
t = log[1/p]n とおき、まず 1<=m<=t の場合を扱う。
0≦x では x^2 e^(-x) <= 4e^(-2)、および 1/(1-p^m) <= 1/(1-p) より
Σ[m=1～floor(t)] e^(-p^m n) p^(2m)n / (1-p^m)
≦Σ[m=1～floor(t)] 4e^(-2) / {n(1-p)}
≦4e^(-2) log[1/p]n / {n(1-p)}
次に、t<m の場合、p^m n < 1 および 1/(1-p^m) <= 1/(1-1/n) = n/(n-1) より
Σ[m=floor(t)+1～∞] e^(-p^m n) p^(2m)n / (1-p^m)
≦Σ[m=floor(t)+1～∞] p^m n / (n-1)
=p^{floor(t)+1} n / {(1-p)(n-1)}
< 1 / {(1-p)(n-1)}
以上をまとめると、
lim[n→∞] {E(n) - A(n)}= lim[n→∞] Σ[m=1～∞] {e^(-p^m n) - (1-p^m)^n}= 0
よって、n が十分大きいとき、E(n) は A(n) で近似できる。
次に、t が十分に大きいとき、D(t) = E(p^(-t)) - t は周期1の周期関数になることを示す。
A(p^(-1)n) = 1 + Σ[m=1～∞] [1 - e^{-p^(m-1) n}]
= 1 + 1 - e^(-n) + Σ[m=1～∞] [1 - e^{-p^m n}]
=1 - e^(-n) + A(n) において、
lim[n→∞] {A(p^(-1)n) - A(n)} = 1
lim[n→∞] {E(p^(-1)n) - E(n)} = 1
lim[t→∞] {E(p^(-t-1)) - E(p^(-t))} = 1
lim[t→∞] [{E(p^(-t-1)) -t-1} - {E(p^(-t)) - t}] = 0
lim[t→∞] {D(t+1) - D(t)} = 0
よって、t が十分に大きいとき、D(t) は周期1の周期関数である。
次に、t が十分に大きいとき、D(t) の1周期分の平均 d[0] は 1/2 + γlog[1/p]e であること
を示す。
d[0] = ∫[t～t+1] D(x) dx
= ∫[t～t+1] {A(p^(-x)) - x} dx
= ∫[t～t+1] [1 - x + Σ[m=1～∞] {1 - e^(-p^(m-x))}] dx
= (t+1) - t - (t+1)^2/2 + t^2/2 + Σ[m=1～∞] ∫[t～t+1] {1 - e^(-p^(m-x))} dx
= 1/2 - t + Σ[m=1～∞] ∫[t-m～t+1-m] {1 - e^(-p^(-x))} dx
= 1/2 - t + ∫[-∞～t] {1 - e^(-p^(-x))} dx
ここで、p^(-x) = y とおくと、x = log[1/p]y, dx = log[1/p]e (dy / y) より
d[0] = 1/2 - log[1/p]n + log[1/p]e ∫[0～n] {1 - e^(-y)} dy/y
= 1/2 - log[1/p]n + log[1/p]e [∫[0～1] {1 - e^(-y)} dy/y + ∫[1～n] {1 - e^(-y)} dy/y]
= 1/2 - log[1/p]n + log[1/p]e {- ∫[0～1] e^(-y) log[e]y dy
+ ∫[1～n] dy/y - ∫[1～n] e^(-y) log[e]y dy}
= 1/2 - log[1/p]n + log[1/p]e {log[e]n - ∫[0～n] e^(-y) log[e]y dy}
= 1/2 - log[1/p]e ∫[0～n] e^(-y) log[e]y dy
ここで、∫[0～∞] e^(-y) log[e]y dy = -γ より、lim[n→∞] d[0] = 1/2 + γlog[1/p]e
最後に、t が十分に大きいとき、D(t) = Σ[k=-∞～∞] d[r] e^(i2rπt) とフーリエ級数を定
義すると、
r≠0 のとき、d[r] = - (i/2rπ) Γ(1 - i2rπ log[1/p]e)
特に、 |d[r]| = √(log[1/p]e / {2r sinh(2rπ^2 log[1/p]e)})
となることを示す。
d[r] = ∫[t～t+1] D(x) e^(-i2rπx) dx
= ∫[t～t+1] {A(p^(-x)) - x} e^(-i2rπx) dx
= ∫[t～t+1] [1 - x + Σ[m=1～∞] {1 - e^(-p^(m-x))}] e^(-i2rπx) dx
= (i/2rπ) {(1 - (t+1)) e^(-i2rπ(t+1)) - (1 - t) e^(-i2rπt)}
+ Σ[m=1～∞] ∫[t～t+1] {1 - e^(-p^(m-x))} e^(-i2rπx) dx
= - (i/2rπ) e^(-i2rπt) + ∫[-∞～t] {1 - e^(-p^(-x))} e^(-i2rπx) dx
= - (i/2rπ) e^(-i2rπt) + (i/2rπ) {1 - e^(-p^(-t))} e^(-i2rπt)
- (i/2rπ) log[e](1/p) ∫[-∞～t] e^(-p^(-x)) p^(-x) e^(-i2rπx) dx
= - (i/2rπ) e^(-p^(-t)) e^(-i2rπt)
- (i/2rπ) log[e](1/p) ∫[-∞～t] e^(-p^(-x)) p^(-x) p^(i2rπx log[1/p]e) dx
ここで、p^(-x) = y とおくと、x = log[1/p]y, dx = log[1/p]e (dy / y) より
d[r] = - (i/2rπ) e^(-n) e^(-i2rπt) - (i/2rπ) ∫[0～n] e^(-y) y^(-i2rπ log[1/p]e) dy
ここで、∫[0～∞] e^(-x) x^z dx = Γ(1+z) より、
lim[n→∞] d[r] = - (i/2rπ) Γ(1 - i2rπ log[1/p]e)
さらに、Γ(1+z) Γ(1-z) = z Γ(z) Γ(1-z) = zπ / sin(zπ) より、
lim[n→∞] |d[r]|^2 = d[r] d[-r] = 1/(2rπ)^2 i2rπ^2 log[1/p]e / sin(i2rπ^2 log[1/p]e)
= log[1/p]e / {2r sinh(2rπ^2 log[1/p]e)}
よって、 lim[n→∞] |d[r]| = √(log[1/p]e / {2r sinh(2rπ^2 log[1/p]e)})
DD++さんからのコメントです。（平成２９年２月１０日付け）
= 1/2 - log[1/p]n + log[1/p]e [∫[0～1] {1 - e^(-y)} dy/y + ∫[1～n] {1 - e^(-y)} dy/y]
= 1/2 - log[1/p]n + log[1/p]e {- ∫[0～1] e^(-y) log[e]y dy
+ ∫[1～n] dy/y - ∫[1～n] e^(-y) log[e]y dy}
の部分は、
= 1/2 - log[1/p]n + log[1/p]e {- ∫[0～1] e^(-y) log[e]y dy
+ ∫[1～n] dy/y - e^(-n) log[e]n - ∫[1～n] e^(-y) log[e]y dy}
でしょうか？どうせ n->∞ の極限で消える項ではありますが。
（というか、部分積分をする上で一度積分区間を1で分ける意味もないような気もします）
lim[n→∞] d[r] = - (i/2rπ) Γ(1 - i2rπ log[1/p]e)
において、u+vi とおくと、
(u+vi) e^(i*2rπt) + (u-vi) e^(-i*2rπt)
= 2u cos(2rπt) - 2v sin(2rπt)
= 2 √(u^2+v^2) sin(2rπt+α)
= 2 √(log[1/p]e / {2r sinh(2rπ^2 log[1/p]e)}) sin(2rπt+α)
で、 2sinh(2rπ^2 log[1/p]e) をだいたい e^(2rπ^2 log[1/p]e) と思った結果がらすかるさん
の推定した式、そのときの謎の誤差が -e^(-2rπ^2 log[1/p]e) の分だったわけですか。そし
て、位相の部分は α = arctan(-2u/2v) = arg[Γ(1 - i2rπ log[1/p]e)] だったわけですね。
つまり、p=1/2 とすれば、
({(理論値(n))-log[2]n}の収束値)
= (γ/log2+1/2)
+ 2sin(2πlog[2]n+arg[Γ(1-2iπ/log2)])/{√(2log2*sinh(2π^2/log2))}
+ 2sin(4πlog[2]n+arg[Γ(1-4iπ/log2)])/{√(4log2*sinh(4π^2/log2))}
+ 2sin(6πlog[2]n+arg[Γ(1-6iπ/log2)])/{√(6log2*sinh(6π^2/log2))}
+ ……
いやあ、すっきり解決、素晴らしいです。
らすかるさんからのコメントです。（平成２９年２月１０日付け）
計算誤差と安易に片づけてしまった妙な誤差は、実は本質的な誤差だったのですね。安易
に考えずにちゃんと検討すれば良かったです。DD++さんと全く同感で、すっきり解決、素晴ら
しいです。
({(理論値(n))-log[2]n}の収束値)の式を変形すると
({(理論値(n))-log[2]n}の収束値)
= (γ/log2+1/2)
+ {Im{Γ(1-2πi/log2)}cos(2πlog[2]n)+Re{Γ(1-2πi/log2)}sin(2πlog[2]n)}/π
+ {Im{Γ(1-4πi/log2)}cos(4πlog[2]n)+Re{Γ(1-4πi/log2)}sin(4πlog[2]n)}/(2π)
+ {Im{Γ(1-6πi/log2)}cos(6πlog[2]n)+Re{Γ(1-6πi/log2)}sin(6πlog[2]n)}/(3π)
+ …
となります。私が最初フーリエ変換した時は、
a[0]+Σ[k=1～∞]{a[k]cos(2πlog[2]n)+b[k]sin(2πlog[2]n)}
という形のa[k]、b[k]を求め、そこから合成してsinだけの式にしたのですが、最初に求めた
a[k]、b[k]の値
a[1]=0.00000120515560610270806152746730845087107713049223…
b[1]=0.00000101115039264735696347935132034009875765141159…
a[2]=0.00000000000042597186715903396553064257655998031101…
b[2]=-0.00000000000059094665021333475733910102118069973679…
a[3]=-0.00000000000000000021303688229650588074674247248521…
b[3]=0.00000000000000000032609113968188936118331382649300…
…
は、実は、
a[1]=Im{Γ(1-2πi/log2)}/π
b[1]=Re{Γ(1-2πi/log2)}/π
a[2]=Im{Γ(1-4πi/log2)}/(2π)
b[2]=Re{Γ(1-4πi/log2)}/(2π)
a[3]=Im{Γ(1-6πi/log2)}/(3π)
b[3]=Re{Γ(1-6πi/log2)}/(3π)
…
だったのですね。値はピタリと一致していました。そして、この式から、n=2^mの場合の収束
値 1.332747382432899225008601098… が、
(A158468)=γ/log2+1/2+Σ[k=1～∞]Im{Γ(1-2kπi/log2)}/(kπ)
と表せることがわかります。
（元々ペントミノさんが計算された式なので私の名前で投稿するのはちょっと気が引けました
が）「A158468」にこの関係式を投稿しました。
# Im{Γ(1-2kπi/log2)} は簡単にならないですよね？
