1

I've been looking frantically for an answer to no avail; I've created a quadtree which is supposed to sort an array of more than 1000 objects declared by the struct:

typedef struct node
{
    char     is_leaf;
    struct Particle *p;
    double    m;

    double center_x;
    double center_y;
    double width;

    struct node *sw;
    struct node *nw;
    struct node *se;
    struct node *ne;

} node;

into a quadtree with the function:

node* quadtree_insert(node *n, struct Particle *p, double center_x, double center_y, double width)
{
    if(n == NULL)
    {
        n = (node*)malloc(sizeof(node));
        n->is_leaf = 1;

        n->p = p;
        //n->m = 0.1;

        n->sw = NULL;
        n->se = NULL;
        n->nw = NULL;
        n->ne = NULL;
        if(width < 1e-300){
            n->width = 1e-300;
            }
        else
            n->width    = width;
        return n;
    }
    else{
        //n = (node*)malloc(sizeof(node));
        double x;
        double y;
        if(width < 1e-300){
            width = 1e-300;
            }
        if(n->is_leaf == 1) //! that is, if the node is not a branch
        {
                        x = (double)n->p->x_pos;
            y = (double)n->p->y_pos;

            if(x <= center_x && y <= center_y) //! first quadrant
            {
                n->sw = quadtree_insert(n->sw, n->p, center_x * 0.5, center_y * 0.5, width * 0.5);
            }
            else if(x <= center_x && y > center_y) //! second quadrant
            {
                n->nw = quadtree_insert(n->nw, n->p, center_x * 0.5, center_y + center_y * 0.5, width * 0.5);
            }
            else if(x > center_x && y <= center_y) //! third quadrant
            {
                n->se = quadtree_insert(n->se, n->p, center_x + center_x * 0.5, center_y * 0.5, width * 0.5);
            }
            else //! fourth quadrant
            {
                n->ne = quadtree_insert(n->ne, n->p, center_x + center_x * 0.5, center_y + center_y * 0.5, width * 0.5);
            }

            n->p = NULL; //! sets branch pointer to nothing...
            n->is_leaf = 0;
                    }
        //}

        x = (double)p->x_pos;
        y = (double)p->y_pos;

        if(x <= center_x && y <= center_y) //! first quadrant
        {
            n->sw = quadtree_insert(n->sw, p, center_x * 0.5, center_y * 0.5, width * 0.5);

        }
        else if(x <= center_x && y > center_y) //! second quadrant
        {
            n->nw = quadtree_insert(n->nw, p, center_x * 0.5, center_y + center_y * 0.5, width * 0.5);

        }
        else if(x > center_x && y <= center_y) //! third quadrant
        {
            n->se = quadtree_insert(n->se, p, center_x + center_x * 0.5, center_y * 0.5, width * 0.5);

        }
        else //! fourth quadrant
        {
            n->ne = quadtree_insert(n->ne, p, center_x + center_x * 0.5, center_y + center_y * 0.5, width * 0.5);            
        }
        return n;
    }
}

This is all done by:

 node *root = NULL;
  root = quadtree_insert(root, &particles[0],0.500,0.5,1);

  for(i = 1; i < nParticles; i++)
  {
    quadtree_insert(root, &particles[i],0.5000,0.5,1);
  }

where "particles" is passed on on the form struct Particle* particles. And a particle is defined such that:

struct Particle {
    double mass;
    double x_pos;
    double y_pos;
    double x_vel;
    double y_vel;
};

typedef struct Particle * Particle_structure;

The root is free'd after every iteration, after the for-loop, and the code works for samples smaller than ~200 and valgrind gives no errors for these. Added to this is a middle function, which performs some arithmetic on the particles:

double quadtree_calculate_forcey(struct Particle *p, node *n, double theta_max, double delta_t, int numP, double epsilon, double G)
{
    if(n != NULL)
    {
            double d_x      = (n->center_x - p->x_pos);
                double d_y      = (n->center_y - p->y_pos);
        double r_2     = d_x * d_x + d_y * d_y;
        r_2 = sqrt(r_2) + epsilon;
        if(theta_max <= (n->width / r_2) && !n->is_leaf){
                double a = 0;
            if(n->sw != NULL)
                        a += quadtree_calculate_forcey(p, n->sw, theta_max, delta_t, numP, epsilon, G);
                    if(n->nw != NULL)
                        a += quadtree_calculate_forcey(p, n->nw, theta_max, delta_t, numP, epsilon, G);
                    if(n->se != NULL)
                        a += quadtree_calculate_forcey(p, n->se, theta_max, delta_t, numP, epsilon, G);
                    if(n->ne != NULL)
                        a += quadtree_calculate_forcey(p, n->ne, theta_max, delta_t, numP, epsilon, G);
                    return a;
        }
        else    
                {
                double fy;
            double mass;
           if(d_x == 0 && d_y == 0){ // could be comparing the same star
                 //printf("RÖÖÖVHATT\n");  
                  return 0;
            }
            else{
            //printf("MASS : %f\n", n->m);
                  mass = n->m;
                  //printf("MASS : %f\n", mass);
                  fy = G * (mass * p->mass/ pow(r_2,3)) * d_y;   
                            //printf("DY:%f   DX:%f   R_2:%f   MASSA:%f\n",d_y, d_x, r_2 - epsilon, mass);          
                     // printf("HIT SKA JAG: %f\n",d_y);
                  return fy;              
                 }
        }
    }
    return 0.0;
}

The ringer is that its rolls for a bit (clearing the root, and redoing it for new positions), so the number of variables in the recursion surely cannot be the problem(?). I'm pretty sure I'm out swimming with the fishes when it comes to some allocation of pointer declaration. Any thoughts/help would put you on the good side of Odin!

Edit: An example from how we find mass-center etc. The node-mass is done the same way.

double set_centerx(node *n)
{
    if(n != NULL)
    {
        if(!(n->is_leaf))
        {
                    double a = set_centerx(n->ne);
                    double b = set_centerx(n->nw);
                    double c = set_centerx(n->se);
                    double d = set_centerx(n->sw);
                    double m1 = get_mass2(n->ne);
                    double m2 = get_mass2(n->nw);
                    double m3 = get_mass2(n->se);
                    double m4 = get_mass2(n->sw);
          n->center_x = (double)(m1*a + m2*b + m3*c + m4*d)/(m1+m2+m3+m4);
          return n->center_x;
        }
                n->center_x =n->p->x_pos;
        return n->p->x_pos;
    }

    return 0;
}
radnaskela
  • 13
  • 3
  • I had similar situations solved by changing my base condition. If you have a stackoverflow in your recursion, it could be because of your base condition. If I were you, I'd check my base condition and debug which conditions the code reaches. – Niklas Rosencrantz Jun 01 '16 at 20:00
  • In the `if(n == NULL)` case you never initiailize `n->center_x` or `n->center_y` or `n->m` in the new node, but your `quadtree_calculate_forcey` function seems to read and compare those values – M.M Jun 01 '16 at 21:06
  • Its done in other functions before the call to forcey, since center is supposed to be mass-center we need the particles to be in the correct places. An example is above. – radnaskela Jun 01 '16 at 21:52
  • Hmmm...did you forget to specify that the positions of the particles are inside the unit square (bottom-left corner at the origin, top-right at (1,1))? That's the trouble with magic numbers; there's nothing to explain what they mean. – Jonathan Leffler Jun 01 '16 at 22:56

2 Answers2

3

Partial analysis

