Changeset 1117


Ignore:
Timestamp:
Dec 29, 2011, 7:36:48 PM (11 years ago)
Author:
sam
Message:

math: move the Remez algorithm implementation to the core.

Location:
trunk
Files:
1 added
1 deleted
10 edited
3 moved

Legend:

Unmodified
Added
Removed

  • trunk/src/Makefile.am

    r1106 r1117  
    33
    44liblol_a_SOURCES = \
    5     core.h matrix.cpp matrix.h tiler.cpp tiler.h dict.cpp dict.h \
     5    core.h matrix.cpp real.cpp tiler.cpp tiler.h dict.cpp dict.h \
    66    audio.cpp audio.h scene.cpp scene.h font.cpp font.h layer.cpp layer.h \
    77    map.cpp map.h entity.cpp entity.h ticker.cpp ticker.h lolgl.h \
    … …  
    1212    worldentity.cpp worldentity.h gradient.cpp gradient.h half.cpp half.h \
    1313    platform.cpp platform.h sprite.cpp sprite.h trig.cpp trig.h \
    14     real.cpp real.h \
    1514    \
    1615    lol/unit.h \
     16    lol/math/matrix.h lol/math/real.h lol/math/remez.h \
    1717    \
    1818    application/application.cpp application/application.h \

  • trunk/src/core.h

    r1106 r1117  
    6565#include "trig.h"
    6666#include "half.h"
    67 #include "real.h"
    68 #include "matrix.h"
     67#include "lol/math/real.h"
     68#include "lol/math/matrix.h"
    6969#include "numeric.h"
    7070#include "timer.h"

  • trunk/src/eglapp.h

    r863 r1117  
    1717#define __LOL_EGLAPP_H__
    1818
    19 #include "matrix.h"
     19#include "lol/math/matrix.h"
    2020
    2121namespace lol

  • trunk/src/image/image.h

    r864 r1117  
    1717#define __LOL_IMAGE_H__
    1818
    19 #include "matrix.h"
     19#include "lol/math/matrix.h"
    2020
    2121namespace lol

  • trunk/src/input.h

    r866 r1117  
    1717#define __LOL_INPUT_H__
    1818
    19 #include "matrix.h"
     19#include "lol/math/matrix.h"
    2020
    2121namespace lol

  • trunk/src/lol/math/matrix.h

    r1116 r1117  
    1414//
    1515
    16 #if !defined __LOL_MATRIX_H__
    17 #define __LOL_MATRIX_H__
     16#if !defined __LOL_MATH_MATRIX_H__
     17#define __LOL_MATH_MATRIX_H__
    1818
    1919#include <cmath>
    … …  
    2424namespace lol
    2525{
     26
     27class half;
     28class real;
    2629
    2730#define VECTOR_TYPES(tname, suffix) \
    … …  
    588591};
    589592
     593/*
     594 * Arbitrarily-sized square matrices; for now this only supports
     595 * naive inversion and is used for the Remez inversion method.
     596 */
     597
     598template<int N, typename T> struct Mat
     599{
     600    inline Mat<N, T>() {}
     601
     602    Mat(T x)
     603    {
     604        for (int j = 0; j < N; j++)
     605            for (int i = 0; i < N; i++)
     606                if (i == j)
     607                    m[i][j] = x;
     608                else
     609                    m[i][j] = 0;
     610    }
     611
     612    /* Naive matrix inversion */
     613    Mat<N, T> inv() const
     614    {
     615        Mat a = *this, b((T)1);
     616
     617        /* Inversion method: iterate through all columns and make sure
     618         * all the terms are 1 on the diagonal and 0 everywhere else */
     619        for (int i = 0; i < N; i++)
     620        {
     621            /* If the expected coefficient is zero, add one of
     622             * the other lines. The first we meet will do. */
     623            if (!a.m[i][i])
     624            {
     625                for (int j = i + 1; j < N; j++)
     626                {
     627                    if (!a.m[i][j])
     628                        continue;
     629                    /* Add row j to row i */
     630                    for (int n = 0; n < N; n++)
     631                    {
     632                        a.m[n][i] += a.m[n][j];
     633                        b.m[n][i] += b.m[n][j];
     634                    }
     635                    break;
     636                }
     637            }
     638
     639            /* Now we know the diagonal term is non-zero. Get its inverse
     640             * and use that to nullify all other terms in the column */
     641            T x = (T)1 / a.m[i][i];
     642            for (int j = 0; j < N; j++)
     643            {
     644                if (j == i)
     645                    continue;
     646                T mul = x * a.m[i][j];
     647                for (int n = 0; n < N; n++)
     648                {
     649                    a.m[n][j] -= mul * a.m[n][i];
     650                    b.m[n][j] -= mul * b.m[n][i];
     651                }
     652            }
     653
     654            /* Finally, ensure the diagonal term is 1 */
     655            for (int n = 0; n < N; n++)
     656            {
     657                a.m[n][i] *= x;
     658                b.m[n][i] *= x;
     659            }
     660        }
     661
     662        return b;
     663    }
     664
     665    T m[N][N];
     666};
     667
    590668} /* namespace lol */
    591669
    592 #endif // __LOL_MATRIX_H__
    593 
     670#endif // __LOL_MATH_MATRIX_H__
     671

  • trunk/src/lol/math/real.h

    • Property svn:keywords set to Id
    r1116 r1117  
    1414//
    1515
    16 #if !defined __LOL_REAL_H__
    17 #define __LOL_REAL_H__
     16#if !defined __LOL_MATH_REAL_H__
     17#define __LOL_MATH_REAL_H__
    1818
    1919#include <stdint.h>
    … …  
    148148} /* namespace lol */
    149149
    150 #endif // __LOL_REAL_H__
     150#endif // __LOL_MATH_REAL_H__
    151151

  • trunk/src/lol/math/remez.h

    r1116 r1117  
    11//
    2 // Lol Engine - Sample math program: Chebyshev polynomials
    3 //
    4 // Copyright: (c) 2005-2011 Sam Hocevar <sam@hocevar.net>
     2// Lol Engine
     3//
     4// Copyright: (c) 2010-2011 Sam Hocevar <sam@hocevar.net>
    55//   This program is free software; you can redistribute it and/or
    66//   modify it under the terms of the Do What The Fuck You Want To
    … …  
    99//
    1010
    11 #if !defined __REMEZ_SOLVER_H__
    12 #define __REMEZ_SOLVER_H__
    13 
    14 template<int ORDER> class RemezSolver
     11//
     12// The Remez class
     13// ---------------
     14//
     15
     16#if !defined __LOL_MATH_REMEZ_H__
     17#define __LOL_MATH_REMEZ_H__
     18
     19namespace lol
     20{
     21
     22template<int ORDER, typename T> class RemezSolver
    1523{
    1624public:
    17     typedef real RealFunc(real const &x);
     25    typedef T RealFunc(T const &x);
    1826
    1927    RemezSolver()
    … …  
    2129    }
    2230
    23     void Run(real a, real b, RealFunc *func, RealFunc *weight, int steps)
     31    void Run(T a, T b, RealFunc *func, RealFunc *weight, int steps)
    2432    {
    2533        m_func = func;
    … …  
    5058    }
    5159
    52     real ChebyEval(real const &x)
    53     {
    54         real ret = 0.0, xn = 1.0;
    55 
    56         for (int i = 0; i < ORDER + 1; i++)
    57         {
    58             real mul = 0;
     60    T ChebyEval(T const &x)
     61    {
     62        T ret = 0.0, xn = 1.0;
     63
     64        for (int i = 0; i < ORDER + 1; i++)
     65        {
     66            T mul = 0;
    5967            for (int j = 0; j < ORDER + 1; j++)
    60                 mul += coeff[j] * (real)Cheby(j, i);
     68                mul += coeff[j] * (T)Cheby(j, i);
    6169            ret += mul * xn;
    6270            xn *= x;
    … …  
    6977    {
    7078        /* Pick up x_i where error will be 0 and compute f(x_i) */
    71         real fxn[ORDER + 1];
    72         for (int i = 0; i < ORDER + 1; i++)
    73         {
    74             zeroes[i] = (real)(2 * i - ORDER) / (real)(ORDER + 1);
     79        T fxn[ORDER + 1];
     80        for (int i = 0; i < ORDER + 1; i++)
     81        {
     82            zeroes[i] = (T)(2 * i - ORDER) / (T)(ORDER + 1);
    7583            fxn[i] = Value(zeroes[i]);
    7684        }
    … …  
    7886        /* We build a matrix of Chebishev evaluations: row i contains the
    7987         * evaluations of x_i for polynomial order n = 0, 1, ... */
    80         Matrix<ORDER + 1> mat;
     88        lol::Mat<ORDER + 1, T> mat;
    8189        for (int i = 0; i < ORDER + 1; i++)
    8290        {
    8391            /* Compute the powers of x_i */
    84             real powers[ORDER + 1];
     92            T powers[ORDER + 1];
    8593            powers[0] = 1.0;
    8694            for (int n = 1; n < ORDER + 1; n++)
    … …  
    9098            for (int n = 0; n < ORDER + 1; n++)
    9199            {
    92                 real sum = 0.0;
     100                T sum = 0.0;
    93101                for (int k = 0; k < ORDER + 1; k++)
    94                     sum += (real)Cheby(n, k) * powers[k];
     102                    sum += (T)Cheby(n, k) * powers[k];
    95103                mat.m[i][n] = sum;
    96104            }
    … …  
    116124        for (int i = 0; i < ORDER + 1; i++)
    117125        {
    118             struct { real value, error; } left, right, mid;
     126            struct { T value, error; } left, right, mid;
    119127
    120128            left.value = control[i];
    … …  
    123131            right.error = ChebyEval(right.value) - Value(right.value);
    124132
    125             static real limit = real::R_1 >> 500;
     133            static T limit = (T)1 >> 500;
     134            static T zero = (T)0;
    126135            while (fabs(left.value - right.value) > limit)
    127136            {
    … …  
    129138                mid.error = ChebyEval(mid.value) - Value(mid.value);
    130139
    131                 if ((left.error < real::R_0 && mid.error < real::R_0)
    132                      || (left.error > real::R_0 && mid.error > real::R_0))
     140                if ((left.error < zero && mid.error < zero)
     141                     || (left.error > zero && mid.error > zero))
    133142                    left = mid;
    134143                else
    … …  
    147156         * compute the relative error, since its extrema are at slightly
    148157         * different locations than the absolute error’s. */
    149         real final = 0;
     158        T final = 0;
    150159
    151160        for (int i = 0; i < ORDER + 2; i++)
    152161        {
    153             real a = -1, b = 1;
     162            T a = -1, b = 1;
    154163            if (i > 0)
    155164                a = zeroes[i - 1];
    … …  
    159168            for (;;)
    160169            {
    161                 real c = a, delta = (b - a) >> 3;
    162                 real maxerror = 0;
    163                 real maxweight = 0;
     170                T c = a, delta = (b - a) >> 3;
     171                T maxerror = 0;
     172                T maxweight = 0;
    164173                int best = -1;
    165174                for (int k = 1; k <= 7; k++)
    166175                {
    167                     real error = ChebyEval(c) - Value(c);
    168                     real weight = Weight(c);
     176                    T error = ChebyEval(c) - Value(c);
     177                    T weight = Weight(c);
    169178                    if (fabs(error * maxweight) >= fabs(maxerror * weight))
    170179                    {
    … …  
    176185                }
    177186
    178                 b = a + (real)(best + 1) * delta;
    179                 a = a + (real)(best - 1) * delta;
    180 
    181                 if (b - a < (real)1e-18)
     187                b = a + (T)(best + 1) * delta;
     188                a = a + (T)(best - 1) * delta;
     189
     190                if (b - a < (T)1e-18)
    182191                {
    183                     real e = maxerror / maxweight;
     192                    T e = maxerror / maxweight;
    184193                    if (e > final)
    185194                        final = e;
    … …  
    197206    {
    198207        /* Pick up x_i where error will be 0 and compute f(x_i) */
    199         real fxn[ORDER + 2];
     208        T fxn[ORDER + 2];
    200209        for (int i = 0; i < ORDER + 2; i++)
    201210            fxn[i] = Value(control[i]);
    … …  
    203212        /* We build a matrix of Chebishev evaluations: row i contains the
    204213         * evaluations of x_i for polynomial order n = 0, 1, ... */
    205         Matrix<ORDER + 2> mat;
     214        lol::Mat<ORDER + 2, T> mat;
    206215        for (int i = 0; i < ORDER + 2; i++)
    207216        {
    208217            /* Compute the powers of x_i */
    209             real powers[ORDER + 1];
     218            T powers[ORDER + 1];
    210219            powers[0] = 1.0;
    211220            for (int n = 1; n < ORDER + 1; n++)
    … …  
    215224            for (int n = 0; n < ORDER + 1; n++)
    216225            {
    217                 real sum = 0.0;
     226                T sum = 0.0;
    218227                for (int k = 0; k < ORDER + 1; k++)
    219                     sum += (real)Cheby(n, k) * powers[k];
     228                    sum += (T)Cheby(n, k) * powers[k];
    220229                mat.m[i][n] = sum;
    221230            }
    … …  
    238247
    239248        /* Compute the error */
    240         real error = 0;
     249        T error = 0;
    241250        for (int i = 0; i < ORDER + 2; i++)
    242251            error += mat.m[ORDER + 1][i] * fxn[i];
    … …  
    265274        /* Transform Chebyshev polynomial weights into powers of X^i
    266275         * in the [-1..1] range. */
    267         real bn[ORDER + 1];
     276        T bn[ORDER + 1];
    268277
    269278        for (int i = 0; i < ORDER + 1; i++)
    … …  
    271280            bn[i] = 0;
    272281            for (int j = 0; j < ORDER + 1; j++)
    273                 bn[i] += coeff[j] * (real)Cheby(j, i);
     282                bn[i] += coeff[j] * (T)Cheby(j, i);
    274283        }
    275284
    276285        /* Transform a polynomial in the [-1..1] range into a polynomial
    277286         * in the [a..b] range. */
    278         real k1p[ORDER + 1], k2p[ORDER + 1];
    279         real an[ORDER + 1];
    280 
    281         for (int i = 0; i < ORDER + 1; i++)
    282         {
    283             k1p[i] = i ? k1p[i - 1] * m_invk1 : real::R_1;
    284             k2p[i] = i ? k2p[i - 1] * m_invk2 : real::R_1;
     287        T k1p[ORDER + 1], k2p[ORDER + 1];
     288        T an[ORDER + 1];
     289
     290        for (int i = 0; i < ORDER + 1; i++)
     291        {
     292            k1p[i] = i ? k1p[i - 1] * m_invk1 : (T)1;
     293            k2p[i] = i ? k2p[i - 1] * m_invk2 : (T)1;
    285294        }
    286295
    … …  
    289298            an[i] = 0;
    290299            for (int j = i; j < ORDER + 1; j++)
    291                 an[i] += (real)Comb(j, i) * k1p[j - i] * bn[j];
     300                an[i] += (T)Comb(j, i) * k1p[j - i] * bn[j];
    292301            an[i] *= k2p[i];
    293302        }
    … …  
    298307            if (j)
    299308                printf("+");
    300             printf("x^%i*", j);
     309            printf("x**%i*", j);
    301310            an[j].print(40);
    302311        }
    … …  
    304313    }
    305314
    306     real Value(real const &x)
     315    T Value(T const &x)
    307316    {
    308317        return m_func(x * m_k2 + m_k1);
    309318    }
    310319
    311     real Weight(real const &x)
     320    T Weight(T const &x)
    312321    {
    313322        return m_weight(x * m_k2 + m_k1);
    … …  
    315324
    316325    /* ORDER + 1 Chebyshev coefficients and 1 error value */
    317     real coeff[ORDER + 2];
     326    T coeff[ORDER + 2];
    318327    /* ORDER + 1 zeroes of the error function */
    319     real zeroes[ORDER + 1];
     328    T zeroes[ORDER + 1];
    320329    /* ORDER + 2 control points */
    321     real control[ORDER + 2];
     330    T control[ORDER + 2];
    322331
    323332private:
    324333    RealFunc *m_func, *m_weight;
    325     real m_k1, m_k2, m_invk1, m_invk2;
     334    T m_k1, m_k2, m_invk1, m_invk2;
    326335};
    327336
    328 #endif /* __REMEZ_SOLVER_H__ */
    329 
     337} /* namespace lol */
     338
     339#endif /* __LOL_MATH_REMEZ_H__ */
     340

  • trunk/src/platform/nacl/naclapp.h

    r1082 r1117  
    1717#define __LOL_NACLAPP_H__
    1818
    19 #include "matrix.h"
     19#include "lol/math/matrix.h"
    2020
    2121namespace lol

  • trunk/src/platform/ps3/ps3app.h

    r1063 r1117  
    1717#define __LOL_PS3APP_H__
    1818
    19 #include "matrix.h"
     19#include "lol/math/matrix.h"
    2020
    2121namespace lol

  • trunk/src/platform/sdl/sdlapp.h

    r1063 r1117  
    1717#define __LOL_SDLAPP_H__
    1818
    19 #include "matrix.h"
     19#include "lol/math/matrix.h"
    2020
    2121namespace lol

  • trunk/test/math/Makefile.am

    r1112 r1117  
    2323poly_DEPENDENCIES = $(top_builddir)/src/liblol.a
    2424
    25 remez_SOURCES = remez.cpp remez-matrix.h remez-solver.h
     25remez_SOURCES = remez.cpp
    2626remez_CPPFLAGS = @LOL_CFLAGS@ @PIPI_CFLAGS@
    2727remez_LDFLAGS = $(top_builddir)/src/liblol.a @LOL_LIBS@ @PIPI_LIBS@

  • trunk/test/math/remez.cpp

    r1091 r1117  
    1313#endif
    1414
    15 #include <cstring>
    1615#include <cstdio>
     16#include <cstdlib>
    1717
    1818#if USE_SDL && defined __APPLE__
    … …  
    2020#endif
    2121
    22 #include "core.h"
     22#include "lol/math/real.h"
     23#include "lol/math/matrix.h"
     24#include "lol/math/remez.h"
    2325
    2426using lol::real;
    25 
    26 #include "remez-matrix.h"
    27 #include "remez-solver.h"
     27using lol::RemezSolver;
    2828
    2929/* The function we want to approximate */
    … …  
    4242int main(int argc, char **argv)
    4343{
    44     RemezSolver<2> solver;
     44    RemezSolver<2, real> solver;
    4545    solver.Run(real::R_1 >> 400, real::R_PI_2 * real::R_PI_2, myfun, myerr, 40);
    4646
Note: See TracChangeset for help on using the changeset viewer.