JavaのModIntを考える

序論

プログラミングコンテストでは「答えの数を1000000007で割った余りを求めよ」という問題が大変よく出題される。1000000007は素数で、ほかにも 998244353 が代わりに出題されることもよくある。これらは真の解は大きくなりすぎて計算できない場合でも、加減乗算をするたびにその素数で割った余りを求めてあげれば解が変わらない。

ModInt

ソースコード中では問題で与えられた大きな素数は定数変数の mod と書くことにする。
a * (b + c) を求めたいときに、演算の度に余りを求めると

ans = a * ((b + c) % mod) % mod;

複雑な式になると %mod を書き忘れてオーバーフローで誤答になる危険があるので何とか記法を工夫したい。このモチベーションはとても大きい。

課題点

演算子オーバーロードがない

+ - * などの二項演算子オーバーロードできるプログラミング言語も存在するが Java はできない。そのためどう工夫しても数式的な中置記法をあきらめなければならない。

new のコスト

(a + b) % mod を求めたいだけなのに new で Object を作るのは実行時間に影響があるかもしれないため、new を使う場合はそのオーバーヘッドを知っておきたい。

結果

適当に加算と乗算を1億回するコードを、mod の行を工夫して実行時間を測る。
テストコードはこれ
ModIntTest · GitHub

for (int i = 0; i < a.length; i++) {
    for (int j = 0; j < b.length; j++) {
        ans = (ans * ((a[i] + b[j]) % mod)) % mod;
    }
}

単位はms。

Test AtCoder Codeforces Ideone
1 simple 439 1898 406
2 non final 803 1978 2053
3 function 444 1944 661
4 function2 486 2048 416
5 mint 582 3843 1317
6 mutable 560 2519 700

なぜか Codeforcesが非常に遅い。Codeforcesはそもそも除算・剰余算が驚くほど遅いため、記法の工夫の以前の問題として、演算回数をモンゴメリ乗算などで減らす方がいい。
AtCoder は new が予想よりも早い。正直、どの実装でもいい。
Ideone が最も予想に近い比率の時間になる。関数にすれば遅くなるし、newを使う回数だけもっと遅くなる。

テスト

Test1 simple

for (int i = 0; i < a.length; i++) {
    for (int j = 0; j < b.length; j++) {
        ans = (ans * ((a[i] + b[j]) % mod)) % mod;
    }
}

modを全て手打ちする。typoをしない限りは簡単かつ高速。
加算は mod の2倍を超えないので「剰余算」の代わりに「条件分岐と減算」でもできるが AtCoder では実行時間はほぼ変わらない。遅くなることも有る。

Test2 non final

// final を書き忘れる
long mod = 1000000007; // NG

final にするだけで速くなるので final は付けよう。もしくは

long mod() { return 1000000007; } // GOOD

のような定数関数にしても高速。オーバーライドしても高速なので 998244353 と共存させることもできる。

Test3 functions

ans = mul(ans, add(a[i], b[j]));

ただしstatic関数を定義する。

class ModIntStatic {
    static final long mod = 1000000007;

    static long add(long x, long y) {
        return (x + y) % mod;
    }

    static long mul(long x, long y) {
        return x * y % mod;
    }
}

関数記法になってしまうが実装が軽く実行速度が速く、可読性もよい。long のままなので関数を通し忘れるとオーバーフローの危険があるが、% mod を書き忘れるよりは防ぎやすい気がする。

Test4 function2

interface ModInt {
    public long mod();

    class ModInt998244353 implements ModInt { public long mod() { return 998244353; } }
    class ModInt1000000007 implements ModInt { public long mod() { return 1000000007; } }

    default long add(long x, long y) {
        return (x + y) % mod();
    }

    default long mul(long x, long y) {
        return x * y % mod();
    }
}

static ではなくなったが同様に高速。オーバーライドで複数の素数に対応できる。

Test5 Mint

ans = ans.mul(a[i].add(b[j]));

ただし、ans, a[i], b[j] は Mint オブジェクト。

class Mint {
    static final long mod = 1000000007;

    final long x;
    Mint(long x) { this.x = x; }

    Mint add(Mint y) {
        return new Mint((x + y.x) % mod);
    }

    Mint mul(Mint y) {
        return new Mint(x * y.x % mod);
    }
}

メソッドにしたので語順が中置記法になった。クラス同士でしか演算できないため、最も安全だと思われる。AtCoder は new がなぜか高速なのでこれも十分良いが、ほかのコンテストサイトでは遅いかもしれない。BigInteger と同じ API にできるのも良い。

Test6 Mutable Mint

tmp.assign(a[i]).addAssign(b[j]);
ans.mulAssign(tmp);

ただし ans, tmp, a[i], b[j] は Mint オブジェクト。

