Skip to content

Instantly share code, notes, and snippets.

@akkie
Created October 14, 2015 10:27
Show Gist options
  • Save akkie/72fd1ce223dfd520db60 to your computer and use it in GitHub Desktop.
Save akkie/72fd1ce223dfd520db60 to your computer and use it in GitHub Desktop.
Best rational approximation
import scala.annotation.tailrec
object Math {
/**
* Calculates the `Best rational approximation` based on the Farey sequence.
*
* Translated from John D. Cook's Python implementation and David Weber's C++ modification.
*
* @see http://www.johndcook.com/blog/2010/10/20/best-rational-approximation/
* @see http://www.johndcook.com/blog/2010/10/20/best-rational-approximation/comment-page-1/#comment-17902
*
* @param x A value to be approximated by two integers.
* @param eps The required precision such that abs(x-a/b) < eps. Eps > 0.
* @param n The maximum size of the numerator allowed.
* @return The best rational approximation for x.
*/
def farey(x: Double, eps: Double, n: Int): (Int, Int) = {
@tailrec
def iterate(a: Int, b: Int, c: Int, d: Int): (Int, Int) = {
if (b <= n && d <= n) {
val mediant = (a + c).toDouble / (b + d).toDouble
if (java.lang.Math.abs(x - mediant) < eps) {
if (b + d <= n) {
(a + c) -> (b + d)
} else if (d > b) {
c -> d
} else {
a -> b
}
} else if (x > mediant) {
iterate(a + c, b + d, c, d)
} else {
iterate(a, b, a + c, b + d)
}
}
else if (b > n) c -> d
else a -> b
}
iterate(0, 1, 1, 0)
}
/**
* Calculates the `Best rational approximation` based on the Farey sequence.
*
* This is an adapted algorithm of the original function which needs slightly
* more iterations in some circumstances but no precision.
*
* @see http://www.johndcook.com/blog/2010/10/20/best-rational-approximation/
*
* @param x A value to be approximated by two integers.
* @param n The maximum size of the numerator allowed.
* @return The best rational approximation for x.
*/
def farey(x: Double, n: Int): (Int, Int) = {
@tailrec
def iterate(a: Int, b: Int, c: Int, d: Int): (Int, Int) = {
if (b <= n && d <= n) {
if (x > (a + c).toDouble / (b + d).toDouble) {
iterate(a + c, b + d, c, d)
} else {
iterate(a, b, a + c, b + d)
}
}
else if (b > n) c -> d
else a -> b
}
iterate(0, 1, 1, 0)
}
}
import play.api.test.PlaySpecification
class MathSpec extends PlaySpecification {
val eps = 0.00005
"The `farey` method" should {
"return the aspect ratio for 3/4" in {
Math.farey(3d / 4d, 10) must be equalTo((3, 4))
Math.farey(3d / 4d, eps, 10) must be equalTo((3, 4))
}
"return the aspect ratio for 4/3" in {
Math.farey(4d / 3d, 10) must be equalTo ((4, 3))
Math.farey(4d / 3d, eps, 10) must be equalTo ((4, 3))
}
"return the aspect ratio for 16/9" in {
Math.farey(16d / 9d, 10) must be equalTo ((16, 9))
Math.farey(16d / 9d, eps, 10) must be equalTo ((16, 9))
}
"return the aspect ratio for 21/9" in {
// http://www.tomshardware.co.uk/answers/id-2038871/resolutions.html
// Reduces to 7/3
Math.farey(21d / 9d, 10) must be equalTo ((7, 3))
Math.farey(21d / 9d, eps, 10) must be equalTo ((7, 3))
}
"return the best rational approximation for 0.127" in {
// http://www.cs.uni.edu/~wallingf/blog/archives/monthly/2010-10.html#e2010-10-25T16_50_29.htm
Math.farey(0.127, 74) must be equalTo((8, 63))
Math.farey(0.127, eps, 74) must be equalTo((8, 63))
}
"return the best rational approximation for 0.367879" in {
// http://www.johndcook.com/blog/2010/10/20/best-rational-approximation/comment-page-1/#comment-17902
Math.farey(0.367879, 10) must be equalTo((3, 8))
Math.farey(0.367879, eps, 10) must be equalTo((3, 8))
}
"return the best rational approximation for 0.36781" in {
// http://www.johndcook.com/blog/2010/10/20/best-rational-approximation/comment-page-1/#comment-17902
Math.farey(0.367879, 100) must be equalTo((32, 87))
Math.farey(0.367879, eps, 100) must be equalTo((32, 87))
}
"return the aspect ratio for resolution 1360x765" in {
// http://stackoverflow.com/a/13466237/2153190
Math.farey(1360d / 765d, 10) must be equalTo ((16, 9))
Math.farey(1360d / 765d, eps, 10) must be equalTo ((16, 9))
}
"return the aspect ratio for resolution 1360x768" in {
// http://stackoverflow.com/a/13466237/2153190
Math.farey(1360d / 768d, 10) must be equalTo ((16, 9))
Math.farey(1360d / 768d, eps, 10) must be equalTo ((16, 9))
}
"return the aspect ratio for resolution 1366x768" in {
// http://stackoverflow.com/a/13466237/2153190
Math.farey(1366d / 768d, 10) must be equalTo ((16, 9))
Math.farey(1366d / 768d, eps, 10) must be equalTo ((16, 9))
}
"return the aspect ratio for resolution 1936x2581" in {
Math.farey(1936d / 2581d, 10) must be equalTo ((3, 4))
Math.farey(1936d / 2581d, eps, 10) must be equalTo ((3, 4))
}
"return the aspect ratio for resolution 542x723" in {
Math.farey(542d / 723d, 10) must be equalTo ((3, 4))
Math.farey(542d / 723d, eps, 10) must be equalTo ((3, 4))
}
}
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment