はじめに

整数 $N,S$ が与えられたときに、$N$ 個の乱数からなり、乱数の総和が $S$ に一致するような整数列を生成せよ。

この記事は、上記の問題を解く方法について考察して実際にプログラムを書いてみた記事です。何に役立つかはわかりませんが、単純に興味が湧いたので考えてみることにしました。

なお、この記事では C++ を使ってコードを書いていますが、難しいテクニックは登場しないので他の言語でも簡単に書き直すことができると思います。

0. C++ で乱数を生成する

本題に入る前に、まずは総和の条件を無視して $L$ 以上 $R$ 以下のランダムな整数を生成するコードを書いてみます。

コード

C++

#include <iostream>
#include <random>
#include <vector>
using namespace std;

int main(void)
{
    long long N, L, R;
    cin >> N >> L >> R;

    //ランダムなシード値を使って乱数生成器を初期化
    random_device rd;
    mt19937_64 mt(rd());

    //生成する乱数の最小値L,最大値Rを指定して初期化
    uniform_int_distribution<long long> uid(L, R);

    vector<long long> array(N);

    for(long long i = 0; i < N; ++i)
    {
        //L以上R以下のランダムな整数を等確率で生成
        array[i] = uid(mt);
    }

    //乱数列を出力
    for(long long i = 0; i < N; ++i)
    {
        cout << array[i] << " ";
    }
    cout << endl;

    return 0;
}

コードの説明

C++ で乱数を生成するには、以下の3つのアイテムを使います¹。いずれも std 内にあり、標準ライブラリ random をインクルードすることで使えます。

クラス名	説明
random_device	予測不能な乱数を生成する。乱数生成器の初期化に必要なシード値を生成するために使う²。
mt19937_64	メルセンヌ・ツイスタを使った疑似乱数生成器。 0〜(64bit 整数型の最大値)の範囲の乱数を生成する。 (32bit 版は mt19937)
uniform_int_distribution	分布生成器(乱数の範囲や分布を調整するもの)の一種。乱数生成器を利用して $L$ 以上 $R$ 以下の整数を等確率で生成する。

まず、random_device rd; と mt19937_64 mt(rd()); で乱数生成器を初期化して、続けて uniform_... uid(L, R); で分布生成器を初期化します。

そして、引数に乱数生成器 mt を指定して uid(mt) を呼び出すことで、$L$ 以上 $R$ 以下のランダムな整数を生成することができます。(uid(mt()) ではないので注意)

入出力例

このプログラムを実行して、試しに $0$ 以上 $100000$ 以下のランダムな整数を $100$ 個生成してみるとこのようになります。

入力例

100 0 100000

出力例

85797 27864 34531 54965 79326 5713 25896 25303 65762 83558 63442 90407 16064 66635 26818 63312 49289 64199 88565 40295 56323 78423 62543 69232 41635 11620 37692 9330 24405 30023 81937 78091 19401 60693 19451 31555 445 97142 21018 1571 74471 81328 47457 56870 33220 48732 55338 84874 42172 85642 21549 91067 14080 60676 92417 72912 56152 90623 34028 84095 7119 61516 75761 31825 36745 20911 73982 47730 46503 76882 57969 8692 73312 50693 28414 24444 121 40532 37494 89090 50845 19376 94481 18526 52143 98558 85295 73337 43206 87846 80052 73900 65891 75649 23696 75494 22712 3735 57814 23155

正しく乱数列が生成できていますね。

個々の乱数を一直線にプロットしたものと、乱数列を箱ひげ図にまとめたものが下の図です。個数が少ないので分布に多少の偏りはありますが、おおよそ均等に分布していて平均値(✕印)も真ん中の $50000$ に近いことが分かります。

1. 【案①】長さ S の区間を N 個の区間に分割する

ではいよいよ、総和が $S$ となる $N$ 個の乱数を生成する方法について考えていきます。

まず思いつく一番簡単で確実な方法は、$0$ 以上 $S$ 以下の区間を $N-1$ 個の乱数で区切って $N$ 個の区間にし、それぞれの区間の幅を求めて数列にする 方法です。

例として $N = 10, S = 100$ の場合を示したのが上の図です。

