数論・素数に関するメソッドを実装する

数論・素数に関するメソッドを C# で一通り実装してまとめました。
通常の数値計算プログラミング用にも、また競技プログラミング用のライブラリとしても使えます。

機能は次の通りです。

Euclid 互除法

  • 最大公約数 (GCD)
  • 最小公倍数 (LCM)
  • ax + by = 1 の解 (Euclid 互除法の拡張)

素数

下の O(x) は計算量のオーダーを表します。

  • 素因数分解
    • O(√n)
  • 約数の列挙
    • O(√n)
  • 素数判定
    • O(√n)
  • n 以下の素数の列挙
    • およそ O(n)
  • m 以上 M 以下の素数の列挙
    • およそ O(√M + (M – m))

ソースコードは以下の通りです。
ここではシンプルな実装とし、例外処理も全体的に省略しています。
さらに無駄な演算処理を除くように改良することで、もう少し速くすることはできるはずです。

using System;
namespace AlgorithmLib.Numerics
{
public static class Euclidean
{
public static int Gcd(int a, int b) { for (int r; (r = a % b) > 0; a = b, b = r) ; return b; }
public static int Lcm(int a, int b) => a / Gcd(a, b) * b;
public static long Gcd(long a, long b) { for (long r; (r = a % b) > 0; a = b, b = r) ; return b; }
public static long Lcm(long a, long b) => a / Gcd(a, b) * b;
// ax + by = 1 の解
// 前提: a と b は互いに素。
// ax + by = GCD(a, b) の解を求める場合、予め GCD(a, b) で割ってからこの関数を利用します。
public static (long x, long y) ExtendedEuclid(long a, long b)
{
if (b == 0) throw new ArgumentOutOfRangeException(nameof(b));
if (b == 1) return (1, 1 - a);
var q = Math.DivRem(a, b, out var r);
var (x, y) = ExtendedEuclid(b, r);
return (y, x - q * y);
}
}
}
view raw Euclidean.cs hosted with ❤ by GitHub
using System;
using System.Collections.Generic;
namespace AlgorithmLib.Numerics
{
public static class Primes
{
/// <summary>
/// 指定された自然数を素因数分解します。O(√n)
/// </summary>
/// <param name="n">自然数。</param>
/// <returns>素因数のコレクション。n = 1 の場合は空の配列。</returns>
public static long[] Factorize(long n)
{
var r = new List<long>();
for (long x = 2; x * x <= n && n > 1; ++x)
while (n % x == 0)
{
r.Add(x);
n /= x;
}
// √n を超える素因数はたかだか 1 個であり、その次数は 1。
if (n > 1) r.Add(n);
return r.ToArray();
}
/// <summary>
/// 指定された自然数の約数をすべて求めます。O(√n)
/// </summary>
/// <param name="n">自然数。</param>
/// <returns>約数のコレクション。</returns>
public static long[] Divisors(long n)
{
var r = new List<long>();
for (long x = 1; x * x <= n; ++x)
if (n % x == 0) r.Add(x);
for (int i = r.Count - 1; i >= 0; --i)
{
var v = n / r[i];
if (v > r[i]) r.Add(v);
}
return r.ToArray();
}
/// <summary>
/// 指定された自然数が素数かどうかを判定します。O(√n)
/// </summary>
/// <param name="n">自然数。</param>
/// <returns>素数である場合に限り true。</returns>
public static bool IsPrime(long n)
{
for (long x = 2; x * x <= n; ++x)
if (n % x == 0) return false;
return n > 1;
}
/// <summary>
/// 指定された自然数以下の素数をすべて求めます。およそ O(n)
/// </summary>
/// <param name="n">自然数。</param>
/// <returns>素数のコレクション。</returns>
public static long[] GetPrimes(long n)
{
var b = new bool[n + 1];
for (long p = 2; p * p <= n; ++p)
if (!b[p])
for (long x = p * p; x <= n; x += p)
b[x] = true;
var r = new List<long>();
for (long x = 2; x <= n; ++x) if (!b[x]) r.Add(x);
return r.ToArray();
}
/// <summary>
/// 指定された範囲内の素数をすべて求めます。およそ O(√M + (M - m))
/// </summary>
/// <param name="m">最小値。</param>
/// <param name="M">最大値。</param>
/// <returns>素数のコレクション。</returns>
public static long[] GetPrimes(long m, long M)
{
var rM = (int)Math.Ceiling(Math.Sqrt(M));
var b1 = new bool[rM + 1];
var b2 = new bool[M - m + 1];
if (m == 1) b2[0] = true;
for (long p = 2; p <= rM; ++p)
if (!b1[p])
{
for (var x = p * p; x <= rM; x += p) b1[x] = true;
for (var x = Math.Max(p, (m + p - 1) / p) * p; x <= M; x += p) b2[x - m] = true;
}
var r = new List<long>();
for (int x = 0; x < b2.Length; ++x) if (!b2[x]) r.Add(m + x);
return r.ToArray();
}
}
}
view raw Primes.cs hosted with ❤ by GitHub

使用例を次に示します。

using System;
using AlgorithmLib.Numerics;
namespace AlgorithmConsole
{
class Program
{
static void Main(string[] args)
{
// 素因数分解
Console.WriteLine(string.Join(" * ", Primes.Factorize(2020)));
// 2 * 2 * 5 * 101
// 約数の列挙
Console.WriteLine(string.Join(" ", Primes.Divisors(2020)));
// 1 2 4 5 10 20 101 202 404 505 1010 2020
// 素数判定
Console.WriteLine(Primes.IsPrime(1000000007));
// True
// n 以下の素数の列挙
Console.WriteLine(string.Join(" ", Primes.GetPrimes(100)));
// 2 3 5 7 11 13 17 19 23 29 31 37 41 43 47 53 59 61 67 71 73 79 83 89 97
// m 以上 M 以下の素数の列挙
Console.WriteLine(string.Join(" ", Primes.GetPrimes(1000000000, 1000000100)));
// 1000000007 1000000009 1000000021 1000000033 1000000087 1000000093 1000000097
}
}
}
view raw Program.cs hosted with ❤ by GitHub

前回: 二分探索のライブラリ化

作成したサンプル

バージョン情報

  • C# 7.0
  • .NET Standard 2.0

参照

カテゴリー: アルゴリズム, 数学. タグ: . 2 Comments »

二分探索のライブラリ化

// Competitive Programming (2) Advent Calendar 2019 の 20 日目の記事です。

競技プログラミングでたびたび利用する二分探索について考察してみます。
といっても、とくに競技プログラミングに限らず活用できると思います。言語は C# です。

まず二分探索の基本的な例として、昇順にソートされた整数のリスト a に対して「値 v をリストに挿入するときのインデックス」を求めてみます。 これは「値 v より大きい値が最初に現れるインデックス」と考え、最後の要素が v 以下のときはリストの長さ a.Count を返すことにすれば、次の IndexForInsert メソッドとして実装できます。

using System;
using System.Collections.Generic;
namespace AlgorithmLib
{
public static class SearchSample0
{
// 指定された値よりも大きい値を持つ最初のインデックスを求めます。
// これは、挿入先のインデックスを意味します。
public static int IndexForInsert(IList<int> a, int v)
{
int l = 0, r = a.Count, m;
while (l < r)
{
m = l + (r - l - 1) / 2;
if (a[m] > v) r = m;
else l = m + 1;
}
return r;
}
// 指定された値以上の値を持つ最初のインデックスを求めます。
// Array.BinarySearch メソッドと同様に、一致する値が存在しない場合はその補数を返します。
// ただし Array.BinarySearch メソッドでは、一致する値が複数存在する場合にその最初のインデックスを返すとは限りません。
public static int IndexOf(IList<int> a, int v)
{
int l = 0, r = a.Count, m;
while (l < r)
{
m = l + (r - l - 1) / 2;
if (a[m] >= v) r = m;
else l = m + 1;
}
return r < a.Count && a[r] == v ? r : ~r;
}
}
}
view raw SearchSample0.cs hosted with ❤ by GitHub

