以前書いたコードで次の様にしていた.

デバッグ版では単位ベクトルにすると誤差が大きい様に感じたが,リリース版ではそれ程でもない?
double Naiseki (const P3& p1,const P3& p2)
{
// → → → →
// a・b=|a||b|cosθ
// a.x*b.x+a.y*b.y+a.z*b.z=sqrt(ax*ax+ay*ay+az*az)*sqrt(bx*bx+by*by+bz*bz)*cosθ
// ... a,bのベクトルが単位ベクトルの時
// sqrt(ax...) と sqrt(bx...) は 1
// a.x*b.x+a.y*b.y+a.z*b.z=1*1*cosθ
// a.x*b.x+a.y*b.y+a.z*b.z=cosθ
// → →
// a・b=ax*bx+ay*by+az*bz
P3 p1u = p1.Uni() ;
P3 p2u = p2.Uni() ;
// return (p1u.x*p2u.x+p1u.y*p2u.y+p1u.z*p2u.z) ;
// 誤差の補正
double dp = p1u.x*p2u.x+p1u.y*p2u.y+p1u.z*p2u.z ;
if (dp<-1 || 1<dp) {
if (dp < -1.) { dp = -1. ; }
if (dp > 1.) { dp = 1. ; }
}
return dp ;
}
P3 Gaiseki (const P3& p1,const P3& p2)
{
// → → → →
// a×b=|a||b|sinθ・c
// → →
// a×b=(ay*bz-az*by,az*bx-ax*bz,ax*by-ay*bx)
P3 pt ;
pt.x = p1.y * p2.z - p1.z * p2.y ;
pt.y = p1.z * p2.x - p1.x * p2.z ;
pt.z = p1.x * p2.y - p1.y * p2.x ;
return pt ;
}