Bayro-Corrochano E., Scheuermann G. Geometric Algebra Computing: in Engineering and Computer Science

Подождите немного. Документ загружается.

488 D. Hildenbrand et al.

Fig. 4 Step A. Translating

by l

and l

in direction of

n to get P

and P

are all in one plane. Every plane can be deﬁned by three points, which here are P

, and the auxiliary point P

H 2

on the x-axis. Another plane, deﬁned by the points

, P

, is orthogonal to plane π

. So the projection of the line through P

and P

onto the plane π

yields L

Proj

, which intersects the sphere S

(with center P

and radius l

), resulting in a point pair. Function pp_get2nd selects link 4 (P

)

from it.

=Sphere(P

=(P

H 2

∧P

∧e

∞

)

∗

|π

P 5P 7

=(P

∧P

∧e

∞

)

∗

Proj

(π

·L

P 5P 7

)

=pp_get2nd



∧L

Proj

)

∗



(5)

See Fig. 5.

Sphere(x, r) generates a sphere around x with the radius r:

Sphere(x, r) =x −

∞

. (6)

The functions pp_get1st and pp_get2nd each pick one point out of a point pair:

pp_get1st(x) =

√

|x ·x|−x

∞

·x

, (7)

pp_get2nd(x) =

−

√

|x ·x|+x

∞

·x

. (8)

Gaalop—High Performance Parallel Computing Based on Conformal Geometric Algebra 489

Fig. 5 Step B. Projection of the line through P

and P

onto the green plane deﬁned by P

and an auxiliary point on the x-axis. Intersection of the sphere around P

with radius l

and the

projected line returns P

6.1.3 Computation of the Position of Link 2

Links 1 and 2 are located on the yz-plane π

, so the intersection of planes π

and

results in a line, with P

and P

on it. The distance between P

and P

is l

;

hence the intersection of the sphere S

around P

with radius l

results in a point

pair, from which P

can be selected

=Sphere(P

=π

∧π

=pp_get1st



∧S

)

∗



(9)

See Fig. 6.

6.1.4 Computation of the Position of Link 3

The intersection of the two spheres S

and S

results in a circle Z

. P

must be

located on circle Z

and on plane π

as well; thus the intersection of Z

and π

490 D. Hildenbrand et al.

Fig. 6 Step C. Intersecting

the sphere around P

with

radius l

with the intersection

of the plane π

and the

yz-plane returns P

Fig. 7 Step D. The

intersection of the spheres

around P

with radius l

and

around P

with radius l

results in the red circle.

Intersecting the circle with

the plane π

returns P

results in a point pair again, from which P

can be selected

=Sphere(P

∧S

=pp_get1st



∧π

)

∗



(10)

See Fig. 7.

Gaalop—High Performance Parallel Computing Based on Conformal Geometric Algebra 491

6.1.5 Compute the Angles of the Links

Since a line is deﬁned by two points, three points are necessary to generate two

intersecting lines and to compute the angle in between. To compute the angle of the

ﬁrst link, a point above the origin (0,0,1) is used as a parameter.

angle(x,y,z)=π −arccos



(x ∧y ∧e

∞

) ·(z ∧y ∧e

∞

)

|x ∧y ∧e

∞

||z ∧y ∧e

∞



. (11)

6.2 Symbolic Optimization of the Kinematic Chain

In this section, we present the CLUCalc code for the just described inverse kinemat-

ics algorithm as an example for the input language of Gaalop. Parts of the generated

C code are presented as an example for a target implementation of Gaalop.

CluCalc models the ﬁrst part of the inverse kinematics algorithms as follows:

PH3 = VecN3(hx-1,hy,hz);

:P1 = VecN3(hx,hy,hz);

:P7 = VecN3(-px,py,pz);

norm1 = FNorm[1]*e1 - FNorm[2]*e2 - FNorm[3]*e3;

?norm2 = 1/abs(norm1);

norm3 = norm1 * norm2;

//generate translator and compute P6

tvec = (norm3*(len[6]))/2;

T1 = 1 - tvec * einf;

?P6 = T1 * P7 * ~T1;

tvec = (norm3*(len[5]+len[6]))/2;

T2 = 1 - tvec * einf;

?P5 = T2 * P7 * ~T2;

S5 = Sphere(P5,len[4]);

Pi3 = *(PH3 ^ e0 ^ P5 ^ einf);

?Pi3a = 1/abs(Pi3);

?Pi3b = Pi3 * Pi3a; // Pi3/abs(Pi3)

