Julia: sin33°、cos33°、tan33°はどんな数？ 再び

sin33°、cos33°、tan33°はどんな数？
https://p-suugaku.blogspot.com/2022/01/sin33cos33tan33.html

については，

Julia/SymPy: sin(33°)，cos(33°)，tan(33°)
https://blog.goo.ne.jp/r-de-r/e/5aa2b18080d95d0e6b9213493d8bceb9

に，すでに書いたけど，tan(33°) の simplify をやった結果も書いておこう。

元ページには写し間違いとは思うが，間違いがある。

手入力は間違いの元，コピペするに限る。検算すればいいだけの話です。

using SymPy
@syms a

##### sin(33°)

sind(33) # 0.5446390350150271

s1 = sind(Sym(33)) |> simplify |> together

s2 = s1(√(5-√Sym(5))=>a)

s3 = s2.coeff(a, 1) |> together |> numer

s4 = 2√(expand(s3^2 * (5-√Sym(5)))/4)/16 + s2.coeff(a, 0)
s4.evalf() # 0.544639035015027
sind(33)   # 0.5446390350150271
s4 |> print
# -sqrt(6)/16 - sqrt(2)/16 + sqrt(10)/16 +
# sqrt(-10*sqrt(3) - 2*sqrt(15) + 4*sqrt(5) + 20)/8 + sqrt(30)/16

(a) 式の3行目10√10 は 10√3
(√6 - √2)/4 * √(10 + 2√5)/4 + (√6 + √2)/4 * (√5 - 1)/4        # 0.544639035015027
√2*((√3 - 1) * √(10 + 2√5) + (√3 + 1) * (√5 - 1))/16          # 0.5446390350150271
# (2*√(20 - 10√10 +4√5  - 2√15) - √2 - √6 + √10 + √30) / 16   # error!
  (2*√(20 - 10√3  + 4√5 - 2√15) - √2 - √6 + √10 + √30) / 16   # 0.5446390350150272

##### cos(33°)

cosd(33) # 0.838670567945424

c1 = cosd(Sym(33)) |> simplify |> together

c2 = c1(√(5-√Sym(5))=>a)

c3 = c2.coeff(a, 1) |> together |> numer

c4 = 2√(expand(c3^2 * (5-√Sym(5)))/4)/16 + c2.coeff(a, 0)

c4.evalf() # 0.838670567945424
cosd(33) # 0.838670567945424
c4 |> print
# -sqrt(30)/16 - sqrt(2)/16 + sqrt(6)/16 + sqrt(10)/16 +
#  sqrt(2*sqrt(15) + 4*sqrt(5) + 10*sqrt(3) + 20)/8

(b) 式の3行目の 10√10 も 10√3 の間違い
(√6 + √2)/4 * √(10 + 2√5)/4 - (√6 - √2)/4 * (√5 - 1)/4      # 0.8386705679454239
√2*((√3 + 1) * √(10 + 2√5) - (√3 - 1) * (√5 - 1))/16        # 0.838670567945424
# (2*√(20 + 10√10 + 4√5 + 2√15) - √2 + √6 + √10 - √30) / 16 # error!
  (2*√(20 + 10√3  + 4√5 + 2√15) - √2 + √6 + √10 - √30) / 16 # 0.8386705679454238

##### tan(33°)

tand(33) # 0.6494075931975105

t1 = tand(Sym(33)) |> simplify |> together

numerator = numer(t1)

denominator = denom(t1)

denominator0 = denominator(√Sym(6)=>-√Sym(6))

denominator2 = expand(denominator * denominator0)

denominator00 = denominator2(8=>-8)

denominator = expand(denominator2 * denominator00) # 256
numerator = expand(numerator * denominator0 * denominator00)

t2 = numerator / denominator |> simplify

t3 = -√(50 - 10√Sym(5))/2 - √(10 - 2√Sym(5)) + √(150 - 30√Sym(5))/4 + 3√(30 - 6√Sym(5))/4t3.evalf()

