Complex Numbers

From Erlang Community

(Difference between revisions)
Revision as of 22:31, 21 November 2006 (edit)
Gordon guthrie (Talk | contribs)

← Previous diff		 Revision as of 18:01, 29 October 2010 (edit) (undo)
ChrisW (Talk | contribs)
(I needed some more complex number functions, so I took the module that was previously here and enhanced it.)
Next diff →
Line 5: Line 5:
== Solution == == Solution ==
-Unfortunately, Erlang does not come with built-in support for complex numbers. Luckily, we can create our own module to do so. +Unfortunately, Erlang does not come with built-in support for complex numbers. Luckily, we can create our own module to do so.
 + 
 +This module contains the basic functions to perform the following operations:
 +1) Arithmetic
 +2) Roots, powers and reciprocal
 +3) Rectangular <-> Polar conversion
 +4) Trigonometric
 +5) Hyperbolic
 + 
<code> <code>
-% Super-primitive Complex Number Implementation+-module(cplx).
-%+-export([make/1, make/2, is_complex/1,
--module(complex).+ real/1, imag/1,
--export([make/2, is_complex/1, add/2, sub/2, mult/2, divide/2,+ add/2, subtract/2, multiply/2, divide/2,
- get_real/1, get_imaginary/1]).+ sqrt/1, root/2, pow/2, reciprocal/1, absolute/1,
 + polar/1, rect/1,
 + cos/1, sin/1, tan/1, cot/1,
 + cosh/1, sinh/1, tanh/1, coth/1]).
 + 
 +%% What does a complex number look like?
 +-record(complex, {r=0.0, i=0.0}).
 + 
 +%% Shorten the complex number validity test by wrapping it as a macro
 +%% This also allows the function to be used a guard.
 +-define(IsComplex(Z), is_record(Z, complex)).
 + 
 +is_complex(Z) -> ?IsComplex(Z).
 + 
 +%% Go make a complex number
 +make(Real, Imag) when is_integer(Real); is_float(Real),
 + is_integer(Imag); is_float(Imag) -> #complex{r=Real, i=Imag}.
 +make(Real) when is_integer(Real); is_float(Real) -> make(Real, 0);
 +make(Real) -> Real.
 + 
 +%% Return the real or imaginary parts of a complex number
 +real(Z) when is_integer(Z); is_float(Z) -> Z;
 +real(Z) when ?IsComplex(Z) -> Z#complex.r.
 + 
 +imag(Z) when ?IsComplex(Z) -> Z#complex.i.
 + 
 +%% ----------------------------------------------------------------------------
 +%% Arithmetic functions
 +%%
 +%% Add complex to complex, or complex to real
 +add(Z1, Z2) when ?IsComplex(Z1), ?IsComplex(Z2) -> #complex{r=Z1#complex.r + Z2#complex.r, i=Z1#complex.i + Z2#complex.i};
 +add(Z, R) when ?IsComplex(Z), is_integer(R); is_float(R) -> #complex{r=Z#complex.r + R, i=Z#complex.i}.
 + 
 +%% Subtract complex Z2 from complex Z1, or real from complex
 +subtract(Z1, Z2) when ?IsComplex(Z1), ?IsComplex(Z2) -> #complex{r=Z1#complex.r - Z2#complex.r, i=Z1#complex.i - Z2#complex.i};
 +subtract(Z, R) when ?IsComplex(Z), is_integer(R); is_float(R) -> #complex{r=Z#complex.r - R, i=Z#complex.i}.
 + 
 +%% Multiply complex by complex, or complex by real
 +multiply(Z1, Z2) when ?IsComplex(Z1), ?IsComplex(Z2) ->
 + #complex{r=Z1#complex.r * Z2#complex.r - Z1#complex.i * Z2#complex.i,
 + i=Z1#complex.r * Z2#complex.i + Z1#complex.i * Z2#complex.r};
 +multiply(Z1, R) when ?IsComplex(Z1),
 + is_integer(R); is_float(R) -> #complex{r=Z1#complex.r * R, i=Z1#complex.i * R}.
 + 
 +%% Divide complex Z1 by complex Z2 or complex by real
 +divide(Z1, Z2) when abs(Z2#complex.r) >= abs(Z2#complex.i) ->
 + P = Z2#complex.i / Z2#complex.r,
 + Q = Z2#complex.r + P * Z2#complex.i,
 + 
 + #complex{r=(Z1#complex.r + P * Z1#complex.i) / Q,
 + i=(Z1#complex.i - P * Z1#complex.r) / Q};
 + 
 +divide(Z1, Z2) when ?IsComplex(Z1), ?IsComplex(Z2) ->
 + P = Z2#complex.r / Z2#complex.i,
 + Q = Z2#complex.i + P * Z2#complex.r,
 + 
 + #complex{r=(Z1#complex.r * P + Z1#complex.i) / Q,
 + i=(Z1#complex.i * P - Z1#complex.r) / Q};
 + 
 +divide(Z,R) when ?IsComplex(Z), is_integer(R); is_float(R) ->
 + #complex{r=Z#complex.r / R, i=Z#complex.i / R}.
 + 
 +%% Absolute value of complex number Z -> |Z|
 +absolute(Z) when Z#complex.r == 0 -> Z#complex.i;
 +absolute(Z) when Z#complex.i == 0 -> Z#complex.r;
 +absolute(Z) when Z#complex.r > Z#complex.i -> abs(Z#complex.r) * math:sqrt(1 + math:pow(abs(Z#complex.i)/abs(Z#complex.r),2));
 +absolute(Z) -> abs(Z#complex.i) * math:sqrt(1 + math:pow(abs(Z#complex.r)/abs(Z#complex.i),2)).
 + 
 + 
 +%% Square root of complex Z
 +%% sqrt(Z) when ?IsComplex(Z) -> root(Z,2).
 +sqrt(Z) when ?IsComplex(Z) ->
 + Modulus = math:pow((Z#complex.r * Z#complex.r) + (Z#complex.i * Z#complex.i), 0.5),
 + #complex{r=math:pow((Modulus + Z#complex.r)/2, 0.5),
 + i=sign(Z#complex.i) * math:pow((Modulus - Z#complex.r)/2, 0.5)}.
 + 
 + 
 +%% Nth root of complex Z
 +root(Z,N) when Z#complex.r >= 0, Z#complex.i == 0 -> #complex{r=math:pow(Z#complex.r, (1/N)), i=0};
 +root(Z,N) when Z#complex.r < 0, Z#complex.i == 0 -> #complex{r=0, i=math:pow(abs(Z#complex.r), (1/N))};
 +root(Z,N) when Z#complex.i /= 0 ->
 + Zpolar = polar(Z),
 + Modulo = math:pow(Zpolar#complex.r, (1/N)),
 + Alpha = (Zpolar#complex.i + 2 * math:pi()) / N,
 +
 + #complex{r=Modulo * math:cos(Alpha), i=Modulo * math:sin(Alpha)}.
--record( complex, {real, imaginary}). 
-is_complex(X) when record(X, complex) -> true;+%% Raise complex Z to the power of N
-is_complex(_) -> false.+pow(Z,N) when Z#complex.i == 0, is_integer(N); is_float(N) -> #complex{r=math:pow(Z#complex.r,N), i=0};
 +pow(Z,N) when ?IsComplex(Z), is_integer(N); is_float(N) ->
 + Zpolar = polar(Z),
 + rect(#complex{r=math:pow(Zpolar#complex.r,N), i=Zpolar#complex.r * N}).
-make(Real, Imaginary) ->+%% Reciprocal of complex Z
- #complex{real = Real, imaginary = Imaginary}.+reciprocal(Z) when ?IsComplex(Z) -> divide(#complex{r=1,i=0},Z).
-make(X) when record(X, complex) -> X; 
-make(Real) -> make(Real, 0)). 
-add(X, Y) ->+%% ----------------------------------------------------------------------------
- A = make(X), B = make(Y),+%% Coordinate conversion functions
- make( A#complex.real + B#complex.real,+%%
- A#complex.imaginary + B#complex.imaginary).+%% Rectangular to polar conversion
 +polar(Z) when Z#complex.r == 0 -> #complex{r=absolute(Z), i=(sign(Z#complex.i) * math:pi() / 2)};
 +polar(Z) when Z#complex.r > 0 -> #complex{r=absolute(Z), i=math:atan(Z#complex.i / Z#complex.r)};
 +polar(Z) when Z#complex.r < 0, Z#complex.i /= 0 -> #complex{r=absolute(Z), i=(sign(Z#complex.i) * math:pi() + math:atan(Z#complex.i / Z#complex.r))};
 +polar(Z) when Z#complex.r < 0, Z#complex.i == 0 -> #complex{r=absolute(Z), i=math:pi()}.
-sub(X, Y) ->+%% Polar to rectangular conversion
- A = make(X), B = make(Y),+rect(Z) when Z#complex.r == 0 -> #complex{r=0, i=0};
- make( A#complex.real - A#complex.real,+rect(Z) when Z#complex.i == 0 -> #complex{r=Z#complex.r, i=0};
- B#complex.imaginary - B#complex.imaginary).+rect(Z) when ?IsComplex(Z) -> #complex{r=Z#complex.r * math:cos(Z#complex.i), i=Z#complex.r * math:sin(Z#complex.i)}.
-mult(X, Y) -> 
- A = make(X), B = make(Y), 
- make( (A#complex.real * B#complex.real) 
- - (A#complex.imaginary * B#complex.imaginary), 
- (A#complex.real * B#complex.imaginary) 
- + (B#complex.real * A#complex.imaginary) ). 
-divide(X,Y) -> 
- A = make(X), B = make(Y), 
- Divisor = math:pow(B#complex.real,2) + math:pow(B#complex.imaginary,2), 
- make( ((A#complex.real * B#complex.real) 
- + (A#complex.imaginary * B#complex.imaginary)) / Divisor, 
- ((A#complex.imaginary * B#complex.real) 
- - (A#complex.real * B#complex.imaginary)) / Divisor). 
-get_real(X) -> X#complex.real.+%% ----------------------------------------------------------------------------
 +%% Trigonometrical functions
 +%%
 +cos(Z) when ?IsComplex(Z) -> #complex{r=math:cos(Z#complex.r) * math:cosh(Z#complex.i), i=math:sin(Z#complex.r) * math:sinh(Z#complex.i) * -1}.
 +sin(Z) when ?IsComplex(Z) -> #complex{r=math:sin(Z#complex.r) * math:cosh(Z#complex.i), i=math:cos(Z#complex.r) * math:sinh(Z#complex.i)}.
 +tan(Z) when ?IsComplex(Z) -> divide(sin(Z), cos(Z)).
 +cot(Z) when ?IsComplex(Z) -> divide(cos(Z), sin(Z)).
-get_imaginary(X) -> X#complex.imaginary.+
 +
 +%% ----------------------------------------------------------------------------
 +%% Hyperbolic functions
 +%%
 +cosh(Z) when ?IsComplex(Z) -> #complex{r=math:cosh(Z#complex.r) * math:cos(Z#complex.i), i=math:sinh(Z#complex.r) * math:sin(Z#complex.i)}.
 +sinh(Z) when ?IsComplex(Z) -> #complex{r=math:sinh(Z#complex.r) * math:cos(Z#complex.i), i=math:cosh(Z#complex.r) * math:sin(Z#complex.i)}.
 +tanh(Z) when ?IsComplex(Z) -> divide(sinh(Z), cosh(Z)).
 +coth(Z) when ?IsComplex(Z) -> divide(cosh(Z), sinh(Z)).
 +
 +
 +
 +
 +
 +%% ----------------------------------------------------------------------------
 +%% Utility functions
 +%%
 +%% Return +1, 0 or -1 depending on comparison of N with 0
 +sign(N) when N > 0 -> +1;
 +sign(N) when N < 0 -> -1;
 +sign(_) -> 0.
</code> </code>
Line 78: Line 188:
-This simple module only scratches the surface of what would be required to implement a full module supporting complex arithmetic. For starters, we are missing trigonometric functions, powers, square root, etc. +Hope this helps...
 +--[[User:ChrisW|ChrisW]] 19:01, 29 October 2010 (BST)
[[Category:CookBook]][[Category:NumberRecipes]] [[Category:CookBook]][[Category:NumberRecipes]]

Revision as of 18:01, 29 October 2010

Problem

You wish to manipulate complex numbers (i.e., numbers with both a real and imaginary component.) This is commonly needed in engineering, science, and mathematics.

Solution

Unfortunately, Erlang does not come with built-in support for complex numbers. Luckily, we can create our own module to do so.

This module contains the basic functions to perform the following operations: 1) Arithmetic 2) Roots, powers and reciprocal 3) Rectangular <-> Polar conversion 4) Trigonometric 5) Hyperbolic 

-module(cplx).
-export([make/1, make/2, is_complex/1,
         real/1, imag/1,
         add/2, subtract/2, multiply/2, divide/2,
         sqrt/1, root/2, pow/2, reciprocal/1, absolute/1,
         polar/1, rect/1,
         cos/1, sin/1, tan/1, cot/1,
         cosh/1, sinh/1, tanh/1, coth/1]).

%% What does a complex number look like?
-record(complex, {r=0.0, i=0.0}).

%% Shorten the complex number validity test by wrapping it as a macro
%% This also allows the function to be used a guard.
-define(IsComplex(Z), is_record(Z, complex)).

is_complex(Z) -> ?IsComplex(Z).

%% Go make a complex number
make(Real, Imag) when is_integer(Real); is_float(Real),
                      is_integer(Imag); is_float(Imag) -> #complex{r=Real, i=Imag}.
make(Real) when is_integer(Real); is_float(Real)       -> make(Real, 0);
make(Real)                                             -> Real.

%% Return the real or imaginary parts of a complex number
real(Z) when is_integer(Z); is_float(Z) -> Z;
real(Z) when ?IsComplex(Z) -> Z#complex.r.

imag(Z) when ?IsComplex(Z) -> Z#complex.i.

%% ----------------------------------------------------------------------------
%% Arithmetic functions
%%
%% Add complex to complex, or complex to real
add(Z1, Z2) when ?IsComplex(Z1), ?IsComplex(Z2)            -> #complex{r=Z1#complex.r + Z2#complex.r, i=Z1#complex.i + Z2#complex.i};
add(Z, R)   when ?IsComplex(Z), is_integer(R); is_float(R) -> #complex{r=Z#complex.r + R, i=Z#complex.i}.

%% Subtract complex Z2 from complex Z1, or real from complex
subtract(Z1, Z2) when ?IsComplex(Z1), ?IsComplex(Z2)            -> #complex{r=Z1#complex.r - Z2#complex.r, i=Z1#complex.i - Z2#complex.i};
subtract(Z, R)   when ?IsComplex(Z), is_integer(R); is_float(R) -> #complex{r=Z#complex.r - R, i=Z#complex.i}.

%% Multiply complex by complex, or complex by real
multiply(Z1, Z2) when ?IsComplex(Z1), ?IsComplex(Z2) ->
  #complex{r=Z1#complex.r * Z2#complex.r - Z1#complex.i * Z2#complex.i,
           i=Z1#complex.r * Z2#complex.i + Z1#complex.i * Z2#complex.r};
multiply(Z1, R)  when ?IsComplex(Z1),
                      is_integer(R); is_float(R)     -> #complex{r=Z1#complex.r * R, i=Z1#complex.i * R}.

%% Divide complex Z1 by complex Z2 or complex by real
divide(Z1, Z2) when abs(Z2#complex.r) >= abs(Z2#complex.i) ->
  P = Z2#complex.i / Z2#complex.r,
  Q = Z2#complex.r + P * Z2#complex.i,

  #complex{r=(Z1#complex.r + P * Z1#complex.i) / Q,
           i=(Z1#complex.i - P * Z1#complex.r) / Q};

divide(Z1, Z2) when ?IsComplex(Z1), ?IsComplex(Z2) ->
  P = Z2#complex.r / Z2#complex.i,
  Q = Z2#complex.i + P * Z2#complex.r,

  #complex{r=(Z1#complex.r * P + Z1#complex.i) / Q,
           i=(Z1#complex.i * P - Z1#complex.r) / Q};

divide(Z,R) when ?IsComplex(Z), is_integer(R); is_float(R) ->
  #complex{r=Z#complex.r / R, i=Z#complex.i / R}.

%% Absolute value of complex number Z -> |Z|
absolute(Z) when Z#complex.r == 0          -> Z#complex.i;
absolute(Z) when Z#complex.i == 0          -> Z#complex.r;
absolute(Z) when Z#complex.r > Z#complex.i -> abs(Z#complex.r) * math:sqrt(1 + math:pow(abs(Z#complex.i)/abs(Z#complex.r),2));
absolute(Z)                                -> abs(Z#complex.i) * math:sqrt(1 + math:pow(abs(Z#complex.r)/abs(Z#complex.i),2)).


%% Square root of complex Z
%% sqrt(Z) when ?IsComplex(Z) -> root(Z,2).
sqrt(Z) when ?IsComplex(Z) ->
  Modulus = math:pow((Z#complex.r * Z#complex.r) + (Z#complex.i * Z#complex.i), 0.5),
  #complex{r=math:pow((Modulus + Z#complex.r)/2, 0.5),
           i=sign(Z#complex.i) * math:pow((Modulus - Z#complex.r)/2, 0.5)}.


%% Nth root of complex Z
root(Z,N) when Z#complex.r >= 0, Z#complex.i == 0 -> #complex{r=math:pow(Z#complex.r, (1/N)), i=0};
root(Z,N) when Z#complex.r  < 0, Z#complex.i == 0 -> #complex{r=0, i=math:pow(abs(Z#complex.r), (1/N))};
root(Z,N) when Z#complex.i /= 0 ->
  Zpolar = polar(Z),
  Modulo = math:pow(Zpolar#complex.r, (1/N)),
  Alpha  = (Zpolar#complex.i + 2 * math:pi()) / N,
  
  #complex{r=Modulo * math:cos(Alpha), i=Modulo * math:sin(Alpha)}.


%% Raise complex Z to the power of N
pow(Z,N) when Z#complex.i == 0, is_integer(N); is_float(N) -> #complex{r=math:pow(Z#complex.r,N), i=0};
pow(Z,N) when ?IsComplex(Z),    is_integer(N); is_float(N) ->
  Zpolar = polar(Z),
  rect(#complex{r=math:pow(Zpolar#complex.r,N), i=Zpolar#complex.r * N}).

%% Reciprocal of complex Z
reciprocal(Z) when ?IsComplex(Z) -> divide(#complex{r=1,i=0},Z).


%% ----------------------------------------------------------------------------
%% Coordinate conversion functions
%%
%% Rectangular to polar conversion
polar(Z) when Z#complex.r == 0                  -> #complex{r=absolute(Z), i=(sign(Z#complex.i) * math:pi() / 2)};
polar(Z) when Z#complex.r > 0                   -> #complex{r=absolute(Z), i=math:atan(Z#complex.i / Z#complex.r)};
polar(Z) when Z#complex.r < 0, Z#complex.i /= 0 -> #complex{r=absolute(Z), i=(sign(Z#complex.i) * math:pi() + math:atan(Z#complex.i / Z#complex.r))};
polar(Z) when Z#complex.r < 0, Z#complex.i == 0 -> #complex{r=absolute(Z), i=math:pi()}.

%% Polar to rectangular conversion
rect(Z) when Z#complex.r == 0 -> #complex{r=0, i=0};
rect(Z) when Z#complex.i == 0 -> #complex{r=Z#complex.r, i=0};
rect(Z) when ?IsComplex(Z)    -> #complex{r=Z#complex.r * math:cos(Z#complex.i), i=Z#complex.r * math:sin(Z#complex.i)}.



%% ----------------------------------------------------------------------------
%% Trigonometrical functions
%%
cos(Z) when ?IsComplex(Z) -> #complex{r=math:cos(Z#complex.r) * math:cosh(Z#complex.i), i=math:sin(Z#complex.r) * math:sinh(Z#complex.i) * -1}.
sin(Z) when ?IsComplex(Z) -> #complex{r=math:sin(Z#complex.r) * math:cosh(Z#complex.i), i=math:cos(Z#complex.r) * math:sinh(Z#complex.i)}.
tan(Z) when ?IsComplex(Z) -> divide(sin(Z), cos(Z)).
cot(Z) when ?IsComplex(Z) -> divide(cos(Z), sin(Z)).

  
  
%% ----------------------------------------------------------------------------
%% Hyperbolic functions
%%
cosh(Z) when ?IsComplex(Z) -> #complex{r=math:cosh(Z#complex.r) * math:cos(Z#complex.i), i=math:sinh(Z#complex.r) * math:sin(Z#complex.i)}.
sinh(Z) when ?IsComplex(Z) -> #complex{r=math:sinh(Z#complex.r) * math:cos(Z#complex.i), i=math:cosh(Z#complex.r) * math:sin(Z#complex.i)}.
tanh(Z) when ?IsComplex(Z) -> divide(sinh(Z), cosh(Z)).
coth(Z) when ?IsComplex(Z) -> divide(cosh(Z), sinh(Z)).
  
  
  
  
  
%% ----------------------------------------------------------------------------
%% Utility functions
%%
%% Return +1, 0 or -1 depending on comparison of N with 0
sign(N) when N > 0 -> +1;
sign(N) when N < 0 -> -1;
sign(_)            -> 0.

Here are some examples of this module's use: 

1> c(complex).
{ok,complex}
2> A = complex:make(15.0e7,3).
{complex,1.50000e+8,3}
3> B = complex:add(1,Aa).
{complex,1.50000e+8,3}
4> io:format("~.16f~n", [complex:get_real(B)])
150000001.0000000000000000
ok
5> complex:is_complex(Aa).
true
6> C = complex:make(14.0e-4,157.2).
{complex,1.40000e-3,157.200}
7> complex:mult(A,C).
{complex,2.09528e+5,2.35800e+10}
8> io:format("~.16f ~.16fi\n",
8>   [complex:get_real(D),complex:get_imaginary(D)]).
209528.3999999999900000 23580000000.0042000000000000i
ok


Hope this helps... --ChrisW 19:01, 29 October 2010 (BST)

Retrieved from "http://www.trapexit.org/Complex_Numbers"
Erlang/OTP Projects
Personal tools