まずは $0$ 以上 $100$ 以下の乱数を $9$ 個生成して、それらに $0$ と $100$ を加えた合計 $11$ 個の数値を小さい順に並べ替えます。上の図では、並び替えた後の数値は順に $0,6,17,25,29,41,56,66,78,91,100$ となっています。

そして、小さい方から順に隣り合う数値の差を求めると $6-0=6,\ 17-6=11,\ \ldots,\ 100-91=9$ となります。これらの和は当然 $100$ ですから、条件を満たす乱数列を得ることができました。

コード

上記の処理を一般化して実装します。

$0$ 以上 $S$ 以下の $N - 1$ 個の乱数を生成して配列 base に追加
さらに $0$ と $S$ を追加してから base を昇順にソート
配列 array の $i$ 番目に base の $i$ 番目と $i + 1$ 番目の要素の差を代入する

C++

#include <iostream>
#include <random>
#include <vector>
#include <algorithm>
using namespace std;

int main(void)
{
    long long N, S;
    cin >> N >> S;

    random_device rd;
    mt19937_64 mt(rd());

    uniform_int_distribution<long long> uid(0, S);

    vector<long long> base(N - 1), array(N);

    for(long long i = 0; i < N - 1; ++i)
    {
        //0以上S以下のランダムな整数を等確率で生成
        base[i] = uid(mt);
    }
    base.push_back(0);
    base.push_back(S);

    //乱数列をソート
    sort(base.begin(), base.end());

    //隣り合う要素の差をとる
    for(long long i = 0; i < N; ++i)
    {
        array[i] = base[i + 1] - base[i];
    }

    //数列とその和を出力
    long long sum = 0;
    for(long long i = 0; i < N; ++i)
    {
        cout << array[i] << " ";
        sum += array[i];
    }
    cout << endl;
    cout << sum << endl;

    return 0;
}

入出力例

先程の図で使用した、総和が $100$ となる $10$ 個の乱数を生成した例も載せておきます。

入力例

10 100

出力例

6 11 8 4 12 15 10 12 13 9 
100

2. 案① の問題点

案① の「長さ $S$ の区間を $N$ 個の区間に分割する」という方法は、簡単ではありますが色々と不都合な点があります。

乱数の範囲を設定できない

案① の方法はその原理上、負の数を含む乱数列を生成することができません。また、「$L$ 以上 $R$ 以下の範囲の乱数列を得る」といったような、乱数の範囲を設定することもできません。

乱数の偏りが大きい

上の図は、1. のコードで $N=100,S=100000$ のときの配列 base および array をそれぞれ箱ひげ図にしたものです。

base は生成された乱数をそのまま要素としているので、もちろんおおよそ均等に分布しています。

一方、base をもとにして生成された array は 乱数の分布が大きく偏っており、乱数の約75%が $0$ 以上 $1500$ 以下の範囲に集中している ことがわかります。乱数列はおおよそ均等な分布であってほしいので、案① で得られた array はあまり使えるものではありません。

ということで、これらの短所を改善した別の案を考える必要がありそうです。

3. 【案②】平均値を基準に乱数の生成範囲を設定する

2. を踏まえて、改めて総和が $S$ となる $L$ 以上 $R$ 以下の乱数 $N$ 個を生成する方法を考えてみます。

ここで使えるのが、0. でも書いたように、$L$ 以上 $R$ 以下の範囲で一様に乱数を生成すると、その平均値は $L$ と $R$ の中間、すなわち $\displaystyle\frac{L+R}{2}$ に近い値になるということです。

そこで、中間値 $\displaystyle\frac{l+r}{2}$ が平均値 $\displaystyle\frac{S}{N}$ に一致するように乱数の生成範囲 $l, r$ を定めることで、乱数分布の均等さを保ちながら、生成した乱数の総和を $S$ に近い値にすることができます。後は、総和が $S$ と一致するように数列全体から均等に誤差を引いて微調整すれば、分布を崩さずに条件を満たす乱数列を得ることができます。

乱数の生成範囲の決め方

乱数の生成範囲 $l, r$ は以下の条件を満たすようにします。

中間値 $\displaystyle\frac{l+r}{2}$ が平均値 $\displaystyle\frac{S}{N}$ に一致する
$L\le l\le r \le R$ を満たす ($L$ 以上 $R$ 以下の範囲をはみ出さない)
上2つの条件を満たすもののうち、生成範囲の幅 $r-l$ が最大となるものを選ぶ

