SX-4を使った円分多項式の計算について
鹿児島大学 理学部 小柴 洋一
１．この計算の目的
この計算は数学・整数論における円分多項式を求めるものです。より詳しいことは文献[１]、文献[2]をご覧下さい。その計算のアルゴリズムの必然性から高速、広領域を必要とする整数の数値計算です。スーパーコンピュータの高速化の発展に応じてその時々においてそれまで不可能であったことが可能に転化しつつ計算可能な範囲が広がりつつあります。
計算公式は３つあるがそのひとつとしてA．GrytczukとB．Tropakは次の計算公式を述べている（文献[3]）。
すなわちｎ－ｔｈ円分多項式を
とするとき
ここでμはMobius関数,ΦはEuler関数、（ｋ,ｄ）はｋとｄの最大公約数をあらわすとする。
その証明は、代数方程式の根の累乗和を同じ根の基本対称式で表す漸化式（文献[4]）、およびRamanujanの和についてのHolderの式から導き出せます（文献[5]）。
この部分については以下のmailing listでの議論が参考になりました。rsnt@math.metro-u.ac.jpおよびtnt@math.metoro-u.ac.jp（98年10月ごろ）
これをプログラム化すると次のようになります。
******************** ファイルの始め ********************
！Ａ．Ｇｒｙｃｚｕｋ ａｎｄ Ｂ．Ｔｒｏｐａｋの計算公式によるプログラム
ｐａｒａｍｅｔｅｒ （ｉｓｉｚｅ=３６４９５３６０）
ｉｎｔｅｇｅｒ*８ Ｔ（ｉｓｉｚｅ）
ｉｎｔｅｇｅｒ*８ ａ（０：ｉｓｉｚｅ）
ｉｎｔｅｇｅｒ*８ ｐｒｉｍｅ（２０）
ｉｎｔｅｇｅｒ*８ ｅｕｌｅｒ，ｒ，ｓ，ｎ，ｓｕｎ，ｍｙｕ，ｃｏｎｓｔ
！
ｎ＝３*５*７*１１*１７*１９
！ ｎ＝３*５*７
ｃａｌｌ ｃｈｅｃｋ（ｎ，ｐｒｉｍｅ，ｓ）
ｍｙｕ＝（‐１）**ｓ
！
ｅｕｌｅｒ＝１
ｄｏ １００ ｊ＝１，ｓ
ｅｕｌｅｒ＝ｅｕｌｅｒ*（ｐｒｉｍｅ（ｊ）‐１）
１００ ｃｏｎｔｉｎｕｅ
！
ｄｏ ５００ ｒ＝１，ｅｕｌｅｒ
Ｔ（ｒ）＝１
５００ ｃｏｎｔｉｎｕｅ
！
ｄｏ ３００ ｊ＝１，ｓ
ｄｏ ３００ ｒ＝１，ｅｕｌｅｒ/ｐｒｉｍｅ（ｊ）
Ｔ（ｐｒｉｍｅ（ｊ）*ｒ）＝Ｔ（ｐｒｉｍｅ（ｊ）*ｒ）*（１－ｐｒｉｍｅ（ｊ））
３００ ｃｏｎｔｉｎｕｅ
！
ａ（０）＝１
ｃｏｎｓｔ＝ｅｕｌｅｒ/１０
ｄｏ ４００ ｊ＝１，ｅｕｌｅｒ
！ ｄｏ ４００ ｊ＝１，ｃｏｎｓｔ
！ ｄｏ ４００ ｊ＝１，１０００００００
ｓｕｍ＝０
ｄｏ ６００ ｍ＝０，ｊ－１
ｓｕｍ＝ｓｕｍ＋ａ（ｍ）*Ｔ（ｊ－ｍ）
６００ ｃｏｎｔｉｎｕｅ
ａ（ｊ）＝－ｍｙｕ*（ｓｕｍ/ｊ）
４００ ｃｏｎｔｉｎｕｅ
ｗｒｉｔｅ（６，*）’ｍａｘ＝’，ＭＡＸＶＡＬ（ａ），’ｍｉｎ＝’，ＭＩＮＶＡＬ（ａ）
ｅｎｄ
******************** ファイルの終り ********************
サブルーチンcheckは本質的でないのでここで掲載することは省略。
時間がかかるのは、端末ラベル４００のdoループです。この計算テーマの問題部分です。
２．SX-4による計算結果
早速ｎ＝３×５×７×１３×１７×１９として計算実行してみた。以下がそのときのジョブ記録である。******************** ファイルの始め ********************
ｃｙｃ_Ｇｒｙｔｃｚｕｋ_Ｔｒｏｐａｋ_６４．９０：
ｆ９０： ｖｅｃ（１）： ｃｙｃ_Ｇｒｙｔｃｚｕｋ_Ｔｒｏｐａｋ_６４．９０， ｌｉｎｅ １３：
Ｖｅｃｔｏｒｉｚｅｄ ｌｏｏｐ
ｆ９０： ｖｅｃ（１）： ｃｙｃ_Ｇｒｙｔｃｚｕｋ_Ｔｒｏｐａｋ_６４．９０， ｌｉｎｅ １７：
Ｖｅｃｔｏｒｉｚｅｄ ｌｏｏｐ
ｆ９０： ｖｅｃ（１）： ｃｙｃ_Ｇｒｙｔｃｚｕｋ_Ｔｒｏｐａｋ_６４．９０， ｌｉｎｅ ２１：
Ｖｅｃｔｏｒｉｚｅｄ ｌｏｏｐ
ｆ９０： ｖｅｃ（１）： ｃｙｃ_Ｇｒｙｔｃｚｕｋ_Ｔｒｏｐａｋ_６４．９０， ｌｉｎｅ ２２：
Ｖｅｃｔｏｒｉｚｅｄ ｌｏｏｐ
ｆ９０： ｖｅｃ（１）： ｃｙｃ_Ｇｒｙｔｃｚｕｋ_Ｔｒｏｐａｋ_６４．９０， ｌｉｎｅ ２９：
Ｖｅｃｔｏｒｉｚｅｄ ｌｏｏｐ
ｆ９０： ｖｅｃ（１）： ｃｙｃ_Ｇｒｙｔｃｚｕｋ_Ｔｒｏｐａｋ_６４．９０， ｌｉｎｅ ３９：
Ｖｅｃｔｏｒｉｚｅｄ ｌｏｏｐ
ｆ９０： ｃｙｃ_Ｇｒｙｔｃｚｕｋ_Ｔｒｏｐａｋ_６４．９０， _ＭＡＩＮ： Ｔｈｅｒｅ ａｒｅ ６ ｄｉａｇｎｏｓｅｓ．
ｃｈｅｃｋ_６４．ｆ９０：
ｆ９０： ｖｅｃ（１）： ｃｙｃ_Ｇｒｙｔｃｚｕｋ_Ｔｒｏｐａｋ_６４．９０， ｌｉｎｅ １４： Ｖｅｃｔｏｒｉｚｅｄ ｌｏｏｐ
ｆ９０： ｃｈｅｃｋ_６４．ｆ９０，ｃｈｅｃｋ： Ｔｈｅｒｅ ｉｓ １ ｄｉａｇｎｏｓｉｓ．
ｍｏｅｂｉｕｓ_６４．f９０：
ｆ９０： ｖｅｃ（３）： ｍｏｅｂｉｕｓ_６４．f９０， ｌｉｎｅ ８ｎｖｅｃｔｏｒｉｚｅｄ ｌｏｏｐ．
ｆ９０： ｍｏｅｂｉｕｓ_６４．９０，ｍｏｅｂｉｕｓ： Ｔｈｅｒｅ ｉｓ １ ｄｉａｇｎｏｓｉｓ．
ｓｔａｒｔ ｔｉｍｅ＝Ｔｕｅ Ｎｏｖ ２ １３：４９：４９ ＪＳＴ １９９９
ｍａｘ＝ ６６９６０６ ｍｉｎ＝ －６５４５８９
．．．．．プログラム 情報．．．．．
経過時間 （秒） ： １４１８．３１９１８１
ユーザ時間 （秒） ： １４０８．７１１８２＝約２３分２９秒
システム時間 （秒） ： ２．０７４６７０
ベクトル命令実行時間 （秒） ： １４０８．１１４８２７
全命令実行時間 ： ８０７１５７０７８１３
ベクトル命令実行数 ： ２１５０９１９２９６６
ベクトル命令実行要素数 ： ５５０５５０７３５８１７８
浮動小数点データ実行要素数 ： １６５９９４１
MOPS 値 ： ３９５０．２１４３３７
MFLOPS 値 ： ０．００１１７８
平均ベクトル長 ： ２５５．９６０６６６
ベクトル演算率 （％） ： ９８．９３６０３７
メモリ使用量 （MB) ： ８０．０３１２５０
ＭＩＰＳ 値 ： ５７．２９７５２８
命令キャッシュミス （秒） ： ０．１２１１７７
オペランドキャッシュミス （秒） ： ０．０５６８１２
バンクコンフリクト時間 （秒） ： ０．１３６７０１
******************** ファイルの終り ********************
係数の最大値として６６９６０６、最小値として-６５４５８９がでていることがわかる。
そこで次の段階としてやるべきことはｎ＝３×５×７×１３×１７×１９×２３の場合である。これは未だ正常実行して計算しきっていません。しかし、部分的にこれらの係数たちを求めることはできました。それが以下。
******************** ファイルの始め ********************
ｃｙｃ_Ｇｒｙｔｃｚｕｋ_Ｔｒｏｐａｋ_６４．９０：
ｆ９０： ｖｅｃ（１）： ｃｙｃ_Ｇｒｙｔｃｚｕｋ_Ｔｒｏｐａｋ_６４．９０， ｌｉｎｅ １４：
Ｖｅｃｔｏｒｉｚｅｄ ｌｏｏｐ
ｆ９０： ｖｅｃ（１）： ｃｙｃ_Ｇｒｙｔｃｚｕｋ_Ｔｒｏｐａｋ_６４．９０， ｌｉｎｅ １８：
Ｖｅｃｔｏｒｉｚｅｄ ｌｏｏｐ
ｆ９０： ｖｅｃ（１）： ｃｙｃ_Ｇｒｙｔｃｚｕｋ_Ｔｒｏｐａｋ_６４．９０， ｌｉｎｅ ２２：
Ｐａｒｔｉａｌｌｙ ｖｅｃｔｏｒｉｚｅｄ ｌｏｏｐ
ｆ９０： ｖｅｃ（１）： ｃｙｃ_Ｇｒｙｔｃｚｕｋ_Ｔｒｏｐａｋ_６４．９０， ｌｉｎｅ ２３：
Ｖｅｃｔｏｒｉｚｅｄ ｌｏｏｐ
ｆ９０： ｖｅｃ（１）： ｃｙｃ_Ｇｒｙｔｃｚｕｋ_Ｔｒｏｐａｋ_６４．９０， ｌｉｎｅ ３２：
Ｖｅｃｔｏｒｉｚｅｄ ｌｏｏｐ
ｆ９０： ｖｅｃ（４）： ｃｙｃ_Ｇｒｙｔｃｚｕｋ_Ｔｒｏｐａｋ_６４．９０， ｌｉｎｅ ３７：
Ｖｅｃｔｏｒｉｚｅｄ ａｒｒａｙ ｅｘｐｒｅｓｓｉｏｎ．
ｆ９０： ｖｅｃ（４）： ｃｙｃ_Ｇｒｙｔｃｚｕｋ_Ｔｒｏｐａｋ_６４．９０， ｌｉｎｅ ３７：
Ｖｅｃｔｏｒｉｚｅｄ ａｒｒａｙ ｅｘｐｒｅｓｓｉｏｎ．
ｆ９０： ｃｙｃ_Ｇｒｙｔｃｚｕｋ_Ｔｒｏｐａｋ_６４．９０， _ＭＡＩＮ： Ｔｈｅｒｅ ａｒｅ ７ ｄｉａｇｎｏｓｅｓ．
ｃｈｅｃｋ_６４．ｆ９０：
ｆ９０： ｖｅｃ（１）： ｃｈｅｃｋ_６４．ｆ９０， ｌｉｎｅ １４： Ｖｅｃｔｏｒｉｚｅｄ ｌｏｏｐ
ｆ９０： ｃｈｅｃｋ_６４．ｆ９０，ｃｈｅｃｋ： Ｔｈｅｒｅ ｉｓ １ ｄｉａｇｎｏｓｉｓ．
ｍｏｅｂｉｕｓ_６４．f９０：
ｆ９０： ｖｅｃ（３）： ｍｏｅｂｉｕｓ_６４．f９０， ｌｉｎｅ ８ｎｖｅｃｔｏｒｉｚｅｄ ｌｏｏｐ．
ｆ９０： ｍｏｅｂｉｕｓ_６４．９０，ｍｏｅｂｉｕｓ： Ｔｈｅｒｅ ｉｓ １ ｄｉａｇｎｏｓｉｓ．
ｓｔａｒｔ ｔｉｍｅ＝Ｔｕｅ Ｎｏｖ ２３ １８：００：１１ＪＳＴ １９９９
ｍａｘ＝ ４０７１７７０３８７ ｍｉｎ＝ －４２４８４５１０８５
．．．．．プログラム 情報．．．．．
経過時間 （秒） ： ２８９４９．７１７４２１
ユーザ時間 （秒） ： ２７０１８．０２０４７８＝約７時間３０分
システム時間 （秒） ： ６４．４６６６９９
ベクトル命令実行時間 （秒） ： ２７０１５．０２８９２５
全命令実行時間 ： １５６１２５３５２９５３０
ベクトル命令実行数 ： ４１６２６８４９９４９９
ベクトル命令実行要素数 ： １０６５６１０１３２７９５８
浮動小数点データ実行要素数 ： ７３００３２２
MOPS 値 ： ３９８６．４５０３９３
MFLOPS 値 ： ０．０００２７０
平均ベクトル長 ： ２５５．９６１０５７
ベクトル演算率 （％） ： ９８．９３６９３５
メモリ使用量 （MB) ： ５６１．０３１２５０
ＭＩＰＳ 値 ： ５７．７８５６３７
命令キャッシュミス （秒） ： ２．７５４３４９
オペランドキャッシュミス （秒） ： １．２２２４６０
バンクコンフリクト時間 （秒） ： ０．６９２６１２
******************** ファイルの終り ********************
const＝euler/５として中途まで計算させました。このときのeulerは３６４９５３６０です。したがってconstは７２９９０７２です。数学書式で説明すると
Φ３×５×７×１１×１３×１７×１９×２３（ｘ）＝ａ７２９９９０７２（１１１５４６４３５９）ｘ７２９９0７２＋・・・
１１１５４６４３５＝３×５×７×１１×１３×１７×１９×２３
最初の項（定数項）からのｘ７２９９0７２項まで求めたことになる。少しばかり関心がひかれるのは係数達の最大値、最小値が４バイト整数としてはオーバーフロウしていることです。このプログラムはオーバーフロウの検証を行っていません。この検証を伴ったプログラムにすると実行時間が絶望的に大きくなるでしょう。この意味では数学者たちを未だ満足させていません。
３．並列プログラムの簡単な例
上記のプログラムの並列化を考えたが、まだ満足すべき並列アルゴリズムを考え出すに至っていない。その過程で次のプログラムを作ったので載せておきます。何かの参考にしてください。******************** ファイルの始め ********************
！並列処理機能の手引 ２０ページ 図３－５ 同時制御の使用例に
加工したもの
ｐｒｏｇｒａｍ ｍａｉｎ
ｕｓｅ ｍｏｄ_ｃｏｍ
！ ｉｎｔｅｇｅｒ ｔ１
ｉｎｔｅｇｅｒ ｔｉｄ，ｔｐａｒａｍ，ｄ１
ｅｘｔｅｒｎａｌ ｓｕｂ
！
ｅｖｅｎｔ
ｃａｌｌ ｐｅａｓｇｎ（ｅｖｅｎｔ）
！
ｔｐａｒａｍ＝１００
ｃａｌｌ ｐｔｆｏｒｋ（ｔｉｄ，ｔｐａｒａｍ，ｓｕｂ）
ｄ１＝５４３２１
ｘ＝ｄ１
ｃａｌｌ ｐｅｐｏｓｔ（ｅｖｅｎｔ）
ｃａｌｌ ｐｔｊｏｉｎ（ｔｉｄ）
！
ｃａｌｌ ｐｅｆｒｅｅ（ｅｖｅｎｔ）
！
！ｃｄｉｒ ｓｅｒｉａｌ
ｗｒｉｔｅ（６，*）’ｘ＝’，ｘ
！ｃｄｉｒ ｅｎｄｓｅｒｉａｌ
ｅｎｄ ｐｒｏｇｒａｍ ｍａｉｎ
******************** ファイルの終り ********************
モジュールmod_comは同部分にあるのと同じ。
また別に以下の例もあります。これは排他制御(lock)の使用例です。
******************** ファイルの始め ********************
ｐｒｏｇｒａｍ ｍａｉｎ
ｕｓｅ ｍｏｄ_ｃｏｍ
ｉｎｔｅｇｅｒ ｔ１，ｔｉｄ
ｅｘｔｅｒｎａｌ ｓｕｂ
！
ｌｏｃｋ＝０
ｃａｌｌ ｐｌａｓｇｎ（ｌｏｃｋ）
！
ｔ１＝２００
ｔｐａｒａｍ＝１２３４５
ｓ＝０
ｃａｌｌ ｐｔｆｏｒｋ（ｔｉｄ，ｔｐａｒａｍ，ｓｕｂ）
ｃａｌｌ ｐｔｊｏｉｎ（ｔｉｄ）
ｓ＝ｓ＋ｔ１
ｃａｌｌ ｐｌｕｎｌｏｃｋ（ｌｏｃｋ）
ｃａｌｌ ｐｔｊｏｉｎ（ｔｉｄ）
！
ｃａｌｌ ｐｌｆｒｅｅ（ｌｏｃｋ）
！
！ｃｄｉｒ ｓｅｒｉａｌ
ｗｒｉｔｅ（６，*）’ｓ＝’，ｓ
！ｃｄｉｒ ｅｎｄｓｅｒｉａｌ
ｅｎｄ ｐｒｏｇｒａｍ ｍａｉｎ
******************** ファイルの終り ********************
参考文献
- 小柴 洋一、円分多項式の係数の計算、大阪大学大型計算機センターニュース、Vol．21 no．2 1191-8，82 51-56p
- 小柴 洋一、円分多項式の係数の計算、九州大学大型計算機センター広報、Vol．30 no．2 1197，141-145p
- A．Grytczuk and B．Tropak，A numerical method for the determination of the cyclotomic polynomial coefficients Computional Number Theory，15-19page，Walter de Gruyter，1991
- ファン・デル・ヴェルデン（銀林 浩訳），現代代数学１，東京図書，110ページ、問題3
- Hardy and Wright，The Theory of Numbers，Oxford，4th editin，1959，238page，定理272