class Mint {
    static final long mod = 1000000007;

    long x;

    Mint(long x) { this.x = x; }

    Mint assign(Mint y) {
        x = y.x;
        return this;
    }

    Mint addAssign(Mint y) {
        x = (x + y.x) % mod;
        return this;
    }

    Mint mulAssign(Mint y) {
        x = x * y.x % mod;
        return this;
    }
}

new がなくなったのでどんな環境でも速いはず。ただし、数式が複雑になるだけ tmp 変数を用意しなければならない。また、変更可能なオブジェクトは関数の引数に渡したときに関数内で書き換えられることを心配する必要があるので、慎重に使わなければならない。

プリミティブタイプはJavaのジェネリクスに入らない

int, long, doubleなどのプリミティブな型を使うことは当然だ。これらに対して同じデータ構造やアルゴリズムを適用したいことがあるのも当然だ。
ジェネリクスがあるのだからint, long, doubleで共通な計算は実装できそうなものだが、できない。

例えば配列の和を求めるクラスを作りたいと思う。

// error
class Sum<T> {
    T sum(T[] a) {
        T ret = 0;
        for (int i = 0; i < a.length; i++) {
            ret += a[i];
        }
        return ret;
    }
}

困難① ジェネリクスにはプリミティブ型を入れられない

そのため `new Sum<long>().sum(array)` のように書けない
では long の代わりにラッパーである Long を使えばいいかと言うとそうではない。

困難② ジェネリクス型は `new` が使えない、`+=` が使えない

Integer型やLong型では演算子が使えるが、ジェネリクスになると利用できない。Javaジェネリクスはそういうもの。

解決

そんな不便なことがあるだろうか?じゃあ`Arrays.sort`はどう実装してるんだろう?JDKソースコードを見てみる。

public static void sort(int[] a)

プリミティブ型の`sort`や`binarySearch`などはそれぞれの型に対して全く同じ実装が(型だけ変えて)繰り返し記述されている。これをマネすれば良さそう。

プログラミングコンテストで使う場合はどうせソースコードをコピー&ペーストする。よって long 型だけ用意しておいて double 型が必要ならエディタで置換するのが最も良い方法だと思う。

疑似乱数メルセンヌツイスタは連続 624 項が分かれば完全に再現できる

標準ライブラリのメルセンヌツイスタは実装と定数が公開されている。すると、出力された疑似乱数列から内部状態をコピーすることができる。

メルセンヌツイスタの実装

英語版の wikipedia の実装を見る。
en.wikipedia.org

定数

標準ライブラリは周期が 2^19937 - 1 になるような定数があらかじめ設定されている。

関数

配列 unsigned mt[624] は長い周期で変化するカウンタの役割をする。mt[624] は実質リングバッファで、index でずれていく。
wikipedia の実装では「初めて乱数を生成するまで初期化しない」「twistは 624 回に 1 回だけまとめて行う」ような効率化がされていた。

乱数を1つ生成する関数の中で tempering (調律)をしている

unsigned y = mt[index];
y ^= (y >> U) & D;
y ^= (y << S) & B;
y ^= (y << T) & C;
y ^= (y >> L);

最終的にこの y を出力するので、tempering の逆関数を作ると配列 mt をコピーできるはず。

逆関数

定数に依存すれば結構簡単な逆関数ができる。

static unsigned inverse_tempering(unsigned y) {
    // inverse y ^= y >> L;
    y ^= y >> L;

    // inverse y ^= (y << T) & C; 
    y ^= (y << T) & C;

    // inverse y ^= (y << S) & B;
    y ^= (y << S) & 0x00001680;
    y ^= (y << S) & 0x000C4000;
    y ^= (y << S) & 0x0D200000;
    y ^= (y << S) & 0x90000000;

    // inverse y ^= (y >> U) & D;
    y ^= (y >> U) & 0x001FFC00;
    y ^= (y >> U) & 0x000003FF;

    return y;
}

コード全体

struct Random {
private:
    static const int W = 32; // word size;
    static const int N = 624; // degree of recurrence;
    static const int M = 397; // shift size;
    static const int R = 31; // mask bits;
    static const unsigned A = 0x9908B0DF; // XOR mask;
    
    // tempering;
    static const int U = 11;
    static const unsigned D = 0xFFFFFFFF;
    static const int S = 7;
    static const unsigned B = 0x9D2C5680;
    static const int T = 15;
    static const unsigned C = 0xEFC60000;
    static const int L = 18;

    static const unsigned F = 1812433253; // initialization;

    static const unsigned LOWER_MASK = (1u << R) - 1u;
    static const unsigned UPPER_MASK = ~0u ^ LOWER_MASK;

    static const unsigned DEFAULT_SEED = 5489;

    unsigned mt[N];
    int index;

    void twist(int i) {
	unsigned x = (mt[i] & UPPER_MASK) | (mt[(i + 1) % N] & LOWER_MASK);
	unsigned xA = x >> 1;
	if (x % 2 != 0) {
	    xA ^= A;
	}
	mt[i] = mt[(i + M) % N] ^ xA;
    }

    void initialize(unsigned seed) {
	mt[0] = seed;
	for (int i=1; i<N; i++) {
	    mt[i] = F * (mt[i-1] ^ (mt[i-1] >> (W - 2))) + i;
	}
	for (int i=0; i<N; i++) {
	    twist(i);
	}
	index = 0;
    }

    unsigned extract_number() {
	if (index >= N) {
	    index = 0;
	}

	unsigned y = mt[index];
	y ^= (y >> U) & D;
	y ^= (y << S) & B;
	y ^= (y << T) & C;
	y ^= (y >> L);

	twist(index);

	index++;
	return y;
    }

    static unsigned inverse_tempering(unsigned y) {
	// inverse y ^= y >> L;
	y ^= y >> L;

	// inverse y ^= (y << T) & C; 
	y ^= (y << T) & C;

	// inverse y ^= (y << S) & B;
	y ^= (y << S) & 0x00001680;
	y ^= (y << S) & 0x000C4000;
	y ^= (y << S) & 0x0D200000;
	y ^= (y << S) & 0x90000000;

	// inverse y ^= (y >> U) & D;
	y ^= (y >> U) & 0x001FFC00;
	y ^= (y >> U) & 0x000003FF;

	return y;
    }

public:
    Random(int seed=DEFAULT_SEED) {
	initialize(seed);
    }

    Random(const vector<unsigned> &prv) {
	assert(prv.size() == N);
	for (int i=0; i<N; i++) {
	    mt[i] = inverse_tempering(prv[i]);
	}
	for (int i=0; i<N; i++) {
	    twist(i);
	}
	index = 0;
    }

    unsigned operator()() {
	return extract_number();
    }
};

void MAIN() {
    mt19937 engine(100);

    vector<unsigned> prv;
    REP (i, 624) {
	prv.push_back(engine());
    }
    Random random(prv);

    assert(engine() == random());
}

Codeforces エイプリルフール 2023 解法

Codeforces April Fools Day Contest 2023 の感想
CF のエイプリルフールコンテストの特徴は

  • タイトルやサンプル入出力がヒント
  • プログラミングコンテストのローカルネタ
  • ミステリーランゲージ (マイナー言語、エソラング)

など
codeforces.com

理不尽に難しいが解けなくもない問題が多い。

A : Are You a Robot?

題:あなたはロボットですか?
答:画像の文字を出力する。YES も NO も間違い

B : Was it Rated?

題:n は rated ですか?
入力:1 ≦ n ≦ 25
答:Codeforces のコンテストのURL通し番号 1 ~ 25 をratedかどうか確認する。15, 20, 21 が unrated。
昨年はAtCoderのコンテスト名前の問題が出題されたのでコンテストの番号は予想できるが、25 がヒントでも何でもないので難しい。

C : Digits

入力:1行目はテストケース数 t。各テストケースは n 個の整数。ただし n は与えられない

input output
3
2
1
4
7
1
2
3
5
8
7
30

答:サンプル入出力を見ると入力の積になっていると気が付くことができる。要素を3個、1個、4個に分割するとうまくいくので円周率の数字毎に区切ればよいと予想できる

D : Trivial Conjecture

問:コラッツの関数 f(n) が書かれている。この関数に k 回繰り返し適用しても1が現れないような整数 n を求めよ。つまり、k項について n, f(n), f(f(n)), ... ≠ 1 となる整数 n を求めよ。
ただし k の上限は秘密
答:現在見つかっているコラッツ関数のチェーンの長さよりはるかに大きな k が来るとどうしようもない。
問題をよく読むと "整数" n を求めるというので、0以下の数を出力するだけ。

E : Not a Geometry Problem

入力:数 x, y, z が与えられる
出力:実数を出力せよ。ただし、相対誤差・もしくは絶対誤差は  10^6 以下でなければならない。つまり参加者のプログラムの出力を  a、審判の出力を  b としたとき、 \frac{|a-b|}{\max(1, |b|)} \le 10^6
答:普通の問題は誤差  10^{-6} で求めるが、本題は  10^6 である。 a = 0 のとき、誤差は必ず 1 以下なので 0 を出力すればよい。

F : Factorization

問:入力の n の最大の素因数を求めよ。
入力は2通りあり、問題文で書かれているので実質 output only