Gaalop optimizes the above CLUCalc code (see Sect. 4.1 for details about the opti-

mization approach) and generates the following C code:

float norm2_opt[32] = {0.0};

norm2_opt[1]=1/sqrt(FNorm[1]*FNorm[1]+FNorm[2]*FNorm[2]

+FNorm[3]*FNorm[3]);

float P6_opt[32] = {0.0};

P6_opt[2]=-px+FNorm[1]*norm2_opt[1]*len[6];

P6_opt[3]=py-FNorm[2]*norm2_opt[1]*len[6];

492 D. Hildenbrand et al.

P6_opt[4]=pz-FNorm[3]*norm2_opt[1]*len[6];

P6_opt[5]=0.5*FNorm[1]*FNorm[1]*norm2_opt[1]*norm2_opt[1]

*len[6]*len[6]+0.5*py*py+0.5*pz*pz-FNorm[1]

*norm2_opt[1]*len[6]*px+0.5*FNorm[3]*FNorm[3]

*norm2_opt[1]*norm2_opt[1]*len[6]*len[6]+0.5

*px*px+0.5*FNorm[2]*FNorm[2]*norm2_opt[1]

*norm2_opt[1]*len[6]*len[6]-FNorm[2]*norm2_opt[1]

*len[6]*py-FNorm[3]*norm2_opt[1]*len[6]*pz;

P6_opt[6]=1;

float P5_opt[32] = {0.0};

P5_opt[2]=-px+FNorm[1]*norm2_opt[1]*len[5]

+FNorm[1]*norm2_opt[1]*len[6];

P5_opt[3]=-FNorm[2]*norm2_opt[1]*len[5]

-FNorm[2]*norm2_opt[1]*len[6]+py;

P5_opt[4]=-FNorm[3]*norm2_opt[1]*len[5]

-FNorm[3]*norm2_opt[1]*len[6]+pz;

P5_opt[5]=0.5*FNorm[2]*FNorm[2]*norm2_opt[1]*norm2_opt[1]

*len[6]*len[6]+0.5*FNorm[3]*FNorm[3]*norm2_opt[1]

*norm2_opt[1]*len[6]*len[6]+0.5*FNorm[1]*FNorm[1]

*norm2_opt[1]*norm2_opt[1]*len[6]*len[6]-FNorm[3]

*norm2_opt[1]*pz*len[5]-FNorm[1]*norm2_opt[1]*px

*len[5]-FNorm[2]*norm2_opt[1]*py*len[5]+0.5*py*py

+0.5*pz*pz+FNorm[2]*FNorm[2]*norm2_opt[1]

*norm2_opt[1]*len[5]*len[6]+FNorm[1]*FNorm[1]

*norm2_opt[1]*norm2_opt[1]*len[5]*len[6]+FNorm[3]

*FNorm[3]*norm2_opt[1]*norm2_opt[1]*len[5]*len[6]

+0.5*px*px-FNorm[1]*norm2_opt[1]*len[6]*px

-FNorm[3]*norm2_opt[1]*len[6]*pz-FNorm[2]

*norm2_opt[1]*len[6]*py+0.5*FNorm[1]*FNorm[1]

*norm2_opt[1]*norm2_opt[1]*len[5]*len[5]

+0.5*FNorm[2]*FNorm[2]*norm2_opt[1]*norm2_opt[1]

*len[5]*len[5]+0.5*FNorm3*FNorm3*norm2_opt[1]

*norm2_opt[1]*len[5]*len[5];

P5_opt[6]=1;

float Pi3a_opt[32] = {0.0};

Pi3a_opt[1]=1/sqrt(hy*hy*P5_opt[4]*P5_opt[4]-2*hy*P5_opt[4]

*hz*P5_opt[3]+hz*hz*P5_opt[3]*P5_opt[3]

+P5_opt[4]*P5_opt[4]*hx*hx-2*P5_opt[4]*hx*hz

*P5_opt[2]-2*P5_opt[4]*P5_opt[4]*hx+hz*hz

*P5_opt[2]*P5_opt[2]+2*hz*P5_opt[2]*P5_opt[4]

+P5_opt[4]*P5_opt[4]+P5_opt[3]*P5_opt[3]*hx*hx

-2*P5_opt[3]*hx*hy*P5_opt[2]-2*P5_opt[3]

*P5_opt[3]*hx+hy*hy*P5_opt[2]*P5_opt[2]+2*hy

*P5_opt[2]*P5_opt[3]+P5_opt[3]*P5_opt[3]);

float Pi3b_opt[32] = {0.0};

