# Diophantine Approximation: find lowest possible denominator to approximate within given precision

## Challenge

Given a number x and a precision e, find the lowest positive integer q such that x can be approximated as a fraction p / q within precision e.

In other words, find the lowest positive integer q such that there exists an integer p such that abs(x - p/q) < e.

## Input

• The pair (x, e) where x is a floating-point number, and e is a positive floating-point number.
• Alternatively, a pair (x, n) where n is a nonnegative integer; then e is implicitly defined as 10**(-n) or 2**(-n), meaning n is the precision in number of digits/bits.

Restricting x to positive floating-point is acceptable.

## Output

The denominator q, which is a positive integer.

## Test cases

• Whenever e > 0.5 ------------------------> 1 because x ≈ an integer
• Whenever x is an integer ----------------> 1 because x ≈ itself
• (3.141592653589793, 0.2) ------------> 1 because x ≈ 3
• (3.141592653589793, 0.0015) --------> 7 because x ≈ 22/7
• (3.141592653589793, 0.0000003) ---> 113 because x ≈ 355/113
• (0.41, 0.01) -------------------------------> 12 for 5/12 or 5 for 2/5, see Rules below

## Rules

• This is code-golf, the shortest code wins!
• Although the input is "a pair", how to encode a pair is unspecified
• The type used for x must allow a reasonable precision
• Floating-point precision errors can be ignored as long as the algorithm is correct. For instance, the output for (0.41, 0.01) should be 12 for 5/12, but the output 5 is acceptable because 0.41-2/5 gives 0.009999999999999953

## Related challenges

• Given the test cases, this OEIS sequence is vaguely related. Commented Sep 4, 2020 at 13:58
• For that kind of challenges, I believe that the consensus is that floating point precision errors can be ignored as long as the algorithm is correct. You may however want to add a note about that. For instance, most answers (including mine) are returning the 5 of 2/5 instead of the 12 of 5/12 for (0.41, 0.01) -- because 0.41-2/5 gives 0.009999999999999953. Commented Sep 4, 2020 at 14:38
• Thanks! I added a note. Also made me realise "Whenever e >= 0.5" should be changed to "Whenever e > 0.5" to be consistent with the strict < in the definition "abs(x - p/q) < e".
– Stef
Commented Sep 4, 2020 at 14:51
• As a mathematician working in Diophantine Approximation, I approve of this challenge.
– A.P.
Commented Sep 5, 2020 at 15:37

# R, Xx bytes

Note: this challenge is quite a good introductory-challenge for R, which is the 'language-of-the-month' for September 2020, so I've blanked-out my answer in the hope of encouraging some other golfers to have a shot at it in R, too...

50 bytes

function(x,e,s=1:e^-1)s[(x-round(x*s)/s)^2<e^2][1]

Try it online!

Function with arguments x & error e. Can handle negative x (even though not required for challenge)

Note 2: dammit! a port of xnor's approach is 6 bytes shorter still:

44 bytes

function(x,e,s=1:e^-1)s[(x+e)%%(1/s)<2*e][1]

Try it online!

# 05AB1E, 13 9 bytes

∞.Δ*Dòα›


Try it online!

Commented:

          # implicit input                    [e, x]
∞         # in the list of natural numbers
.Δ       # find the first that satisfies:    [e, x], q
*      #   multiply                        [e*q, x*q]
#   dump on stack                   e*q, x*q
D    #   duplicate                       e*q, x*q, x*q
ò   #   round to integer                e*q, x*q, round(x*q)
α  #   absolute difference             e*q, abs(x*q - round(x*q))
› #   is this larger?                 e*q > abs(x*q - round(x*q))


# Python, 46 bytes

f=lambda x,e,q=1:(x+e)%(1/q)<e*2or-~f(x,e,q+1)


Try it online!

We want to check that $$\x\$$ is within $$\\pm \epsilon\$$ of a multiple of $$\1/q\$$, that is, it falls within the interval $$\(-\epsilon,\epsilon)\$$ modulo $$\1/q\$$. To do this, we take $$\x+\epsilon\$$, reduce it modulo $$\1/q\$$, and check if the result is at most $$\2 \epsilon\$$.

A same-length alternative using only %1, which might help with porting:

f=lambda x,e,q=1:(x+e)*q%1<e*q*2or-~f(x,e,q+1)


Try it online!

• Did you mean that it falls within the interval...? Commented Sep 5, 2020 at 14:14

# Wolfram Language (Mathematica), 24 bytes

Denominator@*Rationalize


Try it online!

All the credits go to @the default

• Can you just do Denominator@*Rationalize? (I did not post it because I thought it's too long and I do not really like non-winning Mathematica answers) (it's a shame Part doesn't work on Rationals) Commented Sep 4, 2020 at 12:17
• @thedefault. did you change your pronoun? Commented Sep 4, 2020 at 12:31
• Not really, my username has been reset to the default (by the Stack Exchange Community Management Team™), so I set it to the default. Commented Sep 4, 2020 at 12:35

# Python 3, 74 $$\\cdots\$$ 52 50 bytes

Saved a 4 6 bytes thanks to ovs!!!

f=lambda x,e,q=1:not-x*q%1>e*q<x*q%1or-~f(x,e,q+1)


Try it online!

• Instead of testing int(x*q) and int(x*q)+1, you can just use round(x*q) and get rid of the max.
– ovs
Commented Sep 4, 2020 at 14:11
• @ovs Had just realised that but was using int(x*q+.5) which is longer than round - thanks! :-) Commented Sep 4, 2020 at 14:21
• 52 bytes by not using q for the final result.
– ovs
Commented Sep 4, 2020 at 14:43
• @ovs Truely diabolical - thanks! :D Commented Sep 4, 2020 at 14:51
• ... And a slightly different approach for 50 bytes.
– ovs
Commented Sep 4, 2020 at 15:01

