Changeset 996
 Timestamp:
 Sep 28, 2011, 7:05:00 PM (11 years ago)
 Location:
 trunk/test/math
 Files:

 2 edited
Legend:
 Unmodified
 Added
 Removed

trunk/test/math/pi.cpp
r989 r996 19 19 int main(int argc, char **argv) 20 20 { 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(); 24 26 25 /* Degree 240 is required for 512bit mantissa precision */ 26 for (int i = 1; i < 240; i += 2) 27 /* GaussLegendre 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++) 27 32 { 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; 31 39 } 32 40 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); 45 43 sum.print(); 46 44 
trunk/test/math/remez.cpp
r991 r996 20 20 21 21 /* The order of the approximation we're looking for */ 22 static int const ORDER = 4; 22 static int const ORDER = 8; 23 24 static 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 512bit 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 } 23 39 24 40 /* The function we want to approximate */ 25 real myfun(real const &x) 26 { 41 static real myfun(real const &x) 42 { 43 static real const half_pi = compute_pi() * (real)0.5; 27 44 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; 29 49 //return exp(x); 30 50 } … … 47 67 Matrix<N> inv() const 48 68 { 49 Matrix a = *this, b( real(1.0));69 Matrix a = *this, b((real)1.0); 50 70 51 71 /* Inversion method: iterate through all columns and make sure … … 73 93 /* Now we know the diagonal term is nonzero. Get its inverse 74 94 * 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]; 76 96 for (int j = 0; j < N; j++) 77 97 { … … 163 183 for (int i = 0; i < ORDER + 1; i++) 164 184 { 165 zeroes[i] = real(2 * i  ORDER) / real(ORDER + 1);185 zeroes[i] = (real)(2 * i  ORDER) / (real)(ORDER + 1); 166 186 fxn[i] = myfun(zeroes[i]); 167 187 } … … 184 204 for (int k = 0; k < ORDER + 1; k++) 185 205 if (cheby[n][k]) 186 sum += real(cheby[n][k])* powers[k];206 sum += (real)cheby[n][k] * powers[k]; 187 207 mat.m[i][n] = sum; 188 208 } … … 308 328 for (int k = 0; k < ORDER + 1; k++) 309 329 if (cheby[n][k]) 310 sum += real(cheby[n][k])* powers[k];330 sum += (real)cheby[n][k] * powers[k]; 311 331 mat.m[i][n] = sum; 312 332 }
Note: See TracChangeset
for help on using the changeset viewer.