Pi3b_opt[2]=-Pi3a_opt[1]*(-P5_opt[4]*hy+hz*P5_opt[3]);

Pi3b_opt[3]=Pi3a_opt[1]*(-P5_opt[4]*hx+hz*P5_opt[2]

+P5_opt[4]);

Pi3b_opt[4]=-Pi3a_opt[1]*(-P5_opt[3]*hx+hy*P5_opt[2]

+P5_opt[3]);

Gaalop—High Performance Parallel Computing Based on Conformal Geometric Algebra 493

As you can see, the optimized code is very complex in terms of length. Therefore

we only list the CLUCalc code for the second part of the algorithm below.

LP5P7 = *(P5 ^ P7 ^ einf);

LProj = (Pi3b.LP5P7)/Pi3b;

:P4 = pick2nd(*(S5 ^ LProj));

Pi1 = e1;

S1 = Sphere(P1,len(1));

L1 = Pi1 ^ Pi3:Red;

:P2 = pick1st(*(L1 ^ S1));

S2 = Sphere(P2,len(2));

S4 = Sphere(P4,len(3));

C3 = S2 ^ S4;

?P3 = pick1st(*(C3 ^ Pi3));

?angle(VecN3(0,1,0),P1,P2);

?angle(P1,P2,P3);

?angle(P2,P3,P4);

?angle(P3,P4,P5);

?angle(P4,P5,P6);

From the runtime performance point-of-view, our optimized C code achieved re-

sults comparable to the conventional algorithm. This is why the actual speedup that

Gaalop can provide for this inverse kinematics application will result from future

implementations on parallel platforms.

Because of the similarity of this inverse kinematics algorithms to our proof-of-

concept application (see Sect. 2.3 and [9]), we expect for an FPGA implementation

a comparable hardware speedup of about 100 times.

7 Current Status and Future Perspectives

Gaalop is currently able to handle sequential conformal geometric algebra algo-

rithms. The algorithm is currently transformed into C code, CLUCalc code, and

simple L

X formulas. The latest news can be always found on the Gaalop home-

page [7].

We are just extending Gaalop in order to handle control ﬂow with loops, condi-

tions, etc. We are also developing generators for additional output formats like Java

code, CUDA [14], and FPGA descriptions.

One focus will lie on mixed solutions handling reasonable combinations of soft-

ware and hardware implementations.

494 D. Hildenbrand et al.

8Conclusion

Geometric algebra is applicable in many different engineering scenarios and pro-

vides a straightforward and intuitive problem solving approach. With the help of our

Gaalop tool, these algorithms can be automatically transformed into high runtime

performance implementations. With these results, we are convinced that conformal

geometric algebra will be able to become more and more fruitful in a great variety

of engineering applications.

References

1. Abłamowicz, R., Fauser, B.: Clifford/bigebra, a Maple package for Clifford (co)algebra com-

2. Dorst, L., Fontijne, D., Mann, S.: Geometric Algebra for Computer Science, An Object-

Oriented Approach to Geometry. Morgan Kaufman, San Mateo (2007)

3. Fontijne, D.: Efﬁcient implementation of geometric algebra. Ph.D. thesis, University of Ams-

terdam (2007)

4. Franchini, S., Gentile, A., Grimaudo, M., Hung, C.A., Impastato, S., Sorbello, F., Vassallo, G.,

Vitabile, S.: A sliced coprocessor for native Clifford algebra operations. In: Euromico Confer-

ence on Digital System Design, Architectures, Methods and Tools (DSD) (2007)

5. Gentile, A., Segreto, S., Sorbello, F., Vassallo, G., Vitabile, S., Vullo, V.: Cliffosor, an innova-

tive FPGA-based architecture for geometric algebra. In: ERSA 2005, pp. 211–217 (2005)

6. Hildenbrand, D.: Geometric computing in computer graphics using conformal geometric al-

gebra. Comput. Graph. 29(5), 802–810 (2005)

7. Hildenbrand, D., Pitt, J.: The Gaalop homepage. Available at http://www.gaalop.de (2010)

8. Hildenbrand, D., Fontijne, D., Perwass, C., Dorst, L.: Tutorial geometric algebra and its ap-

plication to computer graphics. In: Eurographics Conference Grenoble (2004)

9. Hildenbrand, D., Fontijne, D., Wang, Y., Alexa, M., Dorst, L.: Competitive runtime perfor-

mance for inverse kinematics algorithms using Conformal geometric algebra. In: Eurographics

Conference Vienna (2006)

10. Hildenbrand, D., Lange, H., Stock, F., Koch, A.: Efﬁcient inverse kinematics algorithm

