Consistent cross platform procedural generation

Question

What techniques have people successfully used or can suggest to deal with a consistent cross platform math for procedural world generation? Also, if you have done this, what were the pros and cons of your approach? In particular, I'm trying to decide whether to go with a lightweight fixed point approach, or a use a heavyweight multi-precision library [I went with the lightweight approach].

Context: I'm building a game (in C++) that involves generating a 400 billion star galaxy, for which I'm using procedural generation. The target platforms are Android (Arm-v8a), Linux (x86-64), and probably Windows (x86-64). The issue is that this generation involves quite a bit of math, including exp, pow, sqrt, sin, cos, tan and atan2, and the system is ill-conditioned (tiny differences in the calculations could propagate to large differences in the output, like a star turning up or not, or even a whole region of stars changing). I want the same seed to generate exactly the same galaxy on all platforms in case I ever want a multi-player version of this.

My reading so far is that different processors, compilers or libraries cannot be guaranteed to return the same results for transcendental functions, and even simple floating point math might be a problem (such as Arm using a multiply-add instruction having better intermediate precision than the x86-64 using two instructions).

Things I have considered:

• Using double precision floating point, including the same transcendental function code to generate the same results (which I've had to do anyway so that I can vectorize them), and tweaking the compiler flags so that everything compiles exactly the same way on both processors, and rely on IEE-754 conformance to take care of the rest. [this part answered. It's nearly impossible to get right]

• Using a lightweight fixed point framework such as fixedptc.h so I'm only doing integer math. Pros: will definitely be consistent. Cons: lots of effort to scale everything correctly, particularly as generating a galaxy requires numbers with a large dynamic range that is not normally suitable for fixed point.