これらの条件を考えると、条件を満たす $l, r$ は下の図のように $\displaystyle\frac{S}{N}$ の位置によって $2$ 通りに分かれます。ここで、生成範囲の幅の半分の値 $\displaystyle\frac{r-l}{2}$ について考えると、これは$\displaystyle \frac{S}{N}-L,R-\frac{S}{N}$ のうち小さい方の値に等しいことが分かります。

したがって、$l,r$ の値はそれぞれ次のように求められます。

$\displaystyle l=\frac{S}{N}-\min\biggl(\frac{S}{N}-L,R-\frac{S}{N}\biggr)=\max\biggl(L,2\times\frac{S}{N}-R\biggr)$
$\displaystyle r=\frac{S}{N}+\min\biggl(\frac{S}{N}-L,R-\frac{S}{N}\biggr)=\min\biggl(2\times\frac{S}{N}-L,R\biggr)$

誤差の修正

乱数の生成範囲 $l,r$ を決めて実際に $N$ 個の乱数を生成した後は、総和が $S$ になるように乱数を修正します。

生成された $N$ 個の乱数の和を $S^\prime$ とすると、平均値の誤差は $\displaystyle\frac{S^\prime -S}{N}$ なので、それぞれの乱数から $\displaystyle\frac{S^\prime -S}{N}$ だけ引くことで、乱数の分布を変えずに総和を $S$ にすることができます。

ただし、実際にはそれぞれの乱数は整数である必要があります。そこで、まずは各乱数から$\displaystyle\biggl\lfloor\frac{S^\prime -S}{N}\biggr\rfloor$($\displaystyle\frac{S^\prime -S}{N}$ の整数部分) だけ引いて、総和の誤差を $N$ 未満にします。その後、残りの誤差を $s$ として、$N$ 個の乱数から $s$ 個の乱数を重複なくランダムに選んで $1$ ずつ引くことで、残りの誤差を修正します($S^\prime\gt S$ の場合)。

なお、「$N$ 個の乱数から $s$ 個の乱数を重複なくランダムに選ぶ」方法として、Fisher-Yates shuffle というアルゴリズムを使います。

コード

上記を参考にしてコードを実装します。

なお、与えられた入力において平均値 $\displaystyle\frac{S}{N}$ が $L$ 以上 $R$ 以下の範囲をはみ出している場合は、条件を満たす数列は存在しないのであらかじめ除外します。

C++

#include <iostream>
#include <random>
#include <vector>
using namespace std;

int main(void)
{
    long long N, S, L, R;
    cin >> N >> S >> L >> R;

    random_device rd;
    mt19937_64 mt(rd());

    if(S / N < L || R < S / N)
    {
        //平均が範囲外の場合は生成不可能
        cout << "The average is out of range." << endl;
        return 0;
    }

    //乱数の生成範囲を決める
    long long l = max(L, 2 * S / N - R), r = min(2 * S / N - L, R);

    uniform_int_distribution<long long> uid(l, r);

    vector<long long> array(N);

    long long sum = 0;
    for(long long i = 0; i < N; ++i)
    {
        array[i] = uid(mt);
        sum += array[i];
    }

    //誤差を修正する
    long long err = sum - S;
    for(long long i = 0; i < N; ++i)
    {
        array[i] -= err / N;
    }
    err -= err / N * N;

    //微調整(N個の乱数から重複のないように選ぶ(Fisher-Yates shuffle))
    vector<long long> fisher(N);
    for(long long i = 0; i < N; ++i)
    {
        fisher[i] = i;
    }

    long long idx;
    for(long long upper = N - 1; err > 0; --err, --upper)
    {
        uniform_int_distribution<long long> choose(0, upper);
        idx = choose(mt);
        array[fisher[idx]] -= 1;
        fisher[idx] = fisher[upper];
    }
    for(long long upper = N - 1; err < 0; ++err, --upper)
    {
        uniform_int_distribution<long long> choose(0, upper);
        idx = choose(mt);
        array[fisher[idx]] += 1;
        fisher[idx] = fisher[upper];
    }

    //数列とその和を出力
    sum = 0;
    for(long long i = 0; i < N; ++i)
    {
        cout << array[i] << " ";
        sum += array[i];
    }
    cout << endl;
    cout << sum << endl;

    return 0;
}

