versioninfo()

using BenchmarkTools

N = 2^12

function calc_pi_by_gcd(N)
    s = 0
    for a in 1:N
        for b in 1:N
            s += ifelse(gcd(a,b)==1, 1, 0)
        end
    end
    sqrt(6N^2/s)
end

@show calc_pi_by_gcd(N)
@benchmark calc_pi_by_gcd(N)

# https://github.com/JuliaLang/julia/blob/d386e40c17d43b79fc89d3e579fc04547241787c/base/intfuncs.jl#L28-L49
#
# binary GCD (aka Stein's) algorithm
# about 1.7x (2.1x) faster for random Int64s (Int128s)
@inline function gcd_inline(a::T, b::T) where T<:Union{Int64,UInt64,Int128,UInt128}
    a == 0 && return abs(b)
    b == 0 && return abs(a)
    za = trailing_zeros(a)
    zb = trailing_zeros(b)
    k = min(za, zb)
    u = unsigned(abs(a >> za))
    v = unsigned(abs(b >> zb))
    while u != v
        if u > v
            u, v = v, u
        end
        v -= u
        v >>= trailing_zeros(v)
    end
    r = u << k
    # T(r) would throw InexactError; we want OverflowError instead
    r > typemax(T) && throw(OverflowError())
    r % T
end

function calc_pi_by_gcd_inline(N)
    s = 0
    for a in 1:N
        for b in 1:N
            s += ifelse(gcd_inline(a,b)==1, 1, 0)
        end
    end
    sqrt(6N^2/s)
end

@show calc_pi_by_gcd_inline(N)
@benchmark calc_pi_by_gcd_inline(N)

# Euclidean algorithm
@inline function gcd_euclid(a::T, b::T) where T<:Integer
    while b != 0
        t = b
        b = rem(a, b)
        a = t
    end
    abs(a)
end

function calc_pi_by_gcd_euclid(N)
    s = 0
    for a in 1:N
        for b in 1:N
            s += ifelse(gcd_euclid(a,b)==1, 1, 0)
        end
    end
    sqrt(6N^2/s)
end

@show calc_pi_by_gcd_euclid(N)
@benchmark calc_pi_by_gcd_euclid(N)

@show workers();

function calc_pi_by_gcd_parallel(N)
    s = @parallel (+) for a in 1:N
        s = 0
        for b in 1:N
            s += ifelse(gcd(a,b)==1, 1, 0)
        end
        s
    end
    sqrt(6N^2/s)
end

@show calc_pi_by_gcd_parallel(N)
@benchmark calc_pi_by_gcd_parallel(N)

C_code = raw"""
//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_omp(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);
}

double pi_gcd(long n){
    long count=0;
    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_omp(n));
}
"""

@everywhere file = "calc_pi_by_gcd_gcc_gcc"
filedl = file * "." * Libdl.dlext

open(`gcc -Wall -O3 -march=native -lm -fopenmp -xc -shared -o $filedl -`, "w") do f
     print(f, C_code)
end

run(`ls -l $filedl`)

@everywhere const calc_pi_by_gcd_gcc_lib = file
@everywhere calc_pi_by_gcd_gcc_gcc(n::Int64) = ccall((:pi_gcd, calc_pi_by_gcd_gcc_lib), Float64, (Int64,), n)
@everywhere calc_pi_by_gcd_gcc_gcc_omp(n::Int64) = ccall((:pi_gcd_omp, calc_pi_by_gcd_gcc_lib), Float64, (Int64,), n)
@everywhere gcd_gcc(m::Int64, n::Int64) = ccall((:gcd, calc_pi_by_gcd_gcc_lib), Int64, (Int64, Int64,), m, n)

@show calc_pi_by_gcd_gcc_gcc(N)
@show calc_pi_by_gcd_gcc_gcc_omp(N);

@benchmark calc_pi_by_gcd_gcc_gcc(N)

function calc_pi_by_gcd_gcc_julia(N)
    s = 0
    for a in 1:N
        for b in 1:N
            s += ifelse(gcd_gcc(a,b)==1, 1, 0)
        end
    end
    sqrt(6N^2/s)
end

@show calc_pi_by_gcd_gcc_julia(N)
@benchmark calc_pi_by_gcd_gcc_julia(N)

@benchmark calc_pi_by_gcd_gcc_gcc_omp(N)

function calc_pi_by_gcd_gcc_julia_parallel(N)
    s = @parallel (+) for a in 1:N
        s = 0
        for b in 1:N
            s += ifelse(gcd_gcc(a,b)==1, 1, 0)
        end
        s
    end
    sqrt(6N^2/s)
end

@show calc_pi_by_gcd_gcc_julia_parallel(N)
@benchmark calc_pi_by_gcd_gcc_julia_parallel(N)

# gcc版 (omp parallel を使用**しない**場合)

@btime calc_pi_by_gcd_gcc_gcc(20000)

# julia版

@btime calc_pi_by_gcd(20000)

