次の企画のごあんない
先日グルミクで乙和ちゃんのイベントがあったじゃないですか。ポーカーのやつ*1。
あれのサブミッションみたいな位置付けで「トランプに描かれているキャラクターの組み合わせを揃えると称号が貰える」ってやつがあって、その称号がどうしても欲しかったんですね*2*3。
称号は「乙和と衣舞紀」を揃えて得られる「元気づけてあげたいから」と、「衣舞紀とノアと咲姫」を揃えて得られる「どんな乙和でも、大好きだよ」の2種類がありました。3人役はとくに揃えるのが難しくて、真面目にやってるとめちゃくちゃ時間がかかる。前回のTCGの記事の考察を使うと、どうやら85セット, 425回シャンシャンして称号が得られる期待値だそうです。やってられっかって感じですね。
ポーカーなので、カードが集まる前にもうキャラクターの揃いようがないな、揃う確率が下がったなってなることがあります。そういう時に「カードのリセット」が使えるんですよ。ポーカーを面白くするために付けられたリセット機能（フォールドにあたる）が、称号あつめでも優秀な選択肢になります。
そこで気になるのは「どれくらい出てこなかったらリセットすべきか」という問題です。
最初はこれ、さらっと解けちゃうんじゃなーい？って思ってました。でもナイーブに考えるとヤバいんですこの問題。それぞれの状況について伸るか反るかの二択があって、その全通りを調べようとすると全然間に合わないんですね。シミュレーションにかかる時間は私の一生で間に合わないし、解析を行うには頭が間に合わない。
詰んだと思って諦めかけていたある日、この問題設定ってマルコフ決定過程じゃね？と気付いたんですね。これが第一の進歩でした。MDPなら強化学習で解けるというのも覚えてました。やったことないけど。しかし勉強するのもなあと思い留まっていたら、ちょうど「ゼロから作るDL」の強化学習編が公開レビューをやってるという情報を耳にしました。これは何かの縁、強化学習やるしかない......と読み進めています。そして今4章を終えて、もう称号集め問題いけそうな気がする。なんたって環境（推移確率と報酬）が既知だから、DP組んでまわすだけだもん。むしろ強化学習の深淵をまだ見ていないまである。
DP組んでまわすだけとは言っても教科書のソースコードから対応させる部分は多いし、称号集め問題をどう記述するかが難しい。それでモチベ維持のために、「グルミク称号獲得の最適方策」をここで進捗報告していきます*4。今日はこのへんで！
推しのブロマイドを使ったトレーディングカードゲームが出ると聞いて私がやることは1つだった
何枚引けば推しが揃うかの期待値を求めます。
はじめに
お久しぶりです、こりです
twitterでは今は苗字?みたいなのがついて「いのこり」になっています。
いつ変えたんだかは忘れました。 昨年9月の「遠足のレヴュー@北海道」でお便りを読んでいただけたとき、小山百代さんに「何に居残ってるんだろうね」って言われた覚えがあるのでその前後だと思います。
当時は、大学に年単位の居残りをしていました(あっ......)。春からは働いていますのでご安心を。
では本題に移ります。
推しのブロマイドを使ったトレーディングカードゲームが出るらしいよ〜〜〜！*1
やった〜〜〜！どんな写真が見られるのか楽しみ！！！
カードゲームをしない私でも、推しのくらいは揃えたい！！と思うものです。
するとここでいくつかの問題が出てきます。
1つ目はズバリお金の問題です。トレーディングカードゲーム（以下TCG）はランダムでカードが入っているため、「これだけ買えば必ず入っている」という限度がないです。いっぱい買えばそれだけ推しのカードも出やすいでしょうが資金には限りがあります。この問題の単純な解決法は、TCGの名の通り「トレーディング」です。他プレイヤーとカードを交換/譲渡することで欲しいカードを手に入れることが可能で、それがまた醍醐味とされています（知らんけど）。
2つ目はその「トレーディング」に関わる問題で、「友人が少ない」ことです。TCGをする友人が少なく、また同担の友人が多いため、トレーディングで推しを集めることは困難を極めると想像されます。
従って、最小限の枚数の購入によって、自力で推しを揃えたいと考えるのは自然なことでしょう。
そこで今回は、こちらの問題を定義しました。
「 種類のカードが無限にあって、それぞれの種類は等しい確率で出るものとする。その中のきまった
種類のカードを揃えたい。カードを1枚ずつ引いていくとき、ちょうど
枚目で欲しかった
種類のカードが揃う確率はいくらか？」
この問題に対して私は2通りのアプローチ、すなわち ①ガチャシミュレーションによる統計的計算 ②高校数学による組み合わせと確率の解析的計算 により解を求めました。その結果、両者が一致することも確認しました。
そして動機だった「どれだけ引けば推しを揃えられるか？」という問題に、期待値として答えを与えることができました。その結果、マジで当たり前ですが「誰かと交換したほうがいいわね......」という結論に辿り着きました。*2
本記事は、同じ悩みを抱えるオタクのもとに届けばいいなという想い、そして「あわよくば、誰かがこの記事を見てトレーディングに応じてくださればいいな」という期待のもとに執筆します。
さて以下で計算過程の説明を行っていくのですが、「結局、何枚引けばhrkちゃんが揃うんや」とお急ぎの方は まとめまですっ飛ばしていただいて構いません。*3
問題への取り組み
今回の大目標は「どれだけ頑張れば推しが揃うか」の一点なので、,
,
といった一般の場合はそれなりにして、
,
という特殊なケースのみを考えます。
対象とするTCGは、通常28種, パラレル28種, スペシャル3種の計59種類で構成されています。スペシャルは出現確率が低いので「出ない」ものとして、「残り56種が等確率で出る」との簡略化を行いました*4。そして推しのカードは通常28種中2種類、パラレル28種中2種類の計4種類であると仮定しています。これをパターンAとします。もう一方のパターンBは、「パラレルとは単に通常カードにホログラムが入っただけで、絵柄は一緒（だから4種類も揃える必要はない）」と仮定した場合です。このとき通常カードとパラレルカードは同一視できて、 ,
の場合を考えればよいことになります。
一応「予想」という形で、一般の ,
の場合の表式を提示しますが証明は行いません。
シミュレーションによる解法
ランダムに の整数を重複ありで1つずつ生成し、記録していきます。この整数の組に
すべてが出現した時点で終了し、引いた枚数を記録します。これを1回のシミュレーションとします。シミュレーションを100万回行って、必要枚数の分布を近似的に求め、平均と分散を算出しました。
組合せ計算による解析的解法
種類のうち目当ての
種類がちょうど
枚引いたところで揃う確率を
とします。面倒なので
,
は当たり前に使うこととし、
と書きましょう。このとき
では, 明らかに
種類揃うはずがないので
のとき、欲しい全種類が順番自由で出なくてはならないので
はすぐにわかります。一般の は、漸化的に考えると
となります。このカッコ書きの確率を求めることが本解析のメインテーマとなります。
ここで新たな表記として、すべてのカードの種類の集合 , 推しカードの集合
を定義します。そして推しカードの部分集合を
とし、
枚目までに
は既に引いているが
は引けていない確率を
と表します。たとえば
のとき
は
となります。また空集合は
で表します。
このように表記することで、上のカッコ書きの確率は
と表されます。また、カードの番号付けに意味はないので
としても同じです。なので
の要素数
を用いて、
という表記も使っていきます。
種類揃えることを目指して、順に
種類と集めていきましょう。
まず
は、目当ての4種類以外（= 52種類）が出ていいので
です。
次に は「
枚引いて
は当たったけど
が当たらない」確率であり、
となります。
以降も同様にやっていくと、
となるので、求めたい確率 は
であることがわかりました。
また の表式に二項係数が現れることから、一般の
,
について
は
と表せると予想できます。残念ながらこの予想が正しいかは確認できていないです。
結果を確認する：パターンA
それでは結果をプロットして確認しましょう。
図は、56種類が等しい確率で出る無限枚のカードから、欲しい4種類が出るまで引く枚数の確率を示しています。シミュレーション結果が驚くほど理論値と整合的でした。 100万回もやればそうかもしれん。
この図からピークは 80枚 前後であることがわかります。
また解析解を用いた枚数の平均・分散の計算により、 枚 (
) がもっともらしいとわかりました。
を累積した関数
──実は
に等しい── を調べると、
- 103枚（6パック）引けば 50 % の確率で推しが揃う
- 162枚（9パック）引けば 80 % の確率で推しが揃う
- 333枚（19パック）引けば 99 % の確率で推しが揃う
ことがわかりました。上限である30パック540枚を買うと、99.97 %の確率で推しが揃います。ヤムチャしやがって...
結果を確認する：パターンB
次に示す結果は「パラレルカードと通常カードが同じ絵柄で、片方さえあれば満足できる場合」です。
この場合は、
- 35枚（2パック）引けば 50 % の確率で推しが揃う
- 62枚（4パック）引けば 80 % の確率で推しが揃う
- 146枚（8パック）引けば 99 % の確率で推しが揃う
という結果です。
大前提である「パターンAか」「パターンBか」によって大きく結果が違いますね。
まとめ、兼取り扱いについて
「パラレルカードと通常カードが別の絵柄で、両方揃えたい」場合、
- 103枚（6パック）引けば 50 % の確率で推しが揃う
- 162枚（9パック）引けば 80 % の確率で推しが揃う
- 333枚（19パック）引けば 99 % の確率で推しが揃う
「パラレルカードと通常カードが同じ絵柄で、片方さえあれば満足できる」場合、
- 35枚（2パック）引けば 50 % の確率で推しが揃う
- 62枚（4パック）引けば 80 % の確率で推しが揃う
- 146枚（8パック）引けば 99 % の確率で推しが揃う
ことが予想されました。
これらの結果をどう解釈し、購買行動に結びつけるかは提言しません。 個人的には「誰かと交換したほうがいいわね......」と思いました。 3人で3パックずつくらいかな...
また、これは数学モデルであって、現実の事象を再現するものではないため参考に留めていただきますようお願いします。
*1:文字を推しカラー(黄色)にしたかったのですが、めっちゃ見にくかったのでやめました。緑は推しキャラの髪色です。
*2:世界初のトレーディングカードゲーム「マジック・ザ・ギャザリング」の生みの親は組合せ数学でPh.D.を取った数学者らしく、TCGでこんな計算をしてる自分の「掌で踊らされてる」感がすごかったです。
しぇんゆバッテリー
#ラム乾 ありがとうございました！今週は謎を作れませんでした, すみません. 代わりと言っちゃなんですが第4話「ふたりの射程」をテーマに物理()の問題を作りました！
— いのこり (@ImCallinx2You) 2020年10月29日
雨嘉も神琳もよかったね......#アサルトリリィ #アサルトリリィ夢結 #アサルトリリィ弾丸 #アサルトリリィ木曜日 pic.twitter.com/COv7Zyc4NC
ってことで
ゆーじあさんの弾丸の軌道(第4話)を数値計算しました
誤差5cmで神琳撃てた
— いのこり (@ImCallinx2You) 2020年11月13日
ここまで精度良ければ私が雨嘉でしょ pic.twitter.com/7cnhkFvBQp
なわけないでしょ
神琳から雨嘉方向を見たときの弾道をプロットしました（時間が進むにつれて赤くなります）
— いのこり (@ImCallinx2You) 2020年11月27日
パラメータに依るんですけど風で浮いてからギューンと落ちてきます#アサルトリリィ pic.twitter.com/E98hqUIwoI
きれいだよ神琳
#include <bits/stdc++.h>
#define rep(i,n) for(int i=0; i<(n); i++)
using namespace std;
using vd = vector<double>;
using vv = vector<vd>;
const double l = 1000.0;
const double v0 = 1800.0;
const double g = 9.8;
const double k = 1e-3;
const double dt = 1e-3;
const double w0 = -10;
#define XLEN 6
vv init_f(double theta, double phi){
vv x0(2, vd (XLEN, 0.0));
x0[0][0] = 0;
x0[0][1] = 0;
x0[0][2] = 0;
x0[0][3] = v0*cos(theta);
x0[0][4] = v0*sin(theta)*sin(phi);
x0[0][5] = v0*sin(theta)*cos(phi);
return x0;
}
vd f(vd bx){
vd dxdt(XLEN, 0);
double v2 = pow(bx[3], 2) + pow(bx[4]-w0, 2) + pow(bx[5], 2);
dxdt[0] = bx[3];
dxdt[1] = bx[4];
dxdt[2] = bx[5];
dxdt[3] = -k*sqrt(v2)*bx[3];
dxdt[4] = -k*sqrt(v2)*(bx[4]-w0);
dxdt[5] = -k*sqrt(v2)*bx[5] - g;
return dxdt;
}
vd rk4(vd bx, double dt){
vd ax = bx, k, a = bx, h = {1,2,2,1,1};
rep(j,4){
k = f(a);
rep(i,XLEN){
k[i] *= dt;
a[i] = bx[i] + k[i]/h[j+1];
ax[i] += k[i]*h[j]/6.0;
}
}
return ax;
}
vd shot(double theta, double phi, int n=0){
vd yz(2, 0);
vv x = init_f(theta, phi);
double t = 0.0;
while(x[0][0]<l){
t += dt;
x[1] = rk4(x[0], dt);
if(n){
printf("%lf ", t);
rep(i,XLEN)printf("%.10lf%c", x[1][i], i==XLEN-1?'\n':' ');
}
swap(x[0],x[1]);
}
double a = x[0][0]-l;
double b = l-x[0][1];
rep(i,2)yz[i] = (x[1][i+1]*a+x[0][i+1]*b)/(a+b);
return yz;
}
int main(){
double vy_max = 0.1;
double vy_min = 0.0;
double vz_max = 0.1;
double vz_min = 0.0;
double r;
double vy, vz;
double theta, phi;
vd mark;
while(r>1e-3){
vy = (vy_max+vy_min)/2.0;
vz = (vz_max+vz_min)/2.0;
theta = atan(sqrt(vy*vy + vz*vz));
phi = atan2(vy, vz);
mark = shot(theta, phi);
if(mark[0]>0)vy_max = vy;
else vy_min = vy;
if(mark[1]>0)vz_max = vz;
else vz_min = vz;
r = sqrt(pow(mark[0], 2) + pow(mark[1], 2));
}
printf("\ntheta=%lf[deg], phi=%lf[deg]\n", theta*180/M_PI, phi*180/M_PI);
printf("y=%lf[m], z=%lf[m]\n\n", vy*l, vz*l);
return 0;
}解くべき運動方程式
空気抵抗と風、重力が働くときの弾丸の軌道を計算します。
弾丸のように速度が大きい時、空気抵抗の大きさは二乗に比例します。これは弾丸前後面での空気の圧縮・希薄による圧力差が原因になっています*1。風がある時は、向かい風や追い風（=相対風速 ）をもとに空気抵抗が決まります。
したがって運動方程式は
で表されます。
ただし、毎秒1800 m という弾丸の速度は音速を超えているので二乗則が正常に働いているかどうかは分かりません......（文献を当たりましたが、雨嘉さんに撃てるほど小さな質量でのデータが見当たりませんでした）。
運動方程式の右辺第1項が風場ありの空気抵抗、第2項が重力を表します。空気抵抗と重力があるときの軌道は特殊な例を除いて解析解を持ちません。本問もそうなので、数値的に解を求めましょう。
コードの概説
まずアルゴリズムのフローチャートはこんな感じです。
試し撃ちを行い、着弾位置から判断して発射角度を修正することを繰り返し行います。雨嘉さんは一発で当てちゃうのですごい
ではコードの解説に移ります。
個人的な読みやすさのために以下4つを使います。
rep(i,n): iを0からn-1まで行うforループ
vd: 実数の1次元配列
vv: 実数の2次元配列
XLEN: 6 (= 3次元 * 2)
コードは以下5つの関数から構成されます。
まず
- init_f：初期値の設定
位置と速度の初期値 x0 を与えます。0~2が位置、3~5が速度に対応します*2。雨嘉さんの位置から 方向に速さ v0 で発射します。
vv init_f(double theta, double phi){
vv x0(2, vd (XLEN, 0.0));
x0[0][0] = 0; // x
x0[0][1] = 0; // y
x0[0][2] = 0; // z
x0[0][3] = v0*cos(theta); // v_x
x0[0][4] = v0*sin(theta)*sin(phi); // v_y
x0[0][5] = v0*sin(theta)*cos(phi); // v_z
return x0;
} はx軸と速度方向のなす角、
はy軸方向への回転角を表します.
- f：解くべき運動方程式の指定
運動方程式の右辺を記述します。
vd f(vd bx){
vd dxdt(XLEN, 0);
double v2 = pow(bx[3], 2) + pow(bx[4]-w0, 2) + pow(bx[5], 2);
dxdt[0] = bx[3]; // dx/dt = v_x
dxdt[1] = bx[4]; // dy/dt = v_y
dxdt[2] = bx[5]; // dz/dt = v_z
dxdt[3] = -k*sqrt(v2)*bx[3]; // dv_x/dt = a_x
dxdt[4] = -k*sqrt(v2)*(bx[4]-w0); // dv_y/dt = a_y
dxdt[5] = -k*sqrt(v2)*bx[5] - g; // dv_z/dt = a_z
return dxdt;
}状態 bxを使って微分係数 dxdtを返します。
- rk4：4次のルンゲ・クッタ法を用いた数値積分コア
現在の状態 bxを使って、時間 dt後の状態 axを4次精度で計算します。
vd rk4(vd bx, double dt){
vd ax = bx, k, a = bx, h = {1,2,2,1,1};
rep(j,4){
k = f(a); // k_1~k_4を方程式fに基づいて順次計算
rep(i,XLEN){
k[i] *= dt; // 刻み幅をかける
a[i] = bx[i] + k[i]/h[j+1]; // 刻み幅の係数をかけて次の引数とする
ax[i] += k[i]*h[j]/6.0; // 戻り値に加算
}
}
return ax;
}ルンゲ＝クッタ法 - Wikipedia
めちゃめちゃ短く書けたから見て欲しい。
あのルンゲクッタさんをこんなに可愛くしちゃうなんてすごいことよ(これ5話だった)
- shot：弾丸の軌道計算
から始めて、弾丸が
を通過するまで数値積分を繰り返します。そして
での
,
の値を通過前後の2点から線形補間して返します。
vd shot(double theta, double phi, int n=0){
vd yz(2, 0);
vv x = init_f(theta, phi); // 与えられたtheta, phiで初期値を呼ぶ
double t = 0.0;
while(x[0][0]<l){ // xがlを超えたら終了
t += dt;
x[1] = rk4(x[0], dt); // 次の時刻での位置・速度をルンゲクッタで計算
if(n){ // (軌道出力用オプション)
printf("%lf ", t);
rep(i,XLEN)printf("%.10lf%c", x[1][i], i==XLEN-1?'\n':' ');
}
swap(x[0],x[1]); // 状態の更新
}
double a = x[0][0]-l;
double b = l-x[0][1];
rep(i,2)yz[i] = (x[1][i+1]*a+x[0][i+1]*b)/(a+b); // x=lでの通過位置を求める
return yz;
}- main：神琳さんに（誤差 1 mmで）命中するような初期値を探索する
風と重力の影響は雨嘉さんからみて弾丸を右下に流すよう働くので、いい感じに左上に向けて撃たねば命中しません。そこで射撃シミュレーションを繰り返し行って、初速のy, z成分 (,
) を二分探索します。
での弾丸の通過するy, z座標 (
,
) それぞれについて正負を判定します。
YES (通過位置が左/上) ならもっと初速を右/下にしてよいです。
NO (通過位置が右/下) なら初速は今より左/上でなければいけないです。
このアルゴリズムは、日本語で上手く言えないですが ,
を暗に仮定しています。
つまり通過位置が左だったからと言って最適初速がより右にあるとは限らないのです。
今回はたまたま上手くいったので良いですが、二分探索でこぼしてしまい解に収束しない可能性もある判定法でした。代わりに最急降下法などを使う方が自然です。
神琳さんから誤差 1 mm 以内になれば終了します。1周目での誤差が 100 mくらいあるので、20回程度のシミュレーションで収束することが期待されます。
int main(){
//初速方向の探索範囲
double vy_max = 0.1;
double vy_min = 0.0;
double vz_max = 0.1;
double vz_min = 0.0;
double r;
double vy, vz;
double theta, phi;
vd mark;
while(r>1e-3){ // 1 mm 以内になるまで繰り返す
vy = (vy_max+vy_min)/2.0;
vz = (vz_max+vz_min)/2.0;
// vy, vz から theta, phi に変換
theta = atan(sqrt(vy*vy + vz*vz));
phi = atan2(vy, vz);
// 射撃シミュレーションを行う
mark = shot(theta, phi);
// 探索範囲の更新
if(mark[0]>0)vy_max = vy;
else vy_min = vy;
if(mark[1]>0)vz_max = vz;
else vz_min = vz;
// 神琳さんとの距離を求める
r = sqrt(pow(mark[0], 2) + pow(mark[1], 2));
}
printf("\ntheta=%lf[deg], phi=%lf[deg]\n", theta*180/M_PI, phi*180/M_PI);
printf("y=%lf[m], z=%lf[m]\n\n", vy*l, vz*l);
return 0;
}おわり。パラメータの説明や結果・考察はコピ本に載せました。
これは神雨だ......誰が何を言おうと神雨なんだ......
*1: 抗力 - Wikipedia 圧力抗力の項. 一方速度が小さい時は空気との摩擦が支配的で, 大きさは一乗に比例します(摩擦抗力).
*3:そういえばrk4はt依存に対応していませんが, こちらのアップデートは予定していません. 楽そうだからやってもいいけど.
■
アサルトリリィは西暦何年のできごと？
一柳梨璃さんの誕生日をもとに調べます。
#include <bits/stdc++.h>
using namespace std;
const double syn = 22.0 + 12.0/24.0 + 44.048/1440.0;
const double year = 365.242189;
int main(){
int range;
printf("until: AD");
scanf("%d", &range);
if(range<=0){
printf("error: please input an upcoming year in A.D.");
return 0;
}
int day = 0;
double dday;
int month;
for(int i=2021; i<range; i++){
if(i%400==0)
day += 366;
else if(i%100==0)
day += 365;
else if(i%4==0)
day += 366;
else
day += 365;
if(day%7!=5)continue;//Tuesday?
if(day%60!=48)continue;//Kanoe-Tatsu?
dday = year * (i-2020) + 26.4;
dday /= syn;
dday -= (int)dday;
dday *= syn;
if(dday< 20.2 || dday > 21.2)continue;//moon phase:20.7(\pm 0.5)?
printf("%d %lf\n", i, dday);
}
return 0;
}データの出典はこちら*1
とくに説明はしないですが、アニメ アサルトリリィ BOUQUET 5話*2で示唆された、梨璃ちゃんの誕生日前日6月18日の特徴「①火曜日で、②庚辰で、③月齢20.7」*3をもとに候補を絞りました。これらは太陽と月の運動およびグレゴリオ暦で決まります。月の運動は平均朔望周期を使い、月齢の判定は0.5日まで許すと甘めにしています。
旧暦およびそこから決まる六曜は少し面倒*4なため、CASIO様の旧暦計算サイト*5に任せています。CASIO様のサイトからは「火曜日で旧暦5/22で友引な6月18日は西暦何年？」を一括で調べられないので、このプログラムで絞り込むことが重要です。
計算の結果、次は「西暦2052年」らしいことがわかりました。CASIO様サイトで旧暦も確認したところ合っているようです。その次は私のプログラムによると西暦4261年らしいのですが旧暦は合っていないようです。
またどれだけ未来の結果まで信用できるかは数値的にも物理的にも非自明であり*6、少なくともミランコビッチサイクルが効いてくる1万年スケールではもう使わないほうがいいです。
*1:2020年6月18日の曜日・月齢・干支のデータ：六曜・月齢・旧暦カレンダー (２０２０)令和２年11月(霜月)、 干支：壬辰 - Wikipedia、庚辰 - Wikipedia、 朔望月の自転周期：月 - Wikipedia、 太陽年：太陽年 - Wikipedia
*2:第5話 ヒスイカズラ| STORY | アサルトリリィBOUQUET(ブーケ) アニメ公式サイト
*3:twitterでは友引も調べたと書いたけどやってなかった。すみません
*4:閏月の置き方が面倒。 改暦 - Wikipedia
*6:物理的には明日地球が太陽に落っこちることだってありえる
Atomic String
任意の英小文字列 (
) が与えられます。この文字列からいくつかの小文字を大文字にすることで、元素記号を並べた文字列にできるかを判定し、ないなら"No"、あるならそのような文字列の個数を
で割った余りで答えます。
s = "co" が与えられたとき、元素記号での表記が "CO", "Co" の2通りあるので2と答えます。
s = "a" には対応する表記が無いのでNoと答えます。
以下C++での解答:
#include <iostream>
#include <vector>
#include <set>
#define LMAX 100000
using namespace std;
const int MOD = 1000000007;
int main(){
vector<string> p = {"H", "He", "Li", "Be", "B", "C", "N", "O", "F", "Ne", "Na", "Mg", "Al", "Si", "P", "S", "Cl", "Ar", "K", "Ca", "Sc", "Ti", "V", "Cr", "Mn", "Fe", "Co", "Ni", "Cu", "Zn", "Ga", "Ge", "As", "Se", "Br", "Kr", "Rb", "Sr", "Y", "Zr", "Nb", "Mo", "Tc", "Ru", "Rh", "Pd", "Ag", "Cd", "In", "Sn", "Sb", "Te", "I", "Xe", "Cs", "Ba", "La", "Ce", "Pr", "Nd", "Pm", "Sm", "Eu", "Gd", "Tb", "Dy", "Ho", "Er", "Tm", "Yb", "Lu", "Hf", "Ta", "W", "Re", "Os", "Ir", "Pt", "Au", "Hg", "Tl", "Pb", "Bi", "Po", "At", "Rn", "Fr", "Ra", "Ac", "Th", "Pa", "U", "Np", "Pu", "Am", "Cm", "Bk", "Cf", "Es", "Fm", "Md", "No", "Lr", "Rf", "Db", "Sg", "Bh", "Hs", "Mt", "Ds", "Rg", "Cn", "Nh", "Fl", "Mc", "Lv", "Ts", "Og"};
set<string> atom;
for(auto q: p){
q[0] = tolower(q[0]);
atom.insert(q);
}
string s;
cin >> s;
int l = s.size();
/*
if(l>LMAX){
cout << "s is longer than the limit.\n";
return 0;
}
if(!l){
cout << "s is null.\n";
return 0;
}
*/
vector<long long> dp(l+1, 0);
dp[0] = 1;
for(int i=0; i<l; i++){
string a;
a += s[i];
if(atom.count(a))dp[i+1] += dp[i];
if(!i)continue;
a = s[i-1] + a;
if(atom.count(a))dp[i+1] += dp[i-1];
dp[i+1] %= MOD;
}
if(!dp[l])cout << "No" << endl;
else cout << dp[l] << endl;
return 0;
}元素記号名は高々長さ2なので、「i文字目までが何通りで表せるか」は「i-1文字目~」と「i-2文字目~」を使った漸化式で表せます。上のプログラムではdpにしてるけど変数は3個とっておけばよいです。
凸レーダーチャート問題
問題：
— いのこり (@ImCallinx2You) 2020年4月18日
5段階、6要素のレーダーチャートで凸6角形になるグラフは何通り？
解答例です。言語はC++
#include <bits/stdc++.h>
using namespace std;
bool test(vector <int> a){
bool ok = true;
int ab, ac;
for(int j=0; j<6; j++){
ab = a[(j+1)%6];
ac = a[(j+5)%6];
if(a[j%6]<=ab*ac/(ab+ac)){
ok = false;
break;
}
}
return ok;
}
int all(vector<int> a, int n){
if(n==0)
return test(a);
int k;
for(int i=0; i<5; i++){
a[n-1] = i+1;
k += all(a, n-1);
}
return k;
}
int main(){
vector<int> a(6);
cout << all(a,6) << endl;
return 0;
}
コードの心臓は関数で、すべての要素のスコア
が「隣り合うスコアを結んだ時の線分と要素
の軸の交点のスコア
」より大きければ
を返します。1個でも成り立たなければ
です。
(右の白点がレーダーチャートの中心。ぐるり一周あるうちの3要素だけに注目している。)
は、角の二等分線の性質を使って
*1と求められます。
関数で全てのレーダーチャートの結果について
を行い、
の個数を数えます。
計算の結果、答えは 通りです。
*1:これは二等分すべき角を として、
であるのを
とした特殊なケースです。
一般にM段階N要素のレーダーチャートの場合を考えましょう。
N<3の場合、レーダーは多角形にならないので 0通りです。
N=3の場合、
なる
が 1個までなら三角形になるので 答えは
通りです。
N>3の場合、
として上のコードを一般化し、計算機をブン回せば
で求まるはずです。一般式や対称性を使えばもっと速いのかもしれませんが。
広い広い、Star Diamond。
Star Diamond ──それは、冬の夜空に輝く宝石。
6つの1等星；シリウス、リゲル、プロキオン、ポルックス、カペラ、そしてアルデバランからなる「冬のダイヤモンド」。オリオンや冬の大三角に比べれば知名度は低いが、いざ見つければその形状の美しさと広さに心を奪われることでしょう。
次の画像は、3月中旬の星空のようすです。
もう4月は始まっていますが、夜の早い時間であれば冬のダイヤモンドを西の空に見ることができます。
このダイヤモンドは、正六角形とはいえなくとも、整った形をしていてなんとも引き締まって見えます。まるで、冬の静かな冷気の中で、穏やかに結晶化したかのようです。
星の色も様々で、青白くまばゆいシリウスから穏やかな橙色のアルデバランまで幅広く、見飽きることはありません。
ダイヤモンドを構成する6つの星は別の星座に属しています。
おおいぬ座、こいぬ座、ふたご座、ぎょしゃ座、おうし座、オリオン座。
サーカスのステージのような六角形の上を、個性豊かな演者たちが跳ね回っているようにも見えます。
かくいう私、少女歌劇レヴュースタァライトで「Star Diamond」*1が取り上げられるまで個々の星座の神話も知らなかったし、ダイヤモンドとして認識したことすらありませんでした。
その理由の１つに、冬のダイヤモンドは実際に見上げてみると極めて広く、視界の隅の星々を結ぶ必要があるため、意識しないと見逃してしまうからです。
初めて冬空の六角形をなぞった時のことは私もよく覚えています。
あの頃には戻れない、何も知らなかった日々、胸を刺す衝撃を浴びてしまったから──。
はい。この衝撃を、皆様に定量的にお伝えしたいと思います。
今回のトリビアの種は「冬のダイヤモンドは 満月???個分の広さ」です。
(最後までスキップしてもらって構いません。)
（以下では先達の天文学者達の功績、とくに球面三角法を用いた説明を行いますが、ややこしいことを言わずとも全て高校理系数学程度で理解できます。さらに球面三角形の面積の証明(リンク先にお任せ)は、中学受験の問題のような面白さがあります。）
宇宙は3次元に広がっていますが、夜空の星は「天球」に並んでいるように見えます。この天球面上での面積として求めましょう。天球は仮想的な面なので、半径を定義しません。かわりに天球面上の2点間の距離は角度で表します。
球面は曲がっているので、平面で使っていた面積公式は成り立ちません。球面幾何学 - Wikipediaと呼ばれる数学から求める必要があります。球面幾何の大きな特徴は2つあって、
- 2点を結ぶ直線を球の大円で定義する
- 2直線の角度を大円のなす角度で定義する
ことです。この定義から次の図を見ると、
- 線分ABの長さは、
である
- 角度BACの大きさは、
である
ということになります。
六角形の面積を求めるために、まず図形の基本要素である三角形に分割してみましょう。隣り合っていない頂点と頂点を結べば、六角形は4個の三角形に分割することができます。この1つ1つの三角形の面積を求めて合計すれば、六角形の面積が得られます。
証明は以下のリンクで、ということにして、球面上の三角形の面積
は、球の半径を
, 三角形の内角を
,
,
として
と表されます。つまり内角の和からを引けば三角形の面積になります。
（「内角の和はだから、
になるのでは！？」と思われるかもしれません。実は「内角の和は
」は平らな面上*2でのみ成り立つ定理で、球面のような曲がった面上では一般に成り立ちません。）
証明はこちら。興味がある方はご覧ください。とても視覚的に分かりやすい説明です。
すると冬のダイヤモンドの面積は、
だとわかります。
というわけで六角形の内角を求めていきましょう。
そのためには6つの星の位置に関するデータが必要です。今回は天体アーカイブSIMBAD（http://simbad.u-strasbg.fr/simbad/）から赤経赤緯（地球公転軌道を赤道として天球に張った経度緯度）のデータICRS(ep=J2000)を引っ張ってきました（元論文は van Leeuwen F., A&A 474, 653-664 (2007) https://www.aanda.org/articles/aa/abs/2007/41/aa8357-07/aa8357-07.html）。
このデータを用いて、まず星と星の間の天球面上の距離を求めます。上述の通り距離は 星 観測者
星
のなす角なので、ベクトル
,
の内積で求められます。
星による三角形の辺々の長さ
がわかれば、球面で成り立つ余弦定理：
（球面三角法 - Wikipedia）を適用することで内角の大きさがわかります。
角と同様に六角形の内角を調べ上げて足し、最後に
を引けば冬のダイヤモンドの面積が求まります。
次に、基準となる満月の面積を求めましょう。満月の中心を基準として極座標で積分すれば、
と求まります（は月の半径、
は月までの距離）。
実は、満月の面積にはバラツキ（月だけに）が最大程度あります。
これは月が地球周りを楕円軌道で公転していて、月までの距離が変化するためです。満月の面積を基準にしてしまったことが今回のブログで最大の過ちです*3。
そのことを加味した上で、冬のダイヤモンドの面積と満月の面積
の比を求めると、
「冬のダイヤモンドは 満月個分の広さ」
という結果になりました。ざっと1万倍です。
輝き溢れる広大な冬空のステージを、今宵見上げてはみませんか？
*1:スタァライト九九組 6thシングル「Star Diamond https://revuestarlight.com/music/6thbdcd/、6つの星々をテーマにしたアルバム「ラ・レヴュー・エターナル https://revuestarlight.com/music/revuealbum/、3rd スタァライブ「Starry Diamond」 https://revuestarlight.com/event/4597/。
この記事は「Starry Diamond」ライブBD発売を記念したものであると同時に、岩田陽葵さんのお誕生日（4月3日）をお祝いするものです。...といえるのか...。
*2:このような世界をユークリッド空間、球面のような曲がった世界を非ユークリッド空間と呼びます。
*3:結果の不定要素は数々考えられますが、最大の原因がこの月の楕円軌道 () です。星々の固有運動や年周視差による
の変化は年間
程度で十分無視できます。
