Difference between revisions of "Tc"

From Techniques for computer generated pictures in complex dynamics
Jump to: navigation, search
(Created page with "This is a very useful tool. A tc is a complex number z together with a complex vector dz attached to. The notation dz is probably not the best but it is to mimic physicist no...")
 
m (Operator overloading)
 
(11 intermediate revisions by the same user not shown)
Line 1: Line 1:
This is a very useful tool.
+
A tc is a data structure encoding a complex number and a variation of this number. The operations I define on these numbers are the same as most operations on complex numbers.
 +
The corresponding C++ library of functions that I designed really made writing my programs easier, especially for [[distance estimator|distance estimator methods]].
  
A tc is a complex number z together with a complex vector dz attached to.
+
== Overview ==
 +
 
 +
A tc is a complex number z together with a complex vector dz attached.
 
The notation dz is probably not the best but it is to mimic physicist notation, like
 
The notation dz is probably not the best but it is to mimic physicist notation, like
 
\[(z+dz)\times(w+dw)=w\,z+(w\,dz+z\,dw)+\text{neglected}\]
 
\[(z+dz)\times(w+dw)=w\,z+(w\,dz+z\,dw)+\text{neglected}\]
 
or
 
or
 
\[\cos(z+dz)=\cos(z)-\sin(z)dz.\]
 
\[\cos(z+dz)=\cos(z)-\sin(z)dz.\]
 +
A mathematical reinterpretation in terms of ''jets'' is given below.
  
 
== Operator overloading ==
 
== Operator overloading ==
Line 11: Line 15:
 
