modernc.org/ccgo/v3@v3.16.14/lib/testdata/CompCert-3.6/test/c/almabench.c (about) 1 /* 2 almabench 3 Standard C version 4 17 April 2003 5 6 Written by Scott Robert Ladd. 7 No rights reserved. This is public domain software, for use by anyone. 8 9 A number-crunching benchmark that can be used as a fitness test for 10 evolving optimal compiler options via genetic algorithm. 11 12 This program calculates the daily ephemeris (at noon) for the years 13 2000-2099 using an algorithm developed by J.L. Simon, P. Bretagnon, J. 14 Chapront, M. Chapront-Touze, G. Francou and J. Laskar of the Bureau des 15 Longitudes, Paris, France), as detailed in Astronomy & Astrophysics 16 282, 663 (1994) 17 18 Note that the code herein is design for the purpose of testing 19 computational performance; error handling and other such "niceties" 20 is virtually non-existent. 21 22 Actual benchmark results can be found at: 23 http://www.coyotegulch.com 24 25 Please do not use this information or algorithm in any way that might 26 upset the balance of the universe or otherwise cause planets to impact 27 upon one another. 28 */ 29 30 #include <stdio.h> 31 #include <stdlib.h> 32 #include <math.h> 33 #include <string.h> 34 35 #define PI 3.14159265358979323846 36 #define J2000 2451545.0 37 #define JCENTURY 36525.0 38 #define JMILLENIA 365250.0 39 #define TWOPI (2.0 * PI) 40 #define A2R (PI / 648000.0) 41 #define R2H (12.0 / PI) 42 #define R2D (180.0 / PI) 43 #define GAUSSK 0.01720209895 44 #define TEST_LOOPS 20 45 #define TEST_LENGTH 36525 46 #define sineps 0.3977771559319137 47 #define coseps 0.9174820620691818 48 49 const double amas [8] = { 6023600.0, 408523.5, 328900.5, 3098710.0, 1047.355, 3498.5, 22869.0, 19314.0 }; 50 51 const double a [8][3] = 52 { { 0.3870983098, 0, 0 }, 53 { 0.7233298200, 0, 0 }, 54 { 1.0000010178, 0, 0 }, 55 { 1.5236793419, 3e-10, 0 }, 56 { 5.2026032092, 19132e-10, -39e-10 }, 57 { 9.5549091915, -0.0000213896, 444e-10 }, 58 { 19.2184460618, -3716e-10, 979e-10 }, 59 { 30.1103868694, -16635e-10, 686e-10 } }; 60 61 const double dlm[8][3] = 62 { { 252.25090552, 5381016286.88982, -1.92789 }, 63 { 181.97980085, 2106641364.33548, 0.59381 }, 64 { 100.46645683, 1295977422.83429, -2.04411 }, 65 { 355.43299958, 689050774.93988, 0.94264 }, 66 { 34.35151874, 109256603.77991, -30.60378 }, 67 { 50.07744430, 43996098.55732, 75.61614 }, 68 { 314.05500511, 15424811.93933, -1.75083 }, 69 { 304.34866548, 7865503.20744, 0.21103 } }; 70 71 const double e[8][3] = 72 { { 0.2056317526, 0.0002040653, -28349e-10 }, 73 { 0.0067719164, -0.0004776521, 98127e-10 }, 74 { 0.0167086342, -0.0004203654, -0.0000126734 }, 75 { 0.0934006477, 0.0009048438, -80641e-10 }, 76 { 0.0484979255, 0.0016322542, -0.0000471366 }, 77 { 0.0555481426, -0.0034664062, -0.0000643639 }, 78 { 0.0463812221, -0.0002729293, 0.0000078913 }, 79 { 0.0094557470, 0.0000603263, 0 } }; 80 81 const double pi[8][3] = 82 { { 77.45611904, 5719.11590, -4.83016 }, 83 { 131.56370300, 175.48640, -498.48184 }, 84 { 102.93734808, 11612.35290, 53.27577 }, 85 { 336.06023395, 15980.45908, -62.32800 }, 86 { 14.33120687, 7758.75163, 259.95938 }, 87 { 93.05723748, 20395.49439, 190.25952 }, 88 { 173.00529106, 3215.56238, -34.09288 }, 89 { 48.12027554, 1050.71912, 27.39717 } }; 90 91 const double dinc[8][3] = 92 { { 7.00498625, -214.25629, 0.28977 }, 93 { 3.39466189, -30.84437, -11.67836 }, 94 { 0, 469.97289, -3.35053 }, 95 { 1.84972648, -293.31722, -8.11830 }, 96 { 1.30326698, -71.55890, 11.95297 }, 97 { 2.48887878, 91.85195, -17.66225 }, 98 { 0.77319689, -60.72723, 1.25759 }, 99 { 1.76995259, 8.12333, 0.08135 } }; 100 101 const double omega[8][3] = 102 { { 48.33089304, -4515.21727, -31.79892 }, 103 { 76.67992019, -10008.48154, -51.32614 }, 104 { 174.87317577, -8679.27034, 15.34191 }, 105 { 49.55809321, -10620.90088, -230.57416 }, 106 { 100.46440702, 6362.03561, 326.52178 }, 107 { 113.66550252, -9240.19942, -66.23743 }, 108 { 74.00595701, 2669.15033, 145.93964 }, 109 { 131.78405702, -221.94322, -0.78728 } }; 110 111 const double kp[8][9] = 112 { { 69613.0, 75645.0, 88306.0, 59899.0, 15746.0, 71087.0, 142173.0, 3086.0, 0.0 }, 113 { 21863.0, 32794.0, 26934.0, 10931.0, 26250.0, 43725.0, 53867.0, 28939.0, 0.0 }, 114 { 16002.0, 21863.0, 32004.0, 10931.0, 14529.0, 16368.0, 15318.0, 32794.0, 0.0 }, 115 { 6345.0, 7818.0, 15636.0, 7077.0, 8184.0, 14163.0, 1107.0, 4872.0, 0.0 }, 116 { 1760.0, 1454.0, 1167.0, 880.0, 287.0, 2640.0, 19.0, 2047.0, 1454.0 }, 117 { 574.0, 0.0, 880.0, 287.0, 19.0, 1760.0, 1167.0, 306.0, 574.0 }, 118 { 204.0, 0.0, 177.0, 1265.0, 4.0, 385.0, 200.0, 208.0, 204.0 }, 119 { 0.0, 102.0, 106.0, 4.0, 98.0, 1367.0, 487.0, 204.0, 0.0 } }; 120 121 const double ca[8][9] = 122 { { 4.0, -13.0, 11.0, -9.0, -9.0, -3.0, -1.0, 4.0, 0.0 }, 123 { -156.0, 59.0, -42.0, 6.0, 19.0, -20.0, -10.0, -12.0, 0.0 }, 124 { 64.0, -152.0, 62.0, -8.0, 32.0, -41.0, 19.0, -11.0, 0.0 }, 125 { 124.0, 621.0, -145.0, 208.0, 54.0, -57.0, 30.0, 15.0, 0.0 }, 126 { -23437.0, -2634.0, 6601.0, 6259.0, -1507.0, -1821.0, 2620.0, -2115.0,-1489.0 }, 127 { 62911.0,-119919.0, 79336.0, 17814.0,-24241.0, 12068.0, 8306.0, -4893.0, 8902.0 }, 128 { 389061.0,-262125.0,-44088.0, 8387.0,-22976.0, -2093.0, -615.0, -9720.0, 6633.0 }, 129 { -412235.0,-157046.0,-31430.0, 37817.0, -9740.0, -13.0, -7449.0, 9644.0, 0.0 } }; 130 131 const double sa[8][9] = 132 { { -29.0, -1.0, 9.0, 6.0, -6.0, 5.0, 4.0, 0.0, 0.0 }, 133 { -48.0, -125.0, -26.0, -37.0, 18.0, -13.0, -20.0, -2.0, 0.0 }, 134 { -150.0, -46.0, 68.0, 54.0, 14.0, 24.0, -28.0, 22.0, 0.0 }, 135 { -621.0, 532.0, -694.0, -20.0, 192.0, -94.0, 71.0, -73.0, 0.0 }, 136 { -14614.0,-19828.0, -5869.0, 1881.0, -4372.0, -2255.0, 782.0, 930.0, 913.0 }, 137 { 139737.0, 0.0, 24667.0, 51123.0, -5102.0, 7429.0, -4095.0, -1976.0,-9566.0 }, 138 { -138081.0, 0.0, 37205.0,-49039.0,-41901.0,-33872.0,-27037.0,-12474.0,18797.0 }, 139 { 0.0, 28492.0,133236.0, 69654.0, 52322.0,-49577.0,-26430.0, -3593.0, 0.0 } }; 140 141 const double kq[8][10] = 142 { { 3086.0, 15746.0, 69613.0, 59899.0, 75645.0, 88306.0, 12661.0, 2658.0, 0.0, 0.0 }, 143 { 21863.0, 32794.0, 10931.0, 73.0, 4387.0, 26934.0, 1473.0, 2157.0, 0.0, 0.0 }, 144 { 10.0, 16002.0, 21863.0, 10931.0, 1473.0, 32004.0, 4387.0, 73.0, 0.0, 0.0 }, 145 { 10.0, 6345.0, 7818.0, 1107.0, 15636.0, 7077.0, 8184.0, 532.0, 10.0, 0.0 }, 146 { 19.0, 1760.0, 1454.0, 287.0, 1167.0, 880.0, 574.0, 2640.0, 19.0,1454.0 }, 147 { 19.0, 574.0, 287.0, 306.0, 1760.0, 12.0, 31.0, 38.0, 19.0, 574.0 }, 148 { 4.0, 204.0, 177.0, 8.0, 31.0, 200.0, 1265.0, 102.0, 4.0, 204.0 }, 149 { 4.0, 102.0, 106.0, 8.0, 98.0, 1367.0, 487.0, 204.0, 4.0, 102.0 } }; 150 151 const double cl[8][10] = 152 { { 21.0, -95.0, -157.0, 41.0, -5.0, 42.0, 23.0, 30.0, 0.0, 0.0 }, 153 { -160.0, -313.0, -235.0, 60.0, -74.0, -76.0, -27.0, 34.0, 0.0, 0.0 }, 154 { -325.0, -322.0, -79.0, 232.0, -52.0, 97.0, 55.0, -41.0, 0.0, 0.0 }, 155 { 2268.0, -979.0, 802.0, 602.0, -668.0, -33.0, 345.0, 201.0, -55.0, 0.0 }, 156 { 7610.0, -4997.0,-7689.0,-5841.0,-2617.0, 1115.0, -748.0, -607.0, 6074.0, 354.0 }, 157 { -18549.0, 30125.0,20012.0, -730.0, 824.0, 23.0, 1289.0, -352.0,-14767.0,-2062.0 }, 158 { -135245.0,-14594.0, 4197.0,-4030.0,-5630.0,-2898.0, 2540.0, -306.0, 2939.0, 1986.0 }, 159 { 89948.0, 2103.0, 8963.0, 2695.0, 3682.0, 1648.0, 866.0, -154.0, -1963.0, -283.0 } }; 160 161 const double sl[8][10] = 162 { { -342.0, 136.0, -23.0, 62.0, 66.0, -52.0, -33.0, 17.0, 0.0, 0.0 }, 163 { 524.0, -149.0, -35.0, 117.0, 151.0, 122.0, -71.0, -62.0, 0.0, 0.0 }, 164 { -105.0, -137.0, 258.0, 35.0, -116.0, -88.0, -112.0, -80.0, 0.0, 0.0 }, 165 { 854.0, -205.0, -936.0, -240.0, 140.0, -341.0, -97.0, -232.0, 536.0, 0.0 }, 166 { -56980.0, 8016.0, 1012.0, 1448.0,-3024.0,-3710.0, 318.0, 503.0, 3767.0, 577.0 }, 167 { 138606.0,-13478.0,-4964.0, 1441.0,-1319.0,-1482.0, 427.0, 1236.0, -9167.0,-1918.0 }, 168 { 71234.0,-41116.0, 5334.0,-4935.0,-1848.0, 66.0, 434.0,-1748.0, 3780.0, -701.0 }, 169 { -47645.0, 11647.0, 2166.0, 3194.0, 679.0, 0.0, -244.0, -419.0, -2531.0, 48.0 } }; 170 171 //--------------------------------------------------------------------------- 172 // Normalize angle into the range -pi <= A < +pi. 173 double anpm (double a) 174 { 175 double w = fmod(a,TWOPI); 176 177 if (fabs(w) >= PI) 178 w = w - ((a < 0) ? -TWOPI : TWOPI); 179 180 return w; 181 } 182 183 //--------------------------------------------------------------------------- 184 // The reference frame is equatorial and is with respect to the 185 // mean equator and equinox of epoch j2000. 186 void planetpv (double epoch[2], int np, double pv[2][3]) 187 { 188 // working storage 189 int k; 190 double t, da, dl, de, dp, di, doh, dmu, arga, argl, am; 191 double ae, dae, ae2, at, r, v, si2, xq, xp, tl, xsw; 192 double xcw, xm2, xf, ci2, xms, xmc, xpxq2, x, y, z; 193 194 // time: julian millennia since j2000. 195 t = ((epoch[0] - J2000) + epoch[1]) / JMILLENIA; 196 197 // compute the mean elements. 198 da = a[np][0] + (a[np][1] + a[np][2] * t ) * t; 199 dl = (3600.0 * dlm[np][0] + (dlm[np][1] + dlm[np][2] * t ) * t ) * A2R; 200 de = e[np][0] + (e[np][1] + e[np][2] * t ) * t; 201 dp = anpm((3600.0 * pi[np][0] + (pi[np][1] + pi[np][2] * t ) * t ) * A2R ); 202 di = (3600.0 * dinc[np][0] + (dinc[np][1] + dinc[np][2] * t ) * t ) * A2R; 203 doh = anpm((3600.0 * omega[np][0] + (omega[np][1] + omega[np][2] * t ) * t ) * A2R ); 204 205 // apply the trigonometric terms. 206 dmu = 0.35953620 * t; 207 208 for (k = 0; k < 8; ++k) 209 { 210 arga = kp[np][k] * dmu; 211 argl = kq[np][k] * dmu; 212 da = da + (ca[np][k] * cos(arga) + sa[np][k] * sin(arga)) * 0.0000001; 213 dl = dl + (cl[np][k] * cos(argl) + sl[np][k] * sin(argl)) * 0.0000001; 214 } 215 216 arga = kp[np][8] * dmu; 217 da = da + t * (ca[np][8] * cos(arga) + sa[np][8] * sin(arga)) * 0.0000001; 218 219 for (k = 8; k <= 9; ++k) 220 { 221 argl = kq[np][k] * dmu; 222 dl = dl + t * ( cl[np][k] * cos(argl) + sl[np][k] * sin(argl) ) * 0.0000001; 223 } 224 225 dl = fmod(dl,TWOPI); 226 227 // iterative solution of kepler's equation to get eccentric anomaly. 228 am = dl - dp; 229 ae = am + de * sin(am); 230 k = 0; 231 232 while (1) 233 { 234 dae = (am - ae + de * sin(ae)) / (1.0 - de * cos(ae)); 235 ae = ae + dae; 236 k = k + 1; 237 238 if ((k >= 10) || (fabs(dae) < 1e-12)) 239 break; 240 } 241 242 // true anomaly. 243 ae2 = ae / 2.0; 244 at = 2.0 * atan2(sqrt((1.0 + de)/(1.0 - de)) * sin(ae2), cos(ae2)); 245 246 // distance (au) and speed (radians per day). 247 r = da * (1.0 - de * cos(ae)); 248 v = GAUSSK * sqrt((1.0 + 1.0 / amas[np] ) / (da * da * da)); 249 250 si2 = sin(di / 2.0); 251 xq = si2 * cos(doh); 252 xp = si2 * sin(doh); 253 tl = at + dp; 254 xsw = sin(tl); 255 xcw = cos(tl); 256 xm2 = 2.0 * (xp * xcw - xq * xsw ); 257 xf = da / sqrt(1.0 - de*de); 258 ci2 = cos(di / 2.0); 259 xms = (de * sin(dp) + xsw) * xf; 260 xmc = (de * cos(dp) + xcw) * xf; 261 xpxq2 = 2.0 * xp * xq; 262 263 // position (j2000 ecliptic x,y,z in au). 264 x = r * (xcw - xm2 * xp); 265 y = r * (xsw + xm2 * xq); 266 z = r * (-xm2 * ci2); 267 268 // rotate to equatorial. 269 pv[0][0] = x; 270 pv[0][1] = y * coseps - z * sineps; 271 pv[0][2] = y * sineps + z * coseps; 272 273 // velocity (j2000 ecliptic xdot,ydot,zdot in au/d). 274 x = v * ((-1.0 + 2.0 * xp * xp) * xms + xpxq2 * xmc); 275 y = v * (( 1.0 - 2.0 * xq * xq ) * xmc - xpxq2 * xms); 276 z = v * (2.0 * ci2 * (xp * xms + xq * xmc)); 277 278 // rotate to equatorial. 279 pv[1][0] = x; 280 pv[1][1] = y * coseps - z * sineps; 281 pv[1][2] = y * sineps + z * coseps; 282 } 283 284 //--------------------------------------------------------------------------- 285 // Computes RA, Declination, and distance from a state vector returned by 286 // planetpv. 287 void radecdist(double state[2][3], double rdd[3]) 288 { 289 // distance 290 rdd[2] = sqrt(state[0][0] * state[0][0] + state[0][1] * state[0][1] + state[0][2] * state[0][2]); 291 292 // RA 293 rdd[0] = atan2(state[0][1], state[0][0]) * R2H; 294 if (rdd[0] < 0.0) rdd[0] += 24.0; 295 296 // Declination 297 rdd[1] = asin(state[0][2] / rdd[2]) * R2D; 298 } 299 300 /* Test harness */ 301 302 static void test(void) 303 { 304 int p; 305 double jd[2]; 306 double pv[2][3]; 307 double position[3]; 308 309 jd[0] = J2000; 310 jd[1] = 0.0; 311 for (p = 0; p < 8; ++p) 312 { 313 planetpv(jd,p,pv); 314 radecdist(pv,position); 315 printf("p = %d position[0] = %g position[1] = %g\n", 316 p, position[0], position[1]); 317 } 318 } 319 320 321 static void bench(int nloops) 322 { 323 int i, n, p; 324 double jd[2]; 325 double pv[2][3]; 326 double position[3]; 327 328 for (i = 0; i < nloops; ++i) 329 { 330 jd[0] = J2000; 331 jd[1] = 0.0; 332 333 for (n = 0; n < TEST_LENGTH; ++n) 334 { 335 jd[0] += 1.0; 336 337 for (p = 0; p < 8; ++p) 338 { 339 planetpv(jd,p,pv); 340 radecdist(pv,position); 341 } 342 } 343 } 344 } 345 346 int main(int argc, char ** argv) 347 { 348 test(); 349 bench(1); 350 return 0; 351 }