同じ考え方を利用して、Array.BinarySearch メソッドのように「値が一致するインデックスを求める。一致しない場合は、それより大きな値のインデックスの補数」という仕様にするには、上の IndexOf メソッドのように書けます。

注意点:

  • IList<T> インターフェイスは、 T[]List<T> クラスを含みます。
  • リスト内の値は重複してもかまいません。
  • 例外処理は除いています。
  • (l + r - 1) / 2 とするとオーバーフローすることがあるため、 l + (r - l - 1) / 2 としています。

上記のメソッドは、リストのみ、しかも昇順にソート済みという制約の下で適用することができます。
そこで、さらに汎用的に使えるように、抽象度を上げてライブラリを作ることを考えます。

二分探索により実現できることは、見方を変えると「状態が変化する境界を見つけること」、厳密には「与えられた条件に対し、ある境界が存在して、その前方と後方で真偽が逆転する場合にその境界を見つけること」だと考えます。
この考え方で先ほどのコードから抽象化した部分を抜き出すと、次の First メソッドのようになります。

using System;
namespace AlgorithmLib
{
/// <summary>
/// 二分探索のためのメソッドを提供します。
/// </summary>
public static class BinarySearch
{
/// <summary>
/// 条件 f を満たす最初の値を探索します。
/// [l, x) 上で false、[x, r) 上で true となる x を返します。
/// f(l) が true のとき、l を返します。
/// f(r - 1) が false のとき、r を返します。
/// </summary>
/// <param name="f">半開区間 [l, r) 上で定義される条件。</param>
/// <param name="l">探索範囲の下限。</param>
/// <param name="r">探索範囲の上限。</param>
/// <returns>条件 f を満たす最初の値。</returns>
public static int First(Func<int, bool> f, int l, int r)
{
int m;
while (l < r)
if (f(m = l + (r - l - 1) / 2)) r = m;
else l = m + 1;
return r;
}
}
}
view raw BinarySearch.cs hosted with ❤ by GitHub
using System;
using System.Collections.Generic;
namespace AlgorithmLib
{
public static class SearchSample1
{
// 指定された値よりも大きい値を持つ最初のインデックスを求めます。
// これは、挿入先のインデックスを意味します。
public static int IndexForInsert(IList<int> a, int v) =>
BinarySearch.First(i => a[i] > v, 0, a.Count);
// Array.BinarySearch メソッドと異なる点: 一致する値が複数存在する場合は最初のインデックス。
public static int IndexOf(IList<int> a, int v)
{
var r = BinarySearch.First(i => a[i] >= v, 0, a.Count);
return r < a.Count && a[r] == v ? r : ~r;
}
}
}
view raw SearchSample1.cs hosted with ❤ by GitHub

この形にすれば、降順のリストにも、リスト以外の数値の範囲にも適用できます。
例として、AtCoder の ABC 146 C – Buy an Integer を解いてみます。
先ほどの First メソッドに対して Last メソッドを用意し (コメントは省略)、価格が所持金以下となる最後の整数を求めればよいです。

using System;
class C
{
static void Main()
{
var h = Array.ConvertAll(Console.ReadLine().Split(), long.Parse);
Console.WriteLine(Last(n => h[0] * n + h[1] * n.ToString().Length <= h[2], 0, 1000000000));
}
// 条件 f を満たす最後の値を探索します。
static int Last(Func<int, bool> f, int l, int r)
{
int m;
while (l < r) if (f(m = r - (r - l - 1) / 2)) l = m; else r = m - 1;
return l;
}
}
view raw C.cs hosted with ❤ by GitHub

だいぶ簡潔に書けました。

次に小数向けのメソッドを用意して、二分法で平方根を求めてみます。
誤差を表す桁数を指定できるようにしています。

using System;
namespace AlgorithmLib
{
/// <summary>
/// 二分探索のためのメソッドを提供します。
/// </summary>
public static class BinarySearch
{
/// <summary>
/// 条件 f を満たす最初の値を指定された誤差の範囲内で探索します。
/// (l, x) 上で false、[x, r) 上で true となる x を返します。
/// l 近傍で true のとき、l を返します。
/// r 近傍で false のとき、r を返します。
/// </summary>
/// <param name="f">開区間 (l, r) 上で定義される条件。</param>
/// <param name="l">探索範囲の下限。</param>
/// <param name="r">探索範囲の上限。</param>
/// <param name="digits">誤差を表す小数部の桁数。</param>
/// <returns>条件 f を満たす最初の値。</returns>
public static double First(Func<double, bool> f, double l, double r, int digits = 9)
{
double m;
while (Math.Round(r - l, digits) > 0)
if (f(m = l + (r - l) / 2)) r = m;
else l = m;
return r;
}
}
}
view raw BinarySearch.cs hosted with ❤ by GitHub
using System;
namespace AlgorithmLib
{
public static class SearchSample1
{
public static double Sqrt(double v, int digits = 9) =>
BinarySearch.First(x => x * x >= v, Math.Min(1, v), Math.Max(1, v), digits);
}
}
view raw SearchSample1.cs hosted with ❤ by GitHub

まとめ

というわけで、今後は二分探索の問題を意味論的に「与えられた条件を満たす最初または最後の値を求める」と考えれば、
今回作成したライブラリを使って実装ができます。

その他

  • これまでは二分探索を実装するときに再帰を使っていたのですが、おそらく過去に読んだ解説に再帰を使うものが多かったためだと思います。競技プログラミングを始めてから自然に while (l < r) を使うようになりました。
  • Array.BinarySearch メソッドList.BinarySearch メソッドでは、重複する値を検索する場合、その最初のインデックスが返ってくるとは限りません。

前回: 競技プログラミングでも C# で簡潔に書きたい
次回: 数論・素数に関するメソッドを実装する

参照

カテゴリー: アルゴリズム. タグ: . 2 Comments »

競技プログラミングでも C# で簡潔に書きたい

// Competitive Programming (2) Advent Calendar 2019 の 15 日目の記事です。

競技プログラミングの AtCoder というものを今年の8月に始めて、3~4か月をかけて水色に到達しました。

20191124-Rating

スピード勝負は得意ではないのですが、難しめの問題を少し時間をかけて解くのが向いているようです。
また、簡潔なコードを書くことを心掛けていて、ショートコード C# 部門があるとしたらだいたい優勝していると思います。
(スペースを切り詰めたりはしません。むしろ Visual Studio で既定のフォーマットをしてから提出しています。)

20191201-E

20191201-F

 

さて、12月1日に実施された 三井住友信託銀行プログラミングコンテスト2019 では、別解が多くあり考察を楽しめる問題セットでした。
全体的にだいぶ簡潔に書ける回だったので、その中から問題 B, C, D を解説してみたいと思います。

B – Tax Rate

税込価格で N 円となるときの税抜価格はいくらか。

小数を考慮しなければならないため、第一感ではどんな値が適合するのかもやっとしますが、
そんなときは 1 から N までの全探索でも正解できます。LINQ を使って次のように書けます。

using System;
using System.Linq;
class B
{
static void Main()
{
var n = int.Parse(Console.ReadLine());
var x = Enumerable.Range(1, n).FirstOrDefault(i => (int)(1.08 * i) == n);
Console.WriteLine(x > 0 ? $"{x}" : ":(");
}
}
view raw B.cs hosted with ❤ by GitHub

しっかり考察してみると実は、税抜価格が1だけ離れると税込価格も1以上離れることから税抜価格と税込価格は1対1に対応し、
税込価格 N に対する税抜価格はたかだか1個存在することがわかります。
また、 1.08X ≥ N つまり X ≥ N / 1.08 であることから、X の候補は N / 1.08 の天井 (ceiling) に限られます。
(これより 1 以上大きい整数の税込価格は N+1 以上になってしまう。)
あとはこれの税込価格が実際に N に一致するかどうかを確かめれば OK です。

