*[*][2]*[*][4]*[*][6]*[*][8]*[*][10]*[*][12]*[*][14]*
Vladimir V. Kisil
School of Mathematics
University of Leeds
Leeds LS2 9JT
UK
http://www.maths.leeds.ac.uk/~kisilv/
This is an implementation of the Schwerdtfeger--Fillmore--Springer--Cnops construction (SFSCc) based on the Clifford algebra capacities [->] of the GiNaC computer algebra system. SFSCc linearises the linear-fraction action of the Möbius group. This turns to be very useful in several theoretical and applied fields including engineering.
The core of this realisation of SFSCc is done for an arbitrary dimension, while a subclass for two dimensional cycles add some 2D-specific routines including a visualisation to PostScript files through the MetaPost or Asymptote software.
This library is a backbone of many results published in [->], which serve as illustrations of its usage. It can be ported (with various level of required changes) to other CAS with Clifford algebras capabilities.
[**[*]**][32] The usage of computer algebra system (CAS) in Clifford Algebra research has an established history with the famous ``Green book'' [->] already accompanied by a floppy disk with a REDUCE package. This tradition is very much alive, see for example the recent book [->] accompanied by a software CD. Numerous new packages are developed by various research teams across the world to work with Clifford algebras generally or address specific tasks, see on-line proceedings of the recent IKM-2006 conference [->].
Along this lines the present paper presents an implementation of the Schwerdtfeger--Fillmore--Springer--Cnops construction (SFSCc) along with illustrations of its usage. In the case of circles this technique was already spectacularly developed by H. Schwerdtfeger in 1960-ies, see Schwerdtfeger79a. Unfortunately, that beautiful book was not known to the present author until he accomplished his own works Kisil05b, Kisil12a. SFSCc was rediscovered later several times [Cnops02a] and [Porteous95].
SFSCc linearises the linear-fraction action of the Möbius group in Rn. This has clear advantages in several theoretical and applied fields including engineering. Our implementation is based on the Clifford algebra capacities of the GiNaC computer algebra system[->], which were described in [->].
The core of this realisation of SFSCc is done for an arbitrary dimension of _Rn_with a metric given by an arbitrary bilinear form. We also present a subclass for two dimensional cycles (i.e. circles, parabolas and hyperbolas), which add some 2D specific routines including a visualisation to PostScript files through the MetaPost [->] or Asymptote [->] packages. This software is the backbone of many results published in [->] and we use its application to [->] for the demonstration purpose.
There is a Python wrapper [->] for this library. It is based on BoostPython and pyGiNaC packages. The [wrapper allows to use all functions and methods from the library in][48] Python scripts or Python interactive shell. The drawing of object from [cycle2D][49] may be instantly seen in the interactive mode through the Asymptote. The live DVD supplied with book [->] is based on the library presented in this paper and its Python wrapper.
The present package can be ported (with various level of required changes) to other CAS with Clifford algebras capabilities similar to GiNaC.
[**[*]][54] The [cycle][55] class describes loci of points _xinRn_defined by a quadratic equation [*] kx^2-2lx+m=0, where k,minR, lin**Rn. [The class ][59][cycle][60] correspondingly has member variables k, l, m to describe the equation eq:quadric-def and the Clifford algebra unit to describe the metric of surrounding space. The plenty of methods are supplied for various tasks within SFSCc.
We also define a subclass [cycle2D][63] which has more methods specific to two dimensional environment.
[cycle][66][**[[klzzwxh:0127]][70]**][71]
Here is various constructors for the [cycle][72]s. The first one takes values of k, l, _m_as well as metric supplied directly. Note that _l_is admitted either in form of a lst, [matrix][75] or indexed objects from GiNaC. Similarly metric can be given by an object from either tensor, indexed, [matrix][78] or clifford classes exactly in the same way as metric is provided for a clifford_unit() constructors [->].
<cycle class constructors>= [D]->
public:
cycle(const ex & k1, const ex & l1, const ex & m1,
const ex & metr = -(new tensdelta)->setflag(status_flags::dynallocated));
Defines
cycle(links are to index).
[Constructor for a ][90][cycle][91] eq:quadric-def with k=1_and given l defined by the condition that square of its ``radius'' (which is __det_C, see [->]Defn. [->]]Kisil05a) is r_squared. Note that for the default value of the metric the value of l coincides with the centre of this [cycle][96].
<cycle class constructors>+= <-D->
cycle(const lst & l,
const ex & metr = -(new tensdelta)->setflag(status_flags::dynallocated),
const ex & r_squared = 0, const ex & e = 0,
const ex & sign = (new tensdelta)->setflag(status_flags::dynallocated));
Defines
cycle(links are to index).
[If we want to have a cycle identical to to a given one ][108]C up to a space metric which should be replaced by a new one metr, we can use the next constructor.
<cycle class constructors>+= <-D->
cycle(const cycle & C, const ex & metr);
[To any cycle SFSCc associates a matrix, which is of the][115] form eq:matrix-from-cycle [->] E-eq:matrix-for-cycle]Kisil05a. The following constructor make a [cycle][117] from its matrix representation, i.e. it is the realisation of the inverse of the map Q [->]E-eq:matrix-for-cycle]Kisil05a.
<cycle class constructors>+= <-D
cycle(const matrix & M, const ex & metr, const ex & e = 0, const ex & sign = 0);
[cycle][129]The following set of methods get_*() provide a reading access to the various data in the class.
<accessing the data of a cycle>= [D]->
public:
virtual inline ex get_dim() const { return ex_to<varidx>(unit.op(1)).get_dim(); }
virtual inline ex get_metric() const { return ex_to<clifford>(unit).get_metric(); }
virtual inline ex get_metric(const ex &i0, const ex &i1) const
{ return ex_to<clifford>(unit).get_metric(i0, i1); }
virtual inline ex get_k() const { return k; }
[The member ][142]l can be obtained as the whole by the call get_l(), or its individual component is read, for example, by get_l(1).
<accessing the data of a cycle>+= <-D->
inline ex get_l() const { return l; }
inline ex get_l(const ex & i) const
{ return (l.is_zero()?0:l.subs(l.op(1) == i, subs_options::no_pattern)); }
inline ex get_m() const {return m;}
inline ex get_unit() const {return unit;}
[Methods ][151]nops(), op(), let_op(), is_equal(), subs() are standard for expression in GiNaC and described in the GiNaC tutorial. The first three methods are rarely called by a user. In many cases the method subs() may replaced by more suitable subject_to() [->].
<accessing the data of a cycle>+= <-D->
size_t nops() const {return 4;}
ex op(size_t i) const;
ex & let_op(size_t i);
bool is_equal(const basic & other, bool projectively = true) const;
bool is_zero() const;
cycle subs(const ex & e, unsigned options = 0) const;
inline cycle normal() const
{ return cycle(k.normal(), l.normal(), m.normal(), unit.normal());}
[We also provide a method ][162]the_same_as() which return a GiNaC::lst of identities (i.e. GiNaC::relationals), which defines that two cycles are given by the same point of the projective space P3.
<accessing the data of a cycle>+= <-D
ex the_same_as(const basic & other) const;
[**[*]**][168] Cycles are represented by a points in a projective vector space, thus we wish to have a full set of linear operation on them. The metric is inherited from the first [cycle][169] object. First we define it as an methods of the [cycle][172] class.
<Linear operation as cycle methods>=
virtual cycle add(const cycle & rh) const;
virtual cycle sub(const cycle & rh) const;
virtual cycle exmul(const ex & rh) const;
virtual cycle div(const ex & rh) const;
[After that we overload standard binary operations for ][184][cycle][185].
<Linear operation on cycles>= [D]->
const cycle operator+(const cycle & lh, const cycle & rh);
const cycle operator-(const cycle & lh, const cycle & rh);
const cycle operator(const cycle & lh, const ex & rh);
const cycle operator(const ex & lh, const cycle & rh);
const cycle operator/(const cycle & lh, const ex & rh);
Defines
cycle(links are to index).
[We also define a product of two cycles through their matrix][206] representation eq:matrix-from-cycle.
<Linear operation on cycles>+= <-D
const ex operator*(const cycle & lh, const cycle & rh);
Defines
ex(links are to index).
[cycle][214]We start from some general methods which deal with [cycle][219].
[The next][222] method is needed to get rid of the homogeneous ambiguity in the projective space of cycles. If k_new=0 the [cycle][223] is normalised such that its __det__becomes 1. Otherwise the first non-zero coefficient among k, m, l_0, l_1, ... is set to k_new.
<specific methods of the class cycle>= [D]->
public:
cycle normalize(const ex & k_new = numeric(1), const ex & e = 0) const;
cycle normalize_det(const ex & e = 0,
const ex & sign = (new tensdelta)->setflag(status_flags::dynallocated)) const;
[The method ][234]center() returns a list of components of the cycle centre or the corresponding vector (D matrix) if the dimension is not symbolic. The metric, if not supplied is taken from the cycle.
<specific methods of the class cycle>+= <-D->
virtual ex center(const ex & metr = 0, bool return_matrix = false) const;
[The next method returns the value of the expression][240] -kx^2-2lx+m_for the given cycle and point _x. Obviously it should be __if _x_belongs to the cycle.
<specific methods of the class cycle>+= <-D->
virtual ex val(const ex & y) const;
[Then method ][246]passing() returns a relational defined by the identity kx^2-2lx+m===0, i.e this relational describes incidence of point to a cycle.
<specific methods of the class cycle>+= <-D->
inline ex passing(const ex & y) const {return val(y).numer().normal() == 0;}
[We oftenly need to consider a cycle which satisfies some additional][252] conditions, this can be done by the following method subject_to. Its typical application looks like: C1 = C.subject_to(lst(C.passing(P), C.is_orthogonal(C1))); The second parameters vars specifies which components of the [cycle][253] are considered as unknown. Its default value represents all of them which are symbols.
<specific methods of the class cycle>+= <-D->
virtual cycle subject_to(const ex & condition, const ex & vars = 0) const;
[*] There is a set of specific methods which represent mathematical side of SFSCc.
The next method is the main gateway to the SFSCc, it generates the 2×2_matrix [**[*]][265] [l_i ][266]sigma^i_j[1]e^j --- m
k --- -l_i sigma**^i_j [1]e^j from the cycle kx^2-2lx+m=0. Note, that the Clifford unit [1]e_has an arbitrary metric unrelated to the initial metric stored in the unit member variable.
<specific methods of the class cycle>+= <-D->
virtual matrix to_matrix(const ex & e = 0,
const ex & sign = (new tensdelta)->setflag(status_flags::dynallocated)) const;
[The next method returns the value of determinant of the][273] matrix eq:matrix-from-cycle corresponding to the [cycle][274]. It has explicit geometric meaning, see [->]§ [->]]Kisil05a. Before calculation the cycle is normalised by the condition k==k_norm, if k_norm is zero then no normalisation is done.
<specific methods of the class cycle>+= <-D->
virtual ex det(const ex & e = 0,
const ex & sign = (new tensdelta)->setflag(status_flags::dynallocated),
const ex & k_norm = 0) const;
[The determinant of a k-normalised cycle can be treated as the square of its radius][286]
<specific methods of the class cycle>+= <-D->
virtual inline ex radius_sq(const ex & e = 0,
const ex & sign = (new tensdelta)->setflag(status_flags::dynallocated)) const
{ return this->det(e, sign, numeric(1));}
[The matrix eq:matrix-from-cycle corresponding to a cycle may][293] be multiplied by another matrix, which in turn may be either generated by another cycle or be of a different origin. The next methods multiplies a cycle by another cycle or matrix supplied in C.
<specific methods of the class cycle>+= <-D->
virtual ex mul(const ex & C, const ex & e = 0,
const ex & sign = (new tensdelta)->setflag(status_flags::dynallocated),
const ex & sign1 =0) const;
[Having a matrix ][302]C_which represents a cycle and another matrix _M_we can consider a similar matrix _M^-1CM. The later matrix will correspond to a cycle as well, which may be obtained by the following three methods. In the case then _M_belongs to the __group the next two methods make a proper conversion of _M_into Clifford-valued form.
<specific methods of the class cycle>+= <-D->
virtual cycle sl2_similarity(const ex & a, const ex & b, const ex & c, const ex & d,
const ex & e = 0,
const ex & sign = (new tensdelta)->setflag(status_flags::dynallocated),
bool not_inverse=true,
const ex & sign_inv = (new tensdelta)->setflag(status_flags::dynallocated)) const;
virtual cycle sl2_similarity(const ex & M, const ex & e = 0,
const ex & sign = (new tensdelta)->setflag(status_flags::dynallocated),
bool not_inverse=true,
const ex & sign_inv = (new tensdelta)->setflag(status_flags::dynallocated)) const;
[If ][319]M_is a generic _2×2-matrix of another sort then it is used in the similarity in the unchanged form by the next method.
<specific methods of the class cycle>+= <-D->
virtual cycle matrix_similarity(const ex & M, const ex & e = 0,
const ex & sign = (new tensdelta)->setflag(status_flags::dynallocated),
bool not_inverse=true,
const ex & sign_inv = (new tensdelta)->setflag(status_flags::dynallocated)) const;
[The ][328]2×2-matrix _M= a --- b
c --- d _can be also defined by the collection of its elements.
<specific methods of the class cycle>+= <-D->
virtual cycle matrix_similarity(const ex & a, const ex & b, const ex & c, const ex & d,
const ex & e = 0,
const ex & sign = (new tensdelta)->setflag(status_flags::dynallocated),
bool not_inverse=true,
const ex & sign_inv = (new tensdelta)->setflag(status_flags::dynallocated)) const;
[Finally, we have a method for reflection of a cycle in another cycle][340] C, which is given by the similarity of the representing matrices: C C_1 C, see [->]§ [->]]Kisil05a.
<specific methods of the class cycle>+= <-D->
virtual cycle cycle_similarity(const cycle & C, const ex & e = 0,
const ex & sign = (new tensdelta)->setflag(status_flags::dynallocated),
const ex & sign1 = 0,
const ex & sign_inv = (new tensdelta)->setflag(status_flags::dynallocated)) const;
A cycle in the matrix form eq:matrix-from-cycle naturally defines a Möbius transformations of the points: [**[[klzzwxh:0834]][352]][353] l_i sigma^i_j[1]e^j --- m
k --- -l_i sigma^i_j [1]e^j : x fracl_i sigma^i_j[1]e^j x + m k x -l_i sigma**^i_j [1]e^j The following methods realised this transformations.
<specific methods of the class cycle>+= <-D->
virtual inline ex moebius_map(const ex & P, const ex & e = 0,
const ex & sign = (new tensdelta)->setflag(status_flags::dynallocated)) const
{return clifford_moebius_map(to_matrix(e, sign), P, (e.is_zero()?get_metric():e));}
For two matrices _C_1_and _C_2_obtained from cycles the expression [**[*]**][362] C_1C_2=-(C_1C_2) naturally defines an inner product in the space of cycles. The follwong methods realised it.
<specific methods of the class cycle>+= <-D->
virtual inline ex cycle_product(const cycle & C, const ex & e = 0,
const ex & sign = (new tensdelta)->setflag(status_flags::dynallocated)) const
{return scalar_part(ex_to<matrix>(mul(C, e, sign)).trace());}
[The inner product eq:cycle-inner-product defines an][371] orthogonality relation _C_1C_2===0_in the space of cycles which returned by the method is_orthogonal().
<specific methods of the class cycle>+= <-D->
virtual inline ex is_orthogonal(const cycle & C, const ex & e = 0,
const ex & sign = (new tensdelta)->setflag(status_flags::dynallocated)) const
{return (cycle_product(C, e, sign) == 0);}
[In many cases we need a higher order orthogonal relation between][379] cycles--- so called f-orthogonality, see [->]§ [->]]Kisil05a, which is given by the relation: (s []sss)=0.
<specific methods of the class cycle>+= <-D->
ex is_f_orthogonal(const cycle & C, const ex & e = 0,
const ex & sign = (new tensdelta)->setflag(status_flags::dynallocated),
const ex & sign1 = (new tensdelta)->setflag(status_flags::dynallocated),
const ex & sign_inv = (new tensdelta)->setflag(status_flags::dynallocated)) const;
[The remaining to methods check if a cycle is a liner object and if][391] it is normalised to k=1.
<specific methods of the class cycle>+= <-D
inline ex is_linear() const {return (k == 0);}
inline ex is_normalized() const {return (k == 1);}
[**[*]**][398] Two dimensional cycle [cycle2D][399] is a derived class of [cycle][402]. We need to add only very few specific methods for two dimensions, notably for the visualisation.
[This a specialisation of the constructors from ][405][cycle][406] class to [cycle2D][409]. Here is the main constructor.
<constructors of the class cycle2D>= [D]->
public:
cycle2D(const ex & k1, const ex & l1, const ex & m1,
const ex & metr = -unit_matrix(2));
[Constructor for the ][419][cycle2D][420] from l and square of its radius.
<constructors of the class cycle2D>+= <-D->
cycle2D(const lst & l, const ex & metr = -unit_matrix(2), const ex & r_squared =0,
const ex & e =0, const ex & sign = unit_matrix(2));
Defines
cycle2D(links are to index).
[Construction of ][432][cycle2D][433] from its SFSCc matrix.
<constructors of the class cycle2D>+= <-D->
cycle2D(const matrix & M, const ex & metr, const ex & e = 0, const ex & sign = 0);
[Make a two dimensional cycle out of a general one, if the][444] dimensionality of the space permits. The metric of point space can be replaced as well if a valid metr is supplied.
<constructors of the class cycle2D>+= <-D
cycle2D(const cycle & C, const ex & metr = 0);
[The realisation of 2D cycles through matrices with hypercomplex][450] numbers [->] [->]lead to some important differences with this library using the Clifford algebras. One of them: the determinant of a matrix change since. The next method return the determinant as it will be calculated on those hypercomplex matrices.
<methods specific for class cycle2D>= [D]->
public:
virtual inline ex hdet(const ex & e = 0,
const ex & sign = (new tensdelta)->setflag(status_flags::dynallocated),
const ex & k_norm = 0) const
{return -det(e, sign, k_norm);}
[The method ][459]focus() returns list of the focus coordinates and the focal length is provided by focal_length(). This turns to be meaningful not only for parabolas, see [->].
<methods specific for class cycle2D>+= <-D->
ex focus(const ex & e = diag_matrix(lst(-1, 1)), bool return_matrix = false) const;
inline ex focal_length() const {return (get_l(1)/2/k);} // focal length of the cycle
[The methods ][468]roots() returns values of u(if first = true) such that k(u^2-sigmay^2)-2l_1u-2l_2y+m=0, i.e. solves a quadratic equations. If first = false then values of _v_satisfying to _k(y^2-sigmav^2)-2l_1y-2l_2v+m=0_are returned.
<methods specific for class cycle2D>+= <-D->
lst roots(const ex & y = 0, bool first = true) const;
[The next methods is a generalisation of the previous one: it returns intersection][473] points with the line ax+b.
<methods specific for class cycle2D>+= <-D->
lst line_intersect(const ex & a, const ex & b) const;
The method metapost_draw() outputs to the stream ost MetaPost comands to draw parts of two the [cycle2D][479] within the rectangle with the lower left vertex (xmin, ymin) and upper right (xmax, [ymax][482]). The colour of drawing is specified by color (the default is black) and any additional MetaPost options can be provided in the string more_options. By default each set of the drawing commands is preceded a comment line giving description of the cycle, this can be suppressed by setting with_header = false. The default number of points per arc is reasonable in most cases, however user can override this with supplying a value to points_per_arc. The last parameter is for internal use.
<methods specific for class cycle2D>+= <-D->
void metapost_draw(ostream & ost, const ex & xmin = -5, const ex & xmax = 5,
const ex & ymin = -5, const ex & ymax = 5, const lst & color = lst(),
const string more_options = "",
bool with_header = true, int points_per_arc = 0, bool asymptote = false,
const string picture = "", bool only_path=false, bool is_continuation=false) const;
[The similar method provides a drawing output for][492] Asymptote [->] with the same meaning of parameters. However format of more_options should be adjusted correspondingly. Currently asy_draw() is realised as a wrapper around metapost_draw() but this may be changed.
<methods specific for class cycle2D>+= <-D->
inline void asy_draw(ostream & ost, const string picture,
const ex & xmin = -5, const ex & xmax = 5,
const ex & ymin = -5, const ex & ymax = 5, const lst & color = lst(),
const string more_options = "",
bool with_header = true, int points_per_arc = 0) const
{ metapost_draw(ost, xmin, xmax, ymin, ymax, color, more_options, with_header,
points_per_arc, true, picture); }
inline void asy_draw(ostream & ost = std::cout,
const ex & xmin = -5, const ex & xmax = 5,
const ex & ymin = -5, const ex & ymax = 5, const lst & color = lst(),
const string more_options = "",
bool with_header = true, int points_per_arc = 0) const
{ metapost_draw(ost, xmin, xmax, ymin, ymax, color, more_options, with_header,
points_per_arc, true); }
[Finally, we have a similar method which does not issue drawing][505] command, instead it writes a definition for a (array of) path, which may be manipulated later.
<methods specific for class cycle2D>+= <-D
inline void asy_path(ostream & ost = std::cout,
const ex & xmin = -5, const ex & xmax = 5,
const ex & ymin = -5, const ex & ymax = 5,
int points_per_arc = 0, bool is_continuation = false) const
{ metapost_draw(ost, xmin, xmax, ymin, ymax, lst(), "", false,
points_per_arc, true, "", true, is_continuation); }
A quick illustration of the library usage is the symbolic calculation which proves the Lem. [->] from [->]: We check that a Möbius transformation gin_acts on cycles by similarity _g: C -->gCg^-1. We use the following predefined objects: [cycle2D C(k,lst(l,n),m,e);][516] ex W=lst(u,v);
Firstly we define a [cycle2D][517] C2 by the condition between k, l and m in the generic [cycle2D][520] C that C passes through some point W.
<Moebius transformation of cycles>= [D]->
C2 = C.subject_to(lst(C.passing(W)));
[The point ][525]gW is defined to be the Möbius transform of W by an arbitrary g.
<Moebius transforms of W>=
const matrix gW=ex_to<matrix>(clifford_moebius_map(sl2_clifford(a, b, c, d, e), W, e).subs(sl2_relation1,
subs_options::algebraic | subs_options::no_pattern).normal());
Defines
matrix(links are to index).
[Finally we verify that the new cycle ][530]_gCg^-1_passes through P. This proves Lem. [->] from [->].
<Moebius transformation of cycles>+= <-D
cout << "Conjugation of a cycle comes through Moebius transformation: "
<< C2.sl2_similarity(a, b, c, d, es, S2, true, S2).val(gW).subs(sl2_relation1,
subs_options::algebraic | subs_options::no_pattern).normal().is_zero()
<< endl << endl;
[18]:
[48]:
[71]:
[90]: