Speed test for C and pypy and python to calculate integration using Simpson and Trapezoidal method

integrate.c

#include<stdio.h>
#include <math.h>
#include <time.h>

#define rval (-4.0/9)
float f(float x){
  return sqrt(x)*log(x);
}

int main(){
  int N=100;
  float a=0.0, b=1.0;
  float Sn=0.0, Tn=0.0;
  float errS=0.0, errT=0.0;
  float h;
  clock_t t1,t2;
  t1 = clock();
  /* FILE *fp; */
  /* fp = fopen("simpson.txt","w"); */
  while (N < 100000){		/* Simpson integration*/
    h = (b-a)/N;
    Sn = h*f(b)/6;
    for (int i = 0; i < N; i++){
      if (i == N-1)
	Sn += h/6 * (4*f(a+(i+1.0/2)*h));
      else
	Sn += h/6 * (4*f(a+(i+1.0/2)*h) + 2*f(a+(i+1.0)*h));
    }
    /* errS = fabs(Sn-rval); */
    N += 100;
    /* fprintf (fp,"%d\t%.9f\n",N,errS); */
  }
  /* fclose(fp); */
  /* fp = fopen("trapezoidal.txt","w"); */
  N=1000;
  while (N<100000){		/* Trapezoidal integration */
    h = (b-a)/N;
    Tn = 0.0;
    for (int i = 0; i < N; ++i) {
      if (i!=0)
	Tn += h/2 * (f(a+i*h) + f(a+(i+1)*h));
    }
    /* errT = fabs(Tn-rval); */
    N += 100;
    /* fprintf(fp,"%d\t%.9f\n",N,errT); */
  }
  t2 = clock();
  printf ("%.9f\n",(double)(t2-t1)/CLOCKS_PER_SEC);
  return 0;
}

integrate.py

#!/usr/bin/env pypy
from math import sqrt,log
import time
# pypy doesn't support matplotlib
## import matplotlib.pyplot as plt 
a, b = 0.0, 1.0
f = lambda (x): sqrt(x) * log(x)
real = -4.0/9
# x->0 , f(x) = 0

n = range(100,100000,100)
errS, errT = [],[]
t = time.time()
for N in n:
    h = (b-a)/N
    # Sn
    Sn = h*f(b)/6
    for k in range(0,N,1):
        if k==N-1:
            Sn += h/6 * (4*f(a+(k+1.0/2)*h))
        else:
            Sn += h/6 * (4*f(a+(k+1.0/2)*h) + 2*f(a+(k+1)*h))
    # Tn
    Tn = 0.0
    for k in range(0,N,1):
        if k!=0:
            Tn += h/2 * (f(a+k*h) + f(a+(k+1)*h))
    errS.append(abs(Sn-real))
    errT.append(abs(Tn-real))
    # print "Sn = %f  error(Sn) = %f" %(Sn,errS)
    # print "Tn = %f  error(Tn) = %f" %(Tn,errT)
print time.time() - t
## plt.plot(n,errS,label='Error Sn')
## plt.loglog(n,errT,label='Error Tn')
## plt.legend(loc='upper right')
## plt.title('Relative error with different n')
## plt.show()