source: trunk/test/math/remez-matrix.h @ 1050

Last change on this file since 1050 was 1009, checked in by sam, 11 years ago

test: some refactoring in the Remez solver to prepare multiple function
solving.

  • Property svn:keywords set to Id
File size: 2.6 KB
Line 
1//
2// Lol Engine - Sample math program: Chebyshev polynomials
3//
4// Copyright: (c) 2005-2011 Sam Hocevar <sam@hocevar.net>
5//   This program is free software; you can redistribute it and/or
6//   modify it under the terms of the Do What The Fuck You Want To
7//   Public License, Version 2, as published by Sam Hocevar. See
8//   http://sam.zoy.org/projects/COPYING.WTFPL for more details.
9//
10
11#if !defined __REMEZ_MATRIX_H__
12#define __REMEZ_MATRIX_H__
13
14template<int N> struct Matrix
15{
16    inline Matrix() {}
17
18    Matrix(real x)
19    {
20        for (int j = 0; j < N; j++)
21            for (int i = 0; i < N; i++)
22                if (i == j)
23                    m[i][j] = x;
24                else
25                    m[i][j] = 0;
26    }
27
28    /* Naive matrix inversion */
29    Matrix<N> inv() const
30    {
31        Matrix a = *this, b((real)1.0);
32
33        /* Inversion method: iterate through all columns and make sure
34         * all the terms are 1 on the diagonal and 0 everywhere else */
35        for (int i = 0; i < N; i++)
36        {
37            /* If the expected coefficient is zero, add one of
38             * the other lines. The first we meet will do. */
39            if ((double)a.m[i][i] == 0.0)
40            {
41                for (int j = i + 1; j < N; j++)
42                {
43                    if ((double)a.m[i][j] == 0.0)
44                        continue;
45                    /* Add row j to row i */
46                    for (int n = 0; n < N; n++)
47                    {
48                        a.m[n][i] += a.m[n][j];
49                        b.m[n][i] += b.m[n][j];
50                    }
51                    break;
52                }
53            }
54
55            /* Now we know the diagonal term is non-zero. Get its inverse
56             * and use that to nullify all other terms in the column */
57            real x = (real)1.0 / a.m[i][i];
58            for (int j = 0; j < N; j++)
59            {
60                if (j == i)
61                    continue;
62                real mul = x * a.m[i][j];
63                for (int n = 0; n < N; n++)
64                {
65                    a.m[n][j] -= mul * a.m[n][i];
66                    b.m[n][j] -= mul * b.m[n][i];
67                }
68            }
69
70            /* Finally, ensure the diagonal term is 1 */
71            for (int n = 0; n < N; n++)
72            {
73                a.m[n][i] *= x;
74                b.m[n][i] *= x;
75            }
76        }
77
78        return b;
79    }
80
81    void print() const
82    {
83        for (int j = 0; j < N; j++)
84        {
85            for (int i = 0; i < N; i++)
86                printf("%9.5f ", (double)m[j][i]);
87            printf("\n");
88        }
89    }
90
91    real m[N][N];
92};
93
94#endif /* __REMEZ_MATRIX_H__ */
95
Note: See TracBrowser for help on using the repository browser.