← Back to branch summary

~centralelyon2010/inkscape/imagelinks2

« back to all changes in this revision

Viewing changes to src/2geom/numeric/symmetric-matrix-fs-trace.h

  • Committer: JazzyNico
  • Date: 2011-08-29 20:25:30 UTC
  • Revision ID: nicoduf@yahoo.fr-20110829202530-6deuoz11q90usldv
Code refactoring and merging with trunk (revision 10599).

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/2geom/numeric/symmetric-matrix-fs-trace.h
 
1
/*
 
2
 * SymmetricMatrix trace
 
3
 *
 
4
 * Authors:
 
5
 *      Marco Cecchetti <mrcekets at gmail.com>
 
6
 *
 
7
 * Copyright 2009  authors
 
8
 *
 
9
 * This library is free software; you can redistribute it and/or
 
10
 * modify it either under the terms of the GNU Lesser General Public
 
11
 * License version 2.1 as published by the Free Software Foundation
 
12
 * (the "LGPL") or, at your option, under the terms of the Mozilla
 
13
 * Public License Version 1.1 (the "MPL"). If you do not alter this
 
14
 * notice, a recipient may use your version of this file under either
 
15
 * the MPL or the LGPL.
 
16
 *
 
17
 * You should have received a copy of the LGPL along with this library
 
18
 * in the file COPYING-LGPL-2.1; if not, write to the Free Software
 
19
 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
 
20
 * You should have received a copy of the MPL along with this library
 
21
 * in the file COPYING-MPL-1.1
 
22
 *
 
23
 * The contents of this file are subject to the Mozilla Public License
 
24
 * Version 1.1 (the "License"); you may not use this file except in
 
25
 * compliance with the License. You may obtain a copy of the License at
 
26
 * http://www.mozilla.org/MPL/
 
27
 *
 
28
 * This software is distributed on an "AS IS" basis, WITHOUT WARRANTY
 
29
 * OF ANY KIND, either express or implied. See the LGPL or the MPL for
 
30
 * the specific language governing rights and limitations.
 
31
 */
 
32
 
 
33
 
 
34
#ifndef _NL_TRACE_H_
 
35
#define _NL_TRACE_H_
 
36
 
 
37
 
 
38
#include <2geom/numeric/matrix.h>
 
39
#include <2geom/numeric/symmetric-matrix-fs.h>
 
