20/08/30 11:31:30.73 hfh+x6my.net
円周角の一致でなくて半径と中心が一致することで同一円と判定するようにアルゴリズムを変更
N=7
gr=expand.grid(1:N,1:N)
(node=mapply(function(x,y) x+1i*y,gr[,1],gr[,2]))
tric <- function(A,B,C){ # 複素点３点を通る円の中心と半径を返す
a1=Re(A) ; a2=Im(A)
b1=Re(B) ; b2=Im(B)
c1=Re(C) ; c2=Im(C)
p = (a1^2*(-b2) + a1^2*c2 - a2^2*b2 + a2^2*c2 + a2*b1^2 + a2*b2^2 - a2*c1^2 - a2*c2^2 - b1^2*c2 - b2^2*c2 + b2*c1^2 + b2*c2^2)/(2*(-a1*b2 + a1*c2 + a2*b1 - a2*c1 - b1*c2 + b2*c1))
q = -(a1^2*(-b1) + a1^2*c1 + a1*b1^2 + a1*b2^2 - a1*c1^2 - a1*c2^2 - a2^2*b1 + a2^2*c1 - b1^2*c1 + b1*c1^2 + b1*c2^2 - b2^2*c1)/(2*(-a1*b2 + a1*c2 + a2*b1 - a2*c1 - b1*c2 + b2*c1))
Ce=p+1i*q
r=abs(Ce-A)
c(Center=Ce,Radius=r)
}
onCir <- function(x){ # 中心と半径が一致するかを返す
all(tric(node[x[1]],node[x[2]],node[x[3]])==
tric(node[x[1]],node[x[2]],node[x[4]]))
}
sum(combn(N^2,4,onCir),na.rm=TRUE)
７×７個だと
> sum(combn(N^2,4,onCir),na.rm=TRUE)
[1] 5704