#include <stdio.h>
#include <math.h>
#define N 1000//计算的层数
int main(){
double a=0.,b=0.,c=0.,s=0.,t=0.;
int m, i, j;
for(m=1; m<N+1; m++){
a=(1-(m%2)*2)/(double)m;//*6 面
b=(1-(m%2)*2)/sqrt(3*m*m);//*8 点
c=1/sqrt(2*m*m);//*12 边
for(i=1;i<m;i++){
for(j=0;j<m;j++){
a+=4*(1-((m+i+j)%2)*2)/sqrt(m*m+i*i+j*j);//*24
}
c+=2*(1-(i%2)*2)/sqrt(2*m*m+i*i);
}
t=a*6+b*8+c*12;
if (m%10 == 0){
printf("t[%4d]=%9lf\t\ts=%.15lf\t\ts2=%lf\n",m,t,s+a*3+b+c*3,s+t);
}
s+=t;
}
s=s-t+a*3+b+c*3;
printf("s[%4d]=%.15lf\n", N, s);
return 0;
}

## # Cython重写

• date: 2017-05-05
In [2]:
%%cython
cimport cython
from libc.math cimport sqrt

@cython.cdivision(True)
cdef double face, vertex, edge, E = 0.
cdef int m, i, j
for m in range(1, N + 1):
vertex = (1 - (m % 2) * 2) / sqrt(3 * m * m)
edge = 1 / sqrt(2 * m * m)
face = (1 - (m % 2) * 2) / <double > m
for i in range(1, m):
for j in range(m):
face += 4 * (1 - ((m + i + j) % 2) * 2) / \
sqrt(m * m + i * i + j * j)
edge += 2 * (1 - (i % 2) * 2) / sqrt(2 * m * m + i * i)
E += face * 6 + vertex * 8 + edge * 12
E -= face * 3 + 7 * vertex + 9 * edge
return E
In [3]:
xlabel(r"$\log_2 N$")
ylabel(r"$\log_2 \delta E$");