3 * The Great Computer Language Shootout
4 * http://shootout.alioth.debian.org/
6 * Contributed by Sebastien Loisel
13 double eval_A(int i, int j) { return 1.0/((i+j)*(i+j+1)/2+i+1); }
15 void eval_A_times_u(int N, const double u[], double Au[])
21 for(j=0;j<N;j++) Au[i]+=eval_A(i,j)*u[j];
25 void eval_At_times_u(int N, const double u[], double Au[])
31 Au[i] += eval_A(j,i)*u[j];
35 void eval_AtA_times_u(int N, const double u[], double AtAu[])
38 eval_A_times_u(N,u,v);
39 eval_At_times_u(N,v,AtAu);
42 int main(int argc, char *argv[])
45 int N = ((argc == 2) ? atoi(argv[1]) : 2000);
47 double *u = alloca(sizeof(u[0]) * N);
48 double *v = alloca(sizeof(v[0]) * N);
55 eval_AtA_times_u(N,u,v);
56 eval_AtA_times_u(N,v,u);
66 printf("%0.9f\n", sqrt(vBv/vv));