I'm tolerably certain the trouble is primarily due to the bulk repeated code (near to repeated code) in the primary else clause of quadtree_insert(). I've marked the starts with comments saying Fragment 1A and Fragment 1B — and Fragment 1B is now also framed with #ifdef DO_REPEAT and #endif.

    else
    {
        /* Fragment 1A */
        // n = (node*)malloc(sizeof(node));
        double x;
        double y;
        if (width < 1e-300)
        {
            width = 1e-300;
        }
        if (n->is_leaf == 1) // ! that is, if the node is not a branch
        {
            x = (double)n->p->x_pos;
            y = (double)n->p->y_pos;

            if (x <= center_x && y <= center_y) // ! first quadrant
            {
                n->sw = quadtree_insert(n->sw, n->p, center_x * 0.5, center_y * 0.5, width * 0.5);
            }
            else if (x <= center_x && y > center_y) // ! second quadrant
            {
                n->nw = quadtree_insert(n->nw, n->p, center_x * 0.5, center_y + center_y * 0.5, width * 0.5);
            }
            else if (x > center_x && y <= center_y) // ! third quadrant
            {
                n->se = quadtree_insert(n->se, n->p, center_x + center_x * 0.5, center_y * 0.5, width * 0.5);
            }
            else // ! fourth quadrant
            {
                n->ne = quadtree_insert(n->ne, n->p, center_x + center_x * 0.5, center_y + center_y * 0.5, width * 0.5);
            }

            n->p = NULL; // ! sets branch pointer to nothing...
            n->is_leaf = 0;
        }

#ifdef DO_REPEAT
        /* Fragment 1B */
        x = (double)p->x_pos;
        y = (double)p->y_pos;

        if (x <= center_x && y <= center_y) // ! first quadrant
        {
            n->sw = quadtree_insert(n->sw, p, center_x * 0.5, center_y * 0.5, width * 0.5);
        }
        else if (x <= center_x && y > center_y) // ! second quadrant
        {
            n->nw = quadtree_insert(n->nw, p, center_x * 0.5, center_y + center_y * 0.5, width * 0.5);
        }
        else if (x > center_x && y <= center_y) // ! third quadrant
        {
            n->se = quadtree_insert(n->se, p, center_x + center_x * 0.5, center_y * 0.5, width * 0.5);
        }
        else // ! fourth quadrant
        {
            n->ne = quadtree_insert(n->ne, p, center_x + center_x * 0.5, center_y + center_y * 0.5, width * 0.5);
        }
#endif /* DO_REPEAT */
        return n;
    }

I took your code and reordered some of the fragments, and used this main(). Note that I should not have needed to do this; you should have produced an MCVE (How to create a Minimal, Complete, and Verifiable Example?) or SSCCE (Short, Self-Contained, Correct Example) — two names and links for the same basic idea.

static struct Particle particles[] =
{
    {  19.99,  96.07,  62.79, -99.46,  19.70 },
    {  12.94,   1.43, -33.45,  31.80, -66.08 },
    {   6.49,  16.99, -20.83,  92.51,  35.98 },
    {  17.01, -28.85, -94.10,  42.82,  -1.30 },
    {  14.27,  85.07,  88.21,  11.22,  16.85 },
    {  15.73, -56.37,  46.85,  27.40, -15.15 },
    {   1.53, -49.44, -64.27, -29.45, -38.25 },
    {   8.03,  92.11, -47.50,  63.77, -29.99 },
    {   8.67, -99.81,  73.19,  18.75,  88.66 },
    {  16.36,  66.33,  14.23,  87.65,  40.01 },
};

enum { nParticles = sizeof(particles) / sizeof(particles[0]) };

int main(void)
{
    node *root = NULL;
    printf("Particle 0:\n");
    root = quadtree_insert(root, &particles[0], 0.500, 0.5, 1);

    for (int i = 1; i < nParticles; i++)
    {
        printf("Particle %d:\n", i);
        quadtree_insert(root, &particles[i], 0.5000, 0.5, 1);
    }
    return 0;
}

I used a random number generator for the table of values:

random -n 10 -T '    { %6:2[1:20]f, %6:2[-100:100]f, %6:2[-100:100]f, %6:2[-100:100]f, %6:2[-100:100]f },'

It crashed for me on Particle 6. In outline, valgrind said:

…startup blurb omitted…
Particle 0:
Particle 1:
Particle 2:
Particle 3:
Particle 4:
Particle 5:
Particle 6:
==79528== 
==79528== Process terminating with default action of signal 11 (SIGSEGV)
==79528==  Access not within mapped region at address 0x104003FD0
==79528==    at 0x100008E39: malloc (vg_replace_malloc.c:302)
==79528==  If you believe this happened as a result of a stack
==79528==  overflow in your program's main thread (unlikely but
==79528==  possible), you can try to increase the size of the
==79528==  main thread stack using the --main-stacksize= flag.
==79528==  The main thread stack size used in this run was 8388608.
==79528== 
==79528== HEAP SUMMARY:
==79528==     in use at exit: 6,017,700 bytes in 75,080 blocks
==79528==   total heap usage: 75,162 allocs, 82 frees, 6,023,876 bytes allocated
==79528== 
==79528== LEAK SUMMARY:
==79528==    definitely lost: 4,120 bytes in 2 blocks
==79528==    indirectly lost: 2,288 bytes in 6 blocks
==79528==      possibly lost: 4,880 bytes in 45 blocks
==79528==    still reachable: 5,997,720 bytes in 74,904 blocks
==79528==         suppressed: 8,692 bytes in 123 blocks
==79528== Rerun with --leak-check=full to see details of leaked memory

Even on a Mac where I ran it, that indicates problems; the suppressed stuff is OK, but the rest is mostly not.

When compiled without -DDO_REPEAT (a normal compilation), the sample program runs to completion. It leaks, of course, because there's no code to free the memory.

Particle 0:
Particle 1:
Particle 2:
Particle 3:
Particle 4:
Particle 5:
Particle 6:
Particle 7:
Particle 8:
Particle 9:
==79683== 
==79683== HEAP SUMMARY:
==79683==     in use at exit: 26,580 bytes in 191 blocks
==79683==   total heap usage: 273 allocs, 82 frees, 32,756 bytes allocated
==79683== 
==79683== LEAK SUMMARY:
==79683==    definitely lost: 4,200 bytes in 3 blocks
==79683==    indirectly lost: 2,368 bytes in 7 blocks
==79683==      possibly lost: 4,880 bytes in 45 blocks
==79683==    still reachable: 6,440 bytes in 13 blocks
==79683==         suppressed: 8,692 bytes in 123 blocks
==79683== Rerun with --leak-check=full to see details of leaked memory

Note that there's dramatically less memory in use than before.

If the code eliminated by the absence of DO_REPEAT is in fact crucial, then the bug either lies in that code, or lies in the setup work done in Fragment 1A. However, inserting the same particle twice with the recursive calls seems likely to me to be the source of the trouble.

I also note that the quadtree_calculate_forcey() function isn't used in the code; it is definitely not a part of the MCVE.


Fragment 1B is needed

John Bollinger suggested:

Do note that the "repeated" code is the same in form, but not in detail. That is, fragment 1B uses n->p where fragment 1A uses p. I think this is purposeful: the idea seems to be to force all data (the Particles) out to leaf nodes.

and radnaskela affirmed:

It's exactly like John says, every particle should occupy one end-node. My suspicions say its to do with the movement, and planets aligning so that a ridiculous number of iterations have to be done to open two free squares. I'm unable to see a solution to my problem though, any thoughts?

There is a problem with the logic, though I'm not entirely sure what that problem is. The first step I'd do is to add some instrumentation to the code in quadtree_insert(), like this:

