higher genus AGMは超幾何級数表示を持つか -

a,bの算術幾何平均を、M(a,b)として、FがGauss超幾何関数の時

$\dfrac{1}{M(1,x)} = F(\dfrac{1}{2} , \dfrac{1}{2},1 ; 1-x^2)$

であることは、よく知られてる。

19世紀後半に、Borchardtという人が、種数2のテータ関数(テータ定数)の倍角公式に基づいて、算術幾何平均の一般化を作ったようだ。

[PDF] 算術平均の周辺：Borchardtによる４項算術幾何平均とThomaeの公式
https://www2.tsuda.ac.jp/suukeiken/math/suugakushi/sympo18/18_2shiga.pdf

Borchardtの論文は、オンラインでは見つけられなかった

Über das arithmetisch-geometrische Mittel aus vier Elementen
https://books.google.co.jp/books/about/%C3%9Cber_das_arithmetisch_geometrische_Mitt.html?id=iA0inQAACAAJ&redir_esc=y

通常の算術幾何平均では、楕円曲線 $y^2=4(x-a_1)(x-a_2)(x-a_3)$ ( $a_1 \lt a_2 \lt a_3$ )に対して、2つの周期
$\omega_1 = \dfrac{\pi}{M(\sqrt{a_3-a_1} , \sqrt{a_3-a_2})}$
$\omega_2 = \dfrac{i \pi}{M(\sqrt{a_3-a_1} , \sqrt{a_2-a_1})}$
を計算できるのと同様、Borchardt AGMの場合は、2つの2x2周期行列の行列式が計算できる。BorchardtのAGMは、テータ関数側から作ってるので、超幾何関数に相当するようなものは、調べた限り、明示的に書かれているものはない。時期的には、Borchardtの論文は1876年で、Appelの超幾何級数の論文が1880年、Lauricellaの超幾何級数の論文が1893年(Wikipedia情報)なので、多変数超幾何関数の方が新しい

最近は(?)、他にも、色んな"平均"が構成されていて、これらは、超幾何関数の変換公式を利用している

[PDF] 超幾何関数の関数等式と算術調和平均およびその拡張
https://www.riam.kyushu-u.ac.jp/fluid/meeting/18ME-S5/papers/Article_No_38.pdf

超幾何関数の変換公式と平均反復
http://fe.math.kobe-u.ac.jp/Movies/cm/2009-01-08-matsumoto.pdf

このへんを、何気なく、見てると、とある四項平均が、三変数のLauricella超幾何関数の二乗で書けるという式が書いてあったので、Borchardt AGMも、意外と簡単な形で書けたりするのだろうかと疑問に思った。

頑張れば、超楕円積分をLauricella超幾何関数で書いて解析的に計算できるかもしれないけど、別に重要な問題とも思えないので、そこまでモチベーションがない。Mathematicaなら、超楕円積分を適当な超幾何関数で書くくらい、やってくれないかと期待したけど、無理だったので、数値計算でやってみることにした。

何もなしには出来ないので、いくつか仮定を置く。W(a,b,c,d)をBorchardt AGMとして、
$\dfrac{1}{W(1,x_1,x_2,x_3)} = F(1-x_1^2 , 1-x_2^2 , 1-x_3^2)$
となる冪級数Fがあり、係数は、有理数で書けるものとする。(x1,x2,x3)=(1,1,1)の時、左辺は1なので、F(0,0,0)=1である。収束半径とかは分からないけど、どうせ原点近くでの振る舞いしか見ないので、気にしなくていいだろう

Borchardt AGMの計算は、平方根を使うので、有理数だけで済ますことはできず、倍精度程度の浮動小数点数では全然精度が足りない。ということで、任意精度演算をやるために、pythonのmpmathを使う。やることは簡単なので、別になんでもいいけど

最初は、3変数の内、2つは0とする。
$F(z,0,0) = 1 + c_{100}z + c_{200} z^2 + \cdots$
で、zが十分小さければ、高次の項の寄与は、殆ど無視できるはず

