シェルスクリプトやコマンドラインである数が素数か、以下のように調べたい。
$ gmpxx-isprime 1907 && echo prime || echo composite prime
そうした場合、GNU MP を利用すればそういったコマンドラインプログラムを簡単に作る事ができます。
/*
isprime.cc
Copyright (C) 2015 Taiji Yamada <taiji@aihara.co.jp>
*/
#include <iostream>
#include <string>
#include <libgen.h>
#include <gmpxx.h>
int usage(std::ostream &os, int argc, char *argv[])
{
(void)argc;
os << "Usage:\n\t" << basename(argv[0]) << " [number ...]\n";
return os == std::cerr ? 1 : 0;
}
int main(int argc, char *argv[])
{
int a;
mpz_class n;
if (argc < 2) return usage(std::cerr, argc, argv);
for (a=1; a<argc; ++a) {
n = argv[a];
if (!mpz_probab_prime_p(n.get_mpz_t(), 25))
return 1;
}
return 0;
}
ここでは GNU MP の C++ API である gmpxx.h を利用しています。このように作成しておけば、数列がすべて素数であるかを判定するプログラムになります。
$ gmpxx-isprime 953 1907 && echo primes || echo composites primes
他にも以下のようなさまざまなプログラムを用意しました。接頭辞「gmpxx-」がインストールされた時のコマンド名になります。
isprime 素数判定 iscomposite 合成数判定 primes 素数を範囲で出力 primes_reverse 素数を範囲で逆順に出力 n_primes 素数を個数で出力 n_primes_reverse 素数を個数で逆順に出力 composites 合成数を範囲で出力 factorize 素因数分解 composite 合成数生成 gcd 最大公約数 lcm 最小公倍数 divisors 約数 n_divisors 約数の個数 sum_divisors 約数関数 next_prime 次の素数 previous_prime 前の素数 next_composite 次の合成数 previous_composite 前の合成数 integer_valued_polynomial_progression 整数値多項式の数列生成 prime_polynomial_progression 整数値多項式の素数生成 highlycomposites 高度合成数の数列生成 highlyabundants 高度過剰数の数列生成 superabundants 超過剰数の数列生成 naturals さまざまな自然数や整数の数列生成 totient_function トーシェント関数の生成 cototient_function コトーシェント関数の生成 mebius_function メビウス関数の生成 primes_util-collate C/C++ 用素数表作成及び検証
例えば、以下のように高度合成数を出力できます。
$ gmpxx-highlycomposites -u 1000000 1 2 4 6 12 24 36 48 60 120 180 240 360 720 840 1260 1680 2520 5040 7560 10080 15120 20160 25200 27720 45360 50400 55440 83160 110880 166320 221760 277200 332640 498960 554400 665280 720720
そこで、これを素因数分解するには以下のようにします。
$ gmpxx-factorize $(gmpxx-highlycomposites) 1 1^0 2 2^1 4 2^2 6 2^1 3^1 12 2^2 3^1 24 2^3 3^1 36 2^2 3^2 48 2^4 3^1 60 2^2 3^1 5^1 120 2^3 3^1 5^1 180 2^2 3^2 5^1 240 2^4 3^1 5^1 360 2^3 3^2 5^1 720 2^4 3^2 5^1 840 2^3 3^1 5^1 7^1 1260 2^2 3^2 5^1 7^1 1680 2^4 3^1 5^1 7^1 2520 2^3 3^2 5^1 7^1 5040 2^4 3^2 5^1 7^1 7560 2^3 3^3 5^1 7^1 10080 2^5 3^2 5^1 7^1 15120 2^4 3^3 5^1 7^1 20160 2^6 3^2 5^1 7^1 25200 2^4 3^2 5^2 7^1 27720 2^3 3^2 5^1 7^1 11^1 45360 2^4 3^4 5^1 7^1 50400 2^5 3^2 5^2 7^1 55440 2^4 3^2 5^1 7^1 11^1 83160 2^3 3^3 5^1 7^1 11^1 110880 2^5 3^2 5^1 7^1 11^1 166320 2^4 3^3 5^1 7^1 11^1 221760 2^6 3^2 5^1 7^1 11^1 277200 2^4 3^2 5^2 7^1 11^1 332640 2^5 3^3 5^1 7^1 11^1 498960 2^4 3^4 5^1 7^1 11^1 554400 2^5 3^2 5^2 7^1 11^1 665280 2^6 3^3 5^1 7^1 11^1 720720 2^4 3^2 5^1 7^1 11^1 13^1
また、コマンド「composite」は以下の形式の合成数を生成します。
\[c_1\times c_2\cdots\times v_1^{i_1}\times v_2^{i_2}\cdots\times w_1^{j_1}\times w_2^{j_2}\cdots\]
但しここで、\(0<i_i<m, 0\le j_i\le m\)。
$ gmpxx-composite -h
Usage:
gmpxx-composite [-c number ...] [-v number ...] [-w number ...] [-m times]
$ gmpxx-factorize $(gmpxx-composite -c 2 3 -v 5 7 -w 11 13 -m 2)
210 2^1 3^1 5^1 7^1
1050 2^1 3^1 5^2 7^1
1470 2^1 3^1 5^1 7^2
2310 2^1 3^1 5^1 7^1 11^1
2730 2^1 3^1 5^1 7^1 13^1
7350 2^1 3^1 5^2 7^2
11550 2^1 3^1 5^2 7^1 11^1
13650 2^1 3^1 5^2 7^1 13^1
16170 2^1 3^1 5^1 7^2 11^1
19110 2^1 3^1 5^1 7^2 13^1
30030 2^1 3^1 5^1 7^1 11^1 13^1
80850 2^1 3^1 5^2 7^2 11^1
95550 2^1 3^1 5^2 7^2 13^1
150150 2^1 3^1 5^2 7^1 11^1 13^1
210210 2^1 3^1 5^1 7^2 11^1 13^1
1051050 2^1 3^1 5^2 7^2 11^1 13^1
また、コマンド「totient_function」はオイラーのトーシェント関数\(\phi(n)\):
\[\phi(n) = n\prod_{k=1}^{d}\left(1-\frac{1}{p_k}\right)\]
但しここで、\(p_k\)は\(n\)の互いに素な\(d\)個の素因数。ここで分数があるので有理数の扱いが必要に思われますが、少し工夫すれば素因数分解の過程で割り切れる除算のみで実装できます。
mpz_class totient_function(mpz_class n)
{
mpz_class f(2), r(n);
for (; f <= n; mpz_nextprime(f.get_mpz_t(), f.get_mpz_t())) {
if (n % f == 0) {
r -= r / f;
n /= f;
}
while (n % f == 0)
n /= f;
}
return r;
}
さて、コマンド「naturals」は相当な説明を加えないとなりません。これは「-type name」オプションで以下のようなさまざまな数列を生成します。
natural 自然数 odd 奇数 even 偶数 singly_even 単偶数 doublly_even 複偶数 composite 合成数 prime 素数 sophie_germain_prime ソフィー・ジェルマン素数 interprime インター素数 twin_prime 双子素数 cousin_primes いとこ素数 sexy_primes セクシー素数 square_free 非平方数 cube_free 非立方数 perfect_square 平方数 perfect_power 累乗数 powerful 多冪数 highlycomposite 高度合成数 highlyabundant 高度過剰数 superabundant 超過剰数 abundant 過剰数 primitive_abundant 原始過剰数 primitive_abundant_have_no_abundant_proper_divisor 自身を除く過剰な約数をもたない原始過剰数* non_prime_unusual 非素数異常数* pseudoperfect, semiperfect 擬似完全数 primitive_pseudoperfect, primitive_semiperfect 原始擬似完全数 perfect 完全数 multiperfect 倍積完全数 deficient 不足数 superperfect 超完全数 hemiperfect ヘミ完全数 refactorable リファクタブル数 unitary_perfect ユニタリ完全数 primary_pseudoperfect プライマリ擬似完全数 weird 不思議数 practical 実際数* ruth_aaron_pair ルース＝アーロン・ペア primorial_ruth_aaron_pair 素数階乗ルース＝アーロン・ペア unusual 異常数 arithmetic 算術数* harmonic_divisor 調和数 round_number ラウンド数 lucas_carmichael ルーカス・カーマイケル数 semiprime 擬似素数 sphenic 楔数 amicable 友愛数 betrothed 婚約数、準友愛数 augumented_amicable 拡張友愛数 friendly フレンドリ数 mersenne メルセンヌ数 mersenne_prime メルセンヌ素数 wilson_quotient ウィルソン商 wilson_prime ウィルソン素数 fortunate フォーチューン数 hilbert ヒルベルト数 pronic 矩形数 power_of_two 2の冪数 cullen カレン数 cullen_prime カレン素数 woodall ウッダル数 woodall_prime ウッダル素数 regular 正規数* factorial_prime 階乗素数 euclid ユークリッド数 euclid_prime ユークリッド素数 kummer クンマー数 kummer_prime クンマー素数 primorial_prime 素数階乗素数 triangular 三角数 triangular_pyramidal, tetrahedral 三角錐数、四面体数 pentatope 五胞体数 pentagonal 五角数 square_triangular 平方三角数 squared_triangular 平方化三角数* lazy_caterer 怠惰な仕出し屋数* jacobsthal ジェイコブザル数 jacobsthal_lucas ジェイコブザル・ルカ数 kynea カイネア数 prime_kynea カイネア素数 thabit サビット数 prime_thabit サビット素数 fermat フェルマー数 catalan カタラン数 derangement 完全順列 nontotient 非トーシェント数 noncototient 非コトーシェント数 perfect_totient 完全トーシェント数 highly_totient 高度トーシェント数 highly_cototient 高度コトーシェント数 odd_semifactrial 奇数の半階乗数 even_semifactrial 偶数の半階乗数 fibonacci フィボナッチ数 leonardo レオナルド数 lucas リュカ数 pell ペル数 perrin ペラン数 padovan パドバン数 motzkin モズキン数 schreder_hipparchus シュレーダー・ヒパルカス数 sylvester シルベスター数 telephone 電話数 alternating_factorial 交互階乗数 hyperfactorial 超階乗数 exponential_factorial 指数階乗数
ここでアスタリスク「*」付きの日本語名は筆者の命名です。いくつか実装に工夫が必要なものを紹介しましょう。
調和数(harmonic divisor numbers)は自然数\(n\)のうち\(k\)個すべての約数\(d_i\)の調和平均:
\[\frac{k}{\sum_{i=1}^{k}\frac{1}{d_i}}\]
が整数値をもつ\(n\)のことです。ここで分数があるので有理数の扱いが必要に思われますが、最小公倍数を使って工夫すれば最終的に割り切れるか否かで実装できます。
bool is_harmonic_divisor(const mpz_class n)
{
std::vector<mpz_class> q = divisors(n);
mpz_class d(0), c(q[0]);
size_t i;
for (i=1; i<q.size(); ++i)
mpz_lcm(c.get_mpz_t(), c.get_mpz_t(), q[i].get_mpz_t());
for (i=0; i<q.size(); ++i)
d += c/q[i];
return (mpz_class(q.size())*c) % d == 0;
}
擬似完全数(pseudoperfect numbers)は自然数\(n\)のうち自身を除くいくつかの約数\(d_i\)の和が元の数に等しい数\(n\)のことです。この「いくつかの」というのが厄介に思われますが、ビット演算で約数にフラグを割り当てすべての組み合わせを検証すれば実装できます。
bool is_pseudoperfect(const mpz_class n)
{
std::vector<mpz_class> d = proper_divisors(n);
mpz_class i, c;
size_t j;
mpz_ui_pow_ui(c.get_mpz_t(), 2, d.size());
for (i=1; i<c; ++i) {
mpz_class s(0);
for (j=0; j<d.size(); ++j)
if (mpz_tstbit(i.get_mpz_t(), j))
s += d[j];
if (s == n)
return true;
}
return false;
}
さて、ここまで GNU MP を利用した実装を紹介して来ましたが、10桁程度の数値でも速度が遅く、その程度の桁数の整数なら C/C++ の整数で十分ですし、実行速度が全く違います。そこで C/C++ 用の初等整数コマンドラインツール群を用意しました。
但し、素数については実行速度を実用的なものにするため「/var/tmp/isprime_bit_table.bin」を作成します。作成途中で中断すると不完全な素数表となり、正しく動作しません。よって、初回のみ以下のように作成及び検証をじっくりと行うことをお勧めします。
$ gmpxx-prime_util-collate -u 2000000000 && echo valid || echo invalid
これで「/var/tmp/isprime_bit_table.bin」は 1.2GB ほど、Core i7 で、素数表の作成は５分ほどで終わるのですが、GNU MP による検証で１時間ほど掛かります。
さて、GNU MP には mpz_perfect_square_p、mpz_perfect_power_p が用意されているので、上記の「平方数(perfect square numbers)」「累乗数(perfect power numbers)」はそれを利用すればよかったのです。しかし、C/C++ のみで実装するとなると、これに変わるものを自前で用意しなければなりません。
こちらの場合、素数表があるので素因数分解が使えますし、以下の最大公約数 boost::math::gcd を利用した実装で事足ります。
template <typename T>
bool is_perfect_square(T n)
{
T f(2), c, c0(0);
isprime_bit_table_class isprime_bit_table;
for (; f <= n; f = isprime_bit_table.nextprime(f)) {
c = 0;
while (n % f == 0) {
++c;
n /= f;
}
if (c % 2 == 1)
return false;
if (c0 != 0) {
if (boost::math::gcd(c0, c) % 2 == 1)
return false;
c0 = boost::math::gcd(c0, c);
}
else
c0 = c;
if (!(f < isprime_bit_table.get_upper())) throw std::length_error("isprime_bit_table: over the length");
}
return true;
}
template <typename T>
bool is_perfect_power(T n)
{
T f(2), c, c0(0);
isprime_bit_table_class isprime_bit_table;
for (; f <= n; f = isprime_bit_table.nextprime(f)) {
c = 0;
while (n % f == 0) {
++c;
n /= f;
}
if (c == 1)
return false;
if (c0 != 0) {
if (boost::math::gcd(c0, c) == 1)
return false;
c0 = boost::math::gcd(c0, c);
}
else
c0 = c;
if (!(f < isprime_bit_table.get_upper())) throw std::length_error("isprime_bit_table: over the length");
}
return true;
}