入力1 n = 4167792762229302596005813
入力2 n = 502326648535222453054166634657971818804572580255694785590270206376893052666523759828749572821869200397402455443130219791674914146276480544216264509037323019703862145022909043607926185591029614599889902115472391135622402044979347133959392884686037208893694733655782993294168167973855585231709683012084723677082273198866111120369101303677409522966567521782715484001992772768993119841291702786496058775824381444079748162416745495656333618343487208147794874337933873576016717726298883519261055062303842274145012056670644839715140659887936321934474824687778512706909988484451300384818197143498259061041

答:どちらも2つの素数の積である。一般にこんなに大きな素因数分解ができるはずもない。
入力1は素因数分解できるサイズである。

$ factor 4167792762229302596005813
4167792762229302596005813: 991999999999 4201403994187

ここで小さいほうの素数 991999999999 が見た目が特徴的な数であり、入力2も同じ傾向があると予想できる(できない)。
1 が一つと 9 だけでできる素数で試割りすると解ける。
2種の数字で、一方が1回だけ使われる素数を glitch prime と呼ぶらしい(glitch = 1文字だけエラーで置き換わったという意味?)。

G : Colour Vision

問題文はなく、画像だけ

答:単色になったサイトロゴ。カラーコードを16進数で取ると #01722B。1000番台の数を見ると CF のURLと予想できる。しかもコンテスト通し番号 1722 の問Bはcolorに関する問題で、その問題の解をそのまま提出すればよい。

ブラウザ上でカラーコードを取ると必ずしも同じ値にならないのでダウンロードしてから取得する必要がある。

H : Expected Twist

問:インタラクティブ問題。
ラクルは長さ  n の数列  a_1 ... a_n を隠し持っている。参加者は次の質問をすることができる:

  • 区間  l, r最大値は何か

ラクルは  \max( a_l, ..., a_r ) を返答する。
624 回以内の質問をして数列の最小値を求めよ。
ただし数列の要素はランダムに生成された  0 以上  2^{32} - 1 以下の数。

答:最大値のクエリで最小値を求めるのは無理そう。
題名のTwist、数列は乱数、そして 624 からメルセンヌツイスターを想起できる。メルセンヌツイスターは連続 624 項を見ると内部状態がわかり、次の値が予測できるらしい。よって、最初の 624 項を取得して数列を復元する

I : Mountain Climber

入力:各テストケースは一つの文字列が与えられる。
出力:YES NO で答えよ
入力例:たくさん
答:アルファベットの見た目に関係する問題。
アルファベット(のフォント)は上に大きい文字 Ascender と 下に大きい文字 Descender がある。Ascender を開き括弧、Descender と閉じ括弧、それ以外の文字を無視して括弧の対応が取れていれば YES

これを解くのは不可能では??

J : Unmysterious Language

問:Unmysterious language (謎ではない言語)。ジャッジは参加者の提出をみて、正当かどうか判定する。
答:CFはコンパイルエラーメッセージが見えるので適当に提出してエラーを検索するのが Mysterious Language の定跡。試しに出すと

Checker comment
wrong answer I'm sorry, but I'm not sure what you are asking. If you have any question or need any assistance, feel free to ask and I'll do my best to help you.

そういえば最近 ChatGPT が話題になっている。英語で「 accepted をください」と書けばよい。

USACO 2023 Feb. Problem 1. Hungry Cow 解法

USACO 2023 February Contest, Platinum
Problem 1. Hungry Cow

USACO 2023 の 3rd Contest の復習をする。

全体感想

問題名 略解 目標時間
Hungry Cow セグツリー 2時間
Problem Setting bit dp 20分
Watching Cowflix 木dp 3時間

解説を見てコンテスト4時間で3問満点をだすのは不可能だと思ったが、満点は複数人居る。

Hungry Cow 問題

usaco.org

牧場の牛は毎晩、エサの干し草があれば 1 単位の干し草を食べる、干し草がない場合は何もしない。
最初、牧場に干し草が届く予定はない。
U 個の予定の更新が与えられる。i 番目の予定は:

  •  d_i 日目の朝に  b_i 単位の干し草が届けられる

同じ日の予定がある場合は古い予定は消して、新しいもので上書きする。
牛が干し草を食べられる日にちの総和を出力せよ。(ただし総和は大きい可能性があるので 1000000007 で割った余りを出力せよ)

制約

 1 ≦ U ≦ 100000
 1 ≦ d_i ≦ 10 ^ {14}
 1 ≦ b_i ≦ 10 ^ 9

3
4 3
1 5
1 2

各予定を更新すると、

  • 4 + 5 + 6 = 15
  • 1 + 2 + 3 + 4 + 5 + 6 + 7 + 8 = 36
  • 1 + 2 + 4 + 5 + 6 = 18

