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

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

File:
1 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
Note: See TracChangeset for help on using the changeset viewer.