2019年3月14日木曜日

ディオファントス方程式 $x^3+y^3+z^3=n$(2)

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

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

x^3+y^3+z^3=n が n=9k+4, 9k+5 では解を持たないことは簡単に示すことができる。xを9k+mとして9で割った余りをmとする。m=0,1,... 8である。このときx^3を9で割った余りは,m^3を9で割った余りになるので,(0,1,8,27,64,125,216,343,512) ≡ (0,1,8,0,1,8,0,1,8)。
そこで,この剰余(0,1,8)の3つの和のすべての組み合わせを9で割った余りの可能性を調べれば良い。(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)

0 件のコメント:

コメントを投稿