I've been working more on TOS's Li based pi(x) approximation code.
I've been trying to optimize it in c. It seems that I need someone
more knowledgeable than myself in c to point out some simple mistake I
am making that is preventing the code from giving the correct answer.
I tried copying and pasting the li function below into a Sage cell and
"Pythonizing" it, and I got the expected answer, but as is in a c
file, with commands being piped in the command line, I get unexpected
results in Sage. I believe that the c li function below is correct-I
think I am just making some sort of simple error in moving data
around. I have already asked one of my peers in my number theory
class, and we are both stuck.
Sage:
import subprocess
def fastli(x,y):
thesubprocess=subprocess.Popen(args=["/home/kstueve/Desktop/
riemannpi/li",str(x),str(y)],stdout=subprocess.PIPE)
return float(thesubprocess.communicate()[int(0)])
c:
#include <stdio.h>
#include <math.h>
//This is some rough code to calculate li(x^(0.5+it))+li(x^(0.5-it))
//Authors:
//Tomas Oliveira e Silva: Original Version
//Kevin Stueve: Rewrote and optimized in c
//Released under GPL
const double zeros[]={
14.134725142,21.022039639,25.01085758,30.424876126,32.935061588,
21.022039639,25.01085758,30.424876126,32.935061588,37.586178159,
25.01085758,30.424876126,32.935061588,37.586178159,40.918719012,
30.424876126,32.935061588,37.586178159,40.918719012,43.327073281,
32.935061588,37.586178159,40.918719012,43.327073281,48.005150881,
37.586178159,40.918719012,43.327073281,48.005150881,49.773832478,
40.918719012,43.327073281,48.005150881,49.773832478,52.970321478,
43.327073281,48.005150881,49.773832478,52.970321478,56.446247697,
48.005150881,49.773832478,52.970321478,56.446247697,59.347044003,
49.773832478,52.970321478,56.446247697,59.347044003,60.831778525,
52.970321478,56.446247697,59.347044003,60.831778525,65.112544048,
56.446247697,59.347044003,60.831778525,65.112544048,67.079810529,
59.347044003,60.831778525,65.112544048,67.079810529,69.546401711,
60.831778525,65.112544048,67.079810529,69.546401711,72.067157674,
65.112544048,67.079810529,69.546401711,72.067157674,75.704690699,
};
double li(double x, double t);
int main(int argc,char **argv){
//printf("zero%f",zeros[50000]);
double x=atof(argv[1]);
long int num = atol(argv[2]);
//printf("%f\n",x);
//printf("%ld\n",num);
double result=0;
long int i=0;
for (i=0;i<num;i++){
//printf("%f\n",zeros[i]);
result+=li(x,zeros[i]);
}
printf("%f",result);
return 0;
}
double li(double x, double t) {
double logx = log(x);
double u_real = 0.5 * logx;
double u_imag = t * logx;
double u_norm = u_real * u_real + u_imag * u_imag;
double so_real_num = 2 * sqrt(x) * cos(t * logx);
double so_imag_num = 2 * sqrt(x) * sin(t * logx);
double so_real = (so_real_num * u_real + so_imag_num * u_imag) /
u_norm;
double so_imag = (so_imag_num * u_real - so_real_num * u_imag) /
u_norm;
double tol=pow(10,-9)/(sqrt(so_real * so_real + so_imag *
so_imag));
double tol_squared=tol*tol;
double s1_real = 1;
double s1_imag = 0;
double s2_real = 1;
double s2_imag = 0;
double s2_real_temp;
double s2_imag_temp;
int k = 1;
double norm,lastnorm;
for(;;) {
if (k>100) {
printf("error%f\n",t);
return;
//break;
}
s2_real_temp = k * (s2_real * u_real + s2_imag * u_imag) /
u_norm;
s2_imag_temp = k * (s2_imag * u_real - s2_real * u_imag) /
u_norm;
s2_real = s2_real_temp;
s2_imag = s2_imag_temp;
s1_real += s2_real;
s1_imag += s2_imag;
norm=s2_real * s2_real + s2_imag * s2_imag;
if (norm < tol_squared) {
break;
}
//lastnorm=norm;
k++;
}
return so_real * s1_real - so_imag * s1_imag;
}
--
To post to this group, send an email to sage-devel@googlegroups.com
To unsubscribe from this group, send an email to
sage-devel+unsubscr...@googlegroups.com
For more options, visit this group at http://groups.google.com/group/sage-devel
URL: http://www.sagemath.org
