1 /* Quat - A 3D fractal generation program */
2 /* Copyright (C) 1997-2000 Dirk Meyer */
3 /* (email: dirk.meyer@studserv.uni-stuttgart.de) */
4 /* mail: Dirk Meyer */
5 /* Marbacher Weg 29 */
6 /* D-71334 Waiblingen */
7 /* Germany */
8 /* */
9 /* This program is free software; you can redistribute it and/or */
10 /* modify it under the terms of the GNU General Public License */
11 /* as published by the Free Software Foundation; either version 2 */
12 /* of the License, or (at your option) any later version. */
13 /* */
14 /* This program is distributed in the hope that it will be useful, */
15 /* but WITHOUT ANY WARRANTY; without even the implied warranty of */
16 /* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the */
17 /* GNU General Public License for more details. */
18 /* */
19 /* You should have received a copy of the GNU General Public License */
20 /* along with this program; if not, write to the Free Software */
21 /* Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. */
22
23 #ifdef HAVE_CONFIG_H
24 #include <config.h>
25 #endif
26
27 #include "common.h"
28 #include "qmath.h"
29 #include <math.h>
30
q_mul(point c,const point a,const point b)31 void q_mul(point c, const point a, const point b)
32 {
33 c[0] = a[0]*b[0] - a[1]*b[1] - a[2]*b[2] - a[3]*b[3]; /* e part */
34 c[1] = a[0]*b[1] + a[1]*b[0] + a[2]*b[3] - a[3]*b[2]; /* j part */
35 c[2] = a[0]*b[2] - a[1]*b[3] + a[2]*b[0] + a[3]*b[1]; /* k part */
36 c[3] = a[0]*b[3] + a[1]*b[2] - a[2]*b[1] + a[3]*b[0]; /* l part */
37 return;
38 }
39
q_add(point c,const point a,const point b)40 void q_add(point c, const point a, const point b)
41 {
42 c[0] = a[0] + b[0];
43 c[1] = a[1] + b[1];
44 c[2] = a[2] + b[2];
45 c[3] = a[3] + b[3];
46 return;
47 }
48
q_sub(point c,const point a,const point b)49 void q_sub(point c, const point a, const point b)
50 {
51 c[0] = a[0] - b[0];
52 c[1] = a[1] - b[1];
53 c[2] = a[2] - b[2];
54 c[3] = a[3] - b[3];
55 return;
56 }
57
q_exp(point c,const point q)58 void q_exp(point c, const point q)
59 {
60 static double n, f, ex;
61 n = sqrt(q[1]*q[1] + q[2]*q[2] + q[3]*q[3]);
62 ex = exp(q[0]);
63 f = ex*sin(n);
64 if (n != 0.0) f /= n;
65
66 c[0] = ex*cos(n);
67 c[1] = f*q[1];
68 c[2] = f*q[2];
69 c[3] = f*q[3];
70 }
71
q_log(point c,const point q)72 void q_log(point c, const point q)
73 {
74 static double n, f;
75 n = sqrt(q[1]*q[1] + q[2]*q[2] + q[3]*q[3]);
76 if (n != 0.0)
77 {
78 f = atan2(n, q[0]);
79 f /= n;
80 c[0] = 0.5*log(q[0]*q[0] + n*n);
81 c[1] = f*q[1];
82 c[2] = f*q[2];
83 c[3] = f*q[3];
84 return;
85 }
86 else /* special case: real number */
87 {
88 c[0] = 0.5*log(q[0]*q[0]);
89 c[1] = atan2(0.0,q[0]);
90 c[2] = 0.0;
91 c[3] = 0.0;
92 return;
93 /* comment on asymmetry in i, j and k:
94 e^(i*pi)=-1 && e^(j*pi)=-1 && e^(k*pi)=-1
95 [ and even linear combinations of i, j and k give:
96 e^(a*i+b*j+c*k) = -1 for every a, b, c E R with:
97 sqrt(a^2+b^2+c^2) = pi ]
98 the above means that there is an infinite number of solutions to ln(-1).
99 Practically, one must choose one of them and thus break symmetry.
100 [ Even then there still is an infinite number of solutions due to
101 periodicity of sin and cos... ] */
102 }
103 }
104
q_pow(point c,const point a,const point b)105 void q_pow(point c, const point a, const point b)
106 {
107 /* if a == zero: exp(-inf*b) = 0, if b>0;
108 exp(-inf*b) = inf, if b<0
109 if b isn't real, exp(-inf*b) isn't defined, because
110 lim sin(x) for x->-inf (same with cos) doesn't exists */
111 static double an, bnp;
112 static point p;
113 an = a[0]*a[0] + a[1]*a[1] + a[2]*a[2] + a[3]*a[3];
114 bnp = b[1]*b[1] + b[2]*b[2] + b[3]*b[3];
115 if (fabs(an) < 1E-200 && (b[0] > 0.0 || bnp != 0.0) ) {
116 c[0] = 0.0; c[1] = 0.0; c[2] = 0.0; c[3] = 0.0;
117 return;
118 }
119 q_log(c, a);
120 q_mul(p, c, b);
121 q_exp(c, p);
122 return;
123 }
124