Counting points on a given irregular grid that are inside or outside a circle

Question

Given some points that are represented as (X, Y) coordinates, my task is to count the points that are inside a circle, on the circle, and outside the circle.

The input is given in the following way: first, the number n is provided. Then, n numbers are inputted that represent the SetX. Afterwards, the numbers for m and SetY are given in the same way.

The points that I should consider are Cartesian product of SetX and SetY. For example, given SetX: {0, 3, 4} and SetY: {1, 4}, the points would be: (0, 1), (0, 4), (3, 1), (3, 4), (4, 1), (4, 4).

The maximum element of SetX is guaranteed to be the same as the maximum element in SetY, and must be an even number. The center of previously mentioned circle is determined as (max/2, max/2), with radius max/2.

Additional limitations: n and m are at least 1 and no more than 10⁷. Elements of SetX and SetY are no greater than 10⁷.

My solution is as follows:

#include <iostream>

inline double distance_squared(const unsigned int&, const unsigned int&);

unsigned int center = 0;

int main()
{
    unsigned int n, m, i, inside = 0, on_circut = 0, outside = 0, tmp, distance_sqr;

    std::cin >> n;

    unsigned int set_X[n];
    for(i = 0; i < n; ++i){
        std::cin >> set_X[i]; // reading the set of Xs
        if(set_X[i] > center) center = set_X[i]; // finding maximum
    }
    center = center >> 1; //div by 2

    std::cin >> m;

    while(m--){ // reading elements of Y one by one
        std::cin >> tmp; //
        for(i = 0; i < n; ++i){ // for each pair of tmp (being Y coordinate) and a number from Xs (being coord. X)...
            distance_sqr = distance_squared(set_X[i], tmp); // find the distance^2 and compare it to radius^2
            if(distance_sqr == center*center){ ++on_circut; continue; }
            if(distance_sqr > center*center) ++outside;
            else ++inside;
        }
    }

    std::cout << inside << ' ' << on_circut << ' ' << outside;

    return 0;

}

inline double distance_squared(const unsigned int& x, const unsigned int& y){
    return (x - center)*(x - center) + (y - center)*(y - center);
}

The solution provides correct output, but it appears to be too slow. The online judge gives the "time limit exceeded" on majority of the tests. Please do note that in the original code I do not use std::cin and std::cout, but custom functions that I omitted, since I am certain that they are significantly faster than iostream functions, plus, they are quite long compared to the rest of the program.

Is there any other way I could optimise the code? Is my approach generally bad in terms of time complexity?

Answer 1

Well, you are manually testing every point \$(X_i, Y_j)\$ to see where it falls. You shouldn't, really.
Remember the sage advice: If in doubt, sort.

A better algorithm goes:

1. Read all X-coordinates.
2. Change them to squared distance from center-line.
3. Sort them.
4. Read an Y-coordinate.
5. Calculate the squared-distance from center-line for X-coordinates leading to the point sitting on the circle. \$\Delta x^2 = r^2 - (y - r)^2\$
6. Use binary search to determine how many elements from point 3 are smaller/equal/greater than the target.
7. Repeat at 4 until done.

Your algorithm was \$O(n \cdot m)\$. This one is \$O((n + m) \cdot \log(n))\$.

As an aside, if you simply double the given coordinates before processing them, you don't have to insist on the largest being even.

---

Now let's take a look at your code:

1. unsigned only has a minimal maximum of \$2^{16}-1\$, less than \$10^7\$. Even if it had 32 value-bits as on most desktops, \$2^{32}-1\$ is still below \$10^{14}\$. So, use long long which has a guaranteed 64 value-bits.
2. Why do you use double as the result of distance_squared()? If you did the calculation with double it would be accurate, though slow, but you only convert implicitly on return, which is far too late to help.
3. An unsigned is cheaply copied. Passing it by constant reference is actually a pessimization, luckily the optimizer will inline it so no harm done.
4. return 0; is implicit for main().
5. You should always finish your output with a newline. Not doing so means the prompt after calling the program won't be on a new line, and on ancient systems could have even less desirable outcomes.
6. The one time you use continue; it doesn't actually have any effect.
7. Some people like putting single-statements on the same line as the corresponding if, else if or else. I would desist as it impairs readability.
8. All else being equal, try to define any variable in the smallest scope you can, so a reader doesn't have to trace it all over the place. Concentration is a resource best preserved where possible.
9. inline only has the effect of allowing multiple equivalent definitions left. If it's only used in the current TU, better mark it static to get internal linkage.
10. Also consider defining a function before its use so you don't have to add a separate declaration. Remember, don't repeat yourself (DRY).
11. Restrict line-length so there's no horizontal scrolling. Horizontal scrolling, especially couple with vertical scrolling, is the bane of readability.

