PEをMathematicaで

Project Eulerに挑戦してみよう

Project Euler 234 / 半分割可能数

201-300

整数n(≥4)に対して最大の素数(≤√n)を n の下位素数平方根(lower prime square root)とし，lps(n)であらわす。同様に最小の素数(≥√n)を n の上位素数平方根(upper prime square root)とし，ups(n)であらわす。

lps(4) = 2 = ups(4), lps(1000) = 31, ups(1000) = 37 である。

lps(n) と ups(n) のどちらかが n を割り切るが，両方ではないとき整数n(≥4)を半分割可能(semidivisible)と呼ぶ。

15 を超えない半分割可能な数は8, 10, 12で，それらの合計は 30 である。15 は lps(15) = 3 と ups(15) = 5の両方の倍数なので半分割可能でない。さらに例を挙げると，1000 までの半分割可能な整数 92 個の合計は 34825 である。

999966663333 を超えない半分割可能な数全ての合計を求めよ。

Problem 234 - Project Euler


lps(n)=ups(n) だと「両方ではない」の条件をみたさないので，n は平方数ではありません。 p< \sqrt{n}< q をみたす連続する素数 p, q が存在します。平方すると

p^2< n< q^2

この範囲で
(p の倍数の和) + (q の倍数の和) - (pq の倍数の和)*2
を計算して足しあわせます。

ただし，最後の素数 p についてだけ，n が 999966663333 を超えるかどうかのチェックが必要です。下のコード中で sum1=0 なのですが，一応調べました。


Sum を2次関数であらわに書くと若干早くなります。1秒→0.8秒の微妙な変化ですが。

Project Euler 285 / ピタゴラス・オッズ

201-300

Albert は正整数 k を選び，さらに二つの実数 a, b を区間 [0,1] の一様分布からランダムに選ぶ。さらに和 (ka+1)^2+(kb+1)^2 の平方根を計算し,最も近い整数に丸める。結果が k に等しければ，彼は k 点を得る。それ以外は 0 点である。

k=6, a=0.2, b=0.85 ならば

(ka+1)^2+(kb+1)^2=42.05

42.05 の平方根は 6.484... で, 最も近い整数は 6 になる。この値は k=6 に等しいため, 彼は 6 点を得る。

k=1, k=2, ..., k=10 で 10 回プレイした場合，小数第6位で四捨五入した合計点の期待値は 10.20914 である。

k=1, k=2, k=3, ..., k=10^5 で 10^5 回プレイした場合，小数第6位で四捨五入した合計点の期待値はいくらか？

Problem 285 - Project Euler

確率＝面積比

k 点を得る条件は

\left(k-\dfrac{1}{2}\right)^2 \leqq (ka+1)^2+(kb+1)^2 <\left(k+\dfrac{1}{2}\right)^2

x=ka+1,\, y=kb+1 とおきます。

D_1:\left(k-\dfrac{1}{2}\right)^2 \leqq x^2+y^2 <\left(k+\dfrac{1}{2}\right)^2

点 (x, y) は領域 D_2:[1,\, k+1]\times [1,\, k+1] に一様に分布するので，k 点を得る確率は D_2 と D_1 \cap D_2 の面積比です。

D_2 の面積は k^2 です。 D_1 \cap D_2 の面積を S(k) とすると，得点の期待値 E は次のように表されます。

E=\sum \dfrac{S(k)}{k^2}\cdot k=\sum \dfrac{S(k)}{k}

数値積分で S(k) を求めて，その和をとりました。1分半かかりました。

数値積分→普通の積分

Mathematica は単に Integrate するだけで S(k) の具体形を求められます。これを使うほうが速く，無事1分におさまりました。

Project Euler 218 / 完全な直角三角形

201-300

a=7, b=24, c=25 の直角三角形について考える。この三角形の面積は84であり，完全数6と28で割り切れる。さらに, この三角形は原始(primitive)直角三角形である。つまり，gcd(a,b)=1, gcd(b,c)=1 をみたす。さらに，c は平方数である。

以下の条件をみたす直角三角形を完全(perfect)と呼ぶ。

  • 原始直角三角形である
  • 斜辺が平方数である