例えば、一次の項を決めるのに、以下のようなコードを実行する

from mpmath import *
mp.dps = 500


def M(x,y,z):
   a,b,c,d = 1,x,y,z
   max_iter = 1000
   for n in range(max_iter):
      a2 = (a+b+c+d)/4
      b2 = (sqrt(a*b) + sqrt(c*d))/2
      c2 = (sqrt(a*c) + sqrt(b*d))/2
      d2 = (sqrt(a*d) + sqrt(b*c))/2
      a,b,c,d=a2,b2,c2,d2
   return a


if __name__=="__main__":
  for eps in [1.0e-5 , 1.0e-10, 1.0e-15, 1.0e-20]:
    x1,x2,x3=1-mpf(eps),1,1
    z1,z2,z3=1-x1*x1,1-x2*x2,1-x3*x3
    c = (1/M(x1,x2,x3)-1)/z1
    print(str(c)[:50])

結果は、

0.125001171883984440643884854175753723800930555241
0.125000000011718750000898437927005455667004923082
0.125000000000000117187500000000098949851548046299
0.125000000000000000001171874999999999935735474360

とかなるので、まぁ、0.125=1/8に収束しそう、と分かる。どこに収束するかは勘で判断してるだけので、合ってるかどうか証明とかはない

同様に、 $z_1^2$ の係数は、以下のようなコードで決める

if __name__=="__main__":
  for eps in [1.0e-5 , 1.0e-10, 1.0e-15, 1.0e-20,1.0e-30]:
    x1,x2,x3=1-mpf(eps),1,1
    z1,z2,z3=1-x1*x1,1-x2*x2,1-x3*x3
    c=(1/M(x1,x2,x3)- (1+z1/8))/(z1*z1)
    print(str(c)[:50])

結果は

0.058594492194493161915357892508282237042562173272
0.058593750007421875000699310573188476983608793132
0.058593750000000074218750000000075698227920533509
0.058593750000000000000742187499999999959300436760
0.058593750000000000000000000000074218750000000006

で、0.0589375 = 15/256に収束しそう

続ける。

if __name__=="__main__":
  for eps in [1.0e-5 , 1.0e-10, 1.0e-15, 1.0e-20,1.0e-30]:
    x1,x2,x3=1-mpf(eps),1,1
    z1,z2,z3=1-x1*x1,1-x2*x2,1-x3*x3
    c=(1/M(x1,x2,x3)- (1+z1/8+15*z1*z1/256))/(z1*z1*z1)
    print(str(c)[:50])

で、

0.037109910207646803092967250526967109948749607120
0.037109375005352020264233913616668406369705086927
0.037109375000000053520202636718810362650910452159
0.037109375000000000000535202026367187470651540126
0.037109375000000000000000000000053520202636718754

とかになるので、0.037109375=19/512だろう

だんだん辛くなってきたけど

if __name__=="__main__":
  for eps in [1.0e-5 , 1.0e-10, 1.0e-15, 1.0e-20,1.0e-30]:
    x1,x2,x3=1-mpf(eps),1,1
    z1,z2,z3=1-x1*x1,1-x2*x2,1-x3*x3
    c=(1/M(x1,x2,x3)- (1+z1/8+15*z1*z1/256+19*z1*z1*z1/512))/(z1*z1*z1*z1)
    print(str(c)[:50])

が

0.026760516142735359850255545761188337409031621442
0.026760101322507572174538118990437334607904260568
0.026760101318359416481971740722706058694628176862
0.026760101318359375000414819717407226539753154096
0.026760101318359375000000000000041481971740722659

は、0.026760101318359375=7015/262144らしい。などとやっていくと、2つの変数が0の場合は決まる。

0でない変数が複数ある場合は、面倒くさくなる。とりあえず、二次、三次の係数は、普通にやればいい。 $z_1^3z_2$ なんかは、例えば