入出力例

このプログラムを実行して、総和が $10000$ となる $100$ 以上 $1000$ 以下の乱数 $20$ 個を生成してみるとこのようになります。

入力例1

20 10000 100 1000

出力例1

524 671 468 634 472 479 324 182 406 591 756 114 857 422 517 501 221 579 830 452 
10000

ちゃんと条件に合う乱数列が生成できています！合っているのか気になる人は、実際に電卓などで総和を計算してみてください。

負の数を含む乱数列も生成してみます。総和が $0$ となる $-1000$ 以上 $1000$ 以下の乱数 $100$ 個を生成してみるとこのようになります。

入力例2

100 0 -1000 1000

出力例2

-156 886 595 -758 440 308 246 -940 548 -521 769 255 -855 588 574 792 -865 -115 -169 -433 -176 -201 -782 785 -500 858 529 89 170 676 -196 100 -927 891 916 989 90 -651 -821 -96 -828 731 -453 -853 -132 597 -442 840 9 -708 285 -434 195 -473 -657 -542 601 -212 798 -288 -660 938 -198 15 -928 718 420 -903 -813 -634 -856 428 -298 570 -266 314 -286 109 -314 812 -127 -52 -907 612 -8 -90 173 989 -581 -449 929 -167 614 -946 852 -655 796 -527 464 -54 
0

4. 案② の問題点

3. では案② を用いて適切な乱数列を得ることができましたが、この方法にもまだ問題があります。

それは、乱数を生成した後、誤差を修正するときに一部の乱数が $L$ 以上 $R$ 以下の範囲をはみ出してしまうことがある という問題です。

例えば 3. の入力例1では、以下のように数列に $100$ 未満の数値が含まれてしまうことがあります。

入力例1

20 10000 100 1000

出力例1(失敗例)

769 740 771 552 730 91 90 357 625 707 322 861 220 655 617 114 813 314 529 123
10000

また、3. の入力例2のように乱数の生成範囲 $l,r$ がそれぞれ $L,R$ に一致する場合では、さらに一部の乱数が範囲をはみ出てしまう確率が高くなります。 (実は 3. の出力例2は生成に $20$ 回ほど失敗した後にようやく成功したものです…)

入力例2と同じ入力を用いて乱数列を生成する試行を $10000$ 回おこなうコードを書いて実験したところ、$1$ 回の試行で生成に成功する確率は約26.8%となりました。思ったより低いです。

もちろん成功するまで乱数列の生成を繰り返しおこなえばよいのですが、やはり $1,2$ 回程度で生成に成功するのがベストです。したがって、生成に成功する確率を上げる方法を考える必要があります。

この問題の原因は、乱数生成の幅を $L,R$ の範囲ギリギリまで大きくしている ことです。そのため、誤差を修正する段階で $L$ または $R$ に近い値の乱数が範囲をはみ出てしまいます。

そこで、誤差を修正するときのズレも考慮して、乱数生成の幅をいくらか小さくして隙間を残す必要があります。次の章では、適切な隙間の幅を考えるために、誤差がどの位になるかを調べていきます。

5. 誤差について調べる

案② を改善するために、乱数の平均値 $\displaystyle \frac{S^\prime}{N}$ が目標の平均値 $\displaystyle\frac{S}{N}$ からどの位離れているか、すなわち誤差の絶対値 $\displaystyle\biggl|\frac{S^\prime-S}{N}\biggr|$ について調べます。

特殊なケースとして、総和が $S=0$ でありさらに $L=-R$ であるときは、乱数の生成範囲 $l, r$ は必ず $L,R$ にそれぞれ一致します。そこで、この条件において乱数の個数 $N$ や生成幅 $\mathrm{range}=R-L$ の条件をいくつか用意して、誤差の絶対値 $\displaystyle\biggl|\frac{S^\prime}{N}\biggr|$ を調べてみます。

コード

今回は以下の $4\times 4=16$ 通りを調べます。

乱数の生成幅: $\mathrm{range}=100,1000,10000,100000$
乱数の個数: $N=10,100,1000,10000$

