1 /* Copyright (C) 1992-1998 The Geometry Center
2 * Copyright (C) 1998-2000 Stuart Levy, Tamara Munzner, Mark Phillips
3 *
4 * This file is part of Geomview.
5 *
6 * Geomview is free software; you can redistribute it and/or modify it
7 * under the terms of the GNU Lesser General Public License as published
8 * by the Free Software Foundation; either version 2, or (at your option)
9 * any later version.
10 *
11 * Geomview is distributed in the hope that it will be useful, but
12 * WITHOUT ANY WARRANTY; without even the implied warranty of
13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
14 * Lesser General Public License for more details.
15 *
16 * You should have received a copy of the GNU Lesser General Public
17 * License along with Geomview; see the file COPYING. If not, write
18 * to the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139,
19 * USA, or visit http://www.gnu.org.
20 */
21
22 #if HAVE_CONFIG_H
23 # include "config.h"
24 #endif
25
26 #if 0
27 static char copyright[] = "Copyright (C) 1992-1998 The Geometry Center\n\
28 Copyright (C) 1998-2000 Stuart Levy, Tamara Munzner, Mark Phillips";
29 #endif
30
31 #include <stdio.h>
32 #include "options.h"
33 #include "complex.h"
34
35 complex one = {1.0, 0.0};
36 complex zero = {0.0, 0.0};
37
cplx_plus(z0,z1)38 complex cplx_plus(z0, z1)
39 complex z0, z1;
40 {
41 complex sum;
42
43 sum.real = z0.real + z1.real;
44 sum.imag = z0.imag + z1.imag;
45 return(sum);
46 }
47
48
cplx_minus(z0,z1)49 complex cplx_minus(z0, z1)
50 complex z0, z1;
51 {
52 complex diff;
53
54 diff.real = z0.real - z1.real;
55 diff.imag = z0.imag - z1.imag;
56 return(diff);
57 }
58
59
cplx_div(z0,z1)60 complex cplx_div(z0, z1)
61 complex z0, z1;
62 {
63 double mod_sq;
64 complex quotient;
65
66 mod_sq = z1.real * z1.real + z1.imag * z1.imag;
67 quotient.real = (z0.real * z1.real + z0.imag * z1.imag)/mod_sq;
68 quotient.imag = (z0.imag * z1.real - z0.real * z1.imag)/mod_sq;
69 return(quotient);
70 }
71
72
cplx_mult(z0,z1)73 complex cplx_mult(z0, z1)
74 complex z0, z1;
75 {
76 complex product;
77
78 product.real = z0.real * z1.real - z0.imag * z1.imag;
79 product.imag = z0.real * z1.imag + z0.imag * z1.real;
80 return(product);
81 }
82
83
modulus(z0)84 double modulus(z0)
85 complex z0;
86 {
87 return( sqrt(z0.real * z0.real + z0.imag * z0.imag) );
88 }
89
90
cplx_sqrt(z)91 complex cplx_sqrt(z)
92 complex z;
93 {
94 double mod,
95 arg;
96 complex result;
97
98 mod = sqrt(modulus(z));
99 if (mod == 0.0)
100 return(zero);
101 arg = 0.5 * atan2(z.imag, z.real);
102 result.real = mod * cos(arg);
103 result.imag = mod * sin(arg);
104 return(result);
105 }
106
107
108
sl2c_mult(a,b,product)109 void sl2c_mult(a, b, product)
110 sl2c_matrix a,
111 b,
112 product;
113 {
114 sl2c_matrix temp;
115
116 temp[0][0] = cplx_plus(cplx_mult(a[0][0], b[0][0]), cplx_mult(a[0][1], b[1][0]));
117 temp[0][1] = cplx_plus(cplx_mult(a[0][0], b[0][1]), cplx_mult(a[0][1], b[1][1]));
118 temp[1][0] = cplx_plus(cplx_mult(a[1][0], b[0][0]), cplx_mult(a[1][1], b[1][0]));
119 temp[1][1] = cplx_plus(cplx_mult(a[1][0], b[0][1]), cplx_mult(a[1][1], b[1][1]));
120 product[0][0] = temp[0][0];
121 product[0][1] = temp[0][1];
122 product[1][0] = temp[1][0];
123 product[1][1] = temp[1][1];
124 return;
125 }
126
127
sl2c_copy(a,b)128 void sl2c_copy(a, b)
129 sl2c_matrix a,
130 b;
131 {
132 a[0][0] = b[0][0];
133 a[0][1] = b[0][1];
134 a[1][0] = b[1][0];
135 a[1][1] = b[1][1];
136 return;
137 }
138
139
140 /* normalizes a matrix to have determinant one */
sl2c_normalize(a)141 void sl2c_normalize(a)
142 sl2c_matrix a;
143 {
144 complex det,
145 factor;
146 double arg,
147 mod;
148
149 /* compute determinant */
150 det = cplx_minus(cplx_mult(a[0][0], a[1][1]), cplx_mult(a[0][1], a[1][0]));
151 if (det.real == 0.0 && det.imag == 0.0) {
152 printf("singular sl2c_matrix\n");
153 exit(0);
154 }
155
156 /* convert to polar coordinates */
157 arg = atan2(det.imag, det.real);
158 mod = modulus(det);
159
160 /* take square root */
161 arg *= 0.5;
162 mod = sqrt(mod);
163
164 /* take reciprocal */
165 arg = -arg;
166 mod = 1.0/mod;
167
168 /* return to rectangular coordinates */
169 factor.real = mod * cos(arg);
170 factor.imag = mod * sin(arg);
171
172 /* normalize matrix */
173 a[0][0] = cplx_mult(a[0][0], factor);
174 a[0][1] = cplx_mult(a[0][1], factor);
175 a[1][0] = cplx_mult(a[1][0], factor);
176 a[1][1] = cplx_mult(a[1][1], factor);
177
178 return;
179 }
180
181
182 /* inverts a matrix; assumes determinant is already one */
sl2c_invert(a,a_inv)183 void sl2c_invert(a, a_inv)
184 sl2c_matrix a,
185 a_inv;
186 {
187 complex temp;
188
189 temp = a[0][0];
190 a_inv[0][0] = a[1][1];
191 a_inv[1][1] = temp;
192
193 a_inv[0][1].real = -a[0][1].real;
194 a_inv[0][1].imag = -a[0][1].imag;
195
196 a_inv[1][0].real = -a[1][0].real;
197 a_inv[1][0].imag = -a[1][0].imag;
198
199 return;
200 }
201
202
203 /* computes the square of the norm of a matrix */
204 /* relies on the assumption that the sl2c matrix is stored in memory as 8 consecutive doubles */
205 /* IS THIS RELIABLE? */
sl2c_norm_squared(a)206 double sl2c_norm_squared(a)
207 sl2c_matrix a;
208 {
209 int i;
210 double *p;
211 double sum;
212
213 p = (double *) a;
214 sum = 0.0;
215 for (i=8; --i>=0; ) {
216 sum += *p * *p;
217 p++;
218 }
219 return(sum);
220 }
221
222
sl2c_adjoint(a,ad_a)223 void sl2c_adjoint(a, ad_a)
224 sl2c_matrix a,
225 ad_a;
226 {
227 complex temp;
228
229 /* transpose */
230 temp = a[0][1];
231 ad_a[0][0] = a[0][0];
232 ad_a[0][1] = a[1][0];
233 ad_a[1][0] = temp;
234 ad_a[1][1] = a[1][1];
235
236 /* conjugate */
237 ad_a[0][0].imag = -ad_a[0][0].imag;
238 ad_a[0][1].imag = -ad_a[0][1].imag;
239 ad_a[1][0].imag = -ad_a[1][0].imag;
240 ad_a[1][1].imag = -ad_a[1][1].imag;
241
242 return;
243 }
244
245