By: Mosè Giordano
Re-posted from: http://giordano.github.io/blog/2019-05-26-complex-division/
Multiple dispatch is the ability to dispatch a function based on the number and
the run time types of its arguments. The Julia programming language uses
multiple dispatch as its main
paradigm, and this has been a central design concept of the language from its
inception.
I already talked about multiple dispatch when I showed a simple implementation
of the rock-paper-scissors game.
Multiple dispatch is useful not only to develop games, but also to write
mathematically correct code.
A key idea of multiple dispatch is that methods don’t belong to a single class,
but they are shared by all the arguments. Quoting from Julia
documentation:
does the addition operation in
x + ybelong toxany more than it does to
y? The implementation of a mathematical operator generally depends on the
types of all of its arguments. Even beyond mathematical operations, however,
multiple dispatch ends up being a powerful and convenient paradigm for
structuring and organizing programs.
An interesting example of how multiple dispatch allows one to more easily write
mathematically correct code is the division involving complex numbers. In
particular, the corner case of division by 0. But… is the 0 real or complex?
Indeed the answer is different depending on the nature of the 0. When you
divide a complex finite number by (positive) real 0, you’re computing the
limit:
and it is clear the direction in which you take the limit: the real axis. Thus,
for example, from the division we should expect as a result
, assuming that the limit is taken from the positive part of
the real axis. When the 0 is complex, instead, the direction of the limit in
the complex plane is not well-defined, and so the result of the division by a
complex zero should be undefined. This behaviour is also mandated by the
ISO/IEC 10967 standard (part 3,
section 5.2.5.5 “Fundamental complex floating point arithmetic”).
Comparison between languages
Let’s see how different programming languages behave with this operation.
Python
>>> import numpy as np
>>> np.complex128(1) / 0
__main__:1: RuntimeWarning: divide by zero encountered in cdouble_scalars
__main__:1: RuntimeWarning: invalid value encountered in cdouble_scalars
(inf+nanj) # This is correct
>>> np.float64(1) / np.complex128(0)
__main__:1: RuntimeWarning: divide by zero encountered in true_divide
__main__:1: RuntimeWarning: invalid value encountered in true_divide
(inf+nanj) # Should be (nan+nanj)
>>> np.complex128(1) / np.complex128(0)
(inf+nanj) # Should be (nan+nanj)
>>> np.float64(1) / np.complex128(1)
(1+0j) # Should be (1-0j)
R
> complex(real = 1.0) / 0.0
[1] Inf+NaNi # This is correct
> 1.0 / complex(real = 0.0)
[1] Inf+NaNi # Should be NaN+NaNi
> complex(real = 1.0) / complex(real = 0.0)
[1] Inf+NaNi # Should be NaN+NaNi
> 1.0 / complex(real = 1.0)
[1] 1+0i # Should be 1-0i
GNU Octave
octave:1> (1.0 + 0.0*i) / 0.0
warning: division by zero
ans = Inf
octave:2> 1.0 / (0.0 + 0.0*i)
warning: division by zero
ans = Inf
octave:3> (1.0 + 0.0*i) / (0.0 + 0.0*i)
warning: division by zero
ans = Inf
octave:4> 1.0 / (1.0 + 0.0*i)
ans = 1
C
#include <complex.h>
#include <stdio.h>
int main(void)
{
double complex c_zero = 0.0 + 0.0*I;
double complex c_one = 1.0 + 0.0*I;
double complex w1 = c_one / 0.0;
double complex w2 = 1.0 / c_zero;
double complex w3 = c_one / c_zero;
double complex w4 = 1.0 / c_one;
printf("(%g,%g) # Should be (inf,nan)\n", creal(w1), cimag(w1));
printf("(%g,%g) # Should be (nan,nan)\n", creal(w2), cimag(w2));
printf("(%g,%g) # Should be (nan,nan)\n", creal(w3), cimag(w3));
printf("(%g,%g) # Should be (1,-0)\n", creal(w4), cimag(w4));
return 0;
}
Compiled with GCC
$ gcc -std=c99 test.c&&./a.out
(inf,-nan) # Should be (inf,nan)
(inf,-nan) # Should be (nan,nan)
(inf,-nan) # Should be (nan,nan)
(1,0) # Should be (1,-0)
Compiled with ICC
$ icc -std=c99 test.c&&./a.out
test.c(7): warning #39: division by zero
double complex w1 = c_one / 0.0;
^
(inf,-nan) # Should be (inf,nan)
(inf,-nan) # Should be (nan,nan)
(inf,-nan) # Should be (nan,nan)
(1,0) # Should be (1,-0)
C++
#include <complex>
#include <iostream>
int main(void)
{
std::complex<double> c_zero(0.0, 0.0);
std::complex<double> c_one(1.0, 0.0);
std::complex<double> w1 = c_one / 0.0;
std::complex<double> w2 = 1.0 / c_zero;
std::complex<double> w3 = c_one / c_zero;
std::complex<double> w4 = 1.0 / c_one;
std::cout << w1 << " # Should be (inf,nan)" << std::endl;
std::cout << w2 << " # Should be (nan,nan)" << std::endl;
std::cout << w3 << " # Should be (nan,nan)" << std::endl;
std::cout << w4 << " # Should be (1,-0)" << std::endl;
return 0;
}
Compiled with GCC
$ g++ test.cpp&&./a.out
(inf,-nan) # Should be (inf,nan)
(inf,-nan) # Should be (nan,nan)
(inf,-nan) # Should be (nan,nan)
(1,0) # Should be (1,-0)
Compiled with ICC
$ icc test.cpp&&./a.out
(inf,-nan) # Should be (inf,nan)
(inf,-nan) # Should be (nan,nan)
(inf,-nan) # Should be (nan,nan)
(1,0) # Should be (1,-0)
Fortran90
program main
implicit none
complex :: c_zero = (0.0, 0.0)
complex :: c_one = (1.0, 0.0)
double precision :: f_zero = 0.0
double precision :: f_one = 1.0
write (*, '(F0.0,"+i",F0.0,A)') c_one / f_zero, " # Should be Inf+iNaN"
write (*, '(F0.0,"+i",F0.0,A)') f_one / c_zero, " # Correct"
write (*, '(F0.0,"+i",F0.0,A)') c_one / c_zero, " # Correct"
write (*, '(F0.0,"+i",F0.0,A)') f_one / c_one, " # Should be 1.-i0."
endprogram main
Compiled with Gfortran
$ gfortran test.f90&&./a.out
NaN+iNaN # Should be Inf+iNaN
NaN+iNaN # Should be NaN+iNaN
NaN+iNaN # Should be NaN+iNaN
1.+i0. # Should be 1.-i0.
Compiled with ifort
$ ifort test.f90&&./a.out
Inf+iNaN # Should be Inf+iNaN
Inf+iNaN # Should be NaN+iNaN
Inf+iNaN # Should be NaN+iNaN
1.+i0. # Should be 1.-i0.
Julia
julia> complex(1) / 0
Inf + NaN*im # This is correct
julia> 1 / complex(0)
NaN + NaN*im # This is correct
julia> complex(1) / complex(0)
NaN + NaN*im # This is correct
julia> 1 / complex(1)
1.0 - 0.0im # This is correct
Summary
The difference between Julia and the other languages is that the division is
implemented in a different way depending on the type of the two operands. By
using the
@which
macro we can see which are the methods called in the examples above:
julia> @which complex(1) / 0
/(z::Complex, x::Real) in Base at complex.jl:324
julia> @which 1 / complex(0)
/(a::R, z::S) where {R<:Real, S<:Complex} in Base at complex.jl:323
julia> @which complex(1) / complex(0)
/(a::Complex{T}, b::Complex{T}) where T<:Real in Base at complex.jl:327
Thus each combination of real and complex arguments calls a different method.
This directly matches the different mathematical operations that should be
carried out.
The fact that in the other languages reviewed here we get the same result for
any combination of the types of the operands suggests that in those cases the
arguments are always silently converted to or treated as they were all complex,
thus losing the information about their original nature.
Multiple dispatch is nothing completely magic. One could achieve similar
results in Python with a bunch of ifs that check the types of the arguments,
but multiple dispatch keeps the different implementations neatly separated in
different methods, instead of putting the code for the different operations in a
single catch-em-all function. Also in C++ one could define different methods
for different combinations of the type of the arguments, with the difference
that in that case the dispatch would be dynamic only on the first argument (the
class instance, the receiver of the method) and static on the others, whereas
with multiple dispatch the dispatch is based on the runtime type of all the
arguments.
Note: this post is largely inspired by the discussion in issue
#22983 of the Julia
repository.