-is-this-fft-'s blog

By -is-this-fft-, history, 2 months ago, In English

When solving problems in competitive programming, it is almost never a good idea to use the following inbuilt C++ functions. You will be hacked or fail a pretest or worse, a systest.

Why? Because they use floating-point numbers. They are designed to be used with a floating-point input and a floating-point output. The issue is that on a floating-point number, the result may not be exact. Worse, floating-point numbers may not be able to accurately encode the integer.

To calculate $$$\lfloor \frac{a}{b} \rfloor$$$ for positive integers $$$a$$$ and $$$b$$$:
  • Don't use floor((double) a / (double) b) or similar.
  • Do use a / b. It will round down.
  • Warning: be careful with negative numbers. The answer depends on whether we should round down or towards 0.
To calculate $$$\lceil \frac{a}{b} \rceil$$$ for positive integers $$$a$$$ and $$$b$$$:
  • Don't use ceil((double) a / (double) b) or similar.
  • Do use (a + b - 1) / b.
  • Warning: the same caveat goes for negative numbers as before.
To calculate $$$\lfloor \sqrt{a} \rfloor$$$ for a nonnegative integer $$$a$$$:
  • Don't use sqrt(a) or similar.
  • Do use binary search to calculate the square root.
Implementation
To calculate $$$a^b$$$ for nonnegative integers $$$a$$$ and $$$b$$$:
  • Don't use pow(a, b) or similar.
  • Do use the naive algorithm. If you really want to do this and $$$a^b$$$ fits in a long long, then it will take no more than 64 steps, unless $$$a = 0$$$ or $$$a = 1$$$, in which case you can just return $$$a$$$.
To calculate $$$\lfloor \log_2 (a) \rfloor$$$ for a positive integer $$$a$$$:
  • Don't use log2(a) or similar.
  • Do use __lg or an approach based on __builtin_clz or __builtin_clzll. The "clz" stands for "count leading zeroes" and it does exactly what it says — on the binary representation of the number. If you have access to C++20, there is also the bit_width library function.
  • Warning: __builtin_clz(0) and __builtin_clzll(0) are undefined behavior. Also all these functions starting with __ are specific to GCC and you're technically not meant to use them, but it's fine in competitive programming.

Exceptions

The only major exception is if floating-point numbers are an inherent part of your solution. Most commonly, this happens in problems where the output is itself a floating-point number (especially geometry problems, sometimes probability). There are also some algorithms where you need to work with floating-point numbers even if the input and output are integer (FFT, possibly some LP).

But even in these cases, it's always a good idea to do as much as possible with integers and move to floats as late as possible.

I got WA with one compiler, accepted with another

Most likely it has to do with 64-bit versus 32-bit. Also the widths of the various floating-point types vary. It's possible that in some circumstances your solution is actually correct. However, rely on this only if you know what you're doing.

 
 
 
 
  • Vote: I like it
  • +500
  • Vote: I do not like it

»
2 months ago, # |
  Vote: I like it +142 Vote: I do not like it

sqrt is fine just check your answer once after doing it

»
2 months ago, # |
Rev. 2   Vote: I like it +3 Vote: I do not like it

The code given in spoiler is dirty take a look at it -is-this-fft-

»
2 months ago, # |
  Vote: I like it +1 Vote: I do not like it

lol in hindsight me failing pretest 4 on B 3 times and switching from sqrt() to binary search helped me since I didn't fail sys test 11 xd

»
2 months ago, # |
  Vote: I like it +79 Vote: I do not like it

I'm glad that my contest helped raise competitors awareness to this matters. Many people are still using built-ins without caution or understanding the imprecision it could bring. Especially floating-point related!

»
2 months ago, # |
Rev. 2   Vote: I like it -30 Vote: I do not like it

Very usefull, now

Spoiler
»
2 months ago, # |
  Vote: I like it 0 Vote: I do not like it

M here at 69 minutes ago . Helpful content

»
2 months ago, # |
  Vote: I like it 0 Vote: I do not like it

Thanks for your benefit blog .can you make another blog about precision errors but in python since this blog about c++ ?
Thanks in advance .

»
2 months ago, # |
Rev. 2   Vote: I like it +82 Vote: I do not like it

Very helpful blog for the newcomers, thank you! I will definitely link this to my students.

Some thoughts:

  • in the section about $$$a^b$$$, maybe add a special case for $$$-1$$$? I know that you've written about non-negative integers, but with this additional case, it works for negative $$$a$$$ as well.

  • maybe add a section about comparing fractions? Like, if you want to check that $$$\frac{a}{b} = \frac{c}{d}$$$, you shouldn't divide in floating-point numbers; instead either check $$$a*d == b*c$$$ (if it fits in the integer types you use) or normalize the fractions.

  • isn't it possible to estimate $$$sqrt(n)$$$ almost perfectly with the standard function and then search for the answer in something like $$$[sqrt(n)-5, sqrt(n)+5]$$$? I know it's kinda advanced, but so is binary search.

  • »
    »
    2 months ago, # ^ |
    Rev. 2   Vote: I like it +35 Vote: I do not like it

    $$$\pm 5$$$ should be a rough guess and might work fine, but you can still just increment/decrement linearly until $$$R^2 \leq n$$$ and $$$(R+1)^2 > n$$$ are both true. this is guaranteed to work, and should still work fast enough as the error should be not very big.

  • »
    »
    2 months ago, # ^ |
      Vote: I like it +2 Vote: I do not like it

    if u won't mind can u pls tell where do u teach?And how one can enroll in your courses?

    Thank you.

  • »
    »
    2 months ago, # ^ |
      Vote: I like it -7 Vote: I do not like it

    Would you please update the difficulty rating of the last edu round's problems?

»
2 months ago, # |
  Vote: I like it 0 Vote: I do not like it

Reading this makes me feel very thankful that I AC'd Dytechlab's contest question B using normal std::sqrt()

»
2 months ago, # |
Rev. 2   Vote: I like it +73 Vote: I do not like it

sqrt(x) is fine, just do something like this:

int64_t x;
cin >> x;
int64_t r = sqrtl(x) + 1;
while (r * r > x) r--;

sqrtl() use long double instead of double, and might not be necessary.

+1 should be enough, but if you are paranoid then +small would be better.

Generally this should be better than binary search.

If you for some reason what to avoid having any floating point at all, then try to implement the search as $$$r_{i+1} = (r_{i} + x/r_{i}) / 2$$$ instead, it might be faster than just binary search.

»
2 months ago, # |
  Vote: I like it +105 Vote: I do not like it

You can compute $$$\lfloor \sqrt a \rfloor$$$ faster and easier than with binary search:

long long x = sqrt(a) + 2;
while (x * x > a) x--;
  • »
    »
    2 months ago, # ^ |
    Rev. 2   Vote: I like it 0 Vote: I do not like it

    Is there formal proof that $$$R+2 \geq \sqrt{x}$$$ is true given $$$R$$$ as the result given by the standard function? I know that the error will be minimal, but I can't get an exact upper bound for the error.

    • »
      »
      »
      2 months ago, # ^ |
        Vote: I like it +7 Vote: I do not like it

      It depends on the compiler version, and nothing is guaranteed. Read this old blog by misof: The curious case of the pow function.

      Fortunately, for most operations (including sqrt), the magnitude of the relative precision error is equal to the type precision. It should be around $$$10^{-15}$$$ for double. The square root of long long usually doesn't exceed $$$10^9$$$ so the absolute error should be around $$$10^{-6}$$$. It means that sqrt(a) + 0.0001 is more than enough already.

      more info
      more issues
      • »
        »
        »
        »
        2 months ago, # ^ |
          Vote: I like it +31 Vote: I do not like it

        For sqrt there are some guarantees. https://en.cppreference.com/w/cpp/numeric/math/sqrt

        std::sqrt is required by the IEEE standard to be correctly rounded from the infinitely precise result. In particular, the exact result is produced if it can be represented in the floating-point type. The only other operations which require this are the arithmetic operators and the function std::fma. Other functions, including std::pow, are not so constrained.

      • »
        »
        »
        »
        2 months ago, # ^ |
        Rev. 2   Vote: I like it 0 Vote: I do not like it

        If this is the case, isn't sqrt(a) + 1 ok (using the same code snippet, just changing 2 to 1)? Or am I missing something?

        • »
          »
          »
          »
          »
          2 months ago, # ^ |
            Vote: I like it +4 Vote: I do not like it

          Yes. As I said, even sqrt(a)+0.0001 should be enough.

      • »
        »
        »
        »
        2 months ago, # ^ |
        Rev. 2   Vote: I like it +46 Vote: I do not like it

        Just a warning about the blog The curious case of the pow function by misof. Not a single person back then realized what was going on. There are so many erroneous theories in that blog.

        The thing no one figured out is that using 32 bit g++ compiler, every float and every double, is actually a 80 bit long double. Only when a floating point variable is stored in memory will it actually be stored as its corresponding type. Otherwise, all floating point variables will act as long doubles.

        For example, run this program with the input 123456789.

        #include <bits/stdc++.h>
        using namespace std;
        int main () {
            double n;
            cin >> n;
            double x = n * n;
            cout << (long long)x;
        }
        

        You will get the output 15241578750190521 (which should be impossible since 15241578750190521 cant be represented as a double). The reason it is possible is that x is infact a long double.

        Now try running this program

        #include <bits/stdc++.h>
        using namespace std;
        int main () {
            double n;
            cin >> n;
            double x = n * n;
            pow(x, 10);
            cout << (long long)x;
        }
        

        This time around the result is 15241578750190520. This is because the pow(x, 10) call effectively forces x to be stored in memory, which rounds x from 15241578750190521 to 15241578750190520.

        • »
          »
          »
          »
          »
          2 months ago, # ^ |
            Vote: I like it 0 Vote: I do not like it

          Wow, interesting. So the precision can be better than expected and it might cause a misunderstanding of precision. Is there any source where I can read more about this? If I encountered such a difference during a contest, I would just assume that the compiler optimizes something in the first code.

»
2 months ago, # |
Rev. 2   Vote: I like it -15 Vote: I do not like it

For the round-up, I use this.

ans=ceil(a*1.00/b);

I never got an error in this (till now), but is there a possibility that this way can fail?

Thanks!

Why am I being downvoted? (Sorry if I commented something wrong)

»
2 months ago, # |
  Vote: I like it -8 Vote: I do not like it

for $$$log_2$$$ this way may be more useful:
f[0]=0;for(int i=1;i<=N;i++) f[i]=f[i>>1]+1;

However,my friends told me "__builtin_clz()" is as fast as the integer operator "+".

  • »
    »
    2 months ago, # ^ |
      Vote: I like it -14 Vote: I do not like it

    That's a common method used in implementing Sparse Tables for RMQ. However, in my experience, __lg worked very well also.

  • »
    »
    2 months ago, # ^ |
      Vote: I like it +8 Vote: I do not like it

    Memory access is hell of a lot slower than most (series of) arithmetic operations. Use literally anything else: __lg if your compiler supports it, std::bit_width if you're using C++20, even a loop like while(x > 0) will get optimized away.

»
2 months ago, # |
  Vote: I like it +8 Vote: I do not like it

THanks for the blog!

»
2 months ago, # |
  Vote: I like it 0 Vote: I do not like it

-is-this-fft-, could you please fix (again?) the comment in the attached code? For me, it is shown as // works for at least up to (2**31 &mdash; 1)**2

Also, one might want to use binary exponentiation instead of pow. The fun fact is that in this case even $$$0$$$ or $$$1$$$ to some enormous power are calculated as fast as $$$2^{62}$$$ is suggested to be calculated. So it is not obligatory to if them, and if it's done, the code is even faster. Probably this part is "advanced" too...

»
2 months ago, # |
  Vote: I like it +53 Vote: I do not like it

You can also use the Babylonian method to find the square root. It converges quadratically and is very fast.

