Changeset 996 for trunk/test/math


Ignore:
Timestamp:
Sep 28, 2011, 7:05:00 PM (10 years ago)
Author:
sam
Message:

test: minor updates to the Pi program (now almost deprecated) and the
Remez exchange program.

Location:
trunk/test/math
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • trunk/test/math/pi.cpp

    r989 r996  
    1919int main(int argc, char **argv)
    2020{
    21     /* Approximate Pi using Machin's formula: 16*atan(1/5)-4*atan(1/239) */
    22     real sum = 0.0, x0 = 5.0, x1 = 239.0;
    23     real const m0 = -x0 * x0, m1 = -x1 * x1, r16 = 16.0, r4 = 4.0;
     21    printf("Pi: "); real::R_PI.print();
     22    printf("e: "); real::R_E.print();
     23    printf("ln(2): "); real::R_LN2.print();
     24    printf("sqrt(2): "); real::R_SQRT2.print();
     25    printf("sqrt(1/2): "); real::R_SQRT1_2.print();
    2426
    25     /* Degree 240 is required for 512-bit mantissa precision */
    26     for (int i = 1; i < 240; i += 2)
     27    /* Gauss-Legendre computation of Pi -- doesn't work well at all,
     28     * probably because we use finite precision. */
     29    real a = 1.0, b = sqrt((real)0.5), t = 0.25, p = 1.0;
     30
     31    for (int i = 0; i < 3; i++)
    2732    {
    28         sum += r16 / (x0 * (real)i) - r4 / (x1 * (real)i);
    29         x0 *= m0;
    30         x1 *= m1;
     33        real a2 = (a + b) * (real)0.5;
     34        real b2 = sqrt(a * b);
     35        real tmp = a - a2;
     36        real t2 = t - p * tmp * tmp;
     37        real p2 = p + p;
     38        a = a2; b = b2; t = t2; p = p2;
    3139    }
    3240
    33     sum.print();
    34 
    35     /* Bonus: compute e for fun. */
    36     sum = 0.0;
    37     x0 = 1.0;
    38 
    39     for (int i = 1; i < 100; i++)
    40     {
    41         sum += fres(x0);
    42         x0 *= (real)i;
    43     }
    44 
     41    real sum = a + b;
     42    sum = sum * sum / ((real)4 * t);
    4543    sum.print();
    4644
  • trunk/test/math/remez.cpp

    r991 r996  
    2020
    2121/* The order of the approximation we're looking for */
    22 static int const ORDER = 4;
     22static int const ORDER = 8;
     23
     24static real compute_pi()
     25{
     26    /* Approximate Pi using Machin's formula: 16*atan(1/5)-4*atan(1/239) */
     27    real sum = 0.0, x0 = 5.0, x1 = 239.0;
     28    real const m0 = -x0 * x0, m1 = -x1 * x1, r16 = 16.0, r4 = 4.0;
     29
     30    /* Degree 240 is required for 512-bit mantissa precision */
     31    for (int i = 1; i < 240; i += 2)
     32    {
     33        sum += r16 / (x0 * (real)i) - r4 / (x1 * (real)i);
     34        x0 *= m0;
     35        x1 *= m1;
     36    }
     37    return sum;
     38}
    2339
    2440/* The function we want to approximate */
    25 real myfun(real const &x)
    26 {
     41static real myfun(real const &x)
     42{
     43    static real const half_pi = compute_pi() * (real)0.5;
    2744    static real const one = 1.0;
    28     return cos(x) - one;
     45    if (x == (real)0.0 || x == (real)-0.0)
     46        return half_pi;
     47    return sin(x * half_pi) / x;
     48    //return cos(x) - one;
    2949    //return exp(x);
    3050}
     
    4767    Matrix<N> inv() const
    4868    {
    49         Matrix a = *this, b(real(1.0));
     69        Matrix a = *this, b((real)1.0);
    5070
    5171        /* Inversion method: iterate through all columns and make sure
     
    7393            /* Now we know the diagonal term is non-zero. Get its inverse
    7494             * and use that to nullify all other terms in the column */
    75             real x = fres(a.m[i][i]);
     95            real x = (real)1.0 / a.m[i][i];
    7696            for (int j = 0; j < N; j++)
    7797            {
     
    163183    for (int i = 0; i < ORDER + 1; i++)
    164184    {
    165         zeroes[i] = real(2 * i - ORDER) / real(ORDER + 1);
     185        zeroes[i] = (real)(2 * i - ORDER) / (real)(ORDER + 1);
    166186        fxn[i] = myfun(zeroes[i]);
    167187    }
     
    184204            for (int k = 0; k < ORDER + 1; k++)
    185205                if (cheby[n][k])
    186                     sum += real(cheby[n][k]) * powers[k];
     206                    sum += (real)cheby[n][k] * powers[k];
    187207            mat.m[i][n] = sum;
    188208        }
     
    308328            for (int k = 0; k < ORDER + 1; k++)
    309329                if (cheby[n][k])
    310                     sum += real(cheby[n][k]) * powers[k];
     330                    sum += (real)cheby[n][k] * powers[k];
    311331            mat.m[i][n] = sum;
    312332        }
Note: See TracChangeset for help on using the changeset viewer.