def F3(z1,z2,z3):
   v1 = (z1+z2+z3)/8
   v2 = 15*(z1*z1+z2*z2+z3*z3)/256 + 3*(z1*z2+z2*z3+z3*z1)/128
   v3 = 19*(z1*z1*z1+z2*z2*z2+z3*z3*z3)/512 + 3*(z1*z1*z2+z2*z2*z1+z1*z1*z3+z2*z2*z3+z3*z3*z1+z3*z3*z2)/256 + 3*z1*z2*z3/512
   return 1+v1+v2+v3

if __name__=="__main__":
  mp.dps = 1000
  for eps in [1.0e-20,1.0e-40,1.0e-80,1.0e-120]:
    for ratio in [0.9,0.0011,1.1e-30,1.1e-60]:
       x1,x2,x3=1-mpf(eps),1-mpf(eps)*mpf(ratio),1
       z1,z2,z3=1-x1*x1,1-x2*x2,1-x3*x3
       diff = (1/M(x1,x2,x3)- F3(z1,z2,z3) - 7015*(z1**4+z2**4+z3**4)/262144)
       if diff < z1*z1*z1*z2 and diff > z1*z1*z2*z2:
          c = diff/(z1*z1*z1*z2)
          print(str(c)[:50])

とやると

0.007620766506805420299495398378483141958519102149
0.007620766506805419922275006967238197476573541568
0.007614135745958588340065696042324398959715200301
0.007620766506805419922275006967238197472801337654
0.007614135742187500000000000000006621551513671874
0.007614135742187500000377108834006569608349931553
0.007620766506805419922275006967238197472801337654
0.007614135742187500000000000000006621551513671874
0.007614135742187500000000000000000000000000000000

なので、0.0076141357421875=499/65536と求まる

そんな感じで計算してみた表が、以下。合ってるのかどうか分からないけど、分母が、2のべきになったので、多分合ってる。以下が計算値のまとめ

i	j	k	c(i,j,k)
1	0	0	1/8
2	0	0	15/256
1	1	0	3/128
3	0	0	19/512
2	1	0	3/256
1	1	1	3/512
4	0	0	7015/262144
3	1	0	499/65536
2	2	0	789/131072
2	1	1	201/65536
5	0	0	43487/2097152
4	1	0	11685/(2^21)
3	2	0	4161/(2^20)
3	1	1	1059/(2^19)
2	2	1	15123/(2^21)
6	0	0	282387/(2^24)
5	1	0	36591/(2^23)
4	2	0	49185/(2^24)
7	0	0	946293/(2^26)
8	0	0	3323953575/(2^38)

意外と単純な係数が出たけど、超幾何関数そのものではないし、超幾何関数の二乗とかいう可能性も、なさそう。3変数のままやると辛いので、2つの変数は0にすると
$F(z,0,0) = 1 + z/8 +\dfrac{15}{256} z^2 + \dfrac{19}{512} z^3 + \dfrac{7015}{262144} z^4 + \cdots$
なので、形式的冪級数として、平方根を計算ことができるが、
$\sqrt{F(z,0,0)} = 1 + \dfrac{1}{16} z + \dfrac{7}{16^2} z^2 + \dfrac{69}{16^3} z^3 + \dfrac{6267}{2^19} z^4 + \cdots$
で、だめっぽい。

なんか、それらしいLauricella超幾何級数のいくつかの積になったりしないだろうかと思って、脳死状態で、コンピュータに計算してもらったけど、候補を見つけることはできなかった。まぁ、あんまり単純な式はないのかもしれない

とりあえず、数値計算で、冪級数の係数を予想するというのは初めてやったけど、意外と、うまくいくらしい

種数3以上の超楕円積分に対して、Borchardt AGMの一般化が、以下で与えられている。
Hyperelliptic integrals and generalizedarithmetic–geometric mean
https://link.springer.com/article/10.1007/s11139-011-9353-7