int isqrt(int x)
{
    // Integer square root of x.
    // 
    // Use the Babylonian method, adapted to integers. The next estimate will
    // be less than the current one, as long as the latter is greater than √x.
    // In the end, either the estimates will match, or they will trade places.
    // At this point, return the non-greater estimate.
    // 
    // We round the first estimate up to avoid it becoming zero.
    int a = x, b = (x + 1) / 2;
    while (a > b)
    {
        a = b;
        b = (b + x / b) / 2;
    }
    return a;
}

With a good initial guess the method works even faster.

  • »
    »
    2 months ago, # ^ |
    Rev. 2   Vote: I like it -24 Vote: I do not like it

    How does this compare to using Newton's method? AFAIK, Newton's method has a logarithmic time complexity and should work equally well.

    • »
      »
      »
      2 months ago, # ^ |
        Vote: I like it +3 Vote: I do not like it

      The Babylonian method is a specialization of Newton's method for the square root case.

  • »
    »
    2 months ago, # ^ |
      Vote: I like it -8 Vote: I do not like it

    Can you please mention the complexity of this approach ?

»
2 months ago, # |
  Vote: I like it -10 Vote: I do not like it

We had wasted almost 45 minutes in icpc Gwalior-Pune regional contest, just because sqrt was giving a wrong answer. We had calculated sqrt using binary_search and our solution got accepted. Just because we are using inbuilt sqrt function, we did 2 wrong submition and we were not able to perform as we have planned.

»
2 months ago, # |
Rev. 2   Vote: I like it -23 Vote: I do not like it

So why the imprecise Builtins are there if we can't use them?

»
2 months ago, # |
Rev. 2   Vote: I like it -18 Vote: I do not like it

Meanwhile, people writing ll c = a * b; for $$$a, b \le 10^9$$$.

»
2 months ago, # |
  Vote: I like it +9 Vote: I do not like it

Several more points about make int, not float

Cross product

Cross product calculates $$$(-1)^n \cdot 2S$$$. The sign depends on rotation of points in triangle. It can be used in tasks related to convex hull. Don't calculate angles.

int c_p(int x1, int y1, int x2, int y2, int x3, int y3) {
    return (x2 - x1) * (y3 - y1) - (x3 - x1) * (y2 - y1);
}

Task

Square

If the points of triangle are integers, then the square of it is $$$\frac{int}{2}$$$. We can operate on double square, which are ints. We can get it using vector product, which returns $$$(-1)^n \cdot 2S$$$.

For example, we can safely check if point $$$(x4, y4)$$$ is inside of triangle (including boundaries).

