# Changeset 996Tweet

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

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

Location:
trunk/test/math
Files:
2 edited

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

 r989 int main(int argc, char **argv) { /* Approximate Pi using Machin's formula: 16*atan(1/5)-4*atan(1/239) */ real sum = 0.0, x0 = 5.0, x1 = 239.0; real const m0 = -x0 * x0, m1 = -x1 * x1, r16 = 16.0, r4 = 4.0; printf("Pi: "); real::R_PI.print(); printf("e: "); real::R_E.print(); printf("ln(2): "); real::R_LN2.print(); printf("sqrt(2): "); real::R_SQRT2.print(); printf("sqrt(1/2): "); real::R_SQRT1_2.print(); /* Degree 240 is required for 512-bit mantissa precision */ for (int i = 1; i < 240; i += 2) /* Gauss-Legendre computation of Pi -- doesn't work well at all, * probably because we use finite precision. */ real a = 1.0, b = sqrt((real)0.5), t = 0.25, p = 1.0; for (int i = 0; i < 3; i++) { sum += r16 / (x0 * (real)i) - r4 / (x1 * (real)i); x0 *= m0; x1 *= m1; real a2 = (a + b) * (real)0.5; real b2 = sqrt(a * b); real tmp = a - a2; real t2 = t - p * tmp * tmp; real p2 = p + p; a = a2; b = b2; t = t2; p = p2; } sum.print(); /* Bonus: compute e for fun. */ sum = 0.0; x0 = 1.0; for (int i = 1; i < 100; i++) { sum += fres(x0); x0 *= (real)i; } real sum = a + b; sum = sum * sum / ((real)4 * t); sum.print();
• ## trunk/test/math/remez.cpp

 r991 /* The order of the approximation we're looking for */ static int const ORDER = 4; static int const ORDER = 8; static real compute_pi() { /* Approximate Pi using Machin's formula: 16*atan(1/5)-4*atan(1/239) */ real sum = 0.0, x0 = 5.0, x1 = 239.0; real const m0 = -x0 * x0, m1 = -x1 * x1, r16 = 16.0, r4 = 4.0; /* Degree 240 is required for 512-bit mantissa precision */ for (int i = 1; i < 240; i += 2) { sum += r16 / (x0 * (real)i) - r4 / (x1 * (real)i); x0 *= m0; x1 *= m1; } return sum; } /* The function we want to approximate */ real myfun(real const &x) { static real myfun(real const &x) { static real const half_pi = compute_pi() * (real)0.5; static real const one = 1.0; return cos(x) - one; if (x == (real)0.0 || x == (real)-0.0) return half_pi; return sin(x * half_pi) / x; //return cos(x) - one; //return exp(x); } Matrix inv() const { Matrix a = *this, b(real(1.0)); Matrix a = *this, b((real)1.0); /* Inversion method: iterate through all columns and make sure /* Now we know the diagonal term is non-zero. Get its inverse * and use that to nullify all other terms in the column */ real x = fres(a.m[i][i]); real x = (real)1.0 / a.m[i][i]; for (int j = 0; j < N; j++) { for (int i = 0; i < ORDER + 1; i++) { zeroes[i] = real(2 * i - ORDER) / real(ORDER + 1); zeroes[i] = (real)(2 * i - ORDER) / (real)(ORDER + 1); fxn[i] = myfun(zeroes[i]); } for (int k = 0; k < ORDER + 1; k++) if (cheby[n][k]) sum += real(cheby[n][k]) * powers[k]; sum += (real)cheby[n][k] * powers[k]; mat.m[i][n] = sum; } for (int k = 0; k < ORDER + 1; k++) if (cheby[n][k]) sum += real(cheby[n][k]) * powers[k]; sum += (real)cheby[n][k] * powers[k]; mat.m[i][n] = sum; }
Note: See TracChangeset for help on using the changeset viewer.