/* *********************************************************************
This source file is a part of the standard library of Scol
For the latest info, see http://www.scolring.org
Copyright (c) 2014 Stephane Bisaro aka Iri
Some functions has been originally written to the Openspace3d project.
This program is free software; you can redistribute it and/or modify it under
the terms of the GNU Lesser General Public License as published by the Free Software
Foundation; either version 2 of the License, or (at your option) any later
version.
This program is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more details.
You should have received a copy of the GNU Lesser General Public License along with
this program; if not, write to the Free Software Foundation, Inc., 59 Temple
Place - Suite 330, Boston, MA 02111-1307, USA, or go to
http://www.gnu.org/copyleft/lesser.txt
********************************************************************* */
/*
* Standard functions for complex number
* See http://redmine.scolring.org/projects/tutorials/wiki/Scol_usage
* for more informations
*/
/*! \file complex.pkg
* \author Scol team
* \version 0.1
* \copyright GNU Lesser General Public License 2.0 or later
* \brief Scol Standard Library - Complex number API
*
* \details This API provides some methods to work with the complex numbers.
*
**/
/*
* J'implémenterai les nombres complexes dans le kernel, ce sera beaucoup mieux !
* Lorsque j'aurai un poil plus de temps :)
*/
fun std_ccheck2string (str)=
if (nil == str) || (!strcmp "" str) then
0
else
let 0 -> i in
let 1 -> r in
let "0123456789." -> pattern in
(
while ((i < (strlen str)) && (r == 1)) do
let strfind ctoa nth_char str i pattern 0 -> p in
if p == nil then
set r = 0
else
set i = i+1;
r
);;
fun std_cfroms_real (str)=
let strfind "+" str 0 -> pos in
if pos == nil then
[0.0 nil]
else
let strtrim substr str 0 pos -> szR in
if !std_ccheck2string szR then
[0.0 pos]
else
[atof szR pos];;
fun std_cfroms_img (str, posreal)=
let if nil == posreal then 0 else posreal -> pos in
let strfind "i" str pos -> pos2 in
if pos2 == nil then
0.0
else
let strtrim substr str pos+1 pos -> szI in
if (nil == szI) || (!strcmp szI "") then
1.0
else if !std_ccheck2string szI then
0.0
else
atof szI;;
fun std_cfroms (str)=
let std_cfroms_real str -> [r p] in
[r std_cfroms_img str p];;
fun std_ctos (r, i)=
strcatn (ftoa r)::"+"::(ftoa i)::"i"::nil;;
fun std_carg (a, b, flag)=
if (a == nil) || (b == nil) then
nil
else if (a == 0.0) && (b == 0.0) then
if !flag then 0.0 else nil
else
atan2 a b;;
/*else if (a == 0) && (b <. 0.0) then
(-PIf)/.2.0
else if (a == 0) && (b >. 0.0) then
PIf/.2.0
else if (a <. 0.0) && (b <. 0.0) then
(atan b/.a)-.PIf
else if (a <. 0.0) && (b >=. 0.0) then
(atan b/.a)+.PIf
else // a > 0
atan b/.a;;*/
/*! \struct Complex
*
* \ingroup std_complex
*
* \brief Internal structure. You should not call it directly, use
* API instead !
*
**/
struct Complex = [
cReal : F,
cImg : F
] mkComplex;;
/*! \brief Create a new Complex
*
* \ingroup std_complex
* Prototype: fun [F F] Complex
*
* \param F : the real part
* \param F : the imaginary part
*
* \return Complex : the new complex number
**/
fun std_cNew (fReal, fImg)=
mkComplex [fReal fImg];;
/*! \brief Create a new Complex from a literal string (such as \f$a + bi\f$)
*
* \ingroup std_complex
* Prototype: fun [S] Complex
*
* \param S : a literal complex string.
*
* \return Complex : the new complex number
**/
fun std_cFromS (szC)=
let std_cfroms szC -> [r i] in
mkComplex [r i];;
/*! \brief Create a new zero (0) Complex : \f$0+0i\f$
*
* \ingroup std_complex
* Prototype: fun [] Complex
*
* \return Complex : the new complex number
**/
fun std_cZero ()=
mkComplex [0.0 0.0];;
/*! \brief Check if a complex number is 0.
*
* \ingroup std_complex
* Prototype: fun [Complex] I
*
* \return I : 1 if 0
**/
fun std_cIsZero (c)=
(c.cReal == 0.0) && (c.cImg == 0.0);;
/*! \brief Set the real part of a complex number.
*
* \ingroup std_complex
* Prototype: fun [Complex F] Complex
*
* \param Complex : a complex number.
* \param F : the new real part.
*
* \return Complex : the same complex number structure
**/
fun std_cSetReal (c, fReal)=
set c.cReal = fReal;
c;;
/*! \brief Set the imaginary part of a complex number.
*
* \ingroup std_complex
* Prototype: fun [Complex F] Complex
*
* \param Complex : a complex number.
* \param F : the new imaginary part.
*
* \return Complex : the same complex number structure
**/
fun std_cSetImg (c, fImg)=
set c.cImg = fImg;
c;;
/*! \brief Set a complex number.
*
* \ingroup std_complex
* Prototype: fun [Complex F F] Complex
*
* \param Complex : a complex number.
* \param F : the new real part.
* \param F : the new imaginary part.
*
* \return Complex : the same complex number structure
**/
fun std_cSet (c, fReal, fImg)=
set c.cReal = fReal;
set c.cImg = fImg;
c;;
/*! \brief Get the real part of a complex number.
*
* \ingroup std_complex
* Prototype: fun [Complex] F
*
* \param Complex : a complex number.
*
* \return F : the real part
**/
fun std_cGetReal (c)=
c.cReal;;
/*! \brief Get the imaginary part of a complex number.
*
* \ingroup std_complex
* Prototype: fun [Complex] F
*
* \param Complex : a complex number.
*
* \return F : the imaginary part
**/
fun std_cGetImg (c)=
c.cImg;;
/*! \brief Get a complex number.
*
* \ingroup std_complex
* Prototype: fun [Complex] [F F]
*
* \param Complex : a complex number.
*
* \return [F F] : the real part and the imaginary part
**/
fun std_cGet (c)=
[c.cReal c.cImg];;
/*! \brief Get a complex number to a literal string (like \f$a+bi\f$)
*
* \ingroup std_complex
* Prototype: fun [Complex] S
*
* \param Complex : a complex number.
*
* \return S : the literal string
**/
fun std_cToS (c)=
std_ctos c.cReal c.cImg;;
/*! \brief Create the conjugate of a complex number.
*
* The conjugate of \f$a+bi\f$ is \f$a-bi\f$.
*
* \ingroup std_complex
* Prototype: fun [Complex] Complex
*
* \param Complex : a complex number.
*
* \return Complex : its conjugate (a new Complex number)
**/
fun std_cConjugate (c)=
mkComplex [c.cReal (-.c.cImg)];;
/*! \brief Get the modulus (phasis) of a complex number.
*
* The modulus of \f$a+bi\f$ is \f$|a-bi| = \sqrt{(a²+b²)}\f$
*
* \ingroup std_complex
* Prototype: fun [Complex] F
*
* \param Complex : a complex number.
*
* \return F : its modulus
**/
fun std_cMod (c)=
sqrt ((sqr c.cReal)+.(sqr c.cImg));;
/*! \brief Get the argument of a complex number.
*
* The argument (or phasis) of \f$a+bi\f$ is \f$\arctan {b / a}\f$.
*
* \ingroup std_complex
* Prototype: fun [Complex I] F
*
* \param Complex : a complex number.
* \param I : a flag. In the case where a and b are equals at 0 (zero),
* this function returns 0 if this flag is 0 or nil if this flag has
* another value. Indeed, in mathematics, this value is undefined but
* the many language, like C, returns 0 instead of a 'NaN' (Not a Number).
*
* \return F : its argument, in radians (see above)
**/
fun std_cArg (c, flag)=
std_carg c.cReal c.cImg flag;;
/*! \brief Add two complex numbers.
*
* \ingroup std_complex
* Prototype: fun [Complex Complex] [F F]
*
* \param Complex : a complex number.
* \param Complex : a second complex number.
*
* \return [F F] : the result (real part and imaginary part)
* \see std_cAddNew
**/
fun std_cAdd (c1, c2)=
if (c1 == nil) || (c2 == nil) then nil
else
[c1.cReal+.c2.cReal c1.cImg+.c2.cImg];;
/*! \brief Add two complex numbers. The return is a new complex number.
*
* \ingroup std_complex
* Prototype: fun [Complex Complex] Complex
*
* \param Complex : a complex number.
* \param Complex : a second complex number.
*
* \return Complex : a new Complex number
* \see std_cAdd
**/
fun std_cAddNew (c1, c2)=
let std_cAdd c1 c2 -> [a b] in
mkComplex [a b];;
/*! \brief Substract two complex numbers.
*
* \ingroup std_complex
* Prototype: fun [Complex Complex] [F F]
*
* \param Complex : a complex number.
* \param Complex : a second complex number.
*
* \return [F F] : the result (real part and imaginary part)
* \see std_cSubNew
**/
fun std_cSub (c1, c2)=
if (c1 == nil) || (c2 == nil) then nil
else
[c1.cReal-.c2.cReal c1.cImg-.c2.cImg];;
/*! \brief Substract two complex numbers. The return is a new complex number.
*
* \ingroup std_complex
* Prototype: fun [Complex Complex] Complex
*
* \param Complex : a complex number.
* \param Complex : a second complex number.
*
* \return Complex : a new Complex number
* \see std_cSub
**/
fun std_cSubNew (c1, c2)=
let std_cSub c1 c2 -> [a b] in
mkComplex [a b];;
/*! \brief Multiply two complex numbers.
*
* \ingroup std_complex
* Prototype: fun [Complex Complex] [F F]
*
* \param Complex : a complex number.
* \param Complex : a second complex number.
*
* \return [F F] : the result (real part and imaginary part)
* \see std_cMulNew
**/
fun std_cMul (c1, c2)=
if (c1 == nil) || (c2 == nil) then nil
else
[(c1.cReal*.c2.cReal)-.(c1.cImg*.c2.cImg) (c1.cImg*.c2.cReal)+.(c1.cReal*.c2.cImg)];;
/*! \brief Multiply two complex numbers. The return is a new complex number.
*
* \ingroup std_complex
* Prototype: fun [Complex Complex] Complex
*
* \param Complex : a complex number.
* \param Complex : a second complex number.
*
* \return Complex : a new Complex number
* \see std_cMul
**/
fun std_cMulNew (c1, c2)=
let std_cMul c1 c2 -> [a b] in
mkComplex [a b];;
/*! \brief Divide two complex numbers.
*
* \ingroup std_complex
* Prototype: fun [Complex Complex] [F F]
*
* \param Complex : a complex number.
* \param Complex : a second complex number.
*
* \return [F F] : the result (real part and imaginary part)
* \see std_cDivNew
**/
fun std_cDiv (c1, c2)=
if (c1 == nil) || (c2 == nil) then nil
else
let ((sqr c2.cReal)+.(sqr c2.cImg)) -> div in
[((c1.cReal*.c2.cReal)+.(c1.cImg*.c2.cImg))/.div ((c1.cImg*.c2.cReal)-.(c1.cReal*.c2.cImg))/.div];;
/*! \brief Divide two complex numbers. The return is a new complex number.
*
* \ingroup std_complex
* Prototype: fun [Complex Complex] Complex
*
* \param Complex : a complex number.
* \param Complex : a second complex number.
*
* \return Complex : a new Complex number
* \see std_cDiv
**/
fun std_cDivNew (c1, c2)=
let std_cDiv c1 c2 -> [a b] in
mkComplex [a b];;
/*! \brief Returns the inverse of a complex number.
*
* \ingroup std_complex
* Prototype: fun [Complex] [F F]
*
* \param Complex : a complex number.
*
* \return [F F] : the result (real part and imaginary part), nil if error
* \see std_cInvNew
* \see std_cPow (n = -1)
**/
fun std_cInv (c)=
if (c == nil) || (std_cIsZero c) then nil
else
let (sqr c.cReal)+.(sqr c.cImg) -> d in
[c.cReal/.d (-.c.cImg)/.d];;
/*! \brief Returns the inverse of a complex number. This is a new complex number.
*
* \ingroup std_complex
* Prototype: fun [Complex] Complex
*
* \param Complex : a complex number.
*
* \return Complex : a new Complex number or nil if error
* \see std_cInv
* \see std_cPowNew (n = -1)
**/
fun std_cInvNew (c)=
let std_cInv c -> [a b] in
mkComplex [a b];;
/*! \brief The square of a complex number.
*
* \ingroup std_complex
* Prototype: fun [Complex] [F F]
*
* \param Complex : a complex number.
*
* \return [F F] : the result (real part and imaginary part)
* \see std_cSqrNew
**/
fun std_cSqr (c)=
if c == nil then nil
else
[(sqr c.cReal)-.(sqr c.cImg) 2.0*.c.cReal*.c.cImg];;
/*! \brief Square of a complex number. The return is a new complex number.
*
* \ingroup std_complex
* Prototype: fun [Complex] Complex
*
* \param Complex : a complex number.
*
* \return Complex : a new Complex number
* \see std_cSqr
**/
fun std_cSqrNew (c)=
let std_cSqr c -> [a b] in
mkComplex [a b];;
fun std_csqrt (c)=
let std_cMod c -> m in
let if c.cImg <. 0.0 then -.1.0 else if c.cImg == 0.0 then 0.0 else 1.0 -> sign in
let sqrt ((m+.c.cReal)/.2.0) -> A in
let sign*.(sqrt ((m-.c.cReal)/.2.0)) -> B in
[[A B] [(-.A) (-.B)]];;
/*! \brief Square root of a complex number.
*
* \ingroup std_complex
* Prototype: fun [Complex] [[F F] [F F]]
*
* \param Complex : a complex number.
*
* \return [[F F] [F F]] : the result (two complex numbers with their real part and imaginary part)
* \see std_cSqrtNew
**/
fun std_cSqrt (c)=
if c == nil then nil
else std_csqrt c;;
/*! \brief Square root of a complex number. The return is two new complex numbers.
*
* \ingroup std_complex
* Prototype: fun [Complex] [Complex Complex]
*
* \param Complex : a complex number.
*
* \return [Complex Complex] : two new Complex number
* \see std_cSqrt
**/
fun std_cSqrtNew (c)=
let std_csqrt c -> [[a1 b1] [a2 b2]] in
[mkComplex [a1 b1] mkComplex [a2 b2]];;
/*! \brief The power of a complex number by an integer.
*
* \ingroup std_complex
* Prototype: fun [Complex I] [F F]
*
* \param Complex : a complex number.
* \param I : an integer
*
* \return [F F] : the result (real part and imaginary part)
* \see std_cPowNew
**/
fun std_cPow (c, i)=
if (c == nil) || (i == nil) then nil
else
let itof i -> f in
let pow std_cMod c f -> mf in
let f *. std_cArg c 0 -> af in
[mf*.cos af mf*.sin af];;
/*! \brief The power of a complex number by an integer The return is a
* new complex number.
*
* \ingroup std_complex
* Prototype: fun [Complex I] Complex
*
* \param Complex : a complex number.
* \param I : an integer
*
* \return Complex : a new Complex number
* \see std_cPow
*/
fun std_cPowNew (c, i)=
let std_cPow c i -> [a b] in
mkComplex [a b];;
/*! \brief Returns the first solution of the n-th root of a complex number.
*
* \ingroup std_complex
* Prototype: fun [Complex I] [F F]
*
* \param Complex : a complex number.
* \param I : an integer
*
* \return [F F] : the simplier result (real part and imaginary part)
* \remark In fact, the n-th root of a complex number is 'multi valued'.
*
* \see std_cRootnNew
* \see std_cRootnK for a particular value
* \see std_cRootnAll for all values
**/
fun std_cRootn (c, i)=
if (c == nil) || (i == nil) then nil
else
let itof i -> f in
let rootn std_cMod c f -> mf in
let (std_cArg c 0 )/.f -> af in
[mf*.cos af mf*.sin af];;
/*! \brief Returns the first solution of the n-th root of a complex number.
* The return is a new complex number.
*
* \ingroup std_complex
* Prototype: fun [Complex I] Complex
*
* \param Complex : a complex number.
* \param I : an integer
*
* \return Complex : a new Complex number (only the first value is returned here).
* \see std_cRootn
* \see std_cRootnK
* \see std_cRootnAll
*/
fun std_cRootnNew (c, i)=
let std_cRootn c i -> [a b] in
mkComplex [a b];;
/*! \brief Returns all n-th roots of a complex number.
*
* The part real is :
*
* (the n-th root of modulus) * cosine ((the argument + 2*k*Pi) / n)
*
* \f$\sqrt[n]{(a²+b²)} * \cos {((\arctan {b / a} + 2*k*\Pi)/n)}\f$
*
* The imaginary part is :
*
* (the n-th root of modulus) * sine ((the argument + 2*k*Pi) / n)
*
* \f$\sqrt[n]{(a²+b²)} * \sin {((\arctan {b / a} + 2*k*\Pi)/n)}\f$
*
* where 'k' is an integer, with \f$0 <= k < n\f$.
*
* \ingroup std_complex
* Prototype: fun [Complex I] [[I F F] r1]
*
* \param Complex : a complex number.
* \param I : an integer (the 'n_th')
*
* \return [[I F F] r1] : a list of all values. The first item of each tuple
* is the indice 'k'. The size of the list is 'n'.
*
* \see std_cRootn for 'k' = 0
* \see std_cRootnNew for 'k' = 0
* \see std_cRootnK for a given 'k'
*/
fun std_cRootnAll (c, i)=
if (c == nil) || (i == nil) then nil
else
let itof i -> f in
let i-1 -> k in
let itof k -> fk in
let rootn std_cMod c f -> mf in
let nil -> l in
(
while (k >= 0) do
let ((std_cArg c 0 )+.(2.0*.fk*.PIf))/.f -> af in
(
set l = [k mf*.cos af mf*.sin af] :: l;
set k = k-1;
);
l
);;
/*! \brief Returns a particular soultion of the n-th root of a complex number.
*
* \ingroup std_complex
* Prototype: fun [Complex I I] [F F]
*
* \param Complex : a complex number.
* \param I : an integer, the 'n'-th root
* \param I : k : a particular solution (see std_cRootnAll for more details)
* k must be positive or nul and strictly lower than n (else, nil is returned)
*
* \return [F F] : the result (real part and imaginary part)
*
* \see std_cRootnNew for 'k' = 0
* \see std_cRootn for 'k' = 0
* \see std_cRootnAll for all 'k'
**/
fun std_cRootnK (c, i, k)=
if (c == nil) || (i == nil) || (k < 0) || (k >= i) then nil
else
let itof i -> f in
let itof k -> fk in
let rootn std_cMod c f -> mf in
let ((std_cArg c 0 )+.(2.0*.fk*.PIf))/.f -> af in
[mf*.cos af mf*.sin af];;
/*! \brief The natural logarithm (base 'e') of a complex number.
*
* \ingroup std_complex
* Prototype: fun [Complex] [F F]
*
* \param Complex : a complex number.
*
* \return [F F] : the result (real part and imaginary part)
*
* \remark Cosine and sine being periodic functions, the natural
* logarithm of a complex number is also periodic. So, only the simple
* value is returned here, i.E. when 'k' = 0.
* To obtain all values, it need to apply this formula :
*
* the real part is : log of modulus : \f$\ln {\sqrt{(a²+b²)}}\f$
*
* the imaginary part is : the argument + 2*k*Pi : \f$\arctan {b / a}+ 2*k*\Pi\f$
*
* where 'k' is an integer (in \f$\mathbb{Z}\f$)
*
* If you want a particular result, add \f$2k\Pi\f$ to the returned imaginary part.
**/
fun std_cLog (c)=
if c == nil then nil
else [log std_cMod c std_cArg c 0];;
/*! \brief Returns the first solution of the natural logarithm (base 'e')
* of a complex number. The return is a new complex number.
*
* \ingroup std_complex
* Prototype: fun [Complex] Complex
*
* \param Complex : a complex number.
*
* \return Complex : a new Complex number (only the first value is returned here).
* \see std_cLog for 'k' = 0 (see std_cLog for more details)
*/
fun std_cLogNew (c)=
let std_cLog c -> [a b] in
mkComplex [a b];;
/*! \brief Exponentiation by the Euler's formula : \f$e^{f*i}\f$ where i is
* the imaginary unit (\f$i² = -1\f$) and f is a real number (here, f is a
* floatting point number).
*
* e power fi = cos f + i sin f
*
* \f$e^{f*i} = \cos {f} + i*\sin{f}\f$
*
* \ingroup std_complex
* Prototype: fun [F] [F F]
*
* \param F : a floating point number.
*
* \return [F F] : the result (real part and imaginary part)
**/
fun std_cEuler (f)=
if f == nil then nil
else [cos f sin f];;
/*! \brief Compare two complex numbers.
*
* \ingroup std_complex
* Prototype: fun [Complex Complex] I
*
* \param Complex : a complex number.
* \param Complex : a second complex number.
*
* \return Complex : 1 if equals
*
* \remark It is not possible to define an order between two complex number.
* So, only the equality can be defined and here 1 is returned in this case.
**/
fun std_cCmp (c1, c2)=
(c1.cReal == c2.cReal) && (c1.cImg == c2.cImg);;