23
\$\begingroup\$

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.

\$\endgroup\$
9
  • 1
    \$\begingroup\$ How about using integer math but multiple scales? Like maybe one scale for Intergalactic space (1 unit = 100 gazillion km) another for Galactic space (1 unit = 100 bazillion mm) another for star system space, etc... ? \$\endgroup\$ Sep 10, 2022 at 23:11
  • 1
    \$\begingroup\$ @AcmeNerdGames It's not a space scaling issue - I'm already using 64 bit integers for the star positions and that is working great, it's generating the distribution and properties of the stars. Will edit question to clarify. \$\endgroup\$
    – MadMan
    Sep 10, 2022 at 23:35
  • 4
    \$\begingroup\$ On a slight tangent, do calculate some kind of checksum (e.g. a cryptographic hash) of your final galaxy data and check it against a precalculated reference checksum. That way you can be confident that the galaxy data is the same for all players without having to store all of it. \$\endgroup\$ Sep 11, 2022 at 21:55
  • 3
    \$\begingroup\$ If you're doing hierarchical generation (which you presumably are), consider calculating a checksum of the data at every level (above some reasonable threshold, at least). If something does go wrong and the checksums don't match, this also helps you pinpoint exactly where the mismatch happened. You may even want to look into things like hash trees, although those solve a slightly different problem. The concept is similar, though. \$\endgroup\$ Sep 11, 2022 at 22:13
  • 7
    \$\begingroup\$ By the way, you should be aware of the faults of games like No Man's Sky or Elite:Dangerous. To paraphrase a famous quote from Jurassic Park: "Your developers were so preoccupied with whether or not they could [procedurally generate a whole realistically scaled galaxy], they didn't stop to think if they should". A more compact game universe with conceivable scope and conscious design often provides much better gameplay value than one with is as vast as an ocean but as deep as a puddle. \$\endgroup\$
    – Philipp
    Sep 12, 2022 at 10:19

7 Answers 7

26
\$\begingroup\$

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.

\$\endgroup\$
7
  • \$\begingroup\$ Thank you! I'd already got pieces of why trying to match different floating point implementations was a really bad idea (and I've been using integers wherever I can), but that article puts it all together and nails it shut. It still leave me with the question of how to deal consistently with the transcendental stuff. \$\endgroup\$
    – MadMan
    Sep 11, 2022 at 0:36
  • 1
    \$\begingroup\$ Can you not implement the transcendental functions in software with fixed point, using a combination of table lookup and iterative refinement? I might be misunderstanding the problem you're pointing to. \$\endgroup\$
    – DMGregory
    Sep 11, 2022 at 0:48
  • \$\begingroup\$ Yes, I can (or even better, find someone else that has done it). The issue I'm running into (I will edit the question again) is that I have dynamic range issues which will mean carefully changing my scale factors all over the place. I know I can do that, or use a multi-precision library if I can find a cross platform consistency guarantee, but in either case it's a lot of work. I'm assuming that I'm not the first person to run into this and am interested in reading what approach others have used. \$\endgroup\$
    – MadMan
    Sep 11, 2022 at 0:59
  • 1
    \$\begingroup\$ I think that's different enough that it's worth posting as its own standalone question. Let this one stand as "Can I guarantee cross-platform consistency in procedural generation code using floating point?" (answer: no), and then ask a separate question "How to implement transcendental functions with a wide dynamic range using fixed point" including some concrete use case examples. (Though the latter might be as good to ask on StackOverflow, as I imagine the answer is the same in general programming as it is for gamedev, but you'll get more eyeballs there) \$\endgroup\$
    – DMGregory
    Sep 11, 2022 at 1:02
  • 1
    \$\begingroup\$ If you'd like high-quality answers to the dynamic range issue, then at minimum I'd recommend including some sample code in your question showing use cases. It seems to me (perhaps naively) that the range should be dictated by what kind of parameter you're generating - star positions might be measured in light years or parsecs, while planet orbits might be measured in AU - and so the appropriate range would be known in advance, letting you explicitly specify which fixed point types and methods to use for each generated parameter, no dynamic range selection needed. Can you show where this fails? \$\endgroup\$
    – DMGregory
    Sep 11, 2022 at 1:14
7
\$\begingroup\$

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.

\$\endgroup\$
5
  • \$\begingroup\$ I'm doing something a little like this, and I made that part integer from the start (there are a whole lot of details in getting this right, particularly with the wide variety of distances that stars can be visible from - that would make a really good paper if I had time to write it and somewhere to put it). I'm going to pass a level around with all the calculations so I can shift out unnecessary bits and get more accurate calculations. \$\endgroup\$
    – MadMan
    Sep 11, 2022 at 19:45
  • 1
    \$\begingroup\$ Using pseudo-random number generation makes this easier. Code each level so that it requires only the RNG seed handed down from the higher level, and you should get consistent results. \$\endgroup\$
    – EvilSnack
    Sep 11, 2022 at 21:18
  • 2
    \$\begingroup\$ @EvilSnack The RNG part is easy. It's all the code to make the stars fit the density function of the galaxy that has a lot of geometry - that's the part that I've just about finished moving to fixed point (and might have caused stars to shift/appear/disappear on different platforms if I'd used floating point). \$\endgroup\$
    – MadMan
    Sep 12, 2022 at 5:36
  • \$\begingroup\$ @MadMan You don't have to generate all stars at one level. You could generate the brighest stars on a higher level, making it easier to find them over a large area. \$\endgroup\$ Sep 12, 2022 at 10:01
  • 2
    \$\begingroup\$ @JamesHollis already on it. I currently divide the stars into 13 luminosity bands for this purpose. That was the bit I got working (with integers) first. \$\endgroup\$
    – MadMan
    Sep 12, 2022 at 10:04
