rationalisering — more issues with floating-point numbers

#	User	Rating
1	tourist	4009
2	jiangly	3831
3	Radewoosh	3646
4	jqdai0815	3620
4	Benq	3620
6	orzdevinwang	3529
7	ecnerwala	3446
8	Um_nik	3396
9	gamegame	3386
10	ksun48	3373

#	User	Contrib.
1	cry	164
1	maomao90	164
3	Um_nik	163
4	atcoder_official	160
5	-is-this-fft-	158
6	awoo	157
7	adamant	156
8	TheScrasse	154
8	nor	154
10	Dominater069	153

I was motivated to write this blog post after reading this blog post about possible precision issues in output and requiring printing a rounded answer exactly instead of printing an answer within a tolerance. This blog post feels somewhat related insofar as that it has an issue with floating point numbers, and also a doubt as to the correctness of the underlying data. Alternatively, someone can find the mistake in my logic.

The problem in question is Rationalization which asks to find a positive rational $$$\frac{A}{B}$$$ that is within a range $$$[C-F, C+F]$$$. Among all such fractions, minimize $$$A$$$, and among all such with minimal $$$A$$$, minimize $$$B$$$. The judge data guarantee that $$$1 \le A, B < 10^6$$$.

My solution path is as follows: If $$$C-F \le \frac{A}{B} \le C+F$$$, then $$$B(C-F) \le A \le B(C+F)$$$. Therefore, we can just check all $$$B$$$ in increasing order and find the smallest $$$B$$$ where $$$[B(C-F), B(C+F)]$$$ contains some integer, and the smallest such integer should be our $$$A$$$.

Now, $$$C$$$ and $$$F$$$ are floating-point numbers, so in order to make sure that we can represent $$$C-F$$$ and $$$C+F$$$ exactly, we scale them up by powers of 10 until they are exactly integers. Fortunately, we're guaranteed that these numbers are written in decimal form.

Some poorly written Python code

c, f = input().split()
if '.' not in c:
  c += ".0"
if '.' not in f:
  f += ".0"
if c[2] != '.':
  c = "0" + c
if f[2] != '.':
  f = "0" + f
while len(c) != len(f):
  if len(c) < len(f):
    c += "0"
  else:
    f += "0"