based on conformal geometric algebra using reconﬁgurable hardware. In: GRAPP Confer-

ence Madeira (2008)

11. Kasprzyk, N., Koch, A.: High-level-language compilation for reconﬁgurable computers. In:

Proceedings International Conference on Reconﬁgurable Communication-centric SoCs (Re-

CoSoC) (2005)

12. Mishra, B., Wilson, P.R.: Color edge detection hardware based on geometric algebra. In: Eu-

ropean Conference on Visual Media Production (CVMP) (2006)

13. Mishra, B. Wilson, P.R.: VLSI implementation of a geometric algebra parallel processing core.

Technical report, Electronic Systems Design Group, University of Southampton, UK (2006)

14. NVIDIA. The CUDA homepage. Available at http://www.nvidia.com/object/cuda_home.html

(2009)

15. Perwass, C.: The CLU homepage. Available at http://www.clucalc.info (2010)

16. Perwass, C., Gebken, C., Sommer, G.: Implementation of a Clifford algebra co-processor de-

sign on a ﬁeld programmable gate array. In: Ablamowicz, R. (ed.) CLIFFORD ALGEBRAS:

Application to Mathematics, Physics, and Engineering. Progress in Mathematical Physics,

pp. 561–575. Birkhäuser, Basel (2003). 6th Int. Conf. on Clifford Algebras and Applications,

Cookeville, TN

17. The RoboCup Federation. Robocup ofﬁcial site. Available at http://www.robocup.org

Some Applications of Gröbner Bases in Robotics

and Engineering

Rafał Abłamowicz

Abstract Gröbner bases in polynomial rings have numerous applications in geom-

etry, applied mathematics, and engineering. We show a few applications of Gröbner

bases in robotics, formulated in the language of Clifford algebras, and in engineer-

ing to the theory of curves, including Fermat and Bézier cubics, and interpolation

functions used in ﬁnite element theory.

1 Introduction

Gröbner bases were introduced in 1965 by Buchberger [5–8]. For an excellent ex-

position on their theory, see [10, 11, 17], while for a basic introduction with ap-

plications, see [9]. These bases gave rise to development of computer algebra sys-

tems like muMath, Maple, Mathematica, Reduce, AXIOM, CoCoCA, Macaulay,

etc. Buchberger’s Algorithm to compute Gröbner bases has been made more efﬁ-

cient [5] or replaced with another approach [12]. Algorithms to compute Gröbner

bases have been implemented, for example, in Maple [26], Singular [15, 24], and

FGb [12]. For a multitude of applications of Gröbner bases, see [8, 14] and, in par-

ticular, an online repository [16]. For a recent new application in geodesy, see [4].

2 Gröbner Basis Theory in Polynomial Rings

We follow presentation and notation from [10]. Let k[x

,...,x

] be a polyno-

mial ring in n indeterminates over a ﬁeld k. Let f

,...,f

be polynomials in

R. Abłamowicz (



)

Department of Mathematics, Tennessee Technological University, Box 5054, Cookeville,

TN 38505, USA

e-mail: rablamowicz@tntech.edu

E. Bayro-Corrochano, G. Scheuermann (eds.), Geometric Algebra Computing,

DOI 10.1007/978-1-84996-108-0_23, © Springer-Verlag London Limited 2010

495

496 R. Abłamowicz

k[x

,...,x

]. Then f

,...,f

 denotes an ideal ﬁnitely generated by the cho-

sen polynomials. We say that the polynomials form a basis of the ideal. Then, by

V(f

,...,f

) we denote an afﬁne variety deﬁned by f

,...,f

, that is, a subset of

, possibly empty, consisting of all common zeros of f

,...,f

, namely

V(f

,...,f

) =



,...,a

) ∈k

,...,a

) =0 for all 1 ≤i ≤s



In order to deﬁne Gröbner bases in the ideals of k[x

,...,x

], we ﬁrst need a con-

cept of a monomial order.

Deﬁnition 1 A monomial order on k[x

,...,x

] is any relation > on Z

≥0

{(α

,...,α

) |α

∈ Z

≥0

}, or equivalently, any relation on the set of monomials

,α∈ Z

≥0

, satisfying: (i) > is a total ordering on Z

≥0

(for any α, β ∈ Z

≥0

α>β,α= β,orβ>α); (ii) If α>βand γ ∈ Z

≥0

, then α +γ>β+γ ; (iii) >

is a well-ordering on Z

≥0

(every nonempty subset has smallest element). We will

say that x

when α>β.Let f =



