裏 RjpWiki

Julia ときどき R, Python によるコンピュータプログラム,コンピュータ・サイエンス,統計学

組み合わせの数, nCr, binomial(n, r)

2021年07月12日 | ブログラミング

n から r を取り出す取り出し方,すなわち nCr ではあるが,n, r によってはとんでもなく大きな数になる。

そもそも,そんな数が必要なことは普通はない。競技プログラミング(略して,競プロ)の世界では,それでは困るので(?),「1000000007 で割った余り(つまり剰余だね)」を答えさせる。そこでテクニックとして,いろいろあるが,以下のような Python プログラムを書くことができる。

from functools import reduce

# フェルマーの小定理を応用した(組み合わせの数 % mod)
def comb(n, r, mod=10 ** 9 + 7):
    numerator = reduce(lambda x, y: x * y % mod, range(n-r+1, n+1))
    denominator = reduce(lambda x, y: x * y % mod, range(1, r+1))
    return numerator * pow(denominator, mod - 2, mod) % mod

comb(2021, 37) # 717019240
comb(2021, 37, 4) # 0 この解は誤り

Julia だと

function comb(n, r, mod=10^9 + 7)
    red(x, y) = x * y % mod
    numerator = reduce(red, n-r+1:n)
    denominator = reduce(red, 1:r)
    numerator * powermod(denominator, mod - 2, mod) % mod
end

comb(2021, 37) # 717019240
comb(2021, 37, 4) # 0 この解は誤り

2021C37 mod 4 は

binomial(big(2021), 37) % 4 # 3 正しい答え

でなければならない。

ちなみに,

binomial(big(2021), 37) % (10^9 + 7) # 717019240

なので,正しい。

要するに,このプログラムは mod が素数であるときのみ正しい答えを返す

素数でないときは,不正な値をしれっと返す。

コメント
  • X
  • Facebookでシェアする
  • はてなブックマークに追加する
  • LINEでシェアする

和が 1 に近くなる単位分数の集合(1/n; n = 2, 3, ...)

2021年07月12日 | ブログラミング

分子が1で相異なる分母を持つ分数の和がほぼ1になる分数の組み合わせを求める。

https://mobile.twitter.com/yugokitajima/status/1388709059776286726

紹介しているサイトでは 分母が 2〜9 の場合としているが,
https://qiita.com/takuma-saito/items/4a0e6cfffb532ebc6774

大元の記事の画像を見ると少なくとも分母は12までありうると思われる(もっとあるかも知れないが)ので12まで考えよう。

Julia には // 演算子があるので好都合。

using Combinatorics
n = 11
all_perm(x, n) = Iterators.product([x for i = 1:n]...)
diff = Vector{Float64}()
rational = Vector{Rational{Int64}}()
set = Vector{Vector{Int64}}()
for i in all_perm([0, 1], n)
    s = 0
    num = []
    for (j, d) in enumerate(i)
        if d != 0
            s += 1//(j+1)
            append!(num, j+1)
        end
    end
    if abs(1 - s) < 0.01 # 0.004
        append!(diff, abs(float(1-s)))
        append!(rational, s)
        append!(set, [num])
    end
end
o = sortperm(diff);
diff = diff[o];
rational = rational[o];
set = set[o];
for (diff0, rational0, set0) in zip(diff, rational, set)
    println("$diff0\t$rational0\t$set0")
end

結果は,以下の通り。

0.0004329004329004329 2311//2310 [2, 6, 7, 10, 11] が誤差最小解。
最初のの要素は誤差,二番目は分数の和の分数表現,三番目が分母の集合。
1//2 + 1//6 + 1//7 + 1//10 + 1//11 = 2311//2310 = 1.0004329004329005

1/10, 1/11, 1/12 を使用可にしただけで,解の範囲が拡がった。

0.003968253968253968 253//252 [3, 4, 6, 7, 9] が 1/2 〜 1/9 だけを使った場合の誤差最小解であるが,それまでにずいぶんたくさんある。

0.0 1//1 [2, 3, 6]
0.0 1//1 [2, 4, 6, 12]

0.0004329004329004329 2311//2310 [2, 6, 7, 10, 11]
0.0004329004329004329 2311//2310 [3, 4, 7, 10, 11, 12]

0.0007575757575757576 1319//1320 [3, 4, 5, 8, 11]
0.0007575757575757576 1319//1320 [2, 5, 8, 11, 12]
0.0007575757575757576 1319//1320 [3, 5, 6, 8, 11, 12]

0.00202020202020202 496//495 [2, 5, 9, 10, 11]
0.00202020202020202 496//495 [3, 5, 6, 9, 10, 11]
0.00202020202020202 496//495 [4, 5, 6, 9, 10, 11, 12]

0.002777777777777778 361//360 [2, 6, 8, 9, 10]

0.002777777777777778 361//360 [3, 4, 8, 9, 10, 12]

0.003210678210678211 27809//27720 [3, 5, 7, 8, 9, 11]

0.003210678210678211 27809//27720 [4, 5, 7, 8, 9, 11, 12]

0.003968253968253968 253//252 [2, 4, 7, 9]

0.003968253968253968 253//252 [3, 4, 6, 7, 9]
0.003968253968253968 253//252 [2, 6, 7, 9, 12]

0.004365079365079365 2509//2520 [4, 5, 6, 7, 8, 9]

0.004365079365079365 2509//2520 [3, 5, 7, 8, 9, 12]

0.005555555555555556 179//180 [3, 4, 5, 9, 10]

0.005555555555555556 179//180 [2, 5, 9, 10, 12]
0.005555555555555556 179//180 [3, 5, 6, 9, 10, 12]

0.006313131313131313 787//792 [2, 6, 8, 9, 11]

0.006313131313131313 787//792 [3, 4, 8, 9, 11, 12]

0.007142857142857143 139//140 [2, 4, 7, 10]

0.007142857142857143 139//140 [3, 4, 6, 7, 10]
0.007142857142857143 139//140 [2, 6, 7, 10, 12]

0.007575757575757576 133//132 [2, 4, 6, 11]

0.007575757575757576 133//132 [2, 3, 11, 12]

0.0079004329004329 9167//9240 [3, 5, 7, 8, 10, 11]

0.0079004329004329 9167//9240 [4, 5, 7, 8, 10, 11, 12]

0.008333333333333333 119//120 [2, 5, 6, 8]

0.008333333333333333 119//120 [3, 4, 5, 8, 12]
0.008333333333333333 121//120 [3, 4, 5, 8, 10]

0.008333333333333333 121//120 [2, 5, 8, 10, 12]
0.008333333333333333 121//120 [3, 5, 6, 8, 10, 12]

0.009523809523809525 106//105 [2, 5, 6, 7]

0.009523809523809525 106//105 [3, 4, 5, 7, 12]

 

コメント
  • X
  • Facebookでシェアする
  • はてなブックマークに追加する
  • LINEでシェアする

PVアクセスランキング にほんブログ村

PVアクセスランキング にほんブログ村