40
 
 
41
 
 
42
 
 
43
 
 
44
 
 
45
namespace Geom { namespace NL {
 
46
 
 
47
 
 
48
namespace detail
 
49
{
 
50
 
 
51
/*
 
52
 *  helper routines
 
53
 */
 
54
 
 
55
inline
 
56
int sgn_prod (int x, int y)
 
57
{
 
58
    if (x == 0 || y == 0)  return 0;
 
59
    if (x == y) return 1;
 
60
    return -1;
 
61
}
 
62
 
 
63
inline
 
64
bool abs_less (double x, double y)
 
65
{
 
66
    return (std::fabs(x) < std::fabs(y));
 
67
}
 
68
 
 
69
 
 
70
/*
 
71
 * trace K-th of symmetric matrix S of order N
 
72
 */
 
73
template <size_t K, size_t N>
 
74
struct trace
 
75
{
 
76
    static
 
77
    double evaluate (const ConstBaseSymmetricMatrix<N> & S)
 
78
    {
 
79
        THROW_NOTIMPLEMENTED();
 
80
        return K;
 
81
    }
 
82
 
 
83
};
 
84
 
 
85
template <size_t N>
 
86
struct trace<1,N>
 
87
{
 
88
    static
 
89
    double evaluate (const ConstBaseSymmetricMatrix<N> & S)
 
90
    {
 
91
        double t = 0;
 
92
        for (size_t i = 0; i < N; ++i)
 
93
        {
 
94
            t += S(i,i);
 
95
        }
 
96
        return t;
 
97
    }
 
98
};
 
99
 
 
100
template <size_t N>
 
101
struct trace<N,N>
 
102
{
 
103
    static
 
104
    double evaluate (const ConstBaseSymmetricMatrix<N> & S)
 
105
    {
 
106
        Matrix M(S);
 
107
        return det(M);
 
108
    }
 
109
};
 
110
 
 
111
/*
 
112
 *  trace for symmetric matrix of order 2
 
113
 */
 
114
template <>
 
115
struct trace<1,2>
 
116
{
 
117
    static
 
118
    double evaluate (const ConstBaseSymmetricMatrix<2> & S)
 
119
    {
 
120
        return (S.get<0,0>() + S.get<1,1>());
 
121
    }
 
122
};
 
123
 
 
124
template <>
 
125
struct trace<2,2>
 
126
{
 
127
    static
 
128
    double evaluate (const ConstBaseSymmetricMatrix<2> & S)
 
129
    {
 
130
        return (S.get<0,0>() * S.get<1,1>() - S.get<0,1>() * S.get<1,0>());
 
131
    }
 
132
};
 
133
 
 
134
 
 
135
/*
 
136
 *  trace for symmetric matrix of order 3
 
137
 */
 
138
template <>
 
139
struct trace<1,3>
 
140
{
 
141
    static
 
142
    double evaluate (const ConstBaseSymmetricMatrix<3> & S)
 
143
    {
 
144
        return (S.get<0,0>() + S.get<1,1>() + S.get<2,2>());
 
145
    }
 
146
};
 
147
 
 
148
template <>
 
149
struct trace<2,3>
 
150
{
 
151
    static
 
152
    double evaluate (const ConstBaseSymmetricMatrix<3> & S)
 
153
    {
 
154
        double a00 = S.get<1,1>() * S.get<2,2>() - S.get<1,2>() * S.get<2,1>();
 
155
        double a11 = S.get<0,0>() * S.get<2,2>() - S.get<0,2>() * S.get<2,0>();
 
156
        double a22 = S.get<0,0>() * S.get<1,1>() - S.get<0,1>() * S.get<1,0>();
 
157
        return (a00 + a11 + a22);
 
158
    }
 
159
};
 
160
 
 
161
template <>
 
162
struct trace<3,3>
 
163
{
 
164
    static
 
165
    double evaluate (const ConstBaseSymmetricMatrix<3> & S)
 
166
    {
 
167
        double d = S.get<0,0>() * S.get<1,1>() * S.get<2,2>();
 
168
        d += (2 * S.get<1,0>() * S.get<2,0>() * S.get<2,1>());
 
169
        d -= (S.get<0,0>() * S.get<2,1>() * S.get<2,1>());
 
170
        d -= (S.get<1,1>() * S.get<2,0>() * S.get<2,0>());
 
171
        d -= (S.get<2,2>() * S.get<1,0>() * S.get<1,0>());
 
172
        return d;
 
173
    }
 
174
};
 
175
 
 
176
 
 
177
/*
 
178
 *  sign of trace K-th
 
179
 */
 
180
template <size_t K, size_t N>
 
181
struct trace_sgn
 
182
{
 
183
    static
 
184
    int evaluate (const ConstBaseSymmetricMatrix<N> & S)
 
185
    {
 
186
        double d = trace<K, N>::evaluate(S);
 
187
        return sgn(d);
 
188
    }
 
189
};
 
190
 
 
191
 
 
192
/*
 
193
 *  sign of trace for symmetric matrix of order 2
 
194
 */
 
195
template <>
 
196
struct trace_sgn<2,2>
 
197
{
 
198
    static
 
199
    int evaluate (const ConstBaseSymmetricMatrix<2> & S)
 
200
    {
 
201
        double m00 = S.get<0,0>();
 
202
        double m10 = S.get<1,0>();
 
203
        double m11 = S.get<1,1>();
 
204
 
 
205
        int sm00 = sgn (m00);
 
206
        int sm10 = sgn (m10);
 
207
        int sm11 = sgn (m11);
 
208
 
 
209
        if (sm10 == 0)
 
210
        {
 
211
            return sgn_prod (sm00, sm11);
 
212
        }
 
213
        else
 
214
        {
 
215
            int sm00m11 = sgn_prod (sm00, sm11);
 
216
            if (sm00m11 == 1)
 
217
            {
 
218
                int e00, e10, e11;
 
219
                double f00 = std::frexp (m00, &e00);
 
220
                double f10 = std::frexp (m10, &e10);
 
221
                double f11 = std::frexp (m11, &e11);
 
222
 
 
223
                int e0011 = e00 + e11;
 
224
                int e1010 = e10 << 1;
 
225
                int ed = e0011 - e1010;
 
226
 
 
227
                if (ed > 1)
 
228
                {
 
229
                    return 1;
 
230
                }
 
231
                else if (ed < -1)
 
232
                {
 
233
                    return -1;
 
234
                }
 
235
                else
 
236
                {
 
237
                    double d = std::ldexp (f00 * f11, ed) - f10 * f10;
 
238
                    //std::cout << "trace_sgn<2,2>: det = " << d << std::endl;
 
239
                    double eps = std::ldexp (1, -50);
 
240
                    if (std::fabs(d) < eps)  return 0;
 
241
                    return sgn (d);
 
242
                }
 
243
            }
 
244
            return -1;
 
245
        }
 
246
    }
 
247
};
 
248
 
 
249
 
 
250
/*
 
251
 *  sign of trace for symmetric matrix of order 3
 
252
 */
 
253
template <>
 
254
struct trace_sgn<2,3>
 
255
{
 
256
    static
 
257
    int evaluate (const ConstBaseSymmetricMatrix<3> & S)
 
258
    {
 
259
        double eps = std::ldexp (1, -50);
 
260
        double t[6];
 
261
 
 
262
        t[0] = S.get<1,1>() * S.get<2,2>();
 
263
        t[1] = - S.get<1,2>() * S.get<2,1>();
 
264
        t[2] = S.get<0,0>() * S.get<2,2>();
 
265
        t[3] = - S.get<0,2>() * S.get<2,0>();
 
266
        t[4] = S.get<0,0>() * S.get<1,1>();
 
267
        t[5] = - S.get<0,1>() * S.get<1,0>();
 
268
 
 
269
 
 
270
        double* maxp = std::max_element (t, t+6, abs_less);
 
271
        int em;
 
272
        std::frexp(*maxp, &em);
 
273
        double d = 0;
 
274
        for (size_t i = 0; i < 6; ++i)
 
275
        {
 
276
            d += t[i];
 
277
        }
 
278
        double r = std::fabs (std::ldexp (d, -em));  // relative error
 
279
        //std::cout << "trace_sgn<2,3>: d = " << d << std::endl;
 
280
        //std::cout << "trace_sgn<2,3>: r = " << r << std::endl;
 
281
        if (r < eps) return 0;
 
282
        if (d > 0) return 1;
 
283
        return -1;
 
284
    }
 
285
};
 
286
 
 
287
template <>
 
288
struct trace_sgn<3,3>
 
289
{
 
290
    static
 
291
    int evaluate (const ConstBaseSymmetricMatrix<3> & S)
 
292
    {
 
293
 
 
294
        double eps = std::ldexp (1, -48);
 
295
        double t[5];
 
296
 
 
297
        t[0] = S.get<0,0>() * S.get<1,1>() * S.get<2,2>();
 
298
        t[1] = 2 * S.get<1,0>() * S.get<2,0>() * S.get<2,1>();
 
299
        t[2] = -(S.get<0,0>() * S.get<2,1>() * S.get<2,1>());
 
300
        t[3] = -(S.get<1,1>() * S.get<2,0>() * S.get<2,0>());
 
301
        t[4] = -(S.get<2,2>() * S.get<1,0>() * S.get<1,0>());
 
302
 
 
303
        double* maxp = std::max_element (t, t+5, abs_less);
 
304
        int em;
 
305
        std::frexp(*maxp, &em);
 
306
        double d = 0;
 
307
        for (size_t i = 0; i < 5; ++i)
 
308
        {
 
309
            d += t[i];
 
310
        }
 
311
        //std::cout << "trace_sgn<3,3>: d = " << d << std::endl;
 
312
        double r = std::fabs (std::ldexp (d, -em));  // relative error
 
313
        //std::cout << "trace_sgn<3,3>: r = " << r << std::endl;
 
314
 
 
315
        if (r < eps)  return 0;
 
316
        if (d > 0) return 1;
 
317
        return -1;
 
318
    }
 
319
}; // end struct trace_sgn<3,3>
 
320
 
 
321
} // end namespace detail
 
322
 
 
323
 
 
324
template <size_t K, size_t N>
 
325
inline
 
326
double trace (const ConstBaseSymmetricMatrix<N> & _matrix)
 
327
{
 
328
    return detail::trace<K, N>::evaluate(_matrix);
 
329
}
 
330
 
 
331
template <size_t N>
 
332
inline
 
333
double trace (const ConstBaseSymmetricMatrix<N> & _matrix)
 
334
{
 
335
    return detail::trace<1, N>::evaluate(_matrix);
 
336
}
 
337
 
 
338
template <size_t N>
 
339
inline
 
340
double det (const ConstBaseSymmetricMatrix<N> & _matrix)
 
341
{
 
342
    return detail::trace<N, N>::evaluate(_matrix);
 
343
}
 
344
 
 
345
 
 
346
template <size_t K, size_t N>
 
347
inline
 
348
int trace_sgn (const ConstBaseSymmetricMatrix<N> & _matrix)
 
349
{
 
350
    return detail::trace_sgn<K, N>::evaluate(_matrix);
 
351
}
 
352
 
 
353
template <size_t N>
 
354
inline
 
355
int trace_sgn (const ConstBaseSymmetricMatrix<N> & _matrix)
 
356
{
 
357
    return detail::trace_sgn<1, N>::evaluate(_matrix);
 
358
}
 
359
 
 
360
template <size_t N>
 
361
inline
 
362
int det_sgn (const ConstBaseSymmetricMatrix<N> & _matrix)
 
363
{
 
364
    return detail::trace_sgn<N, N>::evaluate(_matrix);
 
365
}
 
366
 
 
367
/*
 
368
template <size_t N>
 
369
inline
 
370
size_t rank (const ConstBaseSymmetricMatrix<N> & S)
 
371
{
 
372
    THROW_NOTIMPLEMENTED();
 
373
    return 0;
 
374
}
 
375
 
 
376
template <>
 
377
inline
 
378
size_t rank<2> (const ConstBaseSymmetricMatrix<2> & S)
 
379
{
 
380
    if (S.is_zero())  return 0;
 
381
    double d = S.get<0,0>() * S.get<1,1>() - S.get<0,1>() * S.get<1,0>();
 
382
    if (d != 0)  return 2;
 
383
    return 1;
 
384
}
 
385
 
 
386
template <>
 
387
inline
 
388
size_t rank<3> (const ConstBaseSymmetricMatrix<3> & S)
 
389
{
 
390
    if (S.is_zero())  return 0;
 
391
 
 
392
    double a20 = S.get<0,1>() * S.get<1,2>() - S.get<0,2>() * S.get<1,1>();
 
393
    double a21 = S.get<0,2>() * S.get<1,0>() - S.get<0,0>() * S.get<1,2>();
 
394
    double a22 = S.get<0,0>() * S.get<1,1>() - S.get<0,1>() * S.get<1,0>();
 
395
    double d = a20 * S.get<2,0>() + a21 * S.get<2,1>() + a22 * S.get<2,2>();
 
396
 
 
397
    if (d != 0)  return 3;
 
398
 
 
399
    if (a20 != 0 || a21 != 0 || a22 != 0)  return 2;
 
400
 
 
401
    double a00 = S.get<1,1>() * S.get<2,2>() - S.get<1,2>() * S.get<2,1>();
 
402
    if (a00 != 0)  return 2;
 
403
 
 
404
    double a10 = S.get<0,2>() * S.get<2,1>() - S.get<0,1>() * S.get<2,2>();
 
405
    if (a10 != 0)  return 2;
 
406
 
 
407
    double a11 = S.get<0,0>() * S.get<2,2>() - S.get<0,2>() * S.get<2,0>();
 
408
    if (a11 != 0)  return 2;
 
409
 
 
410
    return 1;
 
411
}
 
412
*/
 
413
 
 
414
} /* end namespace NL*/ } /* end namespace Geom*/
 
415
 
 
416
 
 
417
 
 
418
 
 
419
#endif // _NL_TRACE_H_
 
420
 
 
421
 
 
422
 
 
423
 
 
424
/*
 
425
  Local Variables:
 
426
  mode:c++
 
427
  c-file-style:"stroustrup"
 
428
  c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
 
429
  indent-tabs-mode:nil
 
430
  fill-column:99
 
431
  End:
 
432
*/
 
433
// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:fileencoding=utf-8:textwidth=99 :