/* the same scatter in floating point: the same 256 particles from the same
   generator, positions in cells, linear weights, mass and momentum weighted
   onto the four nodes and added into the grid, run in particle order and in
   each scrambled order (k * a + 13) % 256 for every odd a, each grid compared
   bit for bit with the first.
   cc -O2 -ffp-contract=off -DREAL=float transfer_float.c          (or double)
   -DDIV=65535 divides the positions by 65,535 instead of 65,536
   -DLIST also prints each odd multiplier and how many nodes its order changes */
#include <stdio.h>
#include <string.h>
#ifndef REAL
#define REAL float
#endif
#ifndef DIV
#define DIV 65536
#endif
#define N 256
#define W 8
#define NODES ((W + 1) * (W + 1))
typedef REAL real;
static long long next(long long s) { return (s * 1103515245LL + 12345LL) % 2147483648LL; }
static real px[N], py[N], m[N], vx[N], vy[N];
static void scatter(int a, real *gm, real *gpx, real *gpy) {
    memset(gm, 0, NODES * sizeof(real)); memset(gpx, 0, NODES * sizeof(real)); memset(gpy, 0, NODES * sizeof(real));
    for (int k = 0; k < N; k++) {
        int i = a == 0 ? k : (k * a + 13) % N;
        int cx = (int)px[i], cy = (int)py[i];
        real fx = px[i] - cx, fy = py[i] - cy;
        int n00 = cy * (W + 1) + cx, n10 = n00 + 1, n01 = n00 + W + 1, n11 = n01 + 1;
        real w00 = (1 - fx) * (1 - fy), w10 = fx * (1 - fy), w01 = (1 - fx) * fy, w11 = fx * fy;
        real qx = m[i] * vx[i], qy = m[i] * vy[i];
        gm[n00] += m[i] * w00; gm[n10] += m[i] * w10; gm[n01] += m[i] * w01; gm[n11] += m[i] * w11;
        gpx[n00] += qx * w00; gpx[n10] += qx * w10; gpx[n01] += qx * w01; gpx[n11] += qx * w11;
        gpy[n00] += qy * w00; gpy[n10] += qy * w10; gpy[n01] += qy * w01; gpy[n11] += qy * w11;
    }
}
int main(void) {
    long long s = 2207;
    for (int i = 0; i < N; i++) {
        s = next(s); px[i] = (real)(s % (W * 65536LL)) / (real)DIV;
        s = next(s); py[i] = (real)(s % (W * 65536LL)) / (real)DIV;
        s = next(s); m[i] = (real)(1 + s % 1000);
        s = next(s); vx[i] = (real)(s % 2001 - 1000);
        s = next(s); vy[i] = (real)(s % 2001 - 1000);
    }
    static real gm[NODES], gpx[NODES], gpy[NODES], hm[NODES], hpx[NODES], hpy[NODES];
    scatter(0, hm, hpx, hpy);
    double tm = 0, tx = 0, ty = 0;
    for (int n = 0; n < NODES; n++) { tm += hm[n]; tx += hpx[n]; ty += hpy[n]; }
    printf("particle order: grid mass %.17g, momentum (%.17g, %.17g)\n", tm, tx, ty);
    int worst = 0, differ97 = -1, orders = 0, same = 0;
    for (int a = 1; a < N; a += 2) {
        scatter(a, gm, gpx, gpy);
        int differ = 0;
        for (int n = 0; n < NODES; n++)
            if (memcmp(&gm[n], &hm[n], sizeof(real)) || memcmp(&gpx[n], &hpx[n], sizeof(real)) || memcmp(&gpy[n], &hpy[n], sizeof(real))) differ++;
#ifdef LIST
        printf("%d %d\n", a, differ);
#endif
        if (a == 97) differ97 = differ;
        if (differ > worst) worst = differ;
        if (differ == 0) same++;
        orders++;
    }
    printf("order (k * 97 + 13) %% 256: %d of %d nodes differ from particle order\n", differ97, NODES);
    printf("%d scrambled orders: %d give the same grid as particle order, the worst differs at %d of %d nodes\n", orders, same, worst, NODES);
    return 0;
}