解法

  1. クエリをノードとするようなセグツリーを作る。U個のアップデートをセグメントツリーに分解する
  2. 日にちの区間をノードとするセグツリーを作る。このセグツリーはロールバックに対応する必要がある

クエリのセグツリーをDFSしながら、日にちのセグツリーを更新・ロールバックする

クエリのセグツリー

同じ日の更新が複数あって、特に前回より減るような更新があると難しい。そのため各更新に対して上書きされるまでの期間を求めて、その範囲に対応するセグツリーのノードに加算する。
セグツリーをDFSすると、ツリーの子に移動するときにそのノードの更新を適用させ、親に戻るときには更新をロールバックする。

  • 更新の回数が logU 倍になる
  • 増える更新だけ対応すればよい
  • 更新のロールバックが必要になる (または永続化)

日にちのセグツリー

増える更新だけ対応するとき、重なる区間を併合する必要があるが、一つずつ併合するとロールバックの時に復活して計算量が悪くなるので、一回の更新を logD 時間で行いたい(D = max d + sum b)。
日にちの区間を完全に覆うならば、子ノードに再帰しないようなセグツリーを作る必要がある。
日にちの範囲が広いので座標圧縮する。干し草を食べられる連続した区間の左端は  d_i のいずれかと等しい。右端は必ずしも  d_i と一致するとは限らない。

  • リーフ:  d_i 以上  d_{i+1} 未満の区間に対応して、この区間の prefix m 日は連続して干し草を食べられる ( 0 ≦ m ≦ d_{i+1} - d_i )
  • ノード:区間の間が何日干し草を食べられるか、また、その日にちの総和

コード

#pragma GCC optimize ("O3")
#pragma GCC target ("sse4")
#pragma GCC optimize("unroll-loops")
#include<stdio.h>
#include<iostream>
#include<vector>
#include<algorithm>
#include<string>
#include<string.h>
#include<map>
#include<memory>

#ifdef LOCAL
#define eprintf(...) fprintf(stderr, __VA_ARGS__)
#else
#define NDEBUG
#define eprintf(...) do {} while (0)
#endif
#include<cassert>

using namespace std;

typedef long long LL;
typedef vector<int> VI;

#define REP(i,n) for(int i=0, i##_len=(n); i<i##_len; ++i)
#define EACH(i,c) for(__typeof((c).begin()) i=(c).begin(),i##_end=(c).end();i!=i##_end;++i)

template<class T> inline void amin(T &x, const T &y) { if (y<x) x=y; }
template<class T> inline void amax(T &x, const T &y) { if (x<y) x=y; }
#define rprintf(fmt, begin, end) do { const auto end_rp = (end); auto it_rp = (begin); for (bool sp_rp=0; it_rp!=end_rp; ++it_rp) { if (sp_rp) putchar(' '); else sp_rp = true; printf(fmt, *it_rp); } putchar('\n'); } while(0)

template<class T> void sort_unique(vector<T> &v) {
    sort(v.begin(), v.end());
    v.erase(unique(v.begin(), v.end()), v.end());
}

template<unsigned MOD_> struct ModInt {
    static constexpr unsigned MOD = MOD_;
    unsigned x;
    void undef() { x = (unsigned)-1; }
    bool isnan() const { return x == (unsigned)-1; }
    inline int geti() const { return (int)x; }
    ModInt() { x = 0; }
    ModInt(int y) { if (y<0 || (int)MOD<=y) y %= (int)MOD; if (y<0) y += MOD; x=y; }
    ModInt(unsigned y) { if (MOD<=y) x = y % MOD; else x = y; }
    ModInt(long long y) { if (y<0 || MOD<=y) y %= MOD; if (y<0) y += MOD; x=y; }
    ModInt(unsigned long long y) { if (MOD<=y) x = y % MOD; else x = y; }
    ModInt &operator+=(const ModInt y) { if ((x += y.x) >= MOD) x -= MOD; return *this; }
    ModInt &operator-=(const ModInt y) { if ((x -= y.x) & (1u<<31)) x += MOD; return *this; }
    ModInt &operator*=(const ModInt y) { x = (unsigned long long)x * y.x % MOD; return *this; }
    ModInt &operator/=(const ModInt y) { x = (unsigned long long)x * y.inv().x % MOD; return *this; }
    ModInt operator-() const { return (x ? MOD-x: 0); }

    ModInt inv() const { return pow(MOD-2); }
    ModInt pow(long long y) const {
	ModInt b = *this, r = 1;
	if (y < 0) { b = b.inv(); y = -y; }
	for (; y; y>>=1) {
	    if (y&1) r *= b;
	    b *= b;
	}
	return r;
    }

    friend ModInt operator+(ModInt x, const ModInt y) { return x += y; }
    friend ModInt operator-(ModInt x, const ModInt y) { return x -= y; }
    friend ModInt operator*(ModInt x, const ModInt y) { return x *= y; }
    friend ModInt operator/(ModInt x, const ModInt y) { return x *= y.inv(); }
    friend bool operator<(const ModInt x, const ModInt y) { return x.x < y.x; }
    friend bool operator==(const ModInt x, const ModInt y) { return x.x == y.x; }
    friend bool operator!=(const ModInt x, const ModInt y) { return x.x != y.x; }
};

