GMP function categories

There are several categories of functions in GMP:

    High-level signed integer arithmetic functions (mpz). There are about 140 arithmetic and logic functions in this category.
    High-level rational arithmetic functions (mpq). This category consists of about 35 functions, but all signed integer arithmetic functions can be used too, by applying them to the numerator and denominator separately.
    High-level floating-point arithmetic functions (mpf). This is the GMP function category to use if the C type `double' doesn't give enough precision for an application. There are about 65 functions in this category.
    C++ class based interface to all of the above. (The C functions and types can of course be used directly from C++ too.)

1. Code *someone* writes as a perfect language implementation
2. A library (whether compiled or not) provided by *someone*
3. Non-standard extensions

typedef float v4sf __attribute__((vector_size (16)));

v4sf multiply_add(const v4sf a, const v4sf b, const v4sf c)
{
    return a * b + c;
}

int distances(const unsigned int   neighbors, /* Number of neighbors */
              unsigned int *const  neighbor,  /* Neighbor atom indexes */
              double *const        distance,  /* Consecutive output array: distance */
              double *const        deltax,    /* Consecutive output array: unit vector x components */
              double *const        deltay,    /* Consecutive output array: unit vector y components */
              double *const        deltaz,    /* Consecutive output array: unit vector z components */
              const double *const  x,         /* Atoms' x coordinates */
              const double *const  y,         /* Atoms' y coordinates */
              const double *const  z,         /* Atoms' z coordinates */
              const unsigned int   atom,      /* Index for the one atom */
              const double         limit)     /* Distance limit, "pivot" */
{
    {
        const double    x0 = x[atom],  /* This atom coordinates; */
                        y0 = y[atom],  /* subtract from the other */
                        z0 = z[atom];  /* atom to get delta vector. */
        int             i;

        for (i = 0; i < neighbors; i++) {
            const unsigned int n = neighbor[i];
            const double    dx = x[n] - x0,   /* Delta vector components */
                            dy = y[n] - y0,
                            dz = z[n] - z0;
            const double    dd = dx*dx + dy*dy + dz*dz;   /* Distance squared */
            const double    d = sqrt(dd);     /* Distance */

            /* Save unit vector, */
            deltax[i] = dx / d;
            deltay[i] = dy / d;
            deltaz[i] = dz / d;

            /* and distance. */
            distance[i] = d;
        }
    }

    /* The distance[], deltax[], deltay[], deltaz[] are guaranteed to
     * have at least one extra element, and padded to a multiple of 4 elements.
     * Use a sentinel to avoid ihead < itail checks.
    */
    distance[neighbors] = limit;

    {
        unsigned int    ihead = 0;
        unsigned int    itail = neighbors - 1;

        while (1) {

            /* Keep close neighbors at start of list */
            while (distance[ihead] < limit)
                ihead++;

            /* All under limit? */
            if (ihead >= itail)
                return ihead;

            /* Keep far neighbors at end of list */
            while (distance[itail] >= limit)
                itail--;

            /* Done? */
            if (ihead >= itail)
                return ihead;

            /* Swap. */
            {   double  dx = deltax[ihead],
                        dy = deltay[ihead],
                        dz = deltaz[ihead],
                        dd = distance[ihead];
                unsigned int nn = neighbor[ihead];

                deltax[ihead] = deltax[itail]; deltax[itail] = dx;
                deltay[ihead] = deltay[itail]; deltay[itail] = dy;
                deltaz[ihead] = deltaz[itail]; deltaz[itail] = dz;
                distance[ihead] = distance[itail]; distance[itail] = dd;
                neighbor[ihead] = neighbor[itail]; neighbor[itail] = nn;
            }

            /* This pair is now in correct order. */
            ihead++;
            itail--;
        }
    }
}

typedef double v2df __attribute__((vector_size (16)));

/* TODO: See if using doubles would actually be faster.
*/
static inline void swap_doubles(void *const ptr1, void *const ptr2)
{
    const unsigned long t  = *(unsigned long *)ptr1;
    *(unsigned long *)ptr1 = *(unsigned long *)ptr2;
    *(unsigned long *)ptr2 = t;
}

static inline void swap_uints(void *const ptr1, void *const ptr2)
{
    const unsigned int t  = *(unsigned int *)ptr1;
    *(unsigned int *)ptr1 = *(unsigned int *)ptr2;
    *(unsigned int *)ptr2 = t;
}

/* Update the neighbor list for one atom.
 * Return the number of interacting neighbors.
 * The neighbor and distance lists are padded to a multiple of 4,
 * with at least one unused entry at end.
*/
int distances(const unsigned int   neighbors,
               unsigned int *const neighbor,
               double *const       distance,
               double *const       deltax,
               double *const       deltay,
               double *const       deltaz,
               const double *const x,
               const double *const y,
               const double *const z,
               const unsigned int  atom,
               const double        limit)
{
    /* Note: neighbor[] array is padded with a safe index
     *       to a multiple of four entries. */

    /* Distance calculation part. */
    {
        const v2df    x0 = { x[atom], x[atom] },
                      y0 = { y[atom], y[atom] },
                      z0 = { z[atom], z[atom] };
        unsigned int  i = 0;

        while (i < neighbors) {
            const unsigned int i1 = neighbor[i++];
            const unsigned int i2 = neighbor[i++];
            const unsigned int i3 = neighbor[i++];
            const unsigned int i4 = neighbor[i++];
            v2df               x1 = { x[i1], x[i2] },
                               y1 = { y[i1], y[i2] },
                               z1 = { z[i1], z[i2] },
                               x2 = { x[i3], x[i4] },
                               y2 = { y[i3], y[i4] },
                               z2 = { z[i3], z[i4] },
                               d1, d2, t1, t2;

            /* Substract current atom coordinates. */

            x1 = __builtin_ia32_subpd(x1, x0);
            y1 = __builtin_ia32_subpd(y1, y0);
            z1 = __builtin_ia32_subpd(z1, z0);

            x2 = __builtin_ia32_subpd(x2, x0);
            y2 = __builtin_ia32_subpd(y2, y0);
            z2 = __builtin_ia32_subpd(z2, z0);

            /* Calculate squared distances. */

            d1 = x1;
            d2 = x2;
            d1 = __builtin_ia32_mulpd(d1, d1);
            d2 = __builtin_ia32_mulpd(d2, d2);

            t1 = y1;
            t2 = y2;
            t1 = __builtin_ia32_mulpd(t1, t1);
            t2 = __builtin_ia32_mulpd(t2, t2);

            d1 = __builtin_ia32_addpd(d1, t1);
            d2 = __builtin_ia32_addpd(d2, t2);

            t1 = z1;
            t2 = z2;
            t1 = __builtin_ia32_mulpd(t1, t1);
            t2 = __builtin_ia32_mulpd(t2, t2);

            d1 = __builtin_ia32_addpd(d1, t1);
            d2 = __builtin_ia32_addpd(d2, t2);

            /* Take square root to get distances. */

            d1 = __builtin_ia32_sqrtpd(d1);
            d2 = __builtin_ia32_sqrtpd(d2);

            /* Divide deltas by distances to get unit vectors. */

            x1 = __builtin_ia32_divpd(x1, d1);
            y1 = __builtin_ia32_divpd(y1, d1);
            z1 = __builtin_ia32_divpd(z1, d1);

            x2 = __builtin_ia32_divpd(x2, d2);
            y2 = __builtin_ia32_divpd(y2, d2);
            z2 = __builtin_ia32_divpd(z2, d2);

            /* Save coordinate deltas. */

            __builtin_ia32_storeupd(&deltax[i-4], x1);
            __builtin_ia32_storeupd(&deltay[i-4], y1);
            __builtin_ia32_storeupd(&deltaz[i-4], z1);

            __builtin_ia32_storeupd(&deltax[i-2], x2);
            __builtin_ia32_storeupd(&deltay[i-2], y2);
            __builtin_ia32_storeupd(&deltaz[i-2], z2);

            /* Save distances. */

            __builtin_ia32_storeupd(&distance[i-4], d1);
            __builtin_ia32_storeupd(&distance[i-2], d2);
        }
    }

    /* Set the distance sentinel. */
    distance[neighbors] = limit;

    /* Neighbor sorting part.
     * Move atoms closer than limit to start of list,
     * and atoms at limit or further to end of list.
     * Other than that, the order is not important.
    */
    {
        unsigned int ihead = 0;
        unsigned int itail = neighbors - 1;

        while (1) {

            /* Keep close neighbors at start of list */
            while (distance[ihead] < limit)
                ihead++;

            /* All close? */
            if (ihead >= itail)
                return ihead;

            /* Keep far neighbors at end of list */
            while (distance[itail] >= limit)
                itail--;

            /* Done? */
            if (ihead >= itail)
                return ihead;

            /* Swap. */
            swap_doubles(&deltax[ihead], &deltax[itail]);
            swap_doubles(&deltay[ihead], &deltay[itail]);
            swap_doubles(&deltaz[ihead], &deltaz[itail]);
            swap_doubles(&distance[ihead], &distance[itail]);
            swap_uints(&neighbor[ihead], &neighbor[itail]);

            /* This pair is now in correct order. */
            ihead++;
            itail--;
        }
    }
}