# Julia版 (gcd を gcc 版の gcd に置き換えた場合)

@btime calc_pi_by_gcd_gcc_julia(20000)

# gcc版 (omp parallel を使用**した**場合)

@btime calc_pi_by_gcd_gcc_gcc_omp(20000)

# julia版 (@parallel を使った場合)

@btime calc_pi_by_gcd_parallel(20000)

# julia版 (gcd を gcc 版に置き換えて、@parallel を使った場合)

@btime calc_pi_by_gcd_gcc_julia_parallel(20000)

# https://github.com/JuliaLang/julia/blob/d386e40c17d43b79fc89d3e579fc04547241787c/base/intfuncs.jl#L28-L49
#
# binary GCD (aka Stein's) algorithm
# about 1.7x (2.1x) faster for random Int64s (Int128s)
@inline function gcd_inline_nounsigned(a::T, b::T) where T<:Union{Int64,UInt64,Int128,UInt128}
    a == 0 && return abs(b)
    b == 0 && return abs(a)
    za = trailing_zeros(a)
    zb = trailing_zeros(b)
    k = min(za, zb)
    # u = unsigned(abs(a >> za))
    # v = unsigned(abs(b >> zb))
    u = abs(a >> za)
    v = abs(b >> zb)
    while u != v
        if u > v
            u, v = v, u
        end
        v -= u
        v >>= trailing_zeros(v)
    end
    r = u << k
    # T(r) would throw InexactError; we want OverflowError instead
    r > typemax(T) && throw(OverflowError())
    r % T
end

function calc_pi_by_gcd_inline_nounsigned(N)
    s = 0
    for a in 1:N
        for b in 1:N
            s += ifelse(gcd_inline_nounsigned(a,b)==1, 1, 0)
        end
    end
    sqrt(6N^2/s)
end

@show calc_pi_by_gcd_inline_nounsigned(N)
@benchmark calc_pi_by_gcd_inline_nounsigned(N)

function calc_pi_by_gcd_inline_nounsigned_UInt64(N)
    s = 0
    for a in UnitRange{UInt64}(1:N)
        for b in UnitRange{UInt64}(1:N)
            s += ifelse(gcd_inline_nounsigned(a,b)==1, 1, 0)
        end
    end
    sqrt(6N^2/s)
end

@show calc_pi_by_gcd_inline_nounsigned_UInt64(N)
@benchmark calc_pi_by_gcd_inline_nounsigned_UInt64(N)

# 純粋 Julia 版 (gcd を Int64 で計算する)

@btime calc_pi_by_gcd_inline_nounsigned(20000)

# オーバーフローエラー (これは意図通りの挙動)

gcd(2^63,2^63)

# これはエラーにならず、負の Int64 の値を返す!

@show gcd(2^63,0);

# すべて Int64 で計算するとこうなる.

@show gcd_inline_nounsigned(2^63, 2^63)
@show gcd_inline_nounsigned(2^63, 0);

# Uint64(2^63) は InexactError() になる.

UInt64(2^63)

# 2^63 % UInt64 ならばエラーにならない.

@show a = 2^63 % UInt64
@show gcd(a, a)
@show gcd(a, 0);

# Int128 の場合

@show b = Int128(2)^127

gcd(b, b)

gcd(b, zero(b))

@show c = b % UInt128
@show gcd(c, c)
@show gcd(c, zero(c));

@inline function mygcd(a::T, b::T) where T<:Union{Int64,UInt64,Int128,UInt128}
    # When T is Int64 or Int128, then abs(typemin(T)) throws InexactError().
    # We want OverflowError instead.
    a == typemin(T) < 0 && throw(OverflowError())
    b == typemin(T) < 0 && throw(OverflowError())
    a, b = abs(a), abs(b)
    a == 0 && return b
    b == 0 && return a
    za = trailing_zeros(a)
    zb = trailing_zeros(b)
    k = min(za, zb)
    u = a >> za
    v = b >> zb
    while u != v
        if u > v
            u, v = v, u
        end
        v -= u
        v >>= trailing_zeros(v)
    end
    u << k
end

mygcd(0,2^63)

mygcd(2^63, 0)

mygcd(2^63+1, 2^63+1)

a = 2^63 % UInt64

mygcd(a, zero(a))

mygcd(zero(a), a)

b = Int128(2)^127

mygcd(b, b)

mygcd(b, zero(0))

c = b % UInt128

mygcd(c, c)

mygcd(c, zero(c))

srand(2018)
a = rand(Int64, 10^6)
b = rand(Int64, 10^6);

bm_gcd = @benchmark gcd.(a,b)

bm_mygcd = @benchmark mygcd.(a,b)

bm_gcd.times[1]/bm_mygcd.times[1]

srand(2018)
c = rand(Int128, 10^6)
d = rand(Int128, 10^6);

bm_gcd128 = @benchmark gcd.(c,d)

bm_mygcd128 = @benchmark mygcd.(c,d)

bm_gcd128.times[1]/bm_mygcd128.times[1]