be a nonzero polynomial in

k[x

,...,x

]. Then, the multidegree of f is

multideg(f ) =max



α ∈Z

≥0

: a

=0



where the maximum is taken with respect to >.Theleading coefﬁcient of f is

LC(f ) =a

multideg

(f ); the leading monomial of f is LM(f ) =x

multideg(f )

; and the

leading term of f is LT(f ) =LC(f ) ·LM(f ).

Some of the monomial orders are: (i) lexicographic (lex)

: α>

lex

β if, in the

vector difference α − β ∈ Z

, the left-most nonzero entry is positive; (ii) graded

reverse lex: α>

grevlex

β if either |α| > |β|, or |α|=|β| and in α − β ∈ Z

the

right-most nonzero entry is negative; (iii) graded inverse lex: α>

ginvlex

β if either

|α|> |β|, or |α|=|β| and in α −β ∈Z

the right-most nonzero entry is positive.

In order to divide any polynomial f by a list of polynomials f

,...,f

, we need

a generalized division algorithm.

Theorem 1 (General division algorithm) Fix a monomial order > on Z

≥0

, and

let F = (f

,...,f

) be an ordered s-tuple of polynomials. Then every f ∈

k[x

,...,x

] can be written as

f =a

+···+a

+r, (1)

where a

,r ∈k[x

,...,x

], and either r =0 or r is a linear combination, with coef-

ﬁcients in k, of monomials, none of which is divisible by any of LT(f

),...,LT(f

We call r a remainder of f on division by F. Furthermore, if a

=0, then we have

multideg(f ) ≥multideg(a

In the following, lex(x

,...,x

) will denote the lex order in which x

> ···>x

Some Applications of Gröbner Bases in Robotics and Engineering 497

The remainder r in (1) is not unique as it depends on the order of polynomi-

als in F and on the monomial order. This shortcoming of the Division Algorithm

disappears when we divide polynomials by a Gröbner basis.

Deﬁnition 2 Let I ⊂ k[x

,...,x

] be a nonzero ideal. Then, LT(I ) is the set of

leading terms of elements of I , and LT(I ) is the ideal generated by the elements

of LT(I ).

The ideal LT(I ) is an example of a monomial ideal. Dickson’s Lemma states

that every monomial ideal I ⊂ k[x

,...,x

] has a ﬁnite basis [10]. Since the poly-

nomial ring k[x

,...,x

] is Noetherian, the famous theorem of Hilbert states that

every ideal I ⊂k[x

,...,x

] is ﬁnitely generated.

Theorem 2 (Hilbert basis theorem) Every ideal I ⊂k[x

,...,x

] has a ﬁnite gen-

erating set. That is, I =g

,...,g

 for some g

,...,g

∈I.

Deﬁnition 3 Fix a monomial order. A ﬁnite subset G ={g

,...,g

} of an ideal I

is said to be a Gröbner basis if LT(g

),...,LT(g

)=LT(I ).

As a consequence of Hilbert’s theorem, every ideal I ⊂k[x

,...,x

] other than

{0} has a Gröbner basis once a monomial order has been chosen.

When dividing f by a Gröbner basis, we denote the remainder as r =

. Due

to the uniqueness of r, one gets unique coset representatives for elements in the

quotient ring k[x

,...,x

]/I : The coset representative of [f ]∈k[x

,...,x

]/I will

Gröbner bases are computed using various algorithms. The most famous one is

the Buchberger’s algorithm that uses S-polynomials. One of its many modiﬁcations

is discussed in [10], whereas [12] implements a completely different approach.

Deﬁnition 4 The S-polynomial of f

∈k[x

,...,x

] is deﬁned as S(f

) =

LT(f

)

−

LT(f

)

, where x

=lcm(LM(f

), LM(f

)), and LM(f

) is the leading

monomial of f

w.r.t. some monomial order.

Theorem 3 (Buchberger theorem) A basis {g

,...,g

}⊂I is a Gröbner basis of I

if and only if

S(g

)

=0 for all i<j.

Buchberger’s algorithm for ﬁnding a Gröbner basis implements the above crite-

rion: If F ={f

,...,f

} fails because S(f

)

=0forsomei<j,then add this

remainder to F and try again.

Gröbner bases computed with the Buchberger’s algorithm are usually too large:

A standard way to reduce them is to replace any polynomial f

with its remain-

der on division by {f

,...,f

i−1

i+1

,...,f

}, removing zero remainders, and for

Here, lcm(LM(f

), LM(f

)) denotes the least common multiple of the leading monomials

LM(f

) and LM(f