• Using a heavyweight multi-precision framework such as GNU MPFR or Boost (which I think is built on something like MPFR. Pros: this is easier to use than a lightweight library. Cons: This is a lot of machinery to pull in to solve this, and I haven't found any guarantee that it is completely consistent cross platform.

Example of the sort of places where dynamic range can be a problem (although it turns out that in this case much of it can be precomputed and doesn't need to be in procedural generation):

• absolute star luminosities can range from about 430,000 times the mass of the sun to 1/10,000 the mass of the sun (that's not including planets and moons, they will have to have their own scale), that's a factor of 4,300,000,000, which can be represented in a 32.32 fixed number, but calculations before or after risk overflowing or underflowing.

• positions and distances store fine in a 64 bit integer, but when for instance working out how far from a spiral arm a point is (which affects the star density), it would be helpful to work with distance^2 to avoid expensive square roots, and squaring a number doubles the amount of dynamic range required - this may require a different scale factor for the generating a large elliptical galaxy like IC 1101 than it would to generate something like the Pleiades.

None of this is insurmountable, but it certainly requires a lot more than blindly replacing doubles by m.n fixed integers.

Notes:

• I'm not particularly concerned about performance, or even high accuracy. Most of the heavy math code is not called very often, and I can use a lightweight approach for the critical stuff.

• I'm not trying to get the rendering the same, as that appears to be a fools errand with all the different graphics cards out there, different screen resolutions, etc. I do want the underlying structure to be the same though.

• Generating it on one platform and storing it is not a solution - the galaxy is too big for that.

• This question is not asking about general cross-platform development and frameworks, about which I know there are many opinions (I've already made the major decisions for my case), just about getting consistent math for procedural generation.

• I'm not trying to solve the space scaling issues - I've already solved this (with for instance 64 bit integers for the star positions). I'm dealing with star distributions and properties which involve a bit of math.

RESULT: There is no one true answer, even in a single project. There seem to be four properties to think about when determining how to perform calculations on continuous quantities (things that would be represented mathematically with real numbers): accuracy, (cross platform) safety, speed and dynamic range.

• double has great accuracy, speed and dynamic range, but is not safe.

• float has great speed and dynamic range, but is not accurate or safe.

• fixed point has great speed and is safe, but not so good accuracy or dynamic range. More bits extends the dynamic range but slows things down. Either hand roll or use a lightweight library for this.

• multi-precision libraries like MPFR or crlibm have high accuracy and dynamic range, but are not fast, and I haven't found anything definitive on safety. (crlibm provides correctly rounded results, so it should be safe).

Looking at my project as an example:

• There are a lot of constants that can be pre-generated even without the seed, and the dynamic ranges get truly ugly here (I really don't want to do fixed math on calculations involving the speed of light squared and plank's constant). This needs accuracy and safety. For this I'm going to use doubles, but have a checksum or generate it statically on a known platform for safety.

• For the shape of the galaxy (spiral arms, sub-clusters etc.), all we need is safety and moderate dynamic range. Fixed point here (although multi-precision libraries might be ok), with a lightweight library for transcendental stuff.

• For star placement, we need safety, speed and some dynamic range (the dynamic range of luminosity can be factored out of this step). 64 bit fixed with occasional forays to 128 bits seems to be good here, although care needs to be taken with overflows and underflows. I'm using fixedptc.h for the sqrt stuff and rolling the rest myself.

• For star rendering, we need speed and dynamic range. floats are the way to go here, it's what low end graphics cards work with anyway.

• (later) for planet placement, we need safety and probably dynamic range. This will probably be the time to break out one of the multi-precision libraries.

This problem is now solved as far as I'm concerned (thanks for some really useful suggestions!), and I've done the bulk of the transition to fixed. Any further edits/comments are to leave ideas for anyone coming later who is trying to deal with the same issues.

Postscript: I found there were some places where fixed point had insufficient range (mostly in the setting up phase), so I ended up writing a small library to provide the floating point functionality I needed using integer instructions. See my answer below for details.

Answer 1

The best resource I've found on this topic is Bruce Dawson's blog article "Floating-Point Determinism"

It establishes that yes, in theory, you could get cross-platform determinism out of floating point math, but in practice you won't if your use case is sufficiently general. Some of the thornier problems he points out include:

• If you need to support 32-bit x86 platforms, and you use any library code, that library could change the per-thread rounding, precision, or denormal settings, invisible to you. The only defense is to never call library code on your generator thread (unreasonably onerous), or force the setting to what you expect between every library call (both wasteful and error prone).

  And just hope there's no gremlins showing up like misbehaving printer drivers! 😨

• different compilers on different platforms can make different decisions about the order in which to evaluate expressions, or which sub-expressions to cache and re-use. While with infinite precision real numbers these should all give identical results, with floats they might not be bit identical. As you mention, you'd often have to disable optimization to prevent this, costing performance.

• C++ and IEEE-754 standards both leave it undefined what precision to use for intermediate calculations, so you can get different behaviour depending on the platform/compiler unless you explicitly capture every intermediate to a variable (and on x87, store it to memory 😖) and forbid the compiler from applying optimizations (onerous, error prone, and bad for performance, and costs you precision you might otherwise have had available to you).

• If you target any platforms that lack a fused multiply-add, you have to forgo it on all platforms, again costing performance and precision.

So my takeaway from that article is if I need bit-for-bit exact matches on a wide range of platforms, fixed point and integer math is still the way to go.

Answer 2

What I would be trying to do is to generate on multiple levels, where things at one level could be calculated from the level above with integer math.

For example, at the top level I might draw a bunch of spiral arms on a 256x256 map, where each point would store things like density, age etc. For each point on the top level map, I would make another 256x256 map, using information from the top level and adding detail, maybe using something like a plasma fractal. Then for each point on that map, generate anywhere from zero to a few thousand stars.

Answer 3

As others pointed out, if you want true determinism, then avoid floating point variables and use integer math instead. If you need a large universe, use 64bit integers. Their resolution is good enough to address points in our real-world galaxy with a precision of about 50m.

Also keep in mind that if you want to auto-generate a whole galaxy, then it usually makes little sense to handle the position of every single object relative to the origin of the galaxy. A galaxy with realistic proportions has huge distances of empty space between star systems, which means there won't be interactions between two objects in different star systems. So there is no point in simulating more than one star system within the same 3d space. That means can use a two-tiered hierarchy for your game universe. With the first tier being the stars and the second tiers being the contents of each star system. That means you first generate the position of every star relative to the center of the galaxy (where a couple light-minutes of inaccuracy are no big deal). Then you can generate the positions of everything within each star system relative to the main star.

For deterministic sine and cosine calculation from integer angles, you can use the good old technique of a precomputed sine table. A sine table is an array with the hard-coded sines of the angles from 0° to 360°. If you need the sine or cosine of a value, you just need to look it up. If you need to calculate the angle from a sine or cosine, you can perform a binary search on that table.

Sine tables were popular in the days of old where CPUs didn't have dedicated circuits for calculating sine and cosine. Lookup tables resulted in a considerable performance improvement back then. On modern CPUs they hurt performance more than they help while sacrificing accuracy. But in your particular use-case they could be useful again because they are deterministic.

Answer 4

Rather than designing your PCG around real number arithmetic and then trying to make it deterministic, it may be easier to start from this constraint and design around it. There are a wide variety of different PCG algorithms, many of which don't require calculations involving real numbers (whether represented by floats or fixed-point ints) at all.

One particularly nice example is Wave Function Collapse; if you don't use the Shannon entropy variant of the algorithm, then only integer arithmetic is used. The same author, Maxim Gumin, also has created a programming language called MarkovJunior specifically designed for PCG, which appears to be very powerful and flexible; likewise, it has a few features which use floating point arithmetic, but if you don't use those features then it's integers all the way down.

Many integer-based PCG algorithms use discrete grids, but even if your game doesn't have a fixed grid, you can either make the grid scale small enough that things don't appear to be on a grid, or you can apply a separate post-processing stage to shift things away from their grid cells by small pseudorandom amounts.

I understand that it may not be feasible to start again from scratch if you have already designed your game's PCG system to use floats and transcendental functions, but for others finding this Q&A while still in their planning stage, I think it's worth taking the determinism constraint into account from the very beginning.

Answer 5

Partly an answer, partly a followup:

I found that while I can mostly just use fixed point, there were some circumstances where fixed just did not have the dynamic range, which left me looking for some sort of soft floating point library. For this problem the requirements were:

• at least doubles (in some places floats had insufficient range/precision)
• at least somewhat performant and not too heavyweight
• a guarantee that floating point wasn't used somewhere internally
• a permissive license

I couldn't find anything that met those criteria to my satisfaction, so I ended up creating a bare-bones library myself, which I have been using successfully. It can be found at https://github.com/royward/pseudo-double and it includes basic operators, comparison, conversion, exp, log, pow, sqrt and some trig, and some other functions. There is a C and C++ binding.

To give an idea of how it works, rather than try and follow IEEE 754 which is used by nearly all hardware, I went with another format that simplifies things a bit for software:

bits 0..15: exponent
bits 16..63: signed mantissa

Here is a sample of the add, multiply and convert to and from fixed point using this idea, other functions are similar (some error checking has been stripped out here for simplicity):

#include <stdint.h>

#define PSEUDO_DOUBLE_TOTAL_BITS 64
#define PSEUDO_DOUBLE_EXP_BITS 16
#define EXP_MASK ((1LL<<PSEUDO_DOUBLE_EXP_BITS)-1)
#define EXP_MASK_INV (~((1ULL<<PSEUDO_DOUBLE_EXP_BITS)-1))
#define PSEUDO_DOUBLE_EXP_BIAS (1U<<(PSEUDO_DOUBLE_EXP_BITS-1))
#define PSEUDO_DOUBLE_HALF_ULP ((1ULL<<(PSEUDO_DOUBLE_EXP_BITS-1))-1)
#define clz __builtin_clzll

typedef uint64_t pseudo_double_i;
typedef int64_t signed_pd_internal;
typedef __int128 signed_large_pd_internal;

inline signed_pd_internal shift_left_signed(signed_pd_internal x, int shift) {
    if(shift>=0) {
        return x<<shift;
    }
    return x>>-shift;
}

pseudo_double_i int64fixed2_to_pdi(int64_t d, int32_t e) {
    if(d==0) {
        return 0;
    }
    int negative=(d<0);
    int lead_bits=clz(negative?~d:d);
    return ((shift_left_signed(d,PSEUDO_DOUBLE_TOTAL_BITS+lead_bits-65))&EXP_MASK_INV)
        +PSEUDO_DOUBLE_EXP_BIAS+65-lead_bits+e;
}

int64_t pdi_to_int64fixed2(pseudo_double_i x, int32_t e) {
    if(x==0) {
        return 0;
    }
    signed_pd_internal vx=((signed_pd_internal)x)&EXP_MASK_INV;
    int32_t exponent=(x&EXP_MASK)-PSEUDO_DOUBLE_EXP_BIAS-e;
    int64_t ret=vx>>(PSEUDO_DOUBLE_TOTAL_BITS-exponent);
    return ret;
}

inline pseudo_double_i pdi_add(pseudo_double_i x, pseudo_double_i y) {
    int32_t expx=x&EXP_MASK;
    int32_t expy=y&EXP_MASK;
    int32_t ydiffx=expy-expx;
    if(ydiffx>=(PSEUDO_DOUBLE_TOTAL_BITS-1)) {
        return y;
    }
    if(ydiffx<=-(PSEUDO_DOUBLE_TOTAL_BITS-1)) {
        return x;
    }
    int32_t exp_max;
    signed_pd_internal vx=((signed_pd_internal)(x&EXP_MASK_INV))>>1;
    signed_pd_internal vy=((signed_pd_internal)(y&EXP_MASK_INV))>>1;
    if(ydiffx>=0) {
        exp_max=expy;
        vx>>=ydiffx;
    } else {
        exp_max=expx;
        vy>>=-ydiffx;
    }
    exp_max+=1;
    signed_pd_internal vr=(vx+vy+PSEUDO_DOUBLE_HALF_ULP)&~PSEUDO_DOUBLE_HALF_ULP;
    if(vr==0) {
        // special case - a mantissa of zero will always make the whole word zero.
        // Makes comparisons much easier
        return (pseudo_double_i)0;
    }
    int32_t leading_bits=clz(vr>0?vr:~vr)-1;
    if(leading_bits>exp_max) {
        leading_bits=exp_max;
    }
    vr<<=leading_bits;
    int32_t new_exponent=exp_max-leading_bits;
    return (pseudo_double_i)((vr&EXP_MASK_INV)+new_exponent);
}

inline pseudo_double_i pdi_mult(pseudo_double_i x, pseudo_double_i y) {
    int32_t expx=x&EXP_MASK;
    int32_t expy=y&EXP_MASK;
    signed_pd_internal vx=(signed_pd_internal)(x&EXP_MASK_INV);
    signed_pd_internal vy=(signed_pd_internal)(y&EXP_MASK_INV);
    signed_pd_internal vr=(((signed_large_pd_internal)vx)*vy)>>64;
    if(vr==0) {
        // special case - a mantissa of zero will always make the whole word zero.
        // Makes comparisons much easier
        return (pseudo_double_i)0;
    }
    int32_t leading_bits=clz(vr>0?vr:~vr)-1;
    vr<<=leading_bits;
    int32_t new_exponent=expx+expy-PSEUDO_DOUBLE_EXP_BIAS-leading_bits;
    return (pseudo_double_i)((vr&EXP_MASK_INV)+new_exponent);
}

To add some type safety (preventing accidentally using integer operations on pseudo-doubles which will return garbage), there is a wrapper:

typedef struct {
    pseudo_double_i val;
} pseudo_double;

inline pseudo_double create_pseudo_double_from_internal(pseudo_double_i x) {
    pseudo_double ret;
    ret.val=x;
    return ret;
}

inline pseudo_double int64fixed2_to_pd(int64_t d, int32_t e) {
    return create_pseudo_double_from_internal(int64fixed2_to_pdi(d,e));
}
inline int64_t pd_to_int64fixed2(pseudo_double d, int32_t e) {
    return pdi_to_int64fixed2(d.val,e);
}
inline pseudo_double pd_add(pseudo_double x, pseudo_double y) {
    return create_pseudo_double_from_internal(pdi_add(x.val,y.val));
}
inline pseudo_double pd_mult(pseudo_double x, pseudo_double y) {
    return create_pseudo_double_from_internal(pdi_mult(x.val,y.val));
}

Answer 6

Using a noise-based RNG might be an option? It seems to work well, and generates its numbers via bit shifting, so there's no precision issues to get in the way.

Not sure where the slides are, but the video description lists some benefits:

(unordered access, better reseeding, record/playback, network loss tolerance, lock-free parallelization, etc.) while being smaller, faster, and easier to use.

Answer 7

You can just use regular floating point math (single precision or double precision). It's a myth (that maybe once upon a time had some truth to it) that floating point math gives different results on different platforms or that it is difficult to get consistent results.

Regarding intermediate precision: on all modern 64-bit platforms the intermediate precision is the same (always 64 bits). If you are also targeting 32-bit platforms (i386 or arm) you can use _FPU_SETCW with the value _FPU_DOUBLE to also get identical results here. See here for more info: http://christian-seiler.de/projekte/fpmath/

SIMD/SSE should be fine as well. Also, web assembly is considered a 32-bit platform, but it uses 64 bit intermediate precision you are covered here as well with no effort.

This covers all PCs (windows/mac/linux), all mobile devices, and all modern game consoles.

Regarding compiler optimizations: A long time ago compilers were indeed buggy and you could get inconsistent results. But now this is not a problem. Just don't use --ffast-math and you will get consistent results on all compilers and all optimization levels.

Regarding sqrt: From my experience, the standard sqrt and sqrtf functions from <math.h> will give identical results on all modern platforms. I would test this on your target platforms just to be sure, but I would be very surprised if you find a platform that is different.

Regarding trigonometric functions: This is the one case where I have found differences between platforms/compilers. So you should use your own (fast) sin/cos implementation. Here is a good one: https://stackoverflow.com/a/28050328

In summary, don't just believe something you read on the internet from 15 years ago, try for yourself! floating point inconsistenty is a persistent myth that should perish.