算額（その1457）

算額（その1457）

福島県耶麻郡猪苗代町中小松西浜 猪苗代小平潟天満宮 明治14年(1881)
http://www.wasan.jp/fukusima/kohiragata3.html

街角の数学 Street Wasan　～落書き帳「○△□」～ 329．小平潟天満宮（その２）
http://streetwasan.web.fc2.com/math17.11.1.html
キーワード：中国剰余定理
#Julia, #SymPy, #算額, #和算

今神前で子供が字を書いている。毎日書く総字数 M は決まっている。
　　初日は，紙一枚毎に３字ずつ書いていくと，最後に２字余った。
　　次の日は，紙一枚毎に５字ずつ書いていくと，４字余った。
　　三日目は，紙一枚毎に７字ずつ書いていくと，６字余った。
　　さらに四日目最終日は，11字ずつ書いたら10字余ったという。
　　一日に書く総字数 M を求めよ。

中国剰余定理を使えば一発で答えが出る。

Julia の SymPy の元は Python の sympy である。crt はインポートされていないので，明示的に sympy.ntheory.modular.crt として呼び出す。

using SymPy

#from sympy.ntheory.modular import crt

# 剰余の条件
moduli = [3, 5, 7, 11]   # 除数のセット
remainders = [2, 4, 6, 10]  # 剰余のセット

# 中国剰余定理を使用して解を求める
solution = sympy.ntheory.modular.crt(moduli, remainders)  # 最小解と周期が返される

    (1154, 1155)

文字の数は 1154 文字である。

検算で正しい解であることがわかる。

mod.(1154, moduli) |> println

    [2, 4, 6, 10]

ブルートフォースでやってみる。16 ミリ秒で計算できた。

@time begin
for n = 1:1000000000
    all(mod.(n, [3, 5, 7, 11]) .== [2, 4, 6, 10]) && (println(n); break)
end
end

    1154
      0.014902 seconds (22.13 k allocations: 1.218 MiB, 97.79% compilation time)

 

算額（その1456）

算額（その1456）

五輪教一：黄金比の眠るほこら，日本評論社，2015年7月10日
キーワード：整数方程式
#Julia, #SymPy, #算額, #和算

今，神前に 328 円ある。これで，一等馬，二等馬，三等馬を買い入れたい。1等馬の頭数の 3/4 は二等馬の頭数に等しく，二等馬の頭数の 2/3 は三等馬の頭数に等しい。また，二等馬の単価は一等馬の単価より 2 円安く，三等馬の単価は二等馬の単価より 3 円安い。一等馬，二等馬，三等馬各々の頭数と価格を求めよ。

整数方程式に持ち込み，あれこれやるが，ブルートフォースでは簡単である。

一等馬の頭数を x，単価を y とする。

# Julia の Base だけで
for x in 1:10
    for y in 5:30
        x*y + 3x*(y - 2)//4 + x*(y - 5)//2 == 328 && println("$x, $y")
    end
end

    8, 20

一等馬は 8 頭，単価は 20 円。
二等馬は 8*3/4 = 6 頭，単価は 20 - 2 = 18 円。
三等馬は 6*2/3 = 4 頭，単価は 18 - 3 = 15 円。

検算 20*8 + 18*6 + 15*4 = 328

算額（その1455）

算額（その1455）

福島県南相馬市 旧小高町 小高神社
街角の数学 Street Wasan ～ 落書き帳「○△□」～ 406．小高神社
http://streetwasan.web.fc2.com/math18.3.11.html

福島県田村市船引町文珠 安倍文殊菩薩堂 明治10年(1877)
街角の数学 Street Wasan ～ 落書き帳「○△□」～ 333．参詣人の数
http://streetwasan.web.fc2.com/math17.11.9.html

五輪教一：黄金比の眠るほこら，日本評論社，2015年7月10日

キーワード：整数方程式
#Julia, #SymPy, #算額, #和算

小高神社の祭礼で，騎馬の出陣数の参加人員はわからないが，その 7/9 の下 2 桁は 68，また，5/8 の下 2 桁は 60 である。総人数は何人か。

整数方程式に持ち込み，あれこれやるが，ブルートフォースでは簡単である。

