stars in D

stars2a.d

module main;
import std.stdio, std.algorithm, std.array, std.conv, std.parallelism, core.memory;

immutable mask = 0b11111110000000000000;
alias int[] star;

struct MinDist { //information about minimum distance between stars in some area
    long d2;
    star a, b;
}

MinDist minMD(MinDist m1, MinDist m2) { return m1.d2 < m2.d2 ? m1 : m2; }

void main(string[] argv)
{
    star[][int] cells;	// split space into 128x128x128 cubic blocks (cells) of size 8K x 8K x 8K
    GC.disable();
    foreach(s; File("stars.dat", "rt").byLine()) {	   
        int[] coords = s.split().map!(to!int).array(); //parse coords of a star
        // calculate which cell it belongs to
        int sector = ((coords[0] & mask) << 1) + ((coords[1] & mask) >> 6) + ((coords[2] & mask) >> 13);
        if (sector in cells) //add to stars of its block
            cells[sector] ~= coords;
        else
            cells[sector] = [coords];
    }
    auto nbr = [-1, 0, 1];
    int[] deltas = []; //differences of indices of neighbor blocks
    foreach(x; nbr)	foreach(y; nbr)	foreach(z; nbr)	deltas ~= (x << 14) + (y << 7) + z;
    GC.enable();

    MinDist localmin(int cell) { //for given cell find two closest stars in its area
        long mind = 10^^13;
        star besta=null, bestb=null;
        star[1000] area;
        int n = 0;
        foreach(delta; deltas) { //gather stars of current cell and its neigbors
            int sector = cell + delta;
            if (sector in cells)
                foreach(st; cells[sector])
                    area[n++] = st;
        }
        foreach(ia; 0..n) //find distances between all stars in this area
            foreach(ib; 0..ia) {
                star a = area[ia], b = area[ib]; 
                long d = (a[0]-b[0])^^2 + (a[1]-b[1])^^2 + (a[2]-b[2])^^2;
                if (d < mind) { //find minimal
                    mind = d;
                    besta = a; bestb = b;
                }
            }
        return MinDist(mind, besta, bestb); 
    }
    //for each habitated cell find a pair of closest stars in its local area
    MinDist[] mds = new MinDist[cells.length];
    foreach(i, cell; taskPool.parallel(cells.keys, 10000)) 
        mds[i] = localmin(cell);		
    //find global minimum among local minima
    MinDist md = taskPool.reduce!minMD(mds, 10000);
    writeln(std.math.sqrt(to!real(md.d2)), md.a, md.b);
}