// Potential of cube from 12 // clang a.c -Wmost -O3 // gcc -Wall -O3 a.c -std=c99 -lm #include #include #include #include typedef double R; typedef R U(R, R, R); static R lg(R x){return x==0.?0.:log(x);} static R c(R x, R yz, R r){ return yz*lg(x+r) - 0.5*(x==0?0:x*x*atan(yz/(x*r)));} static R u(R x, R y, R z){R r = sqrt(x*x + y*y + z*z); return c(x, y*z, r) + c(y, z*x, r) + c(z, x*y, r);} static R V(R x, R y, R z){return u(x+1, y+1, z+1) - u(x+1, y+1, z-1) - u(x+1, y-1, z+1) + u(x+1, y-1, z-1) - u(x-1, y+1, z+1) + u(x-1, y+1, z-1) + u(x-1, y-1, z+1) - u(x-1, y-1, z-1);} static void sh(R x, R y, R z) {printf("sh %10.6f %10.6f %10.6f %10.6f\n", x, y, z, V(x, y, z));} static R sq(R x){return x*x;} static void ti(R X, R Y, R Z){R s = 0; // Naive tripple integral for(R x = -0.995; x<1; x += 0.01) for(R y = -0.995; y<1; y += 0.01) for(R z = -0.995; z<1; z += 0.01) s += 1/sqrt(sq(x-X) + sq(y-Y) + sq(z-Z)); printf("ti ≈ %7.4f %7.4f %7.4f %14.11f %14.11f\n", X, Y, Z, s/1.e6, V(X, Y, Z));} static R rz(){return 0x1.p-29*random()-2;} // uniform in [-2, 2] static R poi(U V, R x, R y, R z){const R d = 1e-5; return (6*V(x, y, z) - ( V(x+d, y, z) + V(x-d, y, z) + V(x, y+d, z) + V(x, y-d, z) + V(x, y, z+d) + V(x, y, z-d)))/(d*d);} static void rs (R x, R y, R z){printf("rs %e %e\n", u(x, y, z), u(2*x, 2*y, 2*z)/4);} static void rp(U V, R x, R y, R z){R p = poi(V, x, y, z); printf("rp ≈ %10.7f %10.7f %10.7f %10.7f\n", x, y, z, p);} static R d; static void p1 (R x, R y){rp(u, x, y, d); rp(u, x, y, -d);} static void p2 (R x){p1(x, d); p1(x, -d);} static void p3 (R D){d = D; p2(D); p2(-D);} static R ri(R x, R y, R z){return 1./sqrt(x*x + y*y + z*z);} int main(){ R const pi = 4.*atan(1); int const m = 1; if(m) {rp(ri, 2.1, 0.5, -0.8); rp(ri, .3, .4, .5);} // printf("me %10.7f\n", u(1.2, 1.3, 1.4)); if(m) {R r = 1; sh(r, r, r); sh(-r, r, r); sh(r, -r, r); sh(-r, -r, r); sh(r, r, -r); sh(-r, r, -r); sh(r, -r, -r); sh(-r, -r, -r); sh(r, r, 0); sh(r, 0, 0); sh(0, 0, 0);} if(m) { ti(2, 3, 4); ti(1.4, -0.7, 1.1); ti(-0.3, 0.644, 1.6); ti(0.2135, -0.6846, 0.5544);} if(0) for(int j=0; j<10; ++j) printf("%10.5f ", rz()); if(m) printf("v %10.7f %10.7f\n", u(1.5, -0.5, 1.5), u(1.5, -0.5, -1.5)); if(m) for(int j=0; j<10000; ++j) {R x = rz(), y = rz(), z = rz(); if(0) printf("Fox %e %e %e\n", x, y, z); assert(fabs(u(x, y, z) - u(y, x, z))<6.e-15); assert(fabs(u(x, y, z) - u(y, z, x))<6.e-15); {R v = poi(V, x, y, z); char in = (fabs(x)<1. && fabs(y)<1. && fabs(z)<1.); if(in && fabs(v - 4.*pi)>2.5e-4 || !in && fabs(v)>6.e-4) printf("p %10.7f %10.7f %10.7f %10.7f\n", x, y, z, v);}} if(m) printf("x %e %e %e\n", u(-1.5, 0.01, 0), u(-0.5, 2, 2), u(-1, 2, 2)); if(m) {R d = sqrt (1.1*1.1+1.8*1.8+1.47*1.47); for(R s = 1; s<200; s *= 1.5) printf("s %e %10e %10.7f %10.7f\n", s, u(1.1*s, 1.8*s, 1.47*s), V(1.1*s, 1.8*s, 1.47*s), 8./(s*d));} if(m) {rp(V, -1, -1, -1.1); rp(V, 1, 2, 3); rp(V, .9, -.88, .91); rp(V,0,0,0);} if(m) {p3(0.3); p3(0.7); p3(1.5);} if(m) rp(u, 0.7, 0.8, 0.9); // Scale: rs(0.6, 0.53, 1.07); rs(0.8, 1.3, 1.1); return 0;} /* rp ≈ 2.1000000 0.5000000 -0.8000000 -0.0000044 rp ≈ 0.3000000 0.4000000 0.5000000 0.0000000 sh 1.000000 1.000000 1.000000 4.760155 sh -1.000000 1.000000 1.000000 4.760155 sh 1.000000 -1.000000 1.000000 4.760155 sh -1.000000 -1.000000 1.000000 4.760155 sh 1.000000 1.000000 -1.000000 4.760155 sh -1.000000 1.000000 -1.000000 4.760155 sh 1.000000 -1.000000 -1.000000 4.760155 sh -1.000000 -1.000000 -1.000000 4.760155 sh 1.000000 1.000000 0.000000 5.709041 sh 1.000000 0.000000 0.000000 7.171241 sh 0.000000 0.000000 0.000000 9.520309 ti ≈ 2.0000 3.0000 4.0000 1.48575240614 1.48575240614 ti ≈ 1.4000 -0.7000 1.1000 4.22013614283 4.22013614284 ti ≈ -0.3000 0.6440 1.6000 4.53140155949 4.53140155945 The above are suspiciously accurate. Point below is in the cube. ti ≈ 0.2135 -0.6846 0.5544 7.85824739229 7.85787207136 v -0.1402198 -0.5469983 p 0.2713243 -1.9954981 -1.9771634 0.0006395 p 0.4986326 0.0096051 0.3152501 12.5660904 p -0.6925280 0.9751509 0.1070048 12.5660193 p -1.5339681 -1.0376135 -1.6918845 -0.0007105 p -0.5361179 0.5376815 -0.3307357 12.5666588 p -0.0307070 -0.6594299 -0.4312473 12.5660904 p 1.4207088 1.7934336 1.9953992 0.0008704 p 0.9115348 -0.9624370 0.5901631 12.5666588 p -1.3081769 -1.3225770 -1.7677495 0.0006750 p 1.5285441 1.2910773 1.9872469 0.0006395 p -0.5018324 -0.4083132 -0.5509174 12.5666588 p 1.5346436 -1.4944363 -1.6653814 -0.0006040 p 0.3325511 0.9701911 -0.8742520 12.5660904 p 0.6851972 -0.6466878 -0.9657847 12.5659483 p -0.8224880 0.4946785 0.4876887 12.5660193 p 0.6079195 -0.4159997 0.5332462 12.5666588 p 1.9718304 1.3058176 -1.3954711 0.0006395 p 1.5878554 0.8765205 1.4467353 0.0006395 p 0.9114972 -0.9544526 -0.4289243 12.5660904 p 1.6526913 -1.9771252 0.2812810 -0.0006040 p 0.5773589 0.7975339 -0.9350118 12.5666588 p 1.7956587 -1.8596699 1.5463543 -0.0006040 p -0.8018375 0.1256166 -0.0947983 12.5666588 p -0.3775605 -0.3433481 -0.7868246 12.5660904 x -6.082310e-03 1.130948e+00 -1.914513e+00 s 1.000000e+00 7.070475e+00 3.1217946 3.1114276 s 1.500000e+00 2.160404e+01 2.0755553 2.0742850 s 2.250000e+00 6.142388e+01 1.3830188 1.3828567 s 3.375000e+00 1.670370e+02 0.9219255 0.9219045 s 5.062500e+00 4.407083e+02 0.6146057 0.6146030 s 7.593750e+00 1.137562e+03 0.4097357 0.4097353 s 1.139062e+01 2.887944e+03 0.2731569 0.2731569 s 1.708594e+01 7.236840e+03 0.1821046 0.1821046 s 2.562891e+01 1.794556e+04 0.1214031 0.1214031 s 3.844336e+01 4.411854e+04 0.0809354 0.0809354 s 5.766504e+01 1.076840e+05 0.0539569 0.0539569 s 8.649756e+01 2.612279e+05 0.0359713 0.0359713 s 1.297463e+02 6.303752e+05 0.0239809 0.0239809 s 1.946195e+02 1.514222e+06 0.0159872 0.0159872 rp ≈ -1.0000000 -1.0000000 -1.1000000 -0.0000711 rp ≈ 1.0000000 2.0000000 3.0000000 0.0005329 rp ≈ 0.9000000 -0.8800000 0.9100000 12.5663036 rp ≈ 0.0000000 0.0000000 0.0000000 12.5663036 rp ≈ 0.3000000 0.3000000 0.3000000 -4.4292015 rp ≈ 0.3000000 0.3000000 -0.3000000 0.4292044 rp ≈ 0.3000000 -0.3000000 0.3000000 0.4292042 rp ≈ 0.3000000 -0.3000000 -0.3000000 3.5707926 rp ≈ -0.3000000 0.3000000 0.3000000 0.4292042 rp ≈ -0.3000000 0.3000000 -0.3000000 3.5707903 rp ≈ -0.3000000 -0.3000000 0.3000000 3.5707926 rp ≈ -0.3000000 -0.3000000 -0.3000000 -7.5707929 rp ≈ 0.7000000 0.7000000 0.7000000 -4.4291948 rp ≈ 0.7000000 0.7000000 -0.7000000 0.4292033 rp ≈ 0.7000000 -0.7000000 0.7000000 0.4292033 rp ≈ 0.7000000 -0.7000000 -0.7000000 3.5707970 rp ≈ -0.7000000 0.7000000 0.7000000 0.4292033 rp ≈ -0.7000000 0.7000000 -0.7000000 3.5707970 rp ≈ -0.7000000 -0.7000000 0.7000000 3.5707970 rp ≈ -0.7000000 -0.7000000 -0.7000000 -7.5707884 rp ≈ 1.5000000 1.5000000 1.5000000 -4.4291681 rp ≈ 1.5000000 1.5000000 -1.5000000 0.4292033 rp ≈ 1.5000000 -1.5000000 1.5000000 0.4292033 rp ≈ 1.5000000 -1.5000000 -1.5000000 3.5707792 rp ≈ -1.5000000 1.5000000 1.5000000 0.4292033 rp ≈ -1.5000000 1.5000000 -1.5000000 3.5707792 rp ≈ -1.5000000 -1.5000000 1.5000000 3.5707792 rp ≈ -1.5000000 -1.5000000 -1.5000000 -7.5707618 rp ≈ 0.7000000 0.8000000 0.9000000 -4.3361581 */