calc pi from gcd

pi_gcd.c

//gcc pi_gcd.c -O3 -march=native -lm -fopenmp
#include <stdio.h>
#include <stdlib.h>
#include <math.h>

inline int ctz(long x){
    // GCCの組み込み関数
    return __builtin_ctzl(x);
}

inline int min(int a, int b){
    return a<b ? a : b;
}

long gcd(long a, long b){
    if(a==0) return b;
    if(b==0) return a;
    a = labs(a);
    b = labs(b);
    int ta = ctz(a);
    int tb = ctz(b);
    int shift = min(ta,tb);
    a >>= ta;
    b >>= tb;
    if(a<b){
        long t=a;
        a=b;
        b=t;
    }
    do{
        a-=b;
        a >>= ctz(a);
        if(a<b){
            long t=a;
            a=b;
            b=t;
        }
    }while(b!=0);
    return a << shift;
}

double pi_gcd(long n){
    long count=0;
#pragma omp parallel for reduction(+:count) collapse(2)
    for(long i=1; i<=n; i++){
        for(long j=1; j<=n; j++){
            if(gcd(i,j)==1){
                count++;
            }
        }
    }
    return sqrt(6.0*n*n/count);
}

int main(int argc, const char* argv[]){
    long n=20000;
    if(argc>=2){
        n=atol(argv[1]);
    }
    printf("%.15f\n", pi_gcd(n));
}