C = With MinGW =

**~52 sec.**

C# 3.5 SP1 =

**~36 sec.**

Why is C# faster ?

I use the Mersenne Twister to generate the random number. For the C implementation I use the new Fast Mersenne Twister, which make extensive use of the preprocessor to generate efficient code. With C# I use an implementation written in 2003. So the Mersenne Twister that I use with C# should be slower. In fact the simulation in C# is even faster with their standard RNG, but I use the Mersenne Twister to get a fair comparison.

I'm surprised by this result, I thought C# would be as slow as Java for Monte Carlo simulations, perhaps a little better as the language can be optimized for windows, but I never expected C# to beat C. Of course it's just a single program but still I would really like to understand how a relatively high level language such as C# can be faster than C, which is directly compiled to machine code.

Here's the code for the two simulations;

using System; using CenterSpace.Free; // Mersenne Twister (from c-sharpcorner dot com) namespace WrightFisher { class WF { static void Main(string[] args) { MersenneTwister generator = new MersenneTwister(); int x = 1000; // Number of simulations int t = 1000; // Number of generations/simulations int Ne = 1000; // Effective population size int i0 = 500; // Initial number of allele 'A1' int i; // Used for trials int i1; // Alleles at the end of a generation int[] f = new int[x]; // Vector for the number of allele 'A1' at the end of each simulation // Features of A1; double s = 0.00; // Selection coefficient for 'a' // Store the probabilities; double[] p = new double[Ne + 1]; for (int j = 0; j < Ne + 1; j++) { p[j] = ((1.0 + s) * j) / ((1.0 + s) * j + Ne - j); } // The simulation; for (int sim = 0; sim < x; sim++) { i1 = i0; for (int gen = 1; gen < t + 1; gen++) { i = 0; for (int allele = 0; allele < Ne; allele++) if (generator.NextDouble() < p[i1]) ++i; // NextDouble() uses genrand_real2() i1 = i; } f[sim] = i1; } int ext = 0, fix = 0, pol = 0; for (int j = 0; j < x; j++) { if (f[j] == 0) ++ext; else if (f[j] == Ne) ++fix; else ++pol; } Console.WriteLine("Fixation = " + (double)fix / x); Console.WriteLine("Polymorphism = " + (double)pol / x); Console.WriteLine("Extinction = " + (double)ext / x); } } }

#include <stdio.h> #include <stdlib.h> #include <time.h> // For the initialization of the pseudorandom number generator #define MEXP 19937 #include "RNG/SFMT.c" int main(void) { init_gen_rand((unsigned)time(NULL)); // Initialization of the RNG int x = 1000; // Number of simulations int t = 1000; // Number of generations/simulations int Ne = 1000; // Effective population size int i0 = 500; // Initial number of allele 'A1' int i; // Used for trials int i1; // Alleles at the end of a generation int f[x]; // Vector for the number of allele 'A1' at the end of each simulation // Features of A1; double s = 0.0; // Selection coefficient for 'a' double p[Ne+1]; // Store the probability int j; for(j = 0; j < Ne+1; j++) p[j] = ((1.+s)*j)/((1.+s)*j+Ne-j); int sim, gen, allele; for(sim = 0; sim < x; sim++) { i1 = i0; for(gen = 1; gen < t+1; gen++) { i = 0; for(allele = 0; allele < Ne; allele++) { if(genrand_real2() < p[i1]) ++i; // i += (genrand_real2() < p[i1]) // Faster ??? } i1 = i; } f[sim] = i1; } int ext = 0, fix = 0, pol = 0; for(j = 0; j < x; j++) { if(f[j]==0) ++ext; else if(f[j]==Ne) ++fix; else ++pol; } printf("Fixation = %.4f\n", (double)fix/x); printf("Polymorphism = %.4f\n", (double)pol/x); printf("Extinction = %.4f\n", (double)ext/x); return 0; }