C++

#include <iostream>
#include <random>
#include <vector>
using namespace std;

int main(void)
{
    long long Q; //各条件における試行回数
    cin >> Q;

    //個数Nと生成幅RANGE=R-Lの条件を変えて調べる
    //総和S=0
    vector<long long> RANGE = {100, 1000, 10000, 100000};
    vector<long long> N = {10, 100, 1000, 10000};

    random_device rd;
    mt19937_64 mt(rd());

    vector<vector<double>> err(16, vector<double>(Q, 0.0));

    for(long long i = 0; i < 4; ++i)
    {
        long long range = RANGE[i];

        for(long long j = 0; j < 4; ++j)
        {
            long long n = N[j];

            for(long long q = 0; q < Q; ++q)
            {
                uniform_int_distribution<long long> uid(-range / 2, range / 2);

                vector<long long> array(n);

                long long sum = 0;
                for(long long i = 0; i < n; ++i)
                {
                    array[i] = uid(mt);
                    sum += array[i];
                }

                //誤差の絶対値を記録
                err[i * 4 + j][q] = (double)abs(sum) / n;
            }
        }
    }

    //誤差を出力
    cout.setf(ios::fixed);
    cout.precision(4);
    for(long long i = 0; i < 4; ++i)
    {
        long long range = RANGE[i];

        cout << "range: " << range << endl;
        for(long long j = 0; j < 4; ++j)
        {
            long long n = N[j];

            cout << " n: " << n << endl;
            cout << "  ";
            for(long long q = 0; q < Q; ++q)
            {
                cout << err[i * 4 + j][q] << " ";
            }
            cout << endl;
        }
        cout << endl;
    }
    return 0;
}

入出力例

各条件における試行回数を $Q=10$ 回としたときの出力例がこちらです。

入力例

出力例

range: 100
 n: 10
  0.4000 1.3000 11.6000 6.4000 15.0000 2.8000 7.5000 1.7000 0.3000 0.3000 
 n: 100
  1.8400 2.1500 1.5900 4.8900 1.3600 1.4500 1.8900 2.0400 0.4300 3.7300 
 n: 1000
  0.8610 0.4400 0.8020 0.8450 0.3520 0.3610 0.7140 1.6350 0.1020 0.8060 
 n: 10000
  0.3061 0.2815 0.0030 0.2113 0.2600 0.0061 0.0547 0.1594 0.1108 0.3124 

range: 1000
 n: 10
  78.2000 140.6000 54.9000 47.2000 130.7000 28.1000 47.0000 22.6000 136.6000 113.1000 
 n: 100
  28.1800 4.8700 39.1400 14.0200 5.9000 19.2800 58.1600 5.0900 30.3200 1.3400 
 n: 1000
  0.0980 5.2750 6.5690 6.2850 4.2960 11.2070 4.3890 20.9880 9.2260 2.6580 
 n: 10000
  0.4190 2.3005 0.9116 3.3104 1.5774 7.3781 0.7403 6.1676 2.1706 3.5965 

range: 10000
 n: 10
  1293.4000 632.2000 272.5000 570.0000 819.4000 107.1000 148.8000 45.4000 855.1000 994.3000 
 n: 100
  549.2500 600.8600 194.2500 9.7300 112.9600 122.5800 100.8200 287.5600 208.3800 7.7700 
 n: 1000
  123.6380 4.8100 76.2150 99.3550 188.4390 3.1350 134.4140 120.5430 125.9540 10.8060 
 n: 10000
  35.2774 12.4515 16.7975 2.1413 35.6339 62.2156 31.2161 16.8822 26.3665 45.1429 

range: 100000
 n: 10
  5286.5000 12410.8000 16814.1000 11645.7000 1190.9000 779.2000 5653.8000 6532.5000 5583.4000 7707.0000 
 n: 100
  3149.8700 2422.1400 123.1900 747.1600 2341.7600 924.1400 2895.2700 3931.6300 103.4500 733.2800 
 n: 1000
  78.7000 20.2220 168.1920 1750.2650 604.6670 1396.8930 1094.0760 782.4720 188.2440 1304.4640 
 n: 10000
  57.4669 75.0663 324.3013 111.0573 301.3292 337.2189 145.8594 135.4358 86.5960 419.6086