scale = 10 ** (len(c)-3)
c = int(c[0:2] + c[3:])
f = int(f[0:2] + f[3:])
num = scale
MAXVAL = 10 ** 6
for b in range(1, MAXVAL):
  while num < b * (c-f):
    num += scale
  if num <= b * (c+f):
    assert (num // scale) < MAXVAL
    print(num // scale)
    print(b)
    found = True
    break
assert found

This gets 60 points on Kattis with a WA verdict on the final subtask, so it does print something which is deemed incorrect (as opposed to printing nothing or printing a value of $$$A$$$ that is too large and getting RTE). Since the solution does pass the first two subtasks, I am reasonably confident that it is not completely wrong. I presume that the second subtask is meant to admit solutions that have precision issues by representing values using floating point values (perhaps by scaling each of the floating point values directly by $$$B$$$), though I haven't bothered to experiment much there.

I do have a weak prior that some Kattis problems with floating point numbers do something incorrect — for example, for this problem, all reference solutions say that three points are collinear if and only if the turn would be at most $$$10^{-9}$$$ degrees, as opposed to checking collinearity exactly.

Because of this prior, I'd like to identify the test case to see where I went wrong. Sadly, I can't find the input data for this problem, so I'm stuck trying to find a bug in my solution/logic or conclude the test data are wrong. I'd just like to ask the author for test case 77 to verify it by hand. In the interim, I might be able to get the test case in a few hours if I can find a way to AC the first 76 cases and differentiate that I'm in test case 77. (Kattis only allows 10 submissions within a 10 minute window.)

Alternatively, did I go wrong somewhere?

Update 1: Test case 77 has c = 7.80212 and f = 0.0000000684. I wrote a slow program that seems to generate the same output as what my program prints.

More Python code

from fractions import Fraction
c = Fraction('7.80212')
f = Fraction('0.0000000684')
lhs = c - f
rhs = c + f
for a in range(8, 10 ** 6):
  found = False
  b = max(1, int(a / rhs))
  while b < 10 ** 6:
    cand = Fraction(a, b)
    if cand >= lhs and cand <= rhs:
      print(a)
      print(b)
      found = True
      break
    if cand < lhs:
      break
    b += 1
  if found:
    break

This fraction is not equal to either of $$$C-F$$$ and $$$C+F$$$.

Update 2: jeroenodb was able to AC the problem and with his solution we were able to prove that the test data are indeed incorrect. I have contacted Kattis to try to get this rectified.

#include "bits/stdc++.h" using namespace std; const int oo = 1e9; int main() { double c,f; cin >> c >> f; pair<int,int> ans = {oo,oo}; for(int b=1;b<int(1e6);++b) { auto lo = max(c-f-1e-9,0.); double a = llround(ceil(lo*b)); if(a==0) a=1; if(a/b<=(c+f)) { ans=min(ans,{a,b}); } } cout << ans.first << '\n' << ans.second << '\n'; }

Comments (8)

Write comment?

adamant

17 months ago, # |

← Rev. 2 →

My solution, based on continued fractions also gets WA on test 77:

Code

C, F = map(float, input().split())
Q = 10**9
C = int(C*Q)
F = int(F*Q)
L = C-F
R = C+F

def fraction(p, q):
    a = []
    while q:
        a.append(p // q)
        p, q = q, p % q
    return a

def convergents(a):
    p = [0, 1]
    q = [1, 0]
    for it in a:
        p.append(p[-1]*it + p[-2])
        q.append(q[-1]*it + q[-2])
    return p, q

# check if a < b assuming that a[-1] = b[-1] = infty and a != b
def less(a, b):
    a = [(-1)**i*a[i] for i in range(len(a))]
    b = [(-1)**i*b[i] for i in range(len(b))]
    return a < b

# [a0; a1, ..., ak] -> [a0, a1, ..., ak-1, 1]
def expand(a):
    if a: # empty a = inf
        a[-1] -= 1
        a.append(1)
    return a

# return a-eps, a+eps
def pm_eps(a):
    b = expand(a.copy())
    a.append(float('inf'))
    b.append(float('inf'))
    return (a, b) if less(a, b) else (b, a)

def middle(p0, q0, p1, q1):
    a0 = pm_eps(fraction(p0, q0))[1]
    a1 = pm_eps(fraction(p1, q1))[0]
    a = []
    for i in range(min(len(a0), len(a1))):
        a.append(min(a0[i], a1[i]))
        if a0[i] != a1[i]:
            break
    a[-1] += 1
    p, q = convergents(a)
    return p[-1], q[-1]

p, q = middle(L,Q, R,Q)
print(p)
print(q)

And it also gives $$$\frac{133229}{17076}$$$ on your test case, so something seems weird...

UPD: I also tried few other solutions, including basically brute force. Still WA77 with the same output. Looks like an error in test data?

→ Reply

nicksms

Paging esr6vqa.

jeroenodb

+55

I got accepted with this very bad, float code with an epsilon of $$$10^{-9}$$$, but only on the left edge of the interval. Possibly the author did something similar, so this is recreating the same bug.

On this test

7.80212 0.0000000684

It prints:

131021 16793

xiaowuc1

17 months ago, # ^ |

Thanks! It would seem that this is the only test case among the judge data that is incorrect (my Python code gets AC if I patch exactly that test case).

EMBailey

+11

The Kattis problem Dice and Ladders seems to have the same issue. This code, which is based on one of the original contest's official solutions, passes:

Spoiler

#include <bits/stdc++.h>

using namespace std;

constexpr int MAXN = 64;

typedef double T;

struct matrix {
    T m[MAXN][MAXN]{};

    matrix operator*(const matrix &other) const {
        matrix res = matrix();
        for (int i = 0; i < MAXN; i++)
            for (int j = 0; j < MAXN; j++)
                for (int k = 0; k < MAXN; k++)
                    res.m[i][j] += m[i][k] * other.m[k][j];
        return res;
    }
};

matrix matrix_pow(int n, matrix m) {
    if (n == 1) return m;
    matrix rec = matrix_pow(n / 2, m);
    return n % 2 ? m * rec * rec : rec * rec;
}

int r, c, n;
T p;

int main() {
    cin >> r >> c >> n >> p;
    matrix dice;
    for (int i = 0; i < r*c; i++) {
        for (int j = i+1; j <= i+6; j++)
            dice.m[i][min(j, r*c-1)] += 1 / T(6);
    }

    matrix ladder;
    for (int i = 0; i < r * c; i++)
        ladder.m[i][i] = 1;
    for (int i = 0; i < n; i++) {
        int a, b;
        cin >> a >> b;
        a--;
        b--;
        ladder.m[a][a] = 0;
        ladder.m[a][b] = 1;
    }

    matrix mat = dice * ladder;

    int lo = 1, hi = 1e8;
    while (lo < hi) {
        int m = (lo + hi) / 2;
        if (matrix_pow(m, mat).m[0][r*c-1] >= p)
            hi = m;
        else
            lo = m + 1;
    }

    cout << lo << endl;
}

However, increasing the precision by changing T to long double causes it to WA on test case 16/17.

I also contacted Kattis about this in February of 2022. However, it does not seem to have been fixed.

kelzzin2_fan

Thanks for making this blog post! I had previously made a blog post asking for help on this question and unfortunately it went nowhere.

William_Kraft

+50

Hello! Thanks for this post!

You are indeed correct that there is an error with the testdata.

While implementing the solution to this problem we accidentally used an error margin in the solution. As only one programmer wrote judge solutions to this problem, the error margin was the same in all of them, which caused the test case to have a wrong answer.

None of us checking the problem caught this issue, and it was only found during the Doris competition in which it was originally used.

The issue was promptly fixed in the Github repository. But we did not mention it to Kattis, as we thought they would reinstall from the repository when adding it to open.kattis.com, an error from our part.

We later noticed it having an unusually high difficulty rating, and that that the old problem was used, and have since contacted Kattis hoping to have it corrected.

The extra error margin is 10^-9, that is, add 10^-9 to F to have the correct solution.

We are sorry for the inconvenience, and will make sure to have multiple different people write solutions as well as to be more clear with Kattis when these errors happen in the future.

xiaowuc1's blog