using System;
class B
{
static void Main()
{
var n = int.Parse(Console.ReadLine());
var x = Math.Ceiling(n / 1.08);
Console.WriteLine((int)(1.08 * x) == n ? $"{x}" : ":(");
}
}
view raw B.cs hosted with ❤ by GitHub

 

C – 100 to 105

合計価格がちょうど X 円となるような買い物ができるか。

商品の価格の差が1円単位なので、商品の個数が N のとき、合計価格の候補は 100N ≤ X ≤ 105N の範囲の整数値となります。
そこでこの問題も、やはり安全策の全探索で通ります。

using System;
using System.Linq;
class C
{
static void Main()
{
var x = int.Parse(Console.ReadLine());
Console.WriteLine(Enumerable.Range(1, 1000).Any(i => 100 * i <= x && x <= 105 * i) ? 1 : 0);
}
}
view raw C.cs hosted with ❤ by GitHub

では今度もまたしっかり考察してみます。
商品の個数 N が [X / 100] より大きい場合は、合計価格は X を超えてしまいます。
また、N が [X / 100] より小さい場合に 100N ≤ X ≤ 105N を満たしているとすると、N+1 もこれを満たしているはずです。
したがって N = [X / 100] の場合のみを調べればよく、次のようなコードで書けます。

using System;
class C
{
static void Main()
{
var x = int.Parse(Console.ReadLine());
Console.WriteLine(x <= x / 100 * 105 ? 1 : 0);
}
}
view raw C.cs hosted with ❤ by GitHub

 

D – Lucky PIN

半角数字からなる文字列 S から3文字を取り出してその順に並べたものは何種類あるか。

例によって 000 から 999 までを全探索します。
たぶんこれが一番短いと思います (当社比)

using System;
using System.Linq;
class D
{
static void Main()
{
Console.ReadLine();
var s = Console.ReadLine();
int i;
Console.WriteLine(Enumerable.Range(0, 1000).Select(x => $"{x:D3}").Count(p => (i = s.IndexOf(p[0])) > -1 && (i = s.IndexOf(p[1], i + 1)) > -1 && s.IndexOf(p[2], i + 1) > -1));
}
}
view raw D.cs hosted with ❤ by GitHub

いちおう次に計算量を考慮して別解を考えてみましょう。
まず前処理として、どの文字がどのインデックスに現れるかを取得できるように辞書を作っておきます。
この辞書のキーはたかだか10種類です。
そして1文字目から順に、選択可能なインデックスの範囲から文字を選択していきます。
3文字目を選ぶときは、候補のうち最後のインデックスが2文字目のインデックスより大きければ成立します。

多重の foreach 文を使って探索してもよいのですが、今回は C# のクエリ式を使ってみます。
クエリ式は普段はそれほど使わないのですが、このようにコレクションの入れ子構造を平坦化 (flatten) するときは見やすくなると思います。通常の LINQ では SelectMany メソッドに相当します。

using System;
using System.Linq;
class D
{
static void Main()
{
Console.ReadLine();
var d = Console.ReadLine().Select((c, i) => new { c, i }).ToLookup(_ => _.c, _ => _.i);
var q = from i1 in d.Select(g => g.First())
from i2 in d.Select(g => g.FirstOrDefault(x => x > i1)).Where(x => x > 0)
select d.Count(g => g.Last() > i2);
Console.WriteLine(q.Sum());
}
}
view raw D.cs hosted with ❤ by GitHub

計算量は、1文字目が10、2文字目が N (たかだか N 文字しかない)、3文字目が10のため、全体で O(100・N) です。
また、2文字目のインデックスを探索するところを二分探索にすれば O(1000・logN) となります。
ただし、前処理が少なくとも O(N) であるためそちらのほうが影響するかもしれません。

次回: 二分探索のライブラリ化

以前に投稿した記事