if (abs(c_p(x1, y1, x2, y2, x3, y3) == abs(c_p(x4, y4, x1, y1, x2, y2) + abs(c_p(x4, y4, x1, y1, x3, y3) + abs(c_p(x4, y4, x2, y2, x3, y3)) ...

Task

Dot product

Dot product is needed more rarely. It can be used to check, if the angle is acute, right or obtuse. For example, it can be used to find distance from point to segment on plane. The are two cases: when perpendicular from point to segment's line lies on segment, or not. It can be checked with dot product.

int d_p(int x1, int y1, int x2, int y2, int x3, int y3) {
    return (x2 - x1) * (x3 - x1) + (y2 - y1) * (y3 - y1);
}

Task

Sort by polar angle

Assume, you need to sort point by polar angle relative to vector $$$(1, 0)$$$. You can do $$$a[i] = atan2(y[i], x[i]);$$$ (if you don't know, what is atan2(y, x), discover it, it helpful in reducing ifs). Instead we can (unfortunately, with ifs) sort points in one quadrant by their tan (which is monotone there). But compare not $$$\frac{y[i]}{x[i]} - \frac{y[j]}{x[j]}$$$, but $$$y[i] * x[j] - y[j] * x[i]$$$.

Task Task

Check, if overflow happened

Assume, in problem you have to multiply numbers $$$\le 10^{18}$$$, but you know for sure, that if at some moment the number is $$$\ge 10^{18}$$$, the answer is "NO". How can we check, if multiplication won't give overflow?

Spoiler

Assume, you have to check, if $$$x \cdot y \le C$$$. Instead, check if $$$x \le \frac{C}{y}$$$

  • »
    »
    2 months ago, # ^ |
      Vote: I like it +8 Vote: I do not like it

    To check if overflow happened, there is a this sometimes useful list of GCC functions.

»
2 months ago, # |
  Vote: I like it 0 Vote: I do not like it

Floating-point numbers can accurately represents $$$2^n$$$, so maybe long double is still useful for some problems that require $$$2^n$$$ of your output. It can represent $$$2^n$$$ for any integer $$$n<15000$$$ (in the circumstances where long double have 15 bits for Exponent).

»
2 months ago, # |
  Vote: I like it 0 Vote: I do not like it

Some people are suggesting sqrtl. Can someone write or refer to a convincing analysis to show that sqrtl is always correct?

  • »
    »
    2 months ago, # ^ |
    Rev. 2   Vote: I like it +5 Vote: I do not like it

    On GCC, long double is 80 bits, and its mantissa is 63 bits, which is enough to store all long long integers, so the argument of sqrtl can be represented precisely. IEEE 754 requires that the returned value is represented as precisely as possible, which translates to exact integers when the root is whole and something slightly bigger than (or equal to? not sure if this is possible but it doesn't harm the analysis so why not) the floor of the exact root, which translates to $$$\lfloor\sqrt{n}\rfloor$$$ after truncating. Sounds right?

    • »
      »
      »
      2 months ago, # ^ |
        Vote: I like it 0 Vote: I do not like it

      Is it possible that the square root is not representable and gets rounded up to the next representable float, which happens to be an integer?

      • »
        »
        »
        »
        2 months ago, # ^ |
        Rev. 12   Vote: I like it +42 Vote: I do not like it

        IEEE 754 guarantees that the absolute error of a square root is within 0.5 ulp. Note that GCC documentation does not seem to confirm that as a requirement, so that may fail on non-x86 machines, but sqrtl uses the fsqrt on x86, which seems to adhere to the typical rounding requirements. As the square root of a number less than $$$2^{62}$$$ is less than $$$2^{31}$$$, the maximal possible exponent is 30, and hence 0.5 ulp is at most $$$0.5 \cdot 2^{30} \cdot 2^{-63} = 2^{-34}$$$. This guarantees that

        $$$ \left| \mathrm{sqrtl}(n) - \sqrt{n} \right| \le 2^{-34}, \left| \mathrm{sqrtl}(n+1) - \sqrt{n+1} \right| \le 2^{-34}. $$$

        Suppose $$$\mathrm{sqrtl}(n) = \mathrm{sqrtl}(n + 1)$$$, where $$$n < 2^{62}$$$, then due to the triangle inequality the above implies

        $$$ \sqrt{n+1} - \sqrt{n} \le 2^{-33}, $$$

        but we have

        $$$ \sqrt{n+1} - \sqrt{n} = \frac{1}{\sqrt{n+1} + \sqrt{n}} > \frac{1}{2\sqrt{n}} \ge 2^{-32}, $$$

        a clear contradiction.

        Nota bene: Notice that I used $$$2^{62}$$$ instead of $$$2^{64}$$$ or $$$2^{63}$$$. That is because if we allow integers up to $$$2^{63}$$$, we have a failure at $$$\mathrm{sqrtl}(2^{63}-3) = \mathrm{sqrtl}(2^{63}-4)$$$.

        • »
          »
          »
          »
          »
          2 months ago, # ^ |
            Vote: I like it +18 Vote: I do not like it

          I swear, if I have to edit this comment one more time...

»
2 months ago, # |
Rev. 2   Vote: I like it +8 Vote: I do not like it

C++20 also has https://en.cppreference.com/w/cpp/numeric/countr_zero and https://en.cppreference.com/w/cpp/numeric/countl_zero to replace __builtin_ctz and __builtin_clz. They return the correct answer (the number of bits of the corresponding datatype) when passed 0 as input, so there's no undefined behavior.