static node *quadtree_insert(node *n, struct Particle *p, double center_x, double center_y, double width)
{
    printf("Centre (X,Y) = (%6.2f,%6.2f)\n", center_x, center_y);
    if (n == NULL)
    {
        n = (node *)malloc(sizeof(node));
        n->is_leaf = 1;

        n->p = p;
        // n->m = 0.1;

        n->sw = NULL;
        n->se = NULL;
        n->nw = NULL;
        n->ne = NULL;
        if (width < 1e-300)
        {
            n->width = 1e-300;
        }
        else
            n->width    = width;
        return n;
    }
    else
    {
        // n = (node*)malloc(sizeof(node));
        double x;
        double y;
        if (width < 1e-300)
        {
            width = 1e-300;
        }
        if (n->is_leaf == 1) // ! that is, if the node is not a branch
        {
            x = (double)n->p->x_pos;
            y = (double)n->p->y_pos;

            if (x <= center_x && y <= center_y) // ! first quadrant
            {
                printf("Recurse SW 1: ");
                n->sw = quadtree_insert(n->sw, n->p, center_x * 0.5, center_y * 0.5, width * 0.5);
            }
            else if (x <= center_x && y > center_y) // ! second quadrant
            {
                printf("Recurse NW 1: ");
                n->nw = quadtree_insert(n->nw, n->p, center_x * 0.5, center_y + center_y * 0.5, width * 0.5);
            }
            else if (x > center_x && y <= center_y) // ! third quadrant
            {
                printf("Recurse SE 1: ");
                n->se = quadtree_insert(n->se, n->p, center_x + center_x * 0.5, center_y * 0.5, width * 0.5);
            }
            else // ! fourth quadrant
            {
                printf("Recurse NE 1: ");
                n->ne = quadtree_insert(n->ne, n->p, center_x + center_x * 0.5, center_y + center_y * 0.5, width * 0.5);
            }

            n->p = NULL; // ! sets branch pointer to nothing...
            n->is_leaf = 0;
        }

        x = (double)p->x_pos;
        y = (double)p->y_pos;

        if (x <= center_x && y <= center_y) // ! first quadrant
        {
            printf("Recurse SW 2: ");
            n->sw = quadtree_insert(n->sw, p, center_x * 0.5, center_y * 0.5, width * 0.5);
        }
        else if (x <= center_x && y > center_y) // ! second quadrant
        {
            printf("Recurse NW 2: ");
            n->nw = quadtree_insert(n->nw, p, center_x * 0.5, center_y + center_y * 0.5, width * 0.5);
        }
        else if (x > center_x && y <= center_y) // ! third quadrant
        {
            printf("Recurse SE 2: ");
            n->se = quadtree_insert(n->se, p, center_x + center_x * 0.5, center_y * 0.5, width * 0.5);
        }
        else // ! fourth quadrant
        {
            printf("Recurse NE 2: ");
            n->ne = quadtree_insert(n->ne, p, center_x + center_x * 0.5, center_y + center_y * 0.5, width * 0.5);
        }
        return n;
    }
}

When run with the 10-point data set, the output up to Particle 5 is:

Particle 0: Centre (X,Y) = (  0.50,  0.50)
Particle 1: Centre (X,Y) = (  0.50,  0.50)
Recurse NE 1: Centre (X,Y) = (  0.75,  0.75)
Recurse SE 2: Centre (X,Y) = (  0.75,  0.25)
Particle 2: Centre (X,Y) = (  0.50,  0.50)
Recurse SE 2: Centre (X,Y) = (  0.75,  0.25)
Recurse SE 1: Centre (X,Y) = (  1.12,  0.12)
Recurse SE 2: Centre (X,Y) = (  1.12,  0.12)
Recurse SE 1: Centre (X,Y) = (  1.69,  0.06)
Recurse SE 2: Centre (X,Y) = (  1.69,  0.06)
Recurse SW 1: Centre (X,Y) = (  0.84,  0.03)
Recurse SE 2: Centre (X,Y) = (  2.53,  0.03)
Particle 3: Centre (X,Y) = (  0.50,  0.50)
Recurse SW 2: Centre (X,Y) = (  0.25,  0.25)
Particle 4: Centre (X,Y) = (  0.50,  0.50)
Recurse NE 2: Centre (X,Y) = (  0.75,  0.75)
Recurse NE 1: Centre (X,Y) = (  1.12,  1.12)
Recurse NE 2: Centre (X,Y) = (  1.12,  1.12)
Recurse NE 1: Centre (X,Y) = (  1.69,  1.69)
Recurse NE 2: Centre (X,Y) = (  1.69,  1.69)
Recurse NE 1: Centre (X,Y) = (  2.53,  2.53)
Recurse NE 2: Centre (X,Y) = (  2.53,  2.53)
Recurse NE 1: Centre (X,Y) = (  3.80,  3.80)
Recurse NE 2: Centre (X,Y) = (  3.80,  3.80)
Recurse NE 1: Centre (X,Y) = (  5.70,  5.70)
Recurse NE 2: Centre (X,Y) = (  5.70,  5.70)
Recurse NE 1: Centre (X,Y) = (  8.54,  8.54)
Recurse NE 2: Centre (X,Y) = (  8.54,  8.54)
Recurse NE 1: Centre (X,Y) = ( 12.81, 12.81)
Recurse NE 2: Centre (X,Y) = ( 12.81, 12.81)
Recurse NE 1: Centre (X,Y) = ( 19.22, 19.22)
Recurse NE 2: Centre (X,Y) = ( 19.22, 19.22)
Recurse NE 1: Centre (X,Y) = ( 28.83, 28.83)
Recurse NE 2: Centre (X,Y) = ( 28.83, 28.83)
Recurse NE 1: Centre (X,Y) = ( 43.25, 43.25)
Recurse NE 2: Centre (X,Y) = ( 43.25, 43.25)
Recurse NE 1: Centre (X,Y) = ( 64.87, 64.87)
Recurse NE 2: Centre (X,Y) = ( 64.87, 64.87)
Recurse SE 1: Centre (X,Y) = ( 97.31, 32.44)
Recurse NE 2: Centre (X,Y) = ( 97.31, 97.31)
Particle 5: Centre (X,Y) = (  0.50,  0.50)
Recurse NW 2: Centre (X,Y) = (  0.25,  0.75)

I'm a bit surprised by some of the processing there, and the number of entries for Particle 4. That may be indicative of the problem.

And then it processes Particle 6:

Particle 6: Centre (X,Y) = (  0.50,  0.50)
Recurse SW 2: Centre (X,Y) = (  0.25,  0.25)
Recurse SW 1: Centre (X,Y) = (  0.12,  0.12)
Recurse SW 2: Centre (X,Y) = (  0.12,  0.12)
Recurse SW 1: Centre (X,Y) = (  0.06,  0.06)
Recurse SW 2: Centre (X,Y) = (  0.06,  0.06)
Recurse SW 1: Centre (X,Y) = (  0.03,  0.03)
Recurse SW 2: Centre (X,Y) = (  0.03,  0.03)
Recurse SW 1: Centre (X,Y) = (  0.02,  0.02)
Recurse SW 2: Centre (X,Y) = (  0.02,  0.02)
Recurse SW 1: Centre (X,Y) = (  0.01,  0.01)
Recurse SW 2: Centre (X,Y) = (  0.01,  0.01)
Recurse SW 1: Centre (X,Y) = (  0.00,  0.00)
Recurse SW 2: Centre (X,Y) = (  0.00,  0.00)
Recurse SW 1: Centre (X,Y) = (  0.00,  0.00)
Recurse SW 2: Centre (X,Y) = (  0.00,  0.00)
Recurse SW 1: Centre (X,Y) = (  0.00,  0.00)
Recurse SW 2: Centre (X,Y) = (  0.00,  0.00)

It gets very tedious from there onwards. The question that you need to ask is "what should be happening with that point"? It would probably be a good idea to have a function to dump the quadtree structure, too.

