R中的三角精度答案

【问题标题】：Trigonometric Accuracy in RR中的三角精度
【发布时间】：2013-07-08 15:01:29
【问题描述】：

所以我遇到了以下问题：

“在抛物线上y=x2/k，选择了三个点A(a,a2/k),B(b,b2/k)和C(c,c2/k)。

令 F(K, X) 为整数四元组 (k, a, b, c) 的个数，使得三角形 ABC 的至少一个角为 45 度，且 1 ≤ k ≤ K 且 -X ≤ a

例如，F(1, 10) = 41 和 F(10, 100) = 12492。求 F(106, 109)。”

为了解决这个问题，我利用了点积的几何定义：theta = cos^-1((A dot B)/(|A|*|B|))，其中 A 和 B 是欧几里得向量, |A|表示A的大小，theta是它们之间的夹角。

我已多次阅读我的脚本，据我所知，它导致 FoKX=22 而不是 FoKX=41 的唯一原因是三角精度或从弧度到度的转换存在错误。让我知道是否是这种情况，或者我在某个地方犯了可能导致这种情况的错误。永远感谢您的帮助！

K<-1
X<-10
FoKX<-0
for(l in 1:K){
  for(i in (-X):(X-2)){
    for(j in (i+1):(X-1)){
      for(k in (j+1):X){
        vecAB<-c(j-i,(j^2-i^2)/l)
        vecAC<-c(k-i,(k^2-i^2)/l)
        vecBA<--vecAB
        vecBC<-c(k-j,(k^2-j^2)/l)
        vecCA<--vecAC
        vecCB<--vecBC
        magAB<-sqrt(sum(vecAB^2))
        magAC<-sqrt(sum(vecAC^2))
        magBA<-magAB
        magBC<-sqrt(sum(vecBC^2))
        magCA<-magAC
        magCB<-magBC
        ABdotAC<-sum(vecAB*vecAC)
        BAdotBC<-sum(vecBA*vecBC)
        CAdotCB<-sum(vecCA*vecCB)
        angA<-acos(ABdotAC/(magAB*magAC))
        angB<-acos(BAdotBC/(magBA*magBC))
        angC<-acos(CAdotCB/(magCA*magCB))
        if(angA==pi/4||angB==pi/4||angC==pi/4){
          FoKX<-FoKX+1
        }
      }
    }
  }
}

【问题讨论】：

可能是浮点问题；见stackoverflow.com/questions/9508518/…
是的。感谢您提供此链接。
我不假思索地写了一个蛮力算法，计算 F(10^6,10^9) 需要很长时间。是的，我错过了那些胡萝卜！回到绘图板。

标签： r trigonometry dot-product

【解决方案1】：

不要与浮点数比较是否完全相等。始终包含模糊因素。

        ....
        ....
        if(abs(angA - pi/4) < 1e-9 ||
           abs(angB - pi/4) < 1e-9 || 
           abs(angC - pi/4) < 1e-9){
          FoKX<-FoKX+1
        }
      }
    }
  }
}

FoKX
[1] 41

【讨论】：

更准确的方法：if(isTRUE(all.equal(angA, pi/4))||isTRUE(all.equal(angB, pi/4))||isTRUE(all.equal(angC, pi/4)))。请参阅 Aaron 推荐的链接：stackoverflow.com/questions/9508518/。

【解决方案2】：

我发现我的这个玩具在处理浮点相等时很有帮助：（像啤酒一样免费，随意修改等）

approxeq <- function(x, y, tolerance = .Machine$double.eps ^ 0.5,...) {
  #input validation 
  if (length(x) != length(y)) warning('x,y lengths differ. Will recycle.')
  #don't care about dimensions so long as you're smart about inputs
  checkit <- abs(x-y) < tolerance
  return(invisible(checkit))
}

这会返回一个逻辑向量，不像内置的 all.equal 有自己的用途。

【讨论】：