2
\$\begingroup\$

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.

\$\endgroup\$
6
  • \$\begingroup\$ I'd check your math on that 0.5 mm - I make it 50 m. Then you don't want to use the whole range of that or you can overflow coordinates just by subtracting them, so 100 m. After allowing for wanting to view the galaxy from the outside, satellite galaxies etc, I went for 10km which is still good enough. When doing calculations, I frequently want distance squared, which means 128 bit integers for intermediate results. +1 for lookup tables - I won't go that route as there are fixed point routines that should be deterministic, and speed does matter a little, but they do mostly fit my use case. \$\endgroup\$
    – MadMan
    Sep 12, 2022 at 9:58
  • \$\begingroup\$ @MadMan 10km precisions for stars is still plenty. I elaborated with another paragraph why you don't need positions relative to the galaxy center for any objects which are smaller than stars. \$\endgroup\$
    – Philipp
    Sep 12, 2022 at 10:02
  • \$\begingroup\$ Regarding the calculation: Yes, you are correct. I think I somehow mixed up meters and kilometers. \$\endgroup\$
    – Philipp
    Sep 12, 2022 at 10:06
  • \$\begingroup\$ I'm ahead of you there :-) I also use 128 bit integers to store the viewpoint and viewpoint velocity so there will never be any weird effects when changing scale. \$\endgroup\$
    – MadMan
    Sep 12, 2022 at 10:07
  • \$\begingroup\$ Any tips for handling the luminosity case described in the question? \$\endgroup\$
    – DMGregory
    Sep 12, 2022 at 12:18
2
\$\begingroup\$

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));
}
\$\endgroup\$
1
\$\begingroup\$

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.

\$\endgroup\$
1
\$\begingroup\$

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.

\$\endgroup\$
1
  • \$\begingroup\$ This doesn't answer my question (all the RNGs that I've seen are integer based, so it's something I already have), but is something worth following up anyway. \$\endgroup\$
    – MadMan
    Sep 14, 2022 at 6:09
0
\$\begingroup\$

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.

\$\endgroup\$
6
  • 1
    \$\begingroup\$ What about the other points raised in Bruce Dawson's blog, including 3rd party code changing flags, or optimizing compilers changing order of evaluation? \$\endgroup\$
    – DMGregory
    Sep 11, 2022 at 12:25
  • 10
    \$\begingroup\$ Speaking from experience, it’s not a myth. It’s less prevalent now than it used to be, but for any long complex calculation, minute differences from intermediate precision and rounding semantics can and sometimes do stack up. You are correct that on the most common modern platforms (barring 32-bit x86, but that’s dying fast) you’re generally fine, but going beyond that you need proven stable algorithms (that is, ones which are unaffected by rounding mode) and guarantees about intermediate precision to be guaranteed consistent across platforms. \$\endgroup\$ Sep 11, 2022 at 12:49
  • 1
    \$\begingroup\$ On intermediate results: intel.com/content/dam/develop/external/us/en/documents/pdf/… : "... The compiler may generate a single FMA instruction for a combined multiply and add operation, e.g. a = b*c + d. This leads to faster, slightly more accurate calculations, but results may differ in the last bit from separate multiply and add instructions, if all terms are positive, or by much more, if there is a cancellation." so at a minimum that would need to be turned off everywhere (ARM does that by default). \$\endgroup\$
    – MadMan
    Sep 11, 2022 at 19:33
  • 7
    \$\begingroup\$ For the rest of it, I hope you are right. Things have come a long way since the x87 days, but barring some authoratitive source for both gcc and clang that "these are the flags that guarantee the same behaviour in all cases on all processors", there is always the risk of some wierd corner case looking to bite, and my reading has found much use of the word "almost" or similar. It is unfortunate that it is impossible to reliably test a negative for what might be a very rare issue. \$\endgroup\$
    – MadMan
    Sep 11, 2022 at 19:34
  • 2
    \$\begingroup\$ I can say with complete confidence that this is false. I just got done chasing down a bug triggered by 5^2 being equal to 24.999999999999992894 rather than 25 under certain circumstances. \$\endgroup\$
    – Mark
    Sep 12, 2022 at 23:27

You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .