ディオファントス方程式 $x^3+y^3+z^3=n$（２）

（ディオファントス方程式 $x^3+y^3+z^3=n$（１）からの続き）

というわけで，今日，３月14日は円周率の日なのかホワイトデイなのかアインシュタインの誕生日なのかホーキングの命日なのか巷では協議が続いている。

x^3+y^3+z^3=n が n=9k+4, 9k+5 では解を持たないことは簡単に示すことができる。xを9k+mとして９で割った余りをmとする。m=0,1,... 8である。このときx^3を９で割った余りは，m^3を９で割った余りになるので，(0,1,8,27,64,125,216,343,512) ≡ (0,1,8,0,1,8,0,1,8)。

そこで，この剰余(0,1,8)の３つの和のすべての組み合わせを９で割った余りの可能性を調べれば良い。(0,0,0)≡0, (0,0,1)≡1, (0,0,8)≡8, (0,1,1)≡2, (0,1,8)≡0, (0,8,8)≡7, (1,1,1)≡3, (1,1,8)≡1, (1,8,8)≡8, (8,8,8)≡6であり，ここから4と5は出てこない。

少しだけ前回のプログラムを高速化してみる。k^3の配列の次元をmとして，前回のものはO(m^3)，今回のものはO(m^2)である。

（注）Juliaは0の配列添字を許さないので，この度もx,y,zから0は除いて考えている。0を含めれば，91=0^3+3^3+4^3 など簡単なものがあると思うけれど。

- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
function pre(m,t)
  for l in 1:m
    t[l]=l^3
  end
end

function cinv(x)
  return Int(ceil(cbrt(abs(x))))
end

function dpe(n,m,t)
  l=0
  if(n%9 ==4 || n%9 ==5)
    return (n, n%9, 0)
  end
  for i in 1:m
    d=cinv(n-t[i])
    if(n>=t[i])
      for j in 1:d
        k=cinv(n-t[i]-t[j])
        l=l+1
        if k !=0 && (t[i]+t[j]+t[k]==n || t[i]+t[j]+t[k]==-n)
          return (n,l)
        end
      end
      for j in d+1:i
        k=cinv(n-t[i]-t[j])
        l=l+1
        if k !=0 && (t[i]+t[j]-t[k]==n || t[i]+t[j]-t[k]==-n)
          return (n,l)
        end
      end
    else
      for j in d+1:i
        k=cinv(t[i]-n-t[j])
        l=l+1
        if k !=0 && (t[i]-t[j]+t[k]==n || t[i]-t[j]+t[k]==-n)
          return (n,l)
        end
      end
      for j in 1:d
        k=cinv(t[i]-n-t[j])
        l=l+1
        if k !=0 && (t[i]-t[j]-t[k]==n || t[i]-t[j]-t[k]==-n)
          return (n,l)
        end
      end
    end
  end
  return (m,l)
end

function pst(l,m,t)
  for n in 1:l
    println(dpe(n,m,t))
  end
end

m=10000
t=ones(Int64,m)
pre(m,t)
@time pst(99,m,t)
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
(1, 1)
(2, 26)
(3, 1)
(4, 4, 0)
(5, 5, 0)
(6, 1)
(7, 5491)
(8, 1)
(9, 23487)
(10, 1)
(11, 4)
(12, 61)
(13, 4, 0)
(14, 5, 0)
(15, 1057)
(16, 100)
(17, 1)
(18, 4)
(19, 185)
(20, 1561)
(21, 131)
(22, 4, 0)
(23, 5, 0)
(24, 2)
(25, 2)
(26, 48674)
(27, 1)
(28, 146)
(29, 1)
(10000, 50004997)
(31, 4, 0)
(32, 5, 0)
(10000, 50004997)
(34, 10)
(35, 96)
(36, 1)
(37, 1574)
(38, 364)
(10000, 50004997)
(40, 4, 0)
(41, 5, 0)
(10000, 50004997)
(43, 2)
(44, 30)
(45, 152329)
(46, 422)
(47, 31)
(48, 6)
(49, 4, 0)
(50, 5, 0)
(51, 317009)
(10000, 50004997)
(53, 9)
(54, 70)
(55, 1)
(56, 240)
(57, 726)
(58, 4, 0)
(59, 5, 0)
(60, 9)
(61, 466761)
(62, 2)
(63, 23)
(64, 5)
(65, 6188)
(66, 5)
(67, 4, 0)
(68, 5, 0)
(69, 342)
(70, 219)
(71, 12)
(72, 50)
(73, 5)
(10000, 50004998)
(10000, 50004998)
(76, 4, 0)
(77, 5, 0)
(78, 1509)
(79, 612)
(80, 6)
(81, 3)
(82, 96)
(83, 11)
(10000, 50004994)
(85, 4, 0)
(86, 5, 0)
(87, 9120551)
(88, 13)
(89, 884141)
(90, 7)
(91, 72576)
(92, 1)
(93, 20)
(94, 4, 0)
(95, 5, 0)
(96, 239)
(97, 6)
(98, 108)
(99, 2)
 14.962059 seconds (103.76 k allocations: 4.814 MiB)