Bug in Complex

[Tobias Neumann]

There seems to be a bug in Complex with domains Float or DoubleFloat:

Working:
(68) -> (-1 :: Complex(MachineFloat))^(3/2)

   (68)  - %i

                                      Type: Expression(Complex(MachineFloat))

But for Float it's not working as expected:

(65) -> (-1 :: Complex(Float))^(3/2)

   (65)  %i

                                                         Type: Complex(Float)

(this is fine:)

(64) -> (-1 :: Complex(Float))^(1/2)

   (64)  %i

                                                         Type: Complex(Float)

Neither working for DoubleFloat:

(69) -> (-1 :: Complex(DoubleFloat))^(3/2)

   (69)  %i

                                                   Type: Complex(DoubleFloat)

I looked at the Complex implementation and there is some special case

"if R is Float or R is DoubleFloat then" but I don't see how this would lead to wrong results.

A temporary workaround for me is to convert the exponent to Float first.

Thanks and best wishes,

Tobias

[Waldek Hebisch]

Thanks for the report. Should be fixed now in the trunk.

--
Waldek Hebisch

[Tobias Neumann]

Thank you. As it turns out, complex exponentiation (^ : % -> %) with a Complex(DoubleFloat) exponent is also broken:

Good:

(7) -> (-1.0 :: Complex(Float))^(0.5)

   (7)  0.2022266248_7959507324 E -20 + %i

                                                         Type: Complex(Float)

Broken:

(6) -> (-1.0 :: Complex(DoubleFloat))^(0.5)

   (6)  6.123233995736766e-17 - %i

                                                   Type: Complex(DoubleFloat)

[Ralf Hemmecke]

Being a bit nasty... ;-)

Why do you say it is broken?
Modulo an error bound of 10^(-15) it looks OK.

Do you expect exact values when you compute with floating point numbers?

Ralf

(7) -> a := (-1.0 :: Complex(Float))^(0.5)

(7) 0.2022266248_7959507324 E -20 + %i
Type: Complex(Float)

(8) -> b := (-1.0 :: Complex(DoubleFloat))^(0.5)

(8) 6.123233995736766e-17 - %i
Type: Complex(DoubleFloat)
(9) -> a^2

(9) - 1.0 + 0.4044532497_5919014648 E -20 %i
Type: Complex(Float)
(10) -> b^2

(10) - 1.0 - 1.2246467991473532e-16 %i
Type: Complex(DoubleFloat)

[Tobias Neumann]

Hi Ralf,

presumably you didn't notice the different sign for the imaginary part?

It's just the same issue that was present for the fractional integer exponent.

I think that this should be handled consistently so that the cut is along the negative real axis,

i.e. sqrt(-x + %i*0) -> %i*sqrt(x)

Best wishes,

Tobias

[Ralf Hemmecke]

> presumably you didn't notice the different sign for the imaginary part?

Oh, no no. I did see that. Why would that be a problem?

Yes, I agree, from the name Float and DoubleFloat one might think that
these two domains are somehow related. But when you consider them as
separate domains both results are valid (within the bounds). Why would
you prefer +%i and not -%i?

Only when you consider them together, then you might get a problem.

Just let's forget about floating point numbers and consider Gaussian
integers. If you ask me for the square root of -1 and I give you -%i,
would you claim that I am wrong? Or would you just tell me: "Hey guy you
should give me two solutions and not just one."

> It's just the same issue that was present for the fractional integer
> exponent.
> I think that this should be handled consistently so that the cut is along
> the negative real axis,
> i.e. sqrt(-x + %i*0) -> %i*sqrt(x)

When you click on ^ at
http://fricas.github.io/api/Complex.html
then you see

^: (%, %) -> %
x^y returns x to the power y.

and

^: (%, Fraction Integer) -> %
x ^ y is the rational exponentiation of x by the power y.

Not very precise, I admit, but it doesn't say anything about branch
cuts. So why would you assume that your rule holds?

You can also introduce two systems of Gaussian integer. R[i] and R[j]
with i^2 = j^2 = -1 and i=-j. Of course, these rings are isomorphic.
Maybe Complex(Float) works like R[i] and Complex(DoubleFloat) like R[j].

Ralf

[Waldek Hebisch]

That is consequence of IEEE rules about signed zero and precedence
rules:

(21) -> m1 := -1.0 :: Complex(DoubleFloat)

(21) - 1.0
Type: Complex(DoubleFloat)
(22) -> imag(m1)

(22) -0.0
Type: DoubleFloat
(23) -> m1a := complex(-1, 0)$Complex(DoubleFloat)

(23) - 1.0
Type: Complex(DoubleFloat)
(24) -> imag(m1a)

(24) 0.0
Type: DoubleFloat
(25) -> m1a^(0.5)

(25) 6.123233995736766e-17 + %i
Type: Complex(DoubleFloat)
31) -> m1b := (-1.0) ::Complex(DoubleFloat)

(31) - 1.0
Type: Complex(DoubleFloat)
(32) -> imag(m1b)

(32) 0.0
Type: DoubleFloat

Coercion has very high precendece, so your example is equvalent to

-(1::Complex(DoubleFloat))

which according to IEEE has negative imaginary part...

--
Waldek Hebisch

[Tobias Neumann]

> That is consequence of IEEE rules about signed zero and precedence rules:

Now that is a subtle and interesting phenomenon.

From that point of view everything is consistent indeed. But in practice

I think that one easily runs into this issue. I would argue that you want to be

able to swap Float and DoubleFloat in applications, and this prevents it without

writing a wrapper for the exponentiation that always enforces a positive sign for a zero imaginary part.

Also virtually every CAS out there gives you +%i for (-1.0)^(0.5) for floating point numbers.

[Tobias Neumann]

Maybe this is most convincing:

(3) -> ((-1) :: Complex(DoubleFloat))^(0.5)

   (3)  6.123233995736766e-17 - %i

                                                   Type: Complex(DoubleFloat)

(4) -> ((-1) :: Complex(DoubleFloat))^(1/2)

   (4)  %i

                                                   Type: Complex(DoubleFloat)

I would say that this is undesired behavior and just confuses people.

[Kurt Pagani]

Not at all ;) It shows once more that DF ought to have no place in a CAS. DF is
a source of confusion on its own. I'd use Interval(Float) or Float, if need be ...
BTW 0.5 and 1/2 have different types. The interpreter is sometimes too broadminded.

> --
> You received this message because you are subscribed to the Google Groups
> "FriCAS - computer algebra system" group.
> To unsubscribe from this group and stop receiving emails from it, send an email
> to fricas-devel...@googlegroups.com
> <mailto:fricas-devel...@googlegroups.com>.
> To view this discussion on the web visit
> https://groups.google.com/d/msgid/fricas-devel/da5f99f4-3951-4052-9f73-515bb7ef6c66n%40googlegroups.com
> <https://groups.google.com/d/msgid/fricas-devel/da5f99f4-3951-4052-9f73-515bb7ef6c66n%40googlegroups.com?utm_medium=email&utm_source=footer>.

[Tobias Neumann]

That was the point I was trying to make: whether you use ^: ( %, Fraction(Integer)) -> %

or ^: (%, %) -> % should at least be consistent within numerical precision for Complex(DoubleFloat).

And I expect consistency betwen Float and DoubleFloat, so that one can swap them and get the same

results within floating point errors.

Now that I know the issue I can work around it. But I think it would be better to address ("fix") it directly in FriCAS.

[Waldek Hebisch]

Well, DoubleFloat is hardware float and modern hardware follows IEEE
rules. And there is vocal group insiting that library functions
should follow IEEE rules, even if done in software. Practicaly,
FriCAS uses Lisp complex 'log' and it gives us IEEE behaviour.
I am not sure how big problem this is for you: if imaginary
part comes from computations, then you have bigger problem.
If imaginary part is constant, then 'complex(x, 0)' gives
you positve sign for imaginary part, so there is easy
workaround. If you are really determined you can create
a wrapper domain inheriting most from Complex(DoubleFloat),
but overriding a few funcions that you do not like. However,
I do not think that turning Complex(DoubleFloat) into such
wrapper domain is good idea: I can not satisfy everyone.
And adding a wrapper is more flexible than fighting with
unwanted wrapper.

Concerning case when exponent is from Fraction(Integer):
this is generic code. We may add a special case for
DoubleFloat so that it behaves like DoubleFloat exponent.
OTOH different domains show different behaviour: purely
symbolic domains are subtly different in their handling
of branches compared to numeric domains. So users who
care about branches must be careful.

--
Waldek Hebisch