# Julia の Base だけで
for x in 1:10000
    mod(x*7//9, 100) == 68 && mod(x*5//8, 100) == 60 && (println(x); break)
end

    2016

総勢 2016 人である。

算額（その1454）

算額（その1454）

山形県高畠町亀岡 亀岡文殊（松髙山大聖寺文殊堂）
五輪教一：黄金比の眠るほこら，日本評論社，2015年7月10日
キーワード：中国剰余定理
#Julia, #SymPy, #算額, #和算

ある同好会の一行が亀岡文殊を参拝した。総人数を x としたとき，x を 5 で割れば 1 余り，7 で割れば 2 余り，9 で割れば 3 余り，11 で割れば 4 余り，13 で割れば 5 余り，17 で割れば 6 余る。総人数は何人か。

中国剰余定理を使えば一発で答えが出る。

Julia の SymPy の元は Python の sympy である。crt はインポートされていないので，明示的に sympy.ntheory.modular.crt として呼び出す。

using SymPy

#from sympy.ntheory.modular import crt

# 剰余の条件
moduli = [5, 7, 9, 11, 13, 17]   # 除数のセット
remainders = [1, 2, 3, 4, 5, 6]  # 剰余のセット

# 中国剰余定理を使用して解を求める
solution = sympy.ntheory.modular.crt(moduli, remainders)  # 最小解と周期が返される

    (698196, 765765)

総勢 698196 人である。

検算で正しい解であることがわかる。

mod.(698196, moduli) |> println

    [1, 2, 3, 4, 5, 6]

ブルートフォースでやってみる。0.1 秒で計算できた。

@time begin
for n = 1:1000000000
    all(mod.(n, [5, 7, 9, 11, 13, 17]) .== [1, 2, 3, 4, 5, 6]) && (println(n); break)
end
end

    698196
      0.105789 seconds (4.90 M allocations: 213.974 MiB, 23.61% gc time, 14.33% compilation time)

 

算額（その1453）

算額（その1453）

福島県三春町 田村大元神社 文政12年(1829)
令和5年(2023)2月田村高校生徒により復元
http://www.wasan.jp/fukusima/tamuradaigen2.html

街角の数学 Street Wasan　～落書き帳「○△□」～ 257．花見の座興
http://streetwasan.web.fc2.com/math17.4.25.html

キーワード：中国剰余定理
#Julia, #SymPy, #算額, #和算

満開の桜の下で，発句詩歌の会が花見をしている。咲き誇る花数を文字数に見なして，これを17文字で区切れば3文字余り，七言絶句28文字で区切れば 23 文字不足し（＝5文字余り），歌31文字で区切れば8文字余る，と話している。 さて，会員たちは桜の花数を何個と見たてたのだろう。

中国剰余定理を使えば一発で答えが出る。

Julia の SymPy の元は Python の sympy である。crt はインポートされていないので，明示的に sympy.ntheory.modular.crt として呼び出す。

using SymPy

#from sympy.ntheory.modular import crt

# 剰余の条件
moduli = [17, 28, 31]   # 除数のセット
remainders = [3, 5, 8]  # 剰余のセット

# 中国剰余定理を使用して解を求める
solution = sympy.ntheory.modular.crt(moduli, remainders)  # 最小解と周期が返される

    (7789, 14756)

花の数は 7789 個である。

検算で正しい解であることがわかる。

mod.(7789, moduli) |> println

    [3, 5, 8]

ブルートフォースでやってみる。1 ミリ秒で計算できた。

@time begin
for n = 1:1000000000
    all(mod.(n, [17, 28, 31]) .== [3, 5, 8]) && (println(n); break)
end
end

    7789
      0.015988 seconds (68.58 k allocations: 2.803 MiB, 94.85% compilation time)

 

セルフうどん メリケンや 伏石店

かつて，宇高連絡船のデッキでさぬきうどんを食べることができた。
宇高連絡船が廃止されてからも，高松駅構内で「連絡船うどん」を食べることができた。
高松駅周辺の再開発で，伏石へ移転した。店舗の一角にその名残がある。
頼んだのは，しっぽくうどん

算額（その1452）

算額（その1452）

福島県三春町 龍穏院 明治26年(1893)
http://www.wasan.jp/fukusima/ryuonin.html

街角の数学 Street Wasan　～落書き帳「○△□」～ 216. 二円の縁
http://streetwasan.web.fc2.com/math17.2.1.html
街角の数学 Street Wasan　～落書き帳「○△□」～ 217. 二円の縁（その2）
http://streetwasan.web.fc2.com/math17.2.2.html

キーワード：円3個，菱形，最大値
#Julia, #SymPy, #算額, #和算

大円 2 個が交わっており，中に小円と菱形が入っている。大円の直径が与えられたとき，黒積が最大となるときの菱形の対角線の長い方の長さはいかほどか。

大円の半径と中心座標を r1, (r1 - r2, 0)
小円の半径と中心座標を r2, (0, 0)

黒積は，灰色部分 = 扇形ABC - 三角形OAB - 四分円OCD の4倍である。
∠OAB = θ，θ = asin(sqrt(r2*(2r1 - r2))/r1)

黒積は小円の大きさ（半径）によって変化する。まず黒積が最大になるときの小円の半径を求める。
菱形の対角線の長い方の長さは 4r1 - 2r2 である。

include("julia-source.txt");
# julia-source.txt ソース　https://blog.goo.ne.jp/r-de-r/e/ad3a427b84bb416c4f5b73089ae813cf

using SymPy
@syms r1, r2, y, θ, 円積率
円積率 = 0.79
円周率 = 4円積率
#円周率 = PI
y = sqrt(2r2*r1 - r2^2)  # sqrt(r1^2 - (r1 - r2)^2)
θ = asin(y/r1)
S = 4(円周率*r1^2*(θ/2PI) - y*(r1 - r2)/2 - 円周率*r2^2/4)
S |> println

    6.32*r1^2*asin(sqrt(2*r1*r2 - r2^2)/r1)/pi - 3.16*r2^2 - 2*(r1 - r2)*sqrt(2*r1*r2 - r2^2)

pyplot(size=(300, 150), grid=false, aspectratio=:none, label="")
plot(S(r1 => 1/2), xlims=(0, 0.5), xlabel="r2", ylabel="黒積")
vline!([0.28746844208856126])

g = diff(S, r2)
g |> println

    6.32*r1*(r1 - r2)/(pi*sqrt(1 - (2*r1*r2 - r2^2)/r1^2)*sqrt(2*r1*r2 - r2^2)) - 6.32*r2 - 2*(r1 - r2)^2/sqrt(2*r1*r2 - r2^2) + 2*sqrt(2*r1*r2 - r2^2)

plot(g(r1 => 1/2), xlims=(0, 0.35), xlabel="r2", ylabel="黒積の導関数")
hline!([0])

solve(g(r1 => 1/2), r2)[1] |> println

    0.287468442088561

S(r1 => 1/2, r2 => 0.287468442088561).evalf() |> println

    0.115685731044897

大円の半径が r1 = 1/2 のときに，小円の半径が r2 = 0.287468442088561 のとき，黒積は最大値 0.115685731044897 をとる。

菱長 = (4r1 - 2r2)(r1 => 1/2, r2 => 0.287468442088561) |> println

    1.42506311582288

大円の半径が r1 = 1/2 のときに，小円の半径が r2 = 0.287468442088561 のとき，菱長は 1.42506311582288 である。

術は 2r1/(0.125/0.79^2 + 0.5) で，r1 = 1/2 のとき 1.42798306829882 となるが，違いがどこから来るのかわからない。

(2r1/(0.125/0.79^2 + 0.5))(r1 => 1/2) |> println

    1.42798306829882

function draw(r1, r2, more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
    y = sqrt(r1^2 - (r1 - r2)^2)
    θ = asind(y/r1)
    println("θ = $θ")
    θ1 = 180 - θ:0.1:180
    xa = @. r1*cosd(θ1) + (r1 - r2)
    ya = @. r1*sind(θ1)
    θ2 = 180:-0.1:90
    append!(xa, r2.*cosd.(θ2))
    append!(ya, r2.*sind.(θ2))
    append!(xa, xa[1])
    append!(ya, ya[1])
    plot([2r1 - r2, 0, r2 - 2r1, 0, 2r1 - r2], [0, y, 0, -y, 0], color=:green, lw=0.5)
    plot!(xa, ya, color=:orange, lw=0.5, seriestype=:shape, fillcolor=:gray70)
    circle2(r1 - r2, 0, r1)
    circle(0, 0, r2, :blue)
    segment(0, y, r1 - r2, 0, :gray70)
    if more
        delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
        hline!([0], color=:gray80, lw=0.5)
        vline!([0], color=:gray80, lw=0.5)
        point(0, 0, "O", :blue, :center, delta=-delta)
        point(0, r2, "D", :blue, :center, delta=-delta)
        point(-r2, 0, "C ", :blue, :right, :vcenter)
        point(0, sqrt(r2*(2r1 - r2)), "  B:√(r2*(2r1-r2))", :red, :left, :vcenter)
        point(r1 - r2, 0, " A:r1-r2 ", :red, :left, :vcenter)
        point(2r1 - r2, 0, "2r1-r2  ", :green, :right, :vcenter)
        point(r2, 0, " r2", :blue, :left, :vcenter)
    end
end;

draw(1/2, 0.15, true)

算額（その1451）

算額（その1451）

福島県三春町 津島神社 明治16年(1883)
http://www.wasan.jp/fukusima/tusima.html

和算で遊ぼう！！ 「三春まちなか寺子屋」2017レポート
6 月 津島神社
https://miharukoma.com/wp-content/uploads/2018/01/%E4%B8%89%E6%98%A5%E3%81%BE%E3%81%A1%E3%81%AA%E3%81%8B%E5%AF%BA%E5%AD%90%E5%B1%8B2017%E3%83%AC%E3%83%9D%E3%83%BC%E3%83%88.pdf

キーワード：数列，部分和
#Julia, #SymPy, #算額, #和算

猟船を所有する船主がいる。初日は 1 艘、二日目は 5 艘、三日目は 9 艘。
逐次、1 から始まって一つ置きの奇数艘ずつ、毎日船を出した。こうして、7 日間ですべての船が海 へ出て行ったという。この船主の所有する船は、全部で何艘か。

問題は簡単であり，この場合は和をとりながら列挙していけばよいが，一般解を SymPy を使って求めることにする。

1. SymPy の sequence を使うためには，Julia では sympy.sequence としなければならない。
2. 項は 0 から始まることに注意すること。

まずは，一般項を表す式を定義しなければならない。
今の場合は簡単であるが，ChatGPT に聞くか，「オンライン整数列大辞典 https://oeis.org/?language=japanese」などを参照すればよい。

using SymPy
@syms x, n
a = sympy.sequence(1 + 4x)  # x は 0 から始まることに注意

    [1, 5, 9, 13, ...]

部分和を計算するには summation を使う。
最初と最後を指定するとき，添字は 0 から始まることに注意する。算額の場合は，1 日目から 7 日目までなので，第 0 項から 第 6 項までである。
念のために分かっている最初の数項についての結果を出してみるとよい。

算額の場合の和は 91 である。「答え：91 艘」

s = summation(a.formula, (x, 0, 6)) # x は 0 から始まることに注意
s |> println

    91

一般項は，第 0 項から 第 n - 1 項まで指定すればよい。
summation の場合と不一致のような気がすると思うが，最後の指定は n - 1 ではなく，n までなので注意が必要である。つまり，(x, 0, n) である。
これも，念のために分かっている最初の数項についての結果を比較してみるとよい。

s = summation(a.formula, (x, 0, n)) # x は 0 から始まることに注意
s |> factor |> println

    (n + 1)*(2*n + 1)

f(n) = (n + 1)*(2*n + 1)
f(0), f(1), f(2), f(3), f(4), f(5), f(6), f(7)

    (1, 6, 15, 28, 45, 66, 91, 120)

項が 1 から始まる場合の一般項は n*(2*n - 1) である。
術で述べられている「(7 × 2 − 1) × 7」である。

これはsummation で得られる一般項の n に (n - 1)を代入すれば得られる。

s(n => n - 1) |> factor |> println

    n*(2*n - 1)

f2(n) = n*(2*n - 1)
f2(1), f2(2), f2(3), f2(4), f2(5), f2(6), f2(7), f2(8)

    (1, 6, 15, 28, 45, 66, 91, 120)

以上のように若干ややこしいので，力ずくで解を求める方法（ブルートフォース）も最初の試みとしては有望である。

最近のコンピュータは優れているので，単純なプログラムを書いても短時間で解を求めることができる場合もある。
「10000000日で全部出航した」のような場合でも，以下のプログラムは 1 秒以内で答えを出す。

@time f2(10000000)

      0.000000 seconds

    199999990000000

@time begin
s = 0
days = 0
for i = 1:4:9999999999
    days += 1
    s += i
    if days == 10000000
        println("$days 日目までに $s 艘が出航した")
        break
    end
end
end

    10000000 日目までに 199999990000000 艘が出航した
      0.638406 seconds (30.02 M allocations: 459.081 MiB, 4.23% gc time, 1.75% compilation time)