# JavaScript (ES7), 38 bytes

Expects (x)(e).

A port of @xnor's method, which is significantly shorter than my original approach.

(x,q=0)=>g=e=>(x+e)%(1/++q)<e*2?q:g(e)


Try it online!

# JavaScript (ES7), 46 bytes

Expects (x)(e).

(x,q=0)=>g=e=>((x*++q+.5|0)/q-x)**2<e*e?q:g(e)


Try it online!

We want to avoid using the lengthy Math.round() and Math.abs(). So we look for the lowest $$\q>0\$$ such that:

$$\left(\frac{\left\lfloor xq+\frac{1}{2}\right\rfloor}{q}-x\right)^2

# C (gcc), 63 59 58 bytes

Saved a byte using xnor's idea in his Python answer!!!

i;f(x,e,q)float x,e,q;{for(q=0;fmod(x+e,1/++q)>2*e;);i=q;}


Try it online!

# MATL, 17 5 bytes

2$YQ&  ### Explanation 2$   % The next function will take two inputs
YQ   % (Implicit inputs: x, e). Rational approximation with specified tolerance.
% Gives two outputs: numerator and denominator
&    % The next function will use its alternative default input/output
% configuration
% (Implicit) Display. With the alternative specification, this displays
% only the top of the stack, that is, the denominator


## Manual approach: 17 bytes

GZ}1\@:q@/-|>~}@


### Explanation

       % Do...while
GZ}   %   Push input: array [e, x]. Split into e and x
1\    %   Modulo 1: gives fractional part of x (*)
@:q   %   Push [0, 1, ... , n-1], where n is iteration index
@/    %   Divide by n, element-wise: gives [0, 1/n, ..., (n-1)/n]
-|    %   Absolute difference between (*) and each entry of the above
>~    %   Is e not greater than each absolute difference? (**)
}       % Finally (execute on loop exit)
@     %   Push current iteration index. This is the output
% End (implicit). A new iteration is run if all entries of (**) are true;
% that is, if all absolute differences were greater than or equal to e
% Display (implicit)


# Charcoal, 27 bytes

ＮθＮη≔¹ζＷ›↔⁻∕⌊⁺·⁵×θζζθη≦⊕ζＩζ


Try it online! Link is to verbose version of code. Explanation:

ＮθＮη


Input $$\ x \$$ and $$\ \epsilon \$$.

≔¹ζ


Start off with $$\ q = 1 \$$.

Ｗ›↔⁻∕⌊⁺·⁵×θζζθη


Calculate $$\ p = \lfloor 0.5 + q z \rfloor \$$ and repeat while $$\ | \frac p q - x | > \epsilon \$$...

≦⊕ζ


... increment $$\ q \$$.

Ｉζ


Output $$\ q \$$.

# Scala, 8460 52 bytes

Saved a whopping 24 bytes thanks to @Dominic van Essen!

x=>e=>1 to 9<<30 find(q=>(x-(x*q+.5).floor/q).abs<e)


Try it online!

• Surely since you're testing the floor & ceiling separately, and you know that floor<ceiling, you don't need either of the abs? Something like this? Commented Sep 4, 2020 at 15:22
• @DominicvanEssen You're right, thanks!
– user
Commented Sep 4, 2020 at 15:25
• Or maybe you could add 0.5 and just use the floor, like this? Commented Sep 4, 2020 at 15:25
• Wow, that's even better! You saved me 24 bytes
– user
Commented Sep 4, 2020 at 15:30

# Wolfram Language 89 bytes

f[n_,e_]:=Denominator@Cases[{#,Abs[n-#]}&/@Convergents@n,x_/;x[[2]]<=e][[1,1]]


f[0.41,.01]
(* 5. *)


This uses the convergents as candidates for approximations.

pi = 3.1415926535897932384626433832795028842


The first 8 convergents of pi:

Convergents[pi, 8]
(* {3, 22/7, 333/106, 355/113, 103993/33102, 104348/33215, 208341/66317, 312689/99532}*)

f[pi, 0.01]
(* 7 *)

f[pi, 0.001]
(* 106 *)

f[pi, 0.00001]
(* 113 *)

f[pi, 0.0000001]
(* 33102 *)

f[pi, 0.0000000001]
(* 99532 *)

• 65 bytes Commented Sep 5, 2020 at 21:09
• @J42161217 Post it! Commented Sep 6, 2020 at 0:28

# Java (JDK), 52 bytes

Port of xnor’s method

x->e->{int q=0;for(;(x+e)%(1./++q)>=e*2;);return q;}


Try it online!

# Java (JDK), 83 69 bytes

x->e->{int q=0;for(;Math.abs(x-Math.ceil(x*++q-.5)/q)>=e;);return q;}


Try it online!

• Commented Sep 4, 2020 at 15:59
• @Arnauld That's funny, I just wrote the exact same thing :)
– user
Commented Sep 4, 2020 at 16:02

# Perl 5, 58 bytes

sub f{grep{$p=$_[0]*$_;abs$p-int$p+.5<$_[1]*$_}1..1/$_[1]}


Try it online!