Trigonometric Accuracy in R

Question

So I am presented with the following problem:

"On the parabola y = x2/k, three points A(a, a2/k), B(b, b2/k) and C(c, c2/k) are chosen.

Let F(K, X) be the number of the integer quadruplets (k, a, b, c) such that at least one angle of the triangle ABC is 45-degree, with 1 ≤ k ≤ K and -X ≤ a < b < c ≤ X.

For example, F(1, 10) = 41 and F(10, 100) = 12492. Find F(106, 109)."

In aims to solve it I have exploited the geometric definition of the dot product: theta = cos^-1((A dot B)/(|A|*|B|)), where A and B are Euclidean vectors, |A| represents the magnitude of A, and theta is the angle between them.

I have read over my script multiple times and as far as I can see the only reason it is resulting in FoKX=22 instead of FoKX=41 is that there is an error in the trigonometic accuracy or conversion from radians to degrees. Let me know if this is the case or I have made a mistake somewhere that might account for this. Thanks always for the help!

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
        }
      }
    }
  }
}

It could be a floating point issue; see http://stackoverflow.com/questions/9508518/why-are-these-numbers-not-equal — Aaron left Stack Overflow, Jul 08 '13 at 15:05
And without thinking about it I've written a brute force algorithm that will take too long to compute F(10^6,10^9). Yep I was missing those carots! Back to the drawing board. — fowlslegs, Jul 09 '13 at 18:59

score 4 · Accepted Answer · answered Jul 08 '13 at 15:12

4

Don't compare for exact equality with floating-point. Always include a fuzz factor.

        ....
        ....
        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

answered Jul 08 '13 at 15:12

Hong Ooi

56,353
13
134
187

A more accurate way to do this: `if(isTRUE(all.equal(angA, pi/4))||isTRUE(all.equal(angB, pi/4))||isTRUE(all.equal(angC, pi/4)))`. See the link Aaron recommended: stackoverflow.com/questions/9508518/ . – fowlslegs Jul 09 '13 at 19:00

score 2 · Answer 2 · edited Jul 08 '13 at 17:04

I find this toy of mine helpful when dealing with floatingpoint equality: (free as in beer, free to modify as you like, etc)

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))
}

This returns a logical vector, unlike the builtin all.equal which has its own uses.

Trigonometric Accuracy in R

2 Answers2