上記の入出力例だけだとよくわからないので、試行回数を $Q=1000$ 回に増やして、各条件における誤差の分布を箱ひげ図にまとめてみます。

乱数の生成幅を $4$ 種類試してみましたが、同じ個数で比較すると生成幅に対する誤差の割合はどれもあまり変わらないようです。

また、乱数の個数 $N$ が大きくなるにつれて誤差の範囲は小さくなっており、$N$ が $10$ 倍になるごとにその分布幅は約 $\displaystyle\frac{1}{3}$ 倍になっています。

生成幅に対する誤差の割合を求める

上の箱ひげ図によると、$N=100$ のときに誤差のほとんどが生成幅の10%以内に収まっていますね。これを目安にした上で、生成幅 $\mathrm{range}$ に対する誤差の分布幅の割合 $E$ を、$N$ を $10$ 倍すると $E$ が $\displaystyle\frac{1}{3}$ 倍になるような数式で表してみます³。

\begin{alignat}{3}
E_{0.99}&=\frac{1}{10}\times\biggl(\frac{1}{3}\biggr)^{\log_{10}{N}-2}
 &&=\frac{9}{10}\times\biggl(\frac{1}{3}\biggr)^{\log_{10}{N}} \\
 &=\frac{9}{10}\times\biggl(\frac{1}{3}\biggr)^{\frac{\log_{3}{N}}{\log_{3}{10}}}
 &&=\frac{9}{10}\times\biggl(\frac{1}{3}\biggr)^{\log_{3}{N}\times\log_{10}{3}} \\
 &=\frac{9}{10}\times\biggl(\frac{1}{N}\biggr)^{\log_{10}{3}}
\end{alignat}

さらに $\log_{10}{3}=0.4771\ldots\approx 0.5$ なので、この近似を用いると $\displaystyle E_{0.99}=\frac{9}{10}\times\frac{1}{\sqrt{N}}\ \cdots(1)$ と表すことができます。

一旦 $E$ を求めてみましたが、よく考えると全ての試行において誤差の割合が $E$ 以下になっている必要はなく、試行の一部では $E$ 以下を満たさずとも、その大部分において $E$ 以下になっていれば良さそうです。 (生成に成功する確率が上がれば、もし生成に失敗しても再度生成を試せば良いです)

箱ひげ図によると、$N=10$ のときに誤差の約75%が生成幅の10%以内に収まっています。そこで、生成される誤差の約75%をカバーするような $E$ を $(1)$ と同様に求めてみます。

\begin{alignat}{3}
E_{0.75}&=\frac{1}{10}\times\biggl(\frac{1}{3}\biggr)^{\log_{10}{N}-1}
 &&=\frac{3}{10}\times\biggl(\frac{1}{N}\biggr)^{\log_{10}{3}} \\
 &\approx \frac{3}{10}\times\frac{1}{\sqrt{N}}\ \cdots(2)
\end{alignat}

$E_{0.75}$ も同様に数式で表すことができました。

ところで、$(1),(2)$ のどちらの場合も誤差の割合は $\displaystyle E=(定数)\times\frac{1}{\sqrt{N}}$ の形で表されています。つまり、定数部分をパラメータとみなして値を調整すれば、生成される誤差の何%を $E$ がカバーするかを決められるということです。

6. 案② を改良する

5. では、生成幅 $\mathrm{range}$ に対する誤差のおおまかな割合 $E$ を求めることができました。これを用いて案② を改良していきます。

乱数の生成範囲の決め方(改良版)

まずは案② と同じ方法で乱数の仮の生成範囲 $l^\prime,r^\prime$ を求めて、その生成幅を $\mathrm{range}=r^\prime-l^\prime$ とします。

そして、5. で求めた $E$ を用いて誤差の分布幅 $E\times\mathrm{range}$ を求めてから、真の生成範囲 $l,r$ を $l=l^\prime+E\times\mathrm{range},r=r^\prime-E\times\mathrm{range}$ で定めます。

なお、この方法では $E\times\mathrm{range}$ は真の生成幅ではなく仮の生成幅に対する誤差の分布幅なので、厳密には正しい方法ではありません。しかし、そもそも $E$ の求め方の時点で大雑把でしたし、幅が小さくなればその差は無視できるので問題なしとします。(これ以上式を複雑にしたくないというのもありますが…)