以下の条件をみたす直角三角形を超完全(super-perfect)と呼ぶ。

  • 完全な直角三角形である
  • 面積が完全数である6と28の倍数である

c≤10^16 をみたす完全な直角三角形のうち，超完全でないものはいくつあるか。

Problem 218 - PukiWiki


a, b, c は直角三角形の3辺で，互いに素なので次のようにおけます。m と n は互いに素で偶奇が異なります。

a=m^2-n^2,\, b=2mn,\, c=m^2+n^2

c は平方数なので m^2+n^2=k^2 とおけます。m と n が互いに素であることも考えると，もう一度一般形を導けます。

m=p^2-q^2,\, n=2pq,\, k=p^2+q^2

p と q は互いに素で偶奇が異なります。
次に直角三角形の面積をSとおいて p, q で表します。

S=\dfrac{ab}{2}=(m^2-n^2)mn

=2pq(p+q)(p-q)(p^2-2pq-q^2)(p^2+2pq-q^2)

これが「6と28の倍数」＝「84の倍数」となる条件は 3, 4, 7 で割り切れることです。このうち3と4の条件は簡単に処理できます。

  • p と q の偶奇が異なることから 2pq は4の倍数なので，S も4の倍数
  • pq(p+q)(p-q) が3の倍数であることは合同式で簡単に示せます（大学入試でもたまに出題されます）

7の条件は Mathematica でしらみつぶししました。どの場合も S は7で割り切れるので，完全な直角三角形はすべて超完全です。

この問題の答えは「0」。意外な結果になりました。

Project Euler 158 / 左隣の文字より辞書順で後になる文字がちょうど1個となる文字列の探査

101-200

26文字のアルファベットから3個の異なった文字を取ると，長さ3の文字列を作ることができる。例えば 'abc'，'hat'，'zyx' である。
この3つの例について調べると，'abc'は左隣の文字より辞書順で後になる文字が2個ある。'hat'では左隣の文字より辞書順で後になる文字がちょうど1個であり，'zyx'では0個である。
左隣の文字より辞書順で後になる文字がちょうど1個となるような長さ3の文字列は全部で10400個ある。
n ≤ 26 個の異なるアルファベットからなる文字列について考える。p(n) を左隣の文字より辞書順で後になる文字がちょうど1個あるような長さ n の文字列の個数とする。p(n) の最大値を求めよ。

Problem 158 - PukiWiki


入試問題でたまに見るタイプの問題ですね。アルファベットのかわりに1から26までの数で考えます。

  1. 26 個の数から n 個を選ぶ。C(26, n) 通り。
  2. n 個の文字を2つのグループ A, B にわける。A か B が空集合になる場合は不可なので 2^n-2 通り。
  3. A, B の要素をそれぞれ小さい順に並べる。
  4. (Aの最小値)>(Bの最大値)となる場合は不可。n-1 通り。

まとめると

p(n)={}_{26}\mathrm{C}_n \cdot (2^n-n-1)

これの最大値を求めます。


Mathematica らしくグラフも描いてみました。最大値は p(18) です。

f:id:variee:20170916021750p:plain

Project Euler 141 / 平方数でもある累進数 n の調べ上げ

101-200

正の整数 n を d で割った商と余りをそれぞれ q と r で表す。d, q, r を適当に並び替えたときに等比数列になる場合がある。

たとえば58を6で割ると商が9で余りが4である。4, 6, 9は公比3/2の等比数列になっている。このような n を累進数と呼ぶことにする。

いくつかの累進数9や 10404=102^2 は平方数になっている。100000未満の累進平方数の和は 124657 である。

10^12 未満の累進平方数の総和を答えよ。

Problem 141 - Project Euler

実例を作ってみよう

まずは定義どおりに立式して，累進平方数の具体例を見てみました。等比数列の条件は等比中項で処理しています。

この方法で総和を求めることもできますが，10^8までで約40秒かかります。


関数形を限定

等比数列の条件を

d=kab,\, q=ka^2,\, r=kb^2

として処理します。公比の a/b は既約分数としてよいので，aとbは互いに素です。このとき，

n=dq+r=kb(ka^3+b)

これが平方数になる組を探すと，約2分半で答えが出ます。