constexpr LL MOD = 1000000007;

using Mint = ModInt<MOD>;
const Mint inv2 = Mint(1) / 2;

Mint f(LL start, LL len) {
    return (Mint(start) * 2 + len - 1) * len * inv2;
}


int U;
LL D[100011];
LL B[100011];

int SIZE;
vector<pair<LL, LL> > DB[1<<19];

void STORE(int left, int right, LL d, LL b) {
    left += SIZE;
    right += SIZE;
    for (; left<right; left/=2, right/=2) {
	if (left & 1) { DB[left++].emplace_back(d, b); }
	if (right & 1) { DB[--right].emplace_back(d, b); }
    }
}

struct Node {
    Node *left, *right;
    LL start;
    LL len;
    LL num;
    Mint sum;

    void clear() {
	left = right = nullptr;
	start = len = num = 0;
	sum = 0;
    }
} nodes[5000000];
int nodei;

Node* new_node() {
    assert(nodei < 5000000);
    nodes[nodei].clear();
    return nodes + nodei++;
}

Node* new_node(Node *left, Node *right) {
    Node *x = new_node();
    x->left = left;
    x->right = right;
    x->start = left->start;
    x->len = left->len + right->len;
    x->num = left->num + right->num;
    x->sum = left->sum + right->sum;
    return x;
}

vector<LL> Ds;

Node* build(int l, int r) {
    if (l == r) {
	return nullptr;
    } else if (l + 1 == r) {
	Node *x = new_node();
	x->start = Ds[l];
	x->len = Ds[r] - Ds[l];
	return x;
    } else {
	Node *x = new_node();
	x->start = Ds[l];
	x->len = Ds[r] - Ds[l];
	x->left = build(l, (l+r)/2);
	x->right = build((l+r)/2, r);
	assert(x->start == x->left->start);
	assert(x->len = x->left->len + x->right->len);
	return x;
    }
}

pair<Node*, LL> add_rec(Node *tree, LL start, LL num) {
    if (!tree || num == 0 || tree->start + tree->len <= start || tree->len == tree->num) {
	return make_pair(tree, num);
    } else if (start <= tree->start && tree->len <= num + tree->num) {
	Node *x = new_node();
	x->start = tree->start;
	x->len = tree->len;
	x->num = x->len;
	x->sum = f(tree->start, x->num);
	return make_pair(x, num + tree->num - x->num);
    } else if (!tree->left && !tree->right) {
	assert(tree->len >= tree->num + num);
	Node *x = new_node();
	x->start = tree->start;
	x->len = tree->len;
	x->num = tree->num + num;
	x->sum = f(tree->start, x->num);
	return make_pair(x, 0LL);
    } else {
	auto L = add_rec(tree->left, start, num);
	auto R = add_rec(tree->right, start, L.second);
	Node *x = new_node(L.first, R.first);
	return make_pair(x, R.second);
    }
}


Node* add(Node *tree, LL d, LL b) {
    auto p = add_rec(tree, d, b);
    assert(p.second == 0);
    return p.first;
}

void DFS(Node *tree, int k) {
    int backup = nodei;
    EACH (e, DB[k]) {
	tree = add(tree, e->first, e->second);
    }
    if (SIZE <= k) {
	if (k-SIZE < U) {
	    printf("%d\n", tree->sum.geti());
	}
    } else {
	DFS(tree, k*2);
	DFS(tree, k*2+1);
    }
    nodei = backup;
}

void MAIN() {
    scanf("%d", &U);

    REP (i, U) {
	scanf("%lld%lld", D+i, B+i);
    }

    Ds.reserve(U+1);
    Ds.assign(D, D+U);
    Ds.push_back(1LL<<50);
    sort_unique(Ds);


    SIZE = 1;
    while (SIZE < U) SIZE += SIZE;

    vector<pair<LL, int> > V;

    REP (i, U) {
	V.emplace_back(D[i], i);
    }
    sort(V.begin(), V.end());

    for (int i=0, j=0; i<U; i=j) {
	while (j < U && V[i].first == V[j].first) {
	    j++;
	}

	int right = U;
	for (int k=j-1; k>=i; k--) {
	    int left = V[k].second;
	    STORE(left, right, D[left], B[left]);
	    right = left;
	}
    }

    Node *tree = build(0, Ds.size()-1u);
    DFS(tree, 1);
}

