A complex or simple bug?

exp(Complex(0, pi)) == -1 // Euler's formula: true

the above example (Euler’s identity) on the help file of Complex evaluates to false…

Not a bug at all. The calculations you’re doing are in floating point - floating point numbers can accumulate small differences over multiple operations, because not all numbers can be represented perfectly in float. Because of this, you can’t reliably use == to compare the value of a complex floating point operation to an analytical result like -1.
Here’s a decent write-up about the phenomena: https://blog.penjee.com/what-are-floating-point-errors-answered/

1 Like

Well. Then it ‘was’ a bug because when the example on the help file was written the return value was true… which was a bug but more true theoretically :wink:

ok.
so did look further and it is caused by

sin(pi)

being non-zero: 1.2246467991474e-16

while

cos(pi)

is perfectly -1

So it isn’t so much of ‘accumulation of float errors.’ it happens even before exp() or any complex arithmetics.

exp() over Complex is implemented as:

exp { ^exp(real) * Complex.new(cos(imag), sin(imag)) }

so you could see the error has only to do with sin(pi)

We are using the standard C sin() function (https://github.com/supercollider/supercollider/blob/develop/lang/LangSource/PyrMathOps.cpp#L485), leaving two possible analyses:

  • The C sin function may be buggy…? (A lot of eyes on that, over several decades… it would take some chutzpah to report a bug upstream.)

  • Or, our value for pi is wrong.

The second is guaranteed to be true (it’s part of scztt’s point, which may have been overlooked). As a transcendental number, pi requires infinite numeric precision to be exact. Floating point numbers offer only finite numeric precision. Therefore, no matter what value we define pi to be, it will always be inexact.

There’s an assumption here that sin(pi) will return 0 and cos(pi) will return exactly -1. It’s indeed a curious observation that we are passing the same floating point value to both functions, but only one of them yields the expected result. If we adjusted our value of pi to suit sin(), then cos() would show the rounding error and the complex exp result would still be off*. There seems to be no single floating point value for which both of those are true (well, I haven’t tested bit-differences… maybe there is).

  • Edit: Maybe not… The slope of cos at pi is 0, so there will be a wider range of input values very close to -1. The margin of error for sin is much smaller because of the higher slope.

You could actually go to the C standard library team to argue that there should be one double that satisfies both functions… but that code has been around for a long time. I doubt it would get much traction.

Leading back to scztt’s point: floating point numbers are not exact and it’s a conceptual error to assume that they will behave that way. In the end, you kinda just have to accept it.

hjh

Testing empirically: we can hack the least significant bits and test alternate values of pi very close to the one we are using now.

(
var x = pi,
ms32 = x.high32Bits,
ls32 = x.low32Bits;

(ls32 - 5 .. ls32 + 5).do { |ls|
	var rad = Float.from64Bits(ms32, ls),
	sin = sin(rad), cos = cos(rad);
	[ls - ls32, rad, rad == pi, sin, sin == 0, cos, cos == -1].postln;
};
""
)

[ -5, 3.1415926535898, false, 2.342910729165e-15, false, -1.0, true ]
[ -4, 3.1415926535898, false, 1.898821519315e-15, false, -1.0, true ]
[ -3, 3.1415926535898, false, 1.4547323094649e-15, false, -1.0, true ]
[ -2, 3.1415926535898, false, 1.0106430996149e-15, false, -1.0, true ]
[ -1, 3.1415926535898, false, 5.665538897648e-16, false, -1.0, true ]
[ 0, 3.1415926535898, true, 1.2246467991474e-16, false, -1.0, true ]
[ 1, 3.1415926535898, false, -3.2162452993533e-16, false, -1.0, true ]
[ 2, 3.1415926535898, false, -7.6571373978539e-16, false, -1.0, true ]
[ 3, 3.1415926535898, false, -1.2098029496355e-15, false, -1.0, true ]
[ 4, 3.1415926535898, false, -1.6538921594855e-15, false, -1.0, true ]
[ 5, 3.1415926535898, false, -2.0979813693356e-15, false, -1.0, true ]

The 0 line in the middle is the value of pi defined in SC. -1 and +1 deviate from this by the smallest possible amount, and show a larger degree of error (and the error increases, the further you go).

So SC is already using the most precise possible pi representation, in double precision.

So there’s nothing we can do on our side. We can’t invent a new double value in between two adjacent double values.

hjh

Hey james. Thanks for the thorough test!
I didn’t know you could do those bits operations on/of floats. cool!

Indeed in python

math.sin(math.pi) = 1.2246467991474e-16

I guess the help file need to be rewritten to:

exp(Complex(0, pi)) == (exp(0) * Complex(cos(pi), sin(pi)))