Calculation speedup for complex division

Can I divide complex number such that the output is the imaginary part of the division only?
For example in below function:

# Check if point is inside or outside of the polygon
# Check more about operation...
function [ip]=pointin(m,p)

it is imag(a/c)/imag(b/c).
Can that be simplified or somehow removed the real part from division algorithm,
because the real part is not needed at all.

Test program:

octave:5> m0
m0 =

   0 + 0i
   0 + 3i
   2 + 3i
   2 + 2i
   1 + 2i
   1 + 1i
   2 + 1i
   2 + 3i
   3 + 3i
   3 + 0i
   0 + 0i

octave:6> pointin(m0,2.5+1.5i)
ans =  1

imag(a/c)/imag(b/c) can be simplified to imag(a/b), right?
Obviously the remain case is when c == 0 but then output is Nan.

Regarding your second question no you can not remove the real part because it will show in the final output:

As an example take z = a + b i and w = c + d i, then
z/w = (a+ b i)/(c+d i)=((a+b i)(c-d i))/((c+d i)(c-d i))=(ac+bd + (bc-ad)i)/(cc+dd)

So as you can see both the real and imaginary parts of each of the different elements show in the final result for both the real and imaginary components.

1 Like

Right, this is what I mean, if the real part (or ab+bd) can be removed from the division algorithm completely because its not needed at all.
for example if
(a+bi)/(c+di) = (a c+b d + (b c-a d)i)/(c c+d d)
but if it can be defined as
idiv((a+bi),(c+di)) = (b c-a d)/(c c+d d)
such that calculation for ac+bd is removed from the algorigthm as unnecesary.
And similarly:
rdiv((a+bi),(c+di)) = (a c-b d)/(c c+d d)
Would this improve code speed?
This complex arithmetic might need more operators than there is…
It could further simplify more, if the c c+d d are reduced because they appear in both divisions that are divided.
As if there is
imag((a+bi)/(c+di)) / imag((e+fi)/(c+di)) (I)
(b c-a d)/(c c+d d) / (f c-e d)/(c c+d d) (III)
reduced to
(b c-a d)/(f c-e d) (III)
The question should be if there is ‘single operator’ that gives the result as in (III)
for given 3 complex numbers (a+bi), (c+di) and (e+fi).
Even the equation (I) looks simple, there is internally a lot of stuff with cc+dd, etc that is not
really needed.
The point is that if I make myself the function as in (III) that would be slow compared to system functions and also it might not be able to handle vector arguments.
Also if operator as [a,b]=divcomp(p1) could get both real and imaginary separately parts at the same time… (to avoid writing two commands a=real(p1); b=imag(p1) instead of one).
When i tested with tic-toc the fastest was the one using the complicated complex numbers.
testdet.m (2.1 KB)
The desired function would be:

# Function for 3 complex numbers
function [ic]=comp3(p1,p2,p3)

as system function…

The first question is why would you need that?

The usual quote applies here: “premature optimization is the root of all evil.” (popularized by Donald Knuth).

How much time is it taking so that you need to worry about that… Are repeating this operation millions of times and if you do are you sure that there are not other places in the code that can be better optimized?

So in conclusion, as far as I know the answer to your question is negative. There is no operator/function that does what you need.

1 Like

I could not be sure, because this is hand made.
But it is certain that even if the complex division faster it can not remove the unnecessary calculations (for cc+dd, etc).
Maybe if it need really faster it need the conversion to the c-code with suitable libraries as well.
Basically this function could run billion times or more to check some things.

Also it looks it can ‘plot()’ the complex numbers, but it can not ‘line()’ them:

octave:31> bb=[-2-2i; 2+2i]; plot(bb); line(bb)
error: line: invalid number of PROPERTY / VALUE pairs
error: called from
    __line__ at line 50 column 5
    line at line 79 column 10

I found something faster, and that has no cc+dd:

function [k,l,c]=solvb(p)
  s1=conj(p(2)-p(1)); s2=conj(p(4)-p(3));

But it would have ‘excess’ ab-bd in the realpart that is not needed.
This seems fastest though, but maybe with some error.
Below code would ‘free’ from the burden of the cc+dd:

# Check if point is inside or outside of the polygon
# Check more about operation...
function [ip]=pointin(m,p)
  s1=conj(v(1:end)); s2=conj(p1);

But not totally free of burdens.
Had I the imul(a+bi,c-di)=bc-ad or idiv(a+bi,c+di)=bc-ad available
it might solve this without doing excess calculations.

I still do not understand the issue, with the “excess” of work. :slight_smile:

In order to say why I will give two examples from different fields.

The first example is the “Brachistochrone Problem”

given two points what is the curve with the fastest descent time?

In your analogy (taken freely so please ignore the hyperbole of trying to make my point) we should strive for the shortest path, a straight line connecting both points. In that case we do not have any excess of traveled time.

In another more close topic, where performance is of critical importance there is a talk:
CppCon 2017: Carl Cook “When a Microsecond Is an Eternity: High Performance Trading Systems in C++”

There they found that performance wise it was faster to always repeat the same code, even if they knew forehand that there are dead branches instead of shortcuting the code flow. Basically imagine a cycle where we issue a continue to go to the next cycle.

In this case the code is faster even if we do extra calculations. The issue in this case was related to the CPU cache that you want to keep “hot” and so it is good to be boring and to repeat always the same instructions even if we do not need some of them most of the time.

At the same time when optimizing do not forget about the compromises that we need to do between execution and development times. :slight_smile: Development time is not only the design and coding but also the testing and possibly later maintenance.

1 Like

The condition for the ‘cs’ in regards of ‘t’ might be better as here:

function [ip]=pointin(m,p)
  s1=conj(v(1:end)); s2=conj(p1);

But it could need check also the ‘s’ as well such that the left side and right side would be
handled in same manner.

octave:5> pointin(m0,0+1.5i)
ans =  1
octave:6> pointin(m0,3+1.5i)
ans = 0