$E\times\mathrm{range}$ の計算についてですが、$\displaystyle E=(定数)\times\frac{1}{\sqrt{N}}$ を計算して $E\times\mathrm{range}$ の値を求めるよりも $\displaystyle\frac{1}{E}=(定数)\times\sqrt{N}$ を計算して $\displaystyle\mathrm{range}\div\frac{1}{E}$ の値を求めるほうが精度が良さそうなので、こちらを使うことにします。

また、5. の $(2)$ より $\displaystyle\frac{1}{E_{0.75}}=\frac{10}{3}\times\sqrt{N}$ なので、比例定数は $\displaystyle k=\frac{10}{3}\approx 3.33$ をデフォルトの値として使います。($k$ の値を大きくするほど $E$ は小さくなり乱数の生成幅 $r-l$ は大きくなりますが、生成に失敗する確率は増加します)

生成に失敗したかを判定する

乱数の生成範囲を調整することで生成に失敗する確率は減ったものの、まだ生成に失敗する可能性はあります。そこで今回は、生成に失敗したかを判定して、もし失敗したならば生成をやり直す処理も加えることにします。

具体的には、誤差を修正する処理のうち以下の2箇所のタイミングで、乱数が $L$ 以上 $R$ 以下の範囲をはみ出るか判定します。

各乱数から $\displaystyle\biggl\lfloor\frac{S^\prime -S}{N}\biggr\rfloor$ を引くと一部の乱数が範囲からはみ出るとき
$N$ 個の乱数から範囲ギリギリでない($L,R$ に等しくない) $s$ 個の乱数を重複なく選ぶことができないとき(起こるのは極めて稀です)

また、何回生成しても生成に失敗してしまう場合もあるので、連続で $10$ 回生成に失敗したら生成を打ち切ってエラーを出力するようにします。

コード

C++

#include <iostream>
#include <random>
#include <vector>
#include <cmath>
using namespace std;

int main(void)
{
    long long N, S, L, R;
    cin >> N >> S >> L >> R;

    const double k = 3.33; //比例定数
    const long long retry_max = 10; //生成回数の最大値

    random_device rd;
    mt19937_64 mt(rd());

    if(S / N < L || R < S / N)
    {
        //平均が範囲外の場合は生成不可能
        cout << "The average is out of range." << endl;
        return 0;
    }

    //乱数の生成範囲を決める
    long long l_min = max(L, 2 * S / N - R), r_max = min(2 * S / N - L, R);
    double E_inv = k * sqrt(N);
    long long gap = min((long long)((r_max - l_min) / E_inv), (r_max - l_min) / 2);
    long long l = l_min + gap, r = r_max - gap;

    uniform_int_distribution<long long> uid(l, r);

    vector<long long> array(N);

    long long retry = 0; //生成のやり直し回数
    for(; retry <= retry_max; ++retry)
    {
        long long sum = 0;
        for(long long i = 0; i < N; ++i)
        {
            array[i] = uid(mt);
            sum += array[i];
        }

        //誤差を修正する
        long long err = sum - S;
        bool fail = false; //生成失敗判定
        for(long long i = 0; i < N; ++i)
        {
            array[i] -= err / N;
            if(array[i] < L || R < array[i])
            {
                fail = true;
                break;
            }
        }
        err -= err / N * N;
        if(fail)
            continue;

        //微調整(N個の乱数から重複のないように選ぶ(Fisher-Yates shuffle))
        vector<long long> fisher;
        if(err > 0)
        {
            for(long long i = 0; i < N; ++i)
            {
                if(L < array[i])
                    fisher.push_back(i);
            }
        }
        else if(err < 0)
        {
            for(long long i = 0; i < N; ++i)
            {
                if(array[i] < R)
                    fisher.push_back(i);
            }
        }

        long long fisher_size = fisher.size();
        if(abs(err) > fisher_size) //生成失敗判定
            continue;

        long long idx;
        for(long long upper = fisher_size - 1; err > 0; --err, --upper)
        {
            uniform_int_distribution<long long> choose(0, upper);
            idx = choose(mt);
            array[fisher[idx]] -= 1;
            fisher[idx] = fisher[upper];
        }
        for(long long upper = fisher_size - 1; err < 0; ++err, --upper)
        {
            uniform_int_distribution<long long> choose(0, upper);
            idx = choose(mt);
            array[fisher[idx]] += 1;
            fisher[idx] = fisher[upper];
        }

        break;
    }

    //retry_max回生成に失敗した場合
    if(retry > retry_max)
    {
        cout << "Failed to generate sequence." << endl;
        return 0;
    }

    //数列とその和を出力
    long long sum = 0;
    for(long long i = 0; i < N; ++i)
    {
        cout << array[i] << " ";
        sum += array[i];
    }
    cout << endl;
    cout << sum << endl;

    //生成のやり直し回数を出力
    cout << "retry: " << retry << endl;

    return 0;
}

入出力例

3.の入出力例と同じ入力でもう一度生成してみます。

入力例1

20 10000 100 1000

出力例1

801 189 200 526 611 708 373 771 735 270 231 302 243 529 376 326 858 861 495 595 
10000
retry: 0

入力例2

100 0 -1000 1000

出力例2

-158 906 -494 52 -731 751 28 394 864 -911 749 651 184 156 833 -424 -279 509 866 -507 -430 -918 190 499 531 346 -931 -701 311 187 -274 316 174 -131 713 505 -20 718 -310 808 -862 784 664 -333 266 -592 645 588 -856 -430 811 440 -479 -231 -893 -93 737 248 317 -604 378 299 -880 -690 476 525 -687 -202 732 -197 -519 -97 -897 -681 -5 -675 -851 -341 528 -648 -912 176 -355 -447 -744 -293 -388 -233 -880 290 467 132 -382 807 643 672 475 234 167 -146 
0
retry: 1

どちらも少ない回数で生成に成功しました！

4. と同様に、入力例2と同じ入力で乱数列を生成する試行を $10000$ 回おこなうコードを書いて実験したところ、一発で生成に成功する確率は約85.2%、$1$ 回以内のやり直しで成功する確率は約97.8%で、最悪でも $5$ 回のやり直しで生成に成功しました。4. では一発で成功する確率は約26.8%だったので、生成に成功する確率をかなり上げることができました。

なお、比例定数 $k$ は思ったより大きくても問題なく、$k=6.0$ 程度なら $10$ 回以内のやり直しで生成に成功するようです。

おまけ

せっかくなので、今回のコードをクラスにまとめて利用しやすくしておきます。生成のやり直し回数はデフォルトを $10$ 回、比例定数 $k$ はデフォルトを $6.0$ として、オプションとして値を調整できるようにしておきます。

コードはほぼ同じなので、別サイトに載せたコードへのリンクだけ載せておきます。

https://gist.github.com/luuguas/d97cd2d784a2d81f49bbcf53f7e9066d

Wandbox というウェブ上でプログラムを実行できるサイトにも載せてみました。色んな入力を試したい方はどうぞ。

https://wandbox.org/permlink/IHM5UvUHooBWfHhM

おわりに

データを分析したり数式を計算したりと色々してみましたが、無事に上手く実装できたので良かったです。気が向いたら活用方法についても考えてみようと思います。

記事の感想や「この記事に載っているもの以外にも良いアイデアがあるよ」「こんな場面で活用できるよ」といった情報など、コメント欄で教えていただけると嬉しいです。最後まで記事をご覧いただきありがとうございました。

C++ 以外の言語をいくつか調べたところ、例えば Python では random.randint() をそのまま使うことができ、JavaScript では mt.js というメルセンヌ・ツイスタのライブラリを導入すると使える nextInt() というメソッドで同様の実装ができます。
↩
random_device と mt19937_64 はどちらも乱数を生成しますが、random_device はシード値が設定できない(=再現性がない)ことや乱数分布の均等性が保証されていないことから、乱数生成器としては不適切です。
そのため、random_device はシード値を生成するときだけ使い、そのシード値を用いて mt19937_64 を初期化して乱数を生成します。
↩
$E$ の右下の $0.99$ は生成された誤差の約 99% が $E$ 以下であるという意味ですが、後述の $E_{0.75}$ と区別するためだけのものなので深い意味はありません。
↩