Algorithm or JS graph drawing library that can generate a graph of 100,000+ nodes and edges while minimizing edge crossings

I'm trying to plot a directed graph of $2^{16}$ nodes and $2^{16}$ edges (but not simply a cycle). Ultimately, I need to be able to share an interactive graph (zooming, panning, labels).

Mathematica did a fine job of drawing this graph in a way that minimized edge crossings:

What you're seeing is a mass of nodes (blue) mashed up, totally hiding all the edges. This isn't a viable solution because it (1) requires an installation of Mathematica, (2) takes several minutes to generate, and (3) can't be savedexporting the plot as SVG crashed all SVG viewers I tried.

SigmaJS with random initial positions, then ForceAtlas2

It seems for large graphs, rendering with HTML5 Canvas is the way to go, and SigmaJS is a popular choice.

The first problem with SigmaJS though was that it doesn't automatically place nodes the way Mathematica did. So to apply any force-directed drawing algorithm, first I had to supply all nodes with initial positions.

Randomly dispersing the 65,536 nodes in a square caused such a hopeless tangle of edges that, even after several minutes of running ForceAtlas2, I could only see this:

SigmaJS with ring-adjacent placement, then ForceAtlas2

Well, no big deal. Instead of random placement, I decided to do a naive depth-first search and place nodes in a ring as I found them. This way the majority of nodes would start right next to a neighbor. Here's what the evolution of that looked like with ForceAtlas2 at start, a few minutes in, an hour later, and a few hours later:

But this made it really evident that the results of force-directed graphing algorithms depend heavily on their initial states. I can see each of those radial "islands" being stuck in local optima, never reaching the configuration Mathematica did.

About this particular graph, and investigating Mathematica's algorithm

The graph I'm studying is a pseudo-random number generator of the form

$$x_{next} = 5x_{current}+273 \bmod 65536$$

for the most part. (A quirk in the actual assembly code implementation actually causes shift-by-1s for ~700 of the 65536 edges.) In other words, what I'm graphing is the succession of "random" numbers generated by that formula, e.g.

$$0 \rightarrow 273 \rightarrow 1365 \rightarrow 7098 \rightarrow 35763 \rightarrow 48016 \rightarrow \ldots$$

Eventually this succession yields a number we've already seen, closing the loop and forming one of the 3 cycles of this graph. I know this isn't really about data science or "Big Data" but the technique I'm looking for is certainly developed for those applications, and the solution would help visualize similarly large, sparsely connected graphs.

To see what Mathematica's doing, first I plotted just a single succession for the first 1,000 integers, i.e.

\begin{align}0 \rightarrow 273\\1 \rightarrow 279\\2 \rightarrow 285\\ \ldots \\999 \rightarrow 5268\end{align}

(There are some 3- and 4-length chains due to the quirk mentioned before.) And here's the same for 10,000 integers:

Clearly Mathematica organizes subgraphs in some order to do with the size of each subgraph.

"Life" begins to form around 40,000 nodes, and as edges connect subgraphs at varying midpoints to produce more and more interesting shapes, and we converge toward the graph we began with:

So: Is there a Javascript library that can do this, or known algorithm(s) that I can attempt to implement? It's clearly not just a force-directed process. There's some sort of sorting happening for an initial layout.

Here's the Mathematica implementation if anyone is interested in playing around with it.

  Hex[exp_] := FromDigits[exp, 16];
  LByte[exp_] := BitAnd[exp, Hex@"00ff"];
  HByte[exp_] := BitAnd[exp, Hex@"ff00"]~BitShiftRight~8;
  PRNG[v_] := Module[{L5, H5, v1, v2, carry},
    L5 = LByte@v*5;
    H5 = HByte@v*5;
    v1 = LByte@H5 + HByte@L5 + 1;
    carry = HByte@v1~BitGet~0;
    v2 = BitShiftLeft[LByte@v1, 8] + LByte@L5;
    Mod[v2 + Hex@"0011" + carry, Hex@"ffff" + 1]
    ];
  mappings = # - PRNG@#  /@ Range[0, Hex@"ffff"];
  (* WARNING! GraphPlot takes a long time to generate. *)
  (* GraphPlot[mappings] *)