算額（その58）

算額（その58）

静岡県伊豆市 江川邸 享和2年(1802)9月
http://www.wasan.jp/sizuoka/egawa.html

円内に 9 個の円が含まれている。円の径を求めよ。

外円の半径を 1 として，図のように記号を定め，方程式を解く。

eq2, eq3 は自明で，そうすれば eq1 は筆算でも解ける。

using SymPy
@syms r::positive, x::positive, y::positive;
eq1 = x^2 + y^2 - 4r^2;
eq2 = x - y;
eq3 = x - (1 - r)/2;

res = solve([eq1, eq2, eq3], (r, x, y))

   1-element Vector{Tuple{Sym, Sym, Sym}}:
    (-1/7 + 2*sqrt(2)/7, 4/7 - sqrt(2)/7, 4/7 - sqrt(2)/7)

using Plots

function circle(ox, oy, r, color=:red; beginangle=0, endangle=360)
   θ = beginangle:0.1:endangle
   x = r.*cosd.(θ)
   y = r.*sind.(θ)
   plot!(ox .+ x, oy .+ y, color=color, linewidth=0.5)
end;

function point(x, y, string="", color=:green, position=:left, vertical=:top; mark=true)
   mark && scatter!([x], [y], color=color, markerstrokewidth=0)
   annotate!(x, y, text(string, 10, position, color, vertical))
end;

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   plot()
   (r, x, y) = (-1/7 + 2*sqrt(2)/7, 4/7 - sqrt(2)/7, 4/7 - sqrt(2)/7)
   println("r = $r;  x = $x;  y = $y")
   circle(0, 0, 1, :black)
   circle(0, 0, r, :magenta)
   circle(0, 1-r, r, :blue)
   circle(0, r-1, r, :blue)
   circle(x, y, r)
   circle(x, -y, r)
   circle(-x, y, r)
   circle(-x, -y, r)
   circle(1-r, 0, r, :green)
   circle(r-1, 0, r, :green)
   if more
       point(0, 1-r, " 1-r", :blue)
       point(1-r, 0, " 1-r", :blue)
       point(r, 0, " r", :magenta)
       point(x, y, "(x,y)", :magenta)
       vline!([0], color=:black, lw=0.5)
       hline!([0], color=:black, lw=0.5)
   end
end;

   r = 0.2612038749637415;  x = 0.36939806251812923;  y = 0.36939806251812923

算額（その57）

算額（その57）

静岡県伊豆市 江川邸 享和2年(1802)9月
http://www.wasan.jp/sizuoka/egawa.html

円内に 10 個の円が含まれている。円の径を求めよ。

外円の半径を 1 として，図のように記号を定め，方程式を解く。

using SymPy
@syms r::positive, x1::positive, y1::positive, x2::negative;
x2 = 2x1
eq2 = (x2 - x1)^2 + (y1 - r)^2 - 4r^2
eq3 = x1^2 + y1^2 - (1 - r)^2
eq4 = x2^2 + r^2 - (1 - r)^2;

solve([eq2, eq3, eq4], (r, x1, y1))

   1-element Vector{Tuple{Sym, Sym, Sym}}:
    (-1/7 + 2*sqrt(2)/7, sqrt(9/28 - sqrt(2)/7), -5/14 + 5*sqrt(2)/7)
using Plots

function circle(ox, oy, r, color=:red; beginangle=0, endangle=360)
   θ = beginangle:0.1:endangle
   x = r.*cosd.(θ)
   y = r.*sind.(θ)
   plot!(ox .+ x, oy .+ y, color=color, linewidth=0.5)
end;

function point(x, y, string="", color=:green, position=:left, vertical=:top; mark=true)
   mark && scatter!([x], [y], color=color, markerstrokewidth=0)
   annotate!(x, y, text(string, 10, position, color, vertical))
end;

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   plot()
   r, x1, y1 = (-1/7 + 2*sqrt(2)/7, sqrt(9/28 - sqrt(2)/7), -5/14 + 5*sqrt(2)/7)
   x2 = 2x1
   println("r = $r, x1 = $x1, y1 = $y1, x2 = $x2")
   circle(0, 0, 1, :black)
   circle(0, r, r, :blue)
   circle(0, -r, r, :blue)
   circle(x1, y1, r, :magenta)
   circle(x1, -y1, r, :magenta)
   circle(-x1, y1, r, :magenta)
   circle(-x1, -y1, r, :magenta)
   circle(x2, r, r, :brown)
   circle(x2, -r, r, :brown)
   circle(-x2, r, r, :brown)
   circle(-x2, -r, r, :brown)
   if more
       point(0, r, " r", :blue)
       point(x1, y1, "(x1,y1)", :magenta)
       point(x2, r, "(x2, r)\nx2=2x1", :brown)
       vline!([0], color=:black, lw=0.5)
       hline!([0], color=:black, lw=0.5)
   end
end;

   r = 0.2612038749637415, x1 = 0.3455402473202352, y1 = 0.6530096874093536, x2 = 0.6910804946404704

atan2(0, 0) はどのように評価されるか

IEEE754 - IEEE Standard for Floating-Point Arithmetic では、

atan2( 0, -0) = π
atan2(-0, -0) = -π

と定義されている。

R はこれに従っている。

> options(digits=16)
> atan2(0, 0)
[1] 0
> atan2(0, -0)
[1] 3.141592653589793
> atan2(-0, -0)
[1] -3.141592653589793

Octave でも R と同じである。

>> format long
>> atan2(0, 0)
ans = 0
>> atan2(0, -0)
ans = 3.141592653589793
>> atan2(-0, -0)
ans = -3.141592653589793

AWK/GAWK も R と同じである。

$ awk "BEGIN{print atan2(0, 0)}"
0
$ awk "BEGIN{print atan2(0, -0)}"
3.14159
$ awk "BEGIN{print atan2(-0, -0)}"
-3.14159

Julia は atan2(x, y) は atan(x, y) である。結果はいずれも 0.0 とされる。

julia> atan(0, 0)
0.0

julia> atan(0, -0)
0.0

julia> atan(-0, -0)
0.0

Python でも Julia と同じく，結果はいずれも 0.0 とされる。

>>> import math
>>> math.atan2(0, 0)
0.0
>>> math.atan2(0, -0)
0.0
>>> math.atan2(-0, -0)
0.0

C++ (Apple clang version 14.0.0 (clang-1400.0.29.202)) でも Julia と同じく，結果はいずれも 0.0 とされる。

#include <stdio.h>
#include <math.h>

int main()
{
    printf("%f\n", atan2(0, 0));
    printf("%f\n", atan2(0, -0));
    printf("%f\n", atan2(-0, -0));
    return 0;
}

$ ./a.out
0.000000
0.000000
0.000000

WolframAlpha では未定義とされる。

atan2(0, 0), atan2(0, -0), atan2(-0, -0)

いずれも，結果は "(未定義)" と表示される