This basically suggests that you have not yet analyzed the conditions under which the points are to be added sufficiently. I'm not clear why you multiply the center coordinates by 1.5 when the insertion is in the SE or NE quadrants and by 0.5 when in the SW or NW quadrants (nor why you use the add plus multiply rather than simply multiply).

The test on width less than 1e-300 is mildly worrying. Given that you call the function with the value 1.0, it takes a while to get that small when you halve the width each time you recurse.

Better tracing might well report on the values of x and y as well as the centre coordinate.

Community
  • 1
  • 1
Jonathan Leffler
  • 730,956
  • 141
  • 904
  • 1,278
  • Do note that the "repeated" code is the same in *form*, but not in detail. That is, fragment 1B uses `n->p` where fragment 1A uses `p`. I think this is purposeful: the idea seems to be to force all data (the `Particle`s) out to leaf nodes. – John Bollinger Jun 01 '16 at 21:23
  • @JohnBollinger: you might be right, but then the bug is in that code. I'm not sufficiently sure what the quadtree should look like — there's no information in the question. The two lots of recursion for the same particle looks kinda odd. I've not established enough diagnostic tracing to tell what's going on in the crash case versus the non-crash case. The problem isn't "out of heap space", there's only 6 MiB allocated when the crash occurs. – Jonathan Leffler Jun 01 '16 at 21:28
  • That piece of advice about the code was brilliant, next time ill make sure to follow through! Its exactly like John says, every particle should occupy one end-node. My suspicions say its to do with the movement, and planets aligning so that a ridiculous number of iterations have to be done to open two free squares. I'm unable to see a solution to my problem though, any thoughts? – radnaskela Jun 01 '16 at 21:50
  • @JonathanLeffler, I'm learning quadtree on the fly, but I think I've identified the problem. It is, of course, related to the near-duplicate code, just as you say. – John Bollinger Jun 01 '16 at 22:31
  • @radnaskela: See the updated material. Use my data set; it goes wrong quickly (you said yours takes several hundred entries to go wrong). – Jonathan Leffler Jun 01 '16 at 22:31
  • @JohnBollinger: I've not yet gone to study quadtree theory, so you're a step ahead of me. To a large extent, I think I've identified where the problem is — how it should be resolved is mostly radnaskela's problem. _[…time passes…]_ Wikipedia on [quadtree](https://en.wikipedia.org/wiki/Quadtree) — I think this is a point quadtree, isn't it? – Jonathan Leffler Jun 01 '16 at 22:33
  • Thank you for the help, Jonathan! Im ridiculously impressed by your ability to break down the problem, hope all your future endeavours succeed. I can also add that you two have taught me more than the past month of random reading about C! – radnaskela Jun 01 '16 at 23:24
  • @radnaskela: You're welcome. More than anything else, I hope you learned a bit about how to go about debugging such problems. Granted, `printf()` statements aren't the only way to debug code — another is to use a debugger — but they do have merits. Tracking what is going on by judiciously placed print operations can really help. (See also [C `#define` macro for debug printing](https://stackoverflow.com/questions/1644868).) But since John Bollinger did the grunt-work of getting the code to work — rather than just diagnosing where the problem is — you're absolutely correct to accept his answer. – Jonathan Leffler Jun 01 '16 at 23:27
2

I see several problems with your code.


In the first place, you have a near-duplication of a sizable block of non-trivial code. The two blocks seem to differ only in which Particle they reference; this begs to be factored out into a helper function.


In the second place, consider this fragment:

else if (x <= center_x && y > center_y) // ! second quadrant
{
    n->nw = quadtree_insert(n->nw, n->p, center_x * 0.5, center_y + center_y * 0.5, width * 0.5);
}

I take the input center_x and center_y to be the center of a square patch of the overall area, and width to be the edge length for that patch. In that case, the center of the northwest quadrant of the patch is (center_x - width * 0.25, center_y + width * 0.25). Only in certain special cases will that coincide with the computation in your recursive call. You actually can determine the width of the patch from the coordinates of its center, but not so directly as you're trying to do.

Similar applies to all seven of the other recursive calls.


In the third place, consider what happens if you end up with two particles that have very similar -- or worse, identical -- coordinates. They can be the only two in the simulation.

The first particle initially gets positioned in a node representing the whole simulation area. When the second is inserted, the quadtree is traversed to the same node. The original particle is therefore moved to a newly created quadrant of the original node's area, and the insertion process for second node resumes. But the insertion then again reaches the same node that the original particle occupies, and the process repeats. And maybe it repeats again. And again.

There is therefore no definite bound to your recursion. If you happen to end up with two particles with identical coordinates -- for example, maybe they are pinned to the corner of the simulation area -- then the you will definitely recurse endlessly, eventually producing a stack overflow. Even if no two particles have identical coordinates, however, it is conceivable that two will be close enough together for the recursion to go deep enough for a stack overflow.

I suspect that this problem is aggravated by the previous issue, but it does not depend on that issue.

One way to approach this problem would be to merge particles that come close enough together. That would allow you to place a firm bound on the recursion depth, based on the patch width. Alternatively, you may abort if recursion goes too deep. Or, you could just let it crash in that case, as it already does.


For what it's worth, here's how I think your quadtree_insert() ought to look:

// adjust as needed:
#define MIN_PARTICLE_SEPARATION 1e-300

void quadtree_insert_helper(node *n, Particle *p, double center_x,
        double center_y, double width) {
    width *= 0.5;

    double new_center_x = center_x
            + width * ((p->x_pos <= center_x) ? -0.5 : 0.5);
    double new_center_y = center_y
            + width * ((p->y_pos <= center_y) ? -0.5 : 0.5);
    node **quad;

    if (new_center_x <= center_x && new_center_y <= center_y) {
        //! first quadrant
        quad = &n->sw;
    } else if (new_center_x <= center_x && new_center_y > center_y) {
        //! second quadrant
        quad = &n->nw;
    } else if (new_center_x > center_x && new_center_y <= center_y) {
        //! third quadrant
        quad = &n->se;
    } else {
        //! fourth quadrant
        quad = &n->ne;
    }
    *quad = quadtree_insert(*quad, p, new_center_x, new_center_y, width);
}

node* quadtree_insert(node *n, struct Particle *p, double center_x,
        double center_y, double width) {
    if (n == NULL) {
        n = malloc(sizeof(*n));
        n->is_leaf = 1;
        n->p = p;
        //n->m = 0.1;
        n->sw = NULL;
        n->se = NULL;
        n->nw = NULL;
        n->ne = NULL;
        n->center_x = center_x;
        n->center_y = center_y;
        n->width = width;
    } else if (n->is_leaf && (with <= MIN_PARTICLE_SEPARATION)) {
        // ... merge particles ...
    } else {
        if (n->is_leaf) { //! that is, if the node is not a branch
            quadtree_insert_helper(n, n->p, center_x, center_y, width);
            //! the node is now a branch
            n->p = NULL;
            n->is_leaf = 0;
        }

        quadtree_insert_helper(n, p, center_x, center_y, width);
    }

    return n;
}
John Bollinger
  • 160,171
  • 8
  • 81
  • 157
  • You, my dear sir, will be knighted if I ever get the pleasure to be king of something. The moral of this story is that years of preparation cannot prepare me for the simplest things. I changed the conditions for centers and voilá, like the sweetest rain of tasty barley from clear blue sky, shit just worked! If I ever get the pleasure of meeting you I'll probably kiss you, guess perhaps we'll both rue the day but still. All your points have been taken, and had been considered. Gut feeling told me it was something sillier, just didn't realize how silly. – radnaskela Jun 01 '16 at 23:22