「奇数和平方数」問題

「奇数和平方数」問題
http://riverplus.hatenablog.com/entry/2018/02/25/143351

自然数 n に対し、連続する n 以下の奇数の和が平方数（自然数の 2 乗で表される数）となるものを探します。
例えば n = 10 の場合、以下の 7 通りです。
項が 1 つだけのものもカウントしている点に注意して下さい。
　　1 ＝ 1^2

　　1＋3 ＝ 2^2

　　1＋3＋5 ＝ 3^2

　　9 ＝ 3^2

　　1＋3＋5＋7 ＝ 4^2

　　7＋9 ＝ 4^2

　　1＋3＋5＋7＋9 ＝ 5^2
このような表し方の数を F(n)と定義します。例えば F(10) = 7 です。
同様に、F(20)＝14, F(30)＝23, F(1000)＝1272 となることが確かめられます。
F(10^7) の値を求めて下さい。
=================================================

ピタゴラス数を数えることになる。

euclid = function(m, n)    {
  repeat {
    temp = n %% m
    if (temp == 0) {
      break
    }
    n = m
    m = temp
  }
  return(m)
}

f = function(N) {
  N = ceiling(N / 2)
  count = 0
  for (m in 2:floor(sqrt(N - 1))) {
    for (n in 1:m) {
      if (m %% 2 != n %% 2 && euclid(m, n) == 1) {
        c = m ^ 2 + n ^ 2
        if (c > N) {
          break
        }
        count = count + N %/% c
      }
    }
  }
  cat(N + count * 2)
}

f(10000000) # 27368012, 3.073s

=================================================

Python だと，若干速いが，アルゴリズムが同じなので本質的な差はない。

import scipy as sp
import math

def f(N):
  N = int(sp.ceil(N / 2))
  count = 0
  for m in range(1, int(sp.floor(sp.sqrt(N - 1)))):
    for n in range(m + 1, N):
      if m % 2 == n % 2:
        continue
      if math.gcd(m, n) == 1:
        c = m ** 2 + n ** 2
        if c > N:
          break
        count += N // c
  return N + 2 * count

import time
s = time.time()
print(f(10000000)) # 27368012, 1.421s
print(time.time() - s)

解説・別解
http://riverplus.hatenablog.com/entry/2018/03/19/220705