C++ allows [https://en.wikipedia.org/wiki/Operator_overloading operator overloading]. In other words, you can use instructions like d=c+a*b in your programs, with any kind type of objects for a,b,c,d.
 
C++ allows [https://en.wikipedia.org/wiki/Operator_overloading operator overloading]. In other words, you can use instructions like d=c+a*b in your programs, with any kind type of objects for a,b,c,d.
  
This is why it is much easier to program with complex numbers in C++ than in Java for instance. Compare C++
+
This is why it is much easier to program with complex numbers in C++ than in many other languages. Compare C++
 
<syntaxhighlight lang="cpp">
 
<syntaxhighlight lang="cpp">
 
z=u*z*(z*z+cos(c*z)+d);
 
z=u*z*(z*z+cos(c*z)+d);
 
</syntaxhighlight>
 
</syntaxhighlight>
with Java
+
with Java (no operator overloading)
 
<syntaxhighlight lang="java">
 
<syntaxhighlight lang="java">
 
z=mul(u,mul(z,add(mul(z,z),add(cos(mul(c,z)),d))));
 
z=mul(u,mul(z,add(mul(z,z),add(cos(mul(c,z)),d))));
 
</syntaxhighlight>
 
</syntaxhighlight>
or worse, with C
+
or with C
 
<syntaxhighlight lang="c">
 
<syntaxhighlight lang="c">
 
mul(w,c,z);
 
mul(w,c,z);
Line 32: Line 36:
 
== Operations ==
 
== Operations ==
  
One defines
+
One defines operations on tc objects as follows:
 
\[(z,dz) + (w,dw) = (z+w,dz+dw)\]
 
\[(z,dz) + (w,dw) = (z+w,dz+dw)\]
 
\[(z,dz) \times (w,dw) = (z\,w,w\,dz+z\,dw)\]
 
\[(z,dz) \times (w,dw) = (z\,w,w\,dz+z\,dw)\]
Line 48: Line 52:
 
In other words it is an element in the tangent space TC of the complex numbers field C.
 
In other words it is an element in the tangent space TC of the complex numbers field C.
 
This is why I chose the name tc for the C++ class.
 
This is why I chose the name tc for the C++ class.
It can also be considered as 1-jets (can be generalized to higher degree power series expansions, like $(z,b,c)$ representing $z(t)=z+b\,t+c\,t^2+o(t^2)$).
+
It can also be considered as 1-jets (can be generalized to higher degree power series expansions, like the $(z,b,c)$ being a 2-jet representing $z(t)=z+b\,t+c\,t^2+o(t^2)$).
 +
Jets are differential geometry object, i.e. there are specific formulae for computing how their expression (coordinates) changes when changing variables.
  
 
== Derivatives ==
 
== Derivatives ==
  
It computes derivatives for you.
+
The tc objects compute derivatives for you!
  
Say you defined Z=(z,1) and computed W=cos(ZZ). Then you get W=(a,b) with a=cos(z*z) and b= the derivative of a w.r.t. z: you did not need to determine that cos(z*z)=-2*z*sin(z*z), the class computed b iteratively for you.
+
(Let us be clear: it does not give you a formal formula. I meant it gives you a numeric value of f' at a given numeric value of the variable.)
 +
 
 +
Say you defined Z=(z,1), computed W=cos(Z×Z) and got W=(a,b). Then a=cos(z×z) and b=the derivative ∂a/∂z at z: you did not need to determine that $\partial \cos(z^2)/\partial z=-2z\sin(z^2)$, the class computed b iteratively for you.
 +
 
 +
Maybe it is more convincing with a more complicated example: you need the derivative at a given value of z of the quantity $\frac{\displaystyle \log z+\cos\frac{1+e^\sqrt z}{1-\sin(z)(1+z)}}{\displaystyle \sqrt{z^3-2z+12}}$? Here is a recipe: define the tc object Z=(z,1), compute the quantity above using Z, and take the dz part of the result.
  
 
== Implementation ==
 
== Implementation ==
  
Either I will put code here or a link to download code.
+
The code below is a possible implementation.
 +
Only the functions I most use are implemented (not cosh, not ==, etc... for instance).
 +
It not optimized and there is no specific handling of division by 0 and overflow, which I imagine may be problematic.
 +
 
 +
<syntaxhighlight>
 +
#include <complex>
 +
 
 +
typedef double reel;
 +
typedef std::complex<reel> cplx;
 +
 
 +
class tc {
 +
public : // all the members below are accessible by users of the package
 +
  // data members
 +
  cplx z;
 +
  cplx dz;
 +
  // member functions (declarations only, the code is given below*)
 +
  // *: for this C++ unses the following (ugly) syntax
 +
  //    return-type class-name::function-name(parameters) { code; }
 +
  // constructors
 +
  tc(cplx=cplx(0,0), cplx=cplx(0,0)); // default constructor
 +
  tc(const tc&); // copy constructor
 +
  // assigment operations : unlike math operations and math functions, they /have/ to belong to the class
 +
  tc& operator=(const tc&);
 +
  const reel& operator=(const reel&);
 +
  const cplx& operator=(const cplx&);
 +
};
 +
 
 +
// now we give the code of the member functions
 +
 
 +
// a constructor
 +
tc::tc(cplx point, cplx vecteur)
 +
{
 +
z = point;
 +
dz = vecteur;
 +
}
 +
 
 +
// the copy constructor
 +
tc::tc(const tc& t)
 +
{
 +
z = t.z;
 +
dz = t.dz;
 +
}
 +
 +
// assignment operators
 +
 
 +
tc& tc::operator=(const tc& t)
 +
{
 +
z = t.z;
 +
dz = t.dz;
 +
return *this;
 +
}
 +
 
 +
const reel& tc::operator=(const reel& r)
 +
{
 +
z = r;
 +
dz = 0;
 +
return r;
 +
}
 +
 
 +
const cplx& tc::operator=(const cplx& c)
 +
{
 +
z = c;
 +
dz = 0;
 +
return c;
 +
}
 +
 
 +
// non-member functions : usual math operators and math functions
 +
 
 +
tc operator*(const tc& a, const tc& b) {
 +
return tc(a.z*b.z, a.z*b.dz+a.dz*b.z);
 +
}
 +
 
 +
tc operator*(const reel& a, const tc& b) {
 +
return tc(a*b.z, a*b.dz);
 +
}
 +
 
 +
tc operator*(const tc& a, const reel& b) {
 +
return tc(a.z*b, a.dz*b);
 +
}
 +
 
 +
tc operator+(const tc& a, const tc& b) {
 +
return tc(a.z+b.z,a.dz+b.dz);
 +
}
 +
 
 +
tc operator+(const reel& a, const tc& b) {
 +
return tc(a+b.z,b.dz);
 +
}
 +
 
 +
tc operator+(const tc& a, const reel& b) {
 +
return tc(a.z+b,a.dz);
 +
}
 +
 
 +
tc operator-(const tc& a) {
 +
return tc(-a.z,-a.dz);
 +
}
 +
 
 +
tc operator-(const tc& a, const tc& b) {
 +
return tc(a.z-b.z,a.dz-b.dz);
 +
}
 +
 
 +
tc operator-(const reel& a, const tc& b) {
 +
return tc(a-b.z,-b.dz);
 +
}
 +
 
 +
tc operator-(const tc& a, const reel& b) {
 +
return tc(a.z-b,a.dz);
 +
}
 +
 
 +
tc operator/(const tc& a, const tc& b) {
 +
return tc(a.z/b.z,a.dz/b.z-a.z*b.dz/(b.z*b.z));
 +
}
 +
 
 +
tc operator/(const reel& a, const tc& b) {
 +
return tc(a/b.z,-a*b.dz/(b.z*b.z));
 +
}
 +
 
 +
tc operator/(const tc& a, const reel& b) {
 +
return tc(a.z/b,a.dz/b);
 +
}
 +
 
 +
tc exp(const tc& a) {
 +
cplx aux=exp(a.z);
 +
return tc(aux,aux*a.dz);
 +
}
 +
 
 +
tc sin(const tc& a) {
 +
return((exp(a*tc(0,1))-exp(a*tc(0,-1)))*tc(0,-0.5));
 +
}
 +
 
 +
tc cos(const tc& a) {
 +
return((exp(a*tc(0,1))+exp(a*tc(0,-1)))*0.5);
 +
}
 +
 
 +
tc tan(const tc& a) {
 +
  return(sin(a)/cos(a));
 +
}
 +
 
 +
tc log(const tc& a) {
 +
return tc(log(a.z),a.dz/a.z);
 +
}
 +
 
 +
tc sqrt(const tc& a) {
 +
return exp(0.5*log(a));
 +
}
 +
 
 +
tc conj(const tc& a) { // CAUTION : here when interpreting the dz part, only reel values of 'time' make sense
 +
return tc(conj(a.z), conj(a.dz));
 +
}
 +
 
 +
static const cplx I(0,1);
 +
</syntaxhighlight>

Latest revision as of 20:37, 20 December 2015

A tc is a data structure encoding a complex number and a variation of this number. The operations I define on these numbers are the same as most operations on complex numbers. The corresponding C++ library of functions that I designed really made writing my programs easier, especially for distance estimator methods.

Overview

A tc is a complex number z together with a complex vector dz attached. The notation dz is probably not the best but it is to mimic physicist notation, like \[(z+dz)\times(w+dw)=w\,z+(w\,dz+z\,dw)+\text{neglected}\] or \[\cos(z+dz)=\cos(z)-\sin(z)dz.\] A mathematical reinterpretation in terms of jets is given below.

Operator overloading

C++ allows operator overloading. In other words, you can use instructions like d=c+a*b in your programs, with any kind type of objects for a,b,c,d.

This is why it is much easier to program with complex numbers in C++ than in many other languages. Compare C++

z=u*z*(z*z+cos(c*z)+d);

with Java (no operator overloading)

z=mul(u,mul(z,add(mul(z,z),add(cos(mul(c,z)),d))));

or with C

mul(w,c,z);
cos(w,w);
add(w,w,d);
mul(v,z,z);
add(w,w,v);
mul(w,z,w);
mul(z,u,w);

Operations

One defines operations on tc objects as follows: \[(z,dz) + (w,dw) = (z+w,dz+dw)\] \[(z,dz) \times (w,dw) = (z\,w,w\,dz+z\,dw)\] \[-(z,dz) = (-z,-dz)\] \[1/(z,dz) = (1/z,-dz/zˆ2)\] \[\overline{(z,dz)} = (\overline{z},\overline{dz})\] \[\exp(z,dz)=(\exp(z),\exp(z)dz)\] \[\cos(z,dz)=(\cos(z),-\sin(z)dz)\] \[etc...\]

Tangent space and jets

One can imagine that a tc pair (z,dz) represents a moving point, starting from z and moving at speed dz: $z(t)=z+dz\,t+o(t)$ . In other words it is an element in the tangent space TC of the complex numbers field C. This is why I chose the name tc for the C++ class. It can also be considered as 1-jets (can be generalized to higher degree power series expansions, like the $(z,b,c)$ being a 2-jet representing $z(t)=z+b\,t+c\,t^2+o(t^2)$). Jets are differential geometry object, i.e. there are specific formulae for computing how their expression (coordinates) changes when changing variables.

Derivatives

The tc objects compute derivatives for you!

(Let us be clear: it does not give you a formal formula. I meant it gives you a numeric value of f' at a given numeric value of the variable.)

Say you defined Z=(z,1), computed W=cos(Z×Z) and got W=(a,b). Then a=cos(z×z) and b=the derivative ∂a/∂z at z: you did not need to determine that $\partial \cos(z^2)/\partial z=-2z\sin(z^2)$, the class computed b iteratively for you.

Maybe it is more convincing with a more complicated example: you need the derivative at a given value of z of the quantity $\frac{\displaystyle \log z+\cos\frac{1+e^\sqrt z}{1-\sin(z)(1+z)}}{\displaystyle \sqrt{z^3-2z+12}}$? Here is a recipe: define the tc object Z=(z,1), compute the quantity above using Z, and take the dz part of the result.

Implementation

The code below is a possible implementation. Only the functions I most use are implemented (not cosh, not ==, etc... for instance). It not optimized and there is no specific handling of division by 0 and overflow, which I imagine may be problematic.

#include <complex>

typedef double reel;
typedef std::complex<reel> cplx;

class tc {
 public : // all the members below are accessible by users of the package
  // data members
  cplx z;
  cplx dz;
  // member functions (declarations only, the code is given below*)
  // *: for this C++ unses the following (ugly) syntax 
  //    return-type class-name::function-name(parameters) { code; }
  // constructors
  tc(cplx=cplx(0,0), cplx=cplx(0,0)); // default constructor
  tc(const tc&); // copy constructor
  // assigment operations : unlike math operations and math functions, they /have/ to belong to the class
  tc& operator=(const tc&); 
  const reel& operator=(const reel&);
  const cplx& operator=(const cplx&);
};

// now we give the code of the member functions

// a constructor
tc::tc(cplx point, cplx vecteur)
{
 z = point;
 dz = vecteur;
}

// the copy constructor
tc::tc(const tc& t)
{
 z = t.z;
 dz = t.dz;
}
 
// assignment operators

tc& tc::operator=(const tc& t)
{
 z = t.z;
 dz = t.dz;
 return *this;
}

const reel& tc::operator=(const reel& r)
{
 z = r;
 dz = 0;
 return r;
}

const cplx& tc::operator=(const cplx& c)
{
 z = c;
 dz = 0;
 return c;
}

// non-member functions : usual math operators and math functions

tc operator*(const tc& a, const tc& b) {
 return tc(a.z*b.z, a.z*b.dz+a.dz*b.z);
}

tc operator*(const reel& a, const tc& b) {
 return tc(a*b.z, a*b.dz);
}

tc operator*(const tc& a, const reel& b) {
 return tc(a.z*b, a.dz*b);
}

tc operator+(const tc& a, const tc& b) {
 return tc(a.z+b.z,a.dz+b.dz);
}

tc operator+(const reel& a, const tc& b) {
 return tc(a+b.z,b.dz);
}

tc operator+(const tc& a, const reel& b) {
 return tc(a.z+b,a.dz);
}

tc operator-(const tc& a) {
 return tc(-a.z,-a.dz);
}

tc operator-(const tc& a, const tc& b) {
 return tc(a.z-b.z,a.dz-b.dz);
}

tc operator-(const reel& a, const tc& b) {
 return tc(a-b.z,-b.dz);
}

tc operator-(const tc& a, const reel& b) {
 return tc(a.z-b,a.dz);
}

tc operator/(const tc& a, const tc& b) {
 return tc(a.z/b.z,a.dz/b.z-a.z*b.dz/(b.z*b.z));
}

tc operator/(const reel& a, const tc& b) {
 return tc(a/b.z,-a*b.dz/(b.z*b.z));
}

tc operator/(const tc& a, const reel& b) {
 return tc(a.z/b,a.dz/b);
}

tc exp(const tc& a) {
 cplx aux=exp(a.z);
 return tc(aux,aux*a.dz);
}

tc sin(const tc& a) {
 return((exp(a*tc(0,1))-exp(a*tc(0,-1)))*tc(0,-0.5));
}

tc cos(const tc& a) {
 return((exp(a*tc(0,1))+exp(a*tc(0,-1)))*0.5);
}

tc tan(const tc& a) {
  return(sin(a)/cos(a));
}

tc log(const tc& a) {
 return tc(log(a.z),a.dz/a.z);
}

tc sqrt(const tc& a) {
 return exp(0.5*log(a));
}

tc conj(const tc& a) { // CAUTION : here when interpreting the dz part, only reel values of 'time' make sense
 return tc(conj(a.z), conj(a.dz));
}

static const cplx I(0,1);