読者です 読者をやめる 読者になる 読者になる

「セカンド・イクエイション」問題解答 ( CodeIQ )

CodeIQ コードゴルフ 数学

はじめに

CodeIQというところでプログラミング問題「セカンド・イクエイション」に解答しましたので、ネタバレです。 codeiq.jp

問題

問題および解説は、以下の記事をご覧ください。

codeiq.jp

提出解

以下、提出したgolf解 Ruby(60) です。( mじゃなくてnになってるのは、まあ… )

s=0;2.step(n=gets.to_i,2){|i|i.upto(n){|j|s+=j*j/i-i|1}};p s

さて、基本的な考え方は、解説記事にあるものとほぼ変わりません。

https://codeiq.jp/magazine/wp-content/uploads/2015/11/1-300x288.png

解説記事にあるこの図の通り、b を定めた時の a,c の組は、1\leq a \leq b/2 の範囲で  ([ b^2/4a ] - a )\times 2 +1 を足し合わせたものになります。
これを、各 b に渡って足し合わせた、  \sum_{1\leq a\leq b/2, b\leq m }^{} ( ([ b^2/4a ] - a )\times 2 + 1 ) を計算することで答えが出ます。

m=7 の場合だと、次の(a,b)の範囲で計算を行った結果を足し合わせる形になります。 f:id:ange1:20151117003734p:plain

はい。しかし、ここでコードを短くすることを考えると、a をそのまま使うよりも、偶数である 2a をベースにした方がすっきり書けることに気付きます。
なぜなら、 ( b^2/4a - a )\times 2 + 1 = b^2/2a/2\times2 - 2a + 1 = b^2/2a - 2a | 1 ( / は整数除算, | は bit or ) だからです。
※注: x/2*2+1 は x|1 と同等

という事で、i=2a を外側ループとして、j=b を内側ループとすると、冒頭にあるようなコードになります。m=7 の例での i,j の処理範囲は次のようになります。 f:id:ange1:20151117004828p:plain

高速化

 \sum_{1\leq a \leq b/2, b \leq m}^{} (b^2/4a-2a+1) の形に戻しましょう。 \sum_{b}^{} b^2 の形がある訳なので、数列の和として一発で求められそうな気がします。
…が、整数除算 /4a が付いているおかげで切り捨てられる分があり、そう簡単にはいきません。

しかし、/4a で切り捨てられる分、すなわち mod 4a による剰余分は周期的ですから、内側の bループをある程度減らせるのではないかと期待できます。 例えば a=1, m=6 の場合、b2 mod 4a は b=2 から周期的に 0,1,0,1… を繰り返して2周期半、なので1周期分と半端分を把握すれば、
 ((2^2-0)/4-2+1)+((3^2-1)/4-2+1)+((4^2-0)/4-2+1)+((5^2-1)/4-2+1)+((6^2-0)/4-2+1)  = \frac{(2^2+3^2+4^2+5^2+6^2+7^2)-((0+1)\times 2+1)}{4}-(2-1)\times(6-2+1)
とまとめることができ、実質ループ回数は周期1回分で済みます。( 平方和は、公式  \sum_{k=1}^{n} k^2 = n(n+1)(2n+1)/6 を活用 )

この考え方に基づいたコードがこちら。m=3000 に対して提出解が CodeIQの実行君で0.6秒程度に対し、こちらは0.1秒程度。…まあ、これはループ回数削減以外の要因も大きそうですが。

s=0
n=gets.to_i
m=n/2
x=n*(n+1)*(n*2+1)/6
1.upto(m){|i|
  u=v=0
  q,r=(n+i).divmod(i*2)
  if q>1
    if r<i
      c=-1
      r=i-r-1
    else
      c=1
      r-=i
    end
    1.upto(i-1){|j|
      u+=j*j%(i*4)
      v=u if j==r
    }
    v=v*c+(q-1)*(u*2+i%4*i)
  else
    r-=i
    1.upto(r){|j|
      v+=j*j%(i*4)
    }
  end
  s+=(x-i*(2*i-1)*(4*i-1)/3-v)/(i*4)
}
p s*2-((n*3-m*4)*m+1)*m/3

内(j)・外(i)のループ変数の変動幅をプロットしてみると、提出版では大きな三角形 ( ループ回数約 m2/4 ) だったのに対し、こちらでは小さな三角形 ( ループ回数約 m2/12 ) に範囲が狭まることになります。なお、i が大きくなると、j の範囲が狭まるため、そもそも1周期に達する前にループが終了してしまいます。 f:id:ange1:20151117013149p:plain

終わりに

まあ結局高速化しても、Ο(m2)なのは変わらないのですが。…周期を考える時に、素因数分解を使ってもっと周期を絞ればあるいは…? ( 大変ですが )。