Answer 2

@Deduplicator has given good advice: use binary search to find the X-coordinates that are nearest the circle.

I would just like to point out that unsigned int set_X[n]; attempts to declare a variable-length array. Variable-length arrays are not part of the C++standard. This is a C99 feature that some C++ compilers happen to support. Thus, your code is not standard C++, and could likely fail to compile. Furthermore, a VLA of 10⁷ elements would be abusive, and could crash the stack. Use std::vector instead.

Answer 3

The other answers have good suggestions for algorithmic improvements. I'm going to mostly limit myself to some subtle ways your code is less compiler-friendly than it could be. (Also a potential bug: x - center wraps because the operands are unsigned. See below).

Sorting will be much better than just applying brute force better, but that can be fun :) Your version leaves a lot of performance on the table even without any significant algorithmic changes, depending on compiler version/options and target CPU. (My changes let it auto-vectorize better with -march=native on more CPUs, and maybe significant gains on typical x86 CPUs even when it can't auto-vectorize.) Also, hoisting a lot of computation out of the inner loop is very good.

---

I didn't see anything in the problem statement guaranteeing that SetX and SetY values are non-negative, just that their maximum is 10^7. I'm assuming that was something to justify your use of unsigned.

double is slow and totally unnecessary. With 64-bit integers, you can keep everything integer. (Fun fact, clang3.7 and later optimize away the conversion from 64-bit integer to double and back.)

Passing read-only int args by reference is silly. It doesn't matter here (the call inlines away), but in a case where the compiler decided not to inline, an extra level of indirection in the asm output from passing pointers instead of the values directly will hurt latency and throughput. On 32-bit CPUs, long long takes two registers (or two slots on the stack), but it's typically still better to pass wide integers by value.

A separate declaration for a simple inline helper-function just makes the code harder to read. I had to go find the function definition after my eyes found the declaration at the top first.

This would have been more sensible (but center as a function arg would have been better still, see below):

#include <cstdint>
static int32_t center = 0;   // Can be signed without breaking the right-shift, because we know it won't be negative.

static inline                      // note: signed args
 uint64_t distance_squared(int32_t x, int32_t y) {
    int64_t dx = x - center;
    int64_t dy = y - center;       // 32-bit signed subtraction is fine, and allows better auto-vectorization, and better scalar on 32-bit CPUs
    return dx*dx + dy*dy;          // but not ok for the multiply
}

main()...

Use int32_t for set_X[] too, to make sure it's wide enough without wasting space in your cache footprint if long is 64 bit. But actually, you should make it int64_t or uint64_t so you can loop over it once after center is known and replace each value with (x-center)*(x-center). This reduces some of the O(N*M) work do O(N). The compiler hoists the center*center and dy*dy loop invariants from the inner loop for you, but it doesn't change what's stored in set_X, so it has to redo x - center subtraction, and square that, before adding and comparing.

You could even transform the inequality to dx2 <= (r2 - dy2), so the inner loop is only a compare. (r2-dy2 is a loop-invariant). Note that that's problematic if they're unsigned and dy2 can be above r2, but it's safe with signed integers since it can't cause overflow in this case. (Compilers may not do this for you with signed integers for some reason. I think gcc doesn't like to change the order of operations for signed integers, because of the C standard making signed overflow Undefined Behaviour. This is a missed optimization for targets like x86 where overflow in a temporary doesn't raise an exception, and signed integer math is associative just like unsigned.)

Potential serious bug from unsigned subtraction wrapping around:

x - center wraps around instead of becoming negative, since you used unsigned long. This problem corrects itself after squaring, if you use the same type for the multiply as for the subtract. (An N*N => N-bit product gives the same result whether the inputs are treated as unsigned or two's complement, so squaring will produce the unsigned value that you want. It's only for a N*N => 2N-bit product that the upper half of the full-width result depends on the interpretation of the inputs.)

Unsigned wraparound becomes a problem when you widen after subtraction but before multiply, like I'm suggesting. (Helps the compiler auto-vectorize more efficiently, but does make the source more complicated. Getting the types wrong can break the code. It can also lead to compilers using extra instructions to do sign-extension later when they could have more efficiently done it sooner on a 64-bit machine. Clang does a nice job for scalar 64-bit here (loading set_X[i] with a sign-extending load instruction, but gcc is clunky and waits until after the subtraction.)

Simply using signed types everywhere is the easiest solution, since 2^31 - 1 is plenty large enough to hold 10^7.

If you did need to deal with x and y values that needed the full range of a 32-bit unsigned type, but didn't need 64-bit, you could get the compiler to zero-extend by doing
int64_t dx = x - static_cast<int64_t>(center). Making sure the compiler knows that x and y were originally 32-bit may help it optimize by using a constant zero for the upper half, but you might as well make the function arg types long long or int64_t. (Storing dx2 in an int64_t array makes this a non-issue for performance so you can just do whatever is easiest. It's a big win overall unless using twice as much cache means you bottleneck on memory bandwidth instead of cache, in which case some computation for each load is worth it).

---

Global variables aren't only bad style, they can also hurt performance.

unsigned int center = 0; is global, so the compiler has to assume any call to any non-inlined function can modify it. That includes std::iostream::operator>>. In both your loops that touch center, it's reloaded every iteration instead of staying in a call-preserved register.

e.g. in the first loop that finds center = max(SetX) on the fly, gcc6.3 for x86-64 emits this, with -O3 -fverbose-asm with some manual comments replacing the -fverbose-asm gcc-generated comments:

.L14:   # top of the loop
    mov     r13d, r12d        # i, i
    mov     edi, OFFSET FLAT:std::cin #,
    lea     rsi, [rbx+r13*4]  # tmp164,      ## set up args for the function call
    call    std::basic_istream<char, std::char_traits<char> >& std::basic_istream<char, std::char_traits<char> >::_M_extract<unsigned int>(unsigned int&)     #
    mov     eax, DWORD PTR [rbx+r13*4]        # load the cin>>set_X[i] result into eax
    mov     edx, DWORD PTR center[rip]    # load center
    cmp     eax, edx
    jbe     .L4                           # conditionally skip the store
    mov     DWORD PTR center[rip], eax    # Update the value in memory
    mov     edx, eax                      # and the value in a register for use after the loop.
.L4:

    add     r12d, 1                        # ++i
    cmp     DWORD PTR [rbp-44], r12d  # n, i
    ja      .L14                           # } while(n>i)

From the bottom of the loop, we can see that since the address of n was passed to a previous cin >>, the compiler now has to assume that any call to an unknown function might modify n. (e.g. it can't prove that std::iostream::operator>> didn't store &n in a global somewhere). You could let the compiler keep n in a register by using cin >> tmp; n = tmp;. (It won't make a measurable performance difference, though).

---

center does stay in a register for the O(N*M) part of your code, though, so the extra cost is only O(N).

It's probably trivial compared to the cost of the cin >> integer parsing, even though it recomputes center*center for every y value. (Hoisting that out of the loop in the source would have been good.)

Loads that hit in L1 cache are very cheap, but the store/reload in the first loop that updates center while reading introduces unnecessary extra latency into the dependency chain involving the value of center. Recomputing r_squared isn't great, either.

---

If you had used static unsigned int center = 0;, and made static inline distance_squared, the compiler could optimize away the static storage and always keep it in a register. (gcc and clang don't, but I'm not sure why.)

static means the scope is limited to this compilation unit, so the compiler know that functions in other files can't access it directly.
So the compiler knows that only code within the current compilation unit can see it.

static distance_squared() lets the compiler know that it can't be called from outside the compilation unit. (It's not in a .h, but remember that only the preprocessor knows about headers vs. .cpp. inline double distance_squared declares a function that can be called from anywhere, including (for all the compiler knows) from cin >>), so the value in memory has to be in sync. It does let it avoid reloading in the second outer loop (since distance_squared only reads it), but it does still force a store/reload in the first loop.

Making center a local variable in main and passing it as a 3rd function arg would totally solve that problem, and be better style anyway. Or as mentioned above, precalculate x-center squared so the function is split up.

---

Using if (cond) { ++count; continue; } is clumsy compared a normal if () / else if() chain:

    if (distance_sqr == rsqr)
        ++on_circut;    // mis-spelled variable name?
    else if(distance_sqr > rsqr)
        ++outside;
    else 
        ++inside;

Yes this works even without {} if you want. But it's probably best to include the {} once you have a chain of if / else if. Too easy to get it wrong when changing something. I think this style is nice:

    if (distance_sqr == rsqr) {
        ++on_border;
    } else if(distance_sqr > rsqr) {
        ++outside;
    } else {
        ++inside;
    }

Less compact than jamming the controlled statement on the same line as the if, but more readable.

Even if you normally use a wide text window so hscrolling isn't needed, it's not great to put too much on one line.

---

The 3-way branch to decide which counter to increment is also more work than needed.

on_circut + outside + inside == m * n, so you can calculate one of the three from m*n and the other two at the end of the loop. gcc and clang don't do this for you. :/

gcc also actually branches. It is likely cheaper on current x86 CPUs to do it branchlessly, like clang does. Reducing it to two counters is a bigger win for branchless than for branching.

center = center >> 1; //div by 2
uint64_t rsqr = center * static_cast<uint64_t>(center);

std::cin >> m;

// 64-bit counters help compilers auto-vectorize incrementing them with 64-bit packed compare results
uint64_t on_circle = 0, outside = 0, inside = 0;
while(m--){          // reading elements of Y one by one
    std::cin >> tmp;
    for(unsigned i = 0; i < n; ++i){ // for each pair of tmp (being Y coordinate) and a number from Xs (being coord. X)...
        int64_t distance_sqr = distance_squared(set_X[i], tmp, center); // find the distance^2 and compare it to radius^2
        outside += (distance_sqr > rsqr);   // hint the compiler towards making branchless code
        on_circle += (distance_sqr == rsqr);
//        if(distance_sqr > rsqr) ++outside;
//        else if(distance_sqr == rsqr){ ++on_circle; }
    }
}
// infer the last counter from the other two
// saving 1/3 of the work if the compiler used branchless code for the conditions
inside = n * m - on_circle - outside;

clang does very well with this: clang4.0 -O3 for x86-64 emits:

    ; 2nd half of clang's unrolled loop
    movsxd  rdi, dword ptr [r13 + 4*rcx + 4]    # load an x value (sign-extending from 32 to 64b for free).
    sub     rdi, r15                            # dx = x - center
    imul    rdi, rdi                            # dx*dx
    add     rdi, r8                             # dx*dx + dy2(hoisted by the compiler
    xor     r14d, r14d                          # zero two registers to set up for setcc
    xor     ebx, ebx
    cmp     rdi, r12                            # compare once
    seta    r14b                                # r14 = (distance_sqr > rsqr)  unsigned Above
    sete    bl                                  # rbx = (distance_sqr == rsqr)
    add     r14, rax                           # increment the counters by 0 or 1, depending on the compare results.
    add     rbx, rdx

    add     rcx, 2                              # unrolled loop: advance by 2 elements
    cmp     rcx, rsi
    jb      .LBB0_13                            # while(pointer < end_pointer)

gcc7.1 is a lot more clunky: It sign-extends after the subtract. It fails to reuse the flag results from one cmp. But that's not something you can control from the source. (Except by not encouraging the compiler to make branchless code.) Well maybe you could try assigning both compare results to variables before incrementing, because the add instruction clobbers flags. Anyway, that's a compiler missed-optimization, and unless you're tuning specifically for gcc, you might just ignore it.

---

With -msse4 or -mavx2, gcc7.1 and clang4.0 can auto-vectorize pretty well. (They need signed 32x32 => 64-bit multiply (pmuldq), but SSE2 (baseline for x86-64) only has the unsigned version of that (pmuludq)). This wouldn't be an issue if we'd done the dx*dx transformation once before the nested loop, but let's talk about the code without that optimization.

The careful choice of types (int32_t for the subtract) allows that part to auto-vectorize with 32-bit elements, before unpacking to 64-bit. You might consider making a separate version of the function that uses only 32-bit integers if the input numbers happen to be small enough. That would let it vectorize with twice as many elements per vector, or with just SSE2 if the upper half of the vector can be ignored.

clang and gcc can make a mess when auto-vectorizing if you aren't careful. Making the counters 64-bit actually helps a lot (because that's the same width at the compare). Without that, they spend time packing down to 32-bit inside the inner loop. (At least clang does, didn't check gcc for this).

Have a look on the Godbolt compiler explorer for my version, with minimal algorithmic changes (not even optimizing out the dx*dx), but with a version of this that auto-vectorizes fairly nicely, especially with gcc. Doing a separate set_X[i] = square(set_X[i] - center) would be a big win for the auto-vectorized version, too.

---

Pre-declaring all your variables at the top of a function is bad style. Only do that if you are forced to for compatibility with a C89 compiler or something. Since you're writing in C++, just don't.

---

Putting your main computation in a function not called main is usually a good thing. gcc treats main as a "cold" function and optimizes it somewhat less than other functions.