t4 = - √Sym(5) - 2 + √Sym(15)/2 + 3√Sym(3)/2

y = (5 - √Sym(5))*((-2√Sym(10) + √Sym(30) - 4√Sym(2) + 3√Sym(6))/4)^2

expand(y*4)

t6 = √expand(y*4)/2 + t4

t6.evalf() # 0.649407593197511

tand(33) # 0.6494075931975105

Julia: プログラムを速くするための定石--その6

Julia: プログラムを速くするための定石
https://blog.goo.ne.jp/r-de-r/e/2e21bccd94b2f24cdbc2b54030b13a1a

Julia: プログラムを速くするための定石--その2
https://blog.goo.ne.jp/r-de-r/e/4a815e35e7364368ed963e75a966b595

Julia: プログラムを速くするための定石--その3
https://blog.goo.ne.jp/r-de-r/e/bb84137b4936fb4ab5e76a3cbd7ce9e5

Julia: プログラムを速くするための定石--その4
https://blog.goo.ne.jp/r-de-r/e/6b538afd6188153bf5e7db614b67f5e5

Julia: プログラムを速くするための定石--その5 の続きである
https://blog.goo.ne.jp/r-de-r/e/2aad94a1df626356ffab159eb16f00f4

で，今回は，

素数の末尾の桁の続き
http://commutative.world.coocan.jp/blog5/2021/12/post-507.html

である。

速くするための定石（関数化する，ファイルは読まないなど）は同じである。
しかしその他に，このプログラムにおいては，速くするというのとは直接関係がないが，繰り返し記述を避けるということを挙げておく。
カウントのために n11~n99 という16個の変数を使っているが，その後の 16 段にわたる if-elseif での記述や 16 行の print/println を見ると，二次元配列を使うのがまっとうだとわかる。

なにはともあれ，元のプログラムの出力結果と計算時間は以下のとおりである。

11: 21423878,13: 33479893,17: 33692895,19: 25164853
31: 27517671,33: 20741925,37: 31818060,39: 33687968
71: 28959412,73: 30589282,77: 20731084,79: 33484260
91: 35860558,93: 28954524,97: 27521999,99: 21424245

177.184738 seconds (1.82 G allocations: 33.924 GiB, 1.34% gc time, 0.02% compilation time)

元の 68 行のプログラムを，以下の 16 行のプログラムのように書き換えた。
短いプログラムはそれだけで，見通しが良い。

using Primes, Printf
function g(n=10^10)
  cnt = zeros(Int, 9, 9)
  a = 2
  for p in primes(3, n)
    b = p % 10
    cnt[a, b] += 1
    a = b
  end
  for i in [1, 3, 7, 9]
    for j in [1, 3, 7, 9]
      @printf("  %d%d: %7d", i, j, cnt[i, j])
    end
    println()
  end
end

@time g(10^10)

結果は以下のとおり。

11: 21423878  13: 33479893  17: 33692895  19: 25164853
31: 27517671  33: 20741925  37: 31818060  39: 33687968
71: 28959412  73: 30589282  77: 20731084  79: 33484260
91: 35860558  93: 28954524  97: 27521999  99: 21424245

26.982763 seconds (334 allocations: 5.885 GiB, 0.56% gc time)

177.184738/26.982763 ≒ 6.5 倍速

なお，上の結果は，10^10 までの素数(最初の 455052511 個の素数)，元ページに書かれている結果は最初の 10^8 個の素数についての結果。

@time g(2038074743)

11: 4623041  13: 7429438  17: 7504612  19: 5442344
31: 6010981  33: 4442561  37: 7043695  39: 7502896
71: 6373982  73: 6755195  77: 4439355  79: 7431870
91: 7991431  93: 6372940  97: 6012739  99: 4622916

4.755174 seconds (334 allocations: 1.254 GiB, 3.46% gc time)