Created
October 14, 2015 10:27
-
-
Save akkie/72fd1ce223dfd520db60 to your computer and use it in GitHub Desktop.
Best rational approximation
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
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) | |
} | |
} |
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
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