int main() {
    int TC = 1;
//    scanf("%d", &TC);
    REP (tc, TC) MAIN();
    return 0;
}

任意精度平方根と素数テスト (SRM821 CrossPrimesOut)

独特な要素が多く、結構印象深い問題。出題されたのは1年前。
出題はSRM 821 Div1の3問目。

問題

community.topcoder.com
正整数Xが1つ与えられる。その平方根を取って小数点を無視して、無限の長さの数字の文字列が得られる。

  • 例 : sqrt(47) → 685565460040104...

次に、無限に存在する素数に対して小さい順に、以下の手続きを行う。

  • その素数が文字列の中で最初に現れる位置を見つける
  • その文字を空白で置き換える。(無限に長い文字列の中にその素数が存在しない場合は何もしない)

この無限の手続きを終えると空白区切りのトークン列が得られる。各トークンは数字の文字列になっている。

平方数ではない正整数Xと整数kが与えられるので、k番目のトークンを求めよ。

制約

1 ≦ X ≦ 10000
X は平方数ではない
0 ≦ k ≦ 30

X = 47
k = 1

答 5654600

6855654600401044124935871449084848960460643461001326275485108185678...
68_5654600___04___49_58_1___08484__604606__4__00__2627_____08185__8...

(0オリジンで)k 番目のトークンは 5654600。
同じ素数は最初の一回しか消さない。一度空白に置き換えたい値は別の素数にはならない。

解法

  1. 400桁の平方根を求めるて、十進数で文字列とする。
  2. 1000000 以下の素数について、文字列に存在するなら空白に置き換える。
  3. トークンが7文字の素数を含むなら空白に置き換える。以下8文字、9文字、と順に調べる。

公式解説
docs.google.com

問題の独自性

  1. 任意精度小数を使うというのが珍しい。ただしsqrtを求めるだけなのでニュートン法や二分法だけできればいい。
  2. この問題の制約の範囲では400桁でいいが、実行してみるまで何桁必要か予想できない。
  3. 無限の文字列から見つからないときは何もしないが、見つかるか見つからないか不明な時はどうするのか?答えには影響しないので無視すればいいが不安定。
  4. もし100文字以上のトークンが残ったらこの中に素数がないことを判定しなければならない。実装中はトークンの長さが予想できない。

実装の注意

任意精度小数のある言語を選択するのが良い。
結構大きめの素数を判定する必要があるので多倍長を用意する。(例えば7936145404496840003051、2^64 より大きい素数が出てくる)
文字列から素数を探すので誤って "0000002" を7桁の素数と判定しないように。

コード (python3)

from decimal import *

getcontext().prec = 500

MAX = 1000000

prime = []
tbl = [True] * MAX
tbl[0] = False
tbl[1] = False
for i in range(2, MAX):
    if tbl[i]:
        prime.append(i)
        for j in range(i*i, MAX, i):
            tbl[j] = False


def isSprp(a, s, d, n):
    a %= n
    if a == 0:
        return True

    x = pow(a, d, n)

    if x == 1:
        return True

    for i in range(s):
        if x == n - 1:
            return True
        x = x * x % n

    return False


def isPrime(n):
    if n < MAX:
        return tbl[n]

    if n % 2 == 0 or n % 3 == 0 or n % 5 == 0 or n % 7 == 0:
        return False

    d = n - 1
    s = 0
    while d % 2 == 0:
        s += 1
        d //= 2
    
    bases = [2, 325, 9375, 28178, 450775, 9780504, 1795265022]
    for b in bases:
        if not isSprp(b, s, d, n):
            return False

    return True


class CrossPrimesOut:
    def findTerm(self, X, N):
        s = str(Decimal(X).sqrt()).replace('.', '')

        for p in prime:
            p = str(p)
            s = s.replace(p, ' ', 1)

        s = s.split()

        if not s:
            return []

        maxLen = max([len(e) for e in s])

        for curLen in range(7, maxLen+1):
            while True:
                cand = []
                for i in range(len(s)):
                    if len(s[i]) < curLen:
                        continue
                    for j in range(len(s[i]) - curLen + 1):
                        if s[i][j] == '0':
                            continue
                        num = int(s[i][j:j+curLen])

                        if isPrime(num):
                            cand.append((num, i, j))

                if not cand:
                    break

                num, i, j = min(cand)

                if len(s[i]) == curLen:
                    s.pop(i)
                elif j == 0:
                    s[i] = s[i][curLen:]
                elif len(s[i]) == j + curLen:
                    s[i] = s[i][:j]
                else:
                    left = s[i][:j]
                    right = s[i][j+curLen:]
                    s[i] = right
                    s.insert(i, left)

        return int(s[N]) % 1000000007

コミュニケーションの問題(COCI 2022 COI Mensza)

出典

COCI 2022 COI 問2 Mensza
https://codeforces.com/blog/entry/103280
IOI系 ではコミュニケーションと呼ばれる独特な形式の問題が出題されることがある。

問題

4人の人物がいる

  • Mr. Malnar マルナー先生
  • Alojzije アロイジエ
  • Benjamin ベンジャミン
  • Cecilija セシリア

次の試験を行う。

  1. マルナー先生が異なる正整数 A, B を決める(1 ≦ A, B ≦ 500000)
  2. マルナー先生Aアロイジエに伝える。A を見たアロイジエは長さ 20 以下の数列 X を作ってマルナー先生に渡す
  3. マルナー先生Bベンジャミンに伝える。B を見たベンジャミンは長さ 20 以下の数列 Y を作ってマルナー先生に渡す
  4. マルナー先生は数列 X, Y を連結して数の出現回数を数える。出現回数の列を昇順に並べた数列 Z を作る
  5. マルナー先生セシリアZ を渡す。Z を見たセシリアA, B のどちらが大きいと思うか宣言する。正しい宣言をすれば試験は成功である

例えば X = [10, 30, 50], Y = [25, 30, 30, 46, 50] の場合、出現回数は { 10 → 1, 25 → 1, 30 → 3, 46 → 1, 50 → 2}、よって Z = [1, 1, 1, 2, 3] がセシリアに渡される。
アロイジエ、ベンジャミン、セシリアの3人は協力してこの試験に成功したい。あなたは試験が確実に成功するような3人分のプログラムを作りなさい。

補足

数列は長さと順番に並んだ要素で構成され、文字列にして送信される。例えば数列 [10, 30, 50] はスペース区切りの1行の文字列 "3 10 30 50" としてやり取りされる。
数列の要素は 1 以上 10^9 以下でなければならない。
アロイジエ、ベンジャミン、セシリア、は上記の方法以外で情報をやり取りできない。
実際はマルナー先生がどんな順番でアロイジエ、ベンジャミン、セシリアと会話するかわからない。アロイジエ、ベンジャミンにいろいろな数を聞いた後でセシリアに聞いて試験するかもしれない。
問題の制約は満点の場合を基準にしている。

解答

アロイジエは A だけを知っている。ベンジャミンは B だけを知っている。セシリアは A も B も知らない、復元もできないがどちらが大きいかだけは知っている。という数学パズル。


























考察1

まず、XとYを連結した数列の長さは40以下であるのでセシリアの受け取る数列は 40 以下の分割数のパターンしかない。これは A と B の組み合わせのパターンより少ないのでセシリアは A, B を同時に復元することはできない。
アロイジエとベンジャミンは A,B の一方の値しか知らないのでどちらの値が大きいかを数列 X, Y に書き込むことはできない。
アロイジエとベンジャミンが完全に同じアルゴリズムであると、A と B が入れ替わった時に Z が全く同じ数列になり、セシリアは区別できない。そのためアロイジエとベンジャミンは異なるアルゴリズムのはずである。

考察2

0 ≦ A, B ≦ 1 の場合を考える。

  • アロイジエ:A=0ならば空列、1ならば [1] を作る
  • ベンジャミン:B=0ならば [1]、1ならば空列を作る

すると (A, B) 4通りについて、

  • (0, 0) : Z = [1]
  • (1, 0) : Z = [2]
  • (0, 1) : Z = 空列
  • (1, 1) : Z = [1]

となり、(1, 0) と (0, 1) は区別することができ、必要のない (0, 0) と (1, 1) が同じ数列になって効率よくまとめられている。

解答

A, B を2進数で表記したときに初めて異なる上位ビットを見つけることができればそれが大小関係になる

A = 0b00101011
B = 0b00100110
          ^

このビットより下位ビットは 0 に上書きし、B の 0 を 1 にフリップすれば A と B が一致する。

A → 0b00101000
B → 0b00101000
          ^

この方法で「上位ビットから見て初めて異なるビットが A は 1、B は 0」 の場合のみ値が等しくなる、つまり、Z に 2 が含まれる。A, B は20ビットで表現できるため、数列の長さも 20 以下で条件を満たす。

よって、

  1. アロイジエ:A を2進数で表して、あるビットが 1 ならばそれより下位のビットを0にしてXにpush
  2. ベンジャミン:B を2進数で表して、あるビットが 0 ならばそのビットを1にして、それより下位のビットを0にしてYにpush
  3. セシリア:数列Zに 2 が含まれていれば A>B、そうでなければ A < B

IOI ではこの3つのプログラムを別のファイルに書いて提出するはずだが、COI では一つのファイルに書くことになっていた。