By translating the following program written in C ++
#include <iostream> #include <math.h> using namespace std; inline int gcd (int a, int b) { if(a==0) return b; return gcd(b%a,a); } double sol(const int n) { const double p=n/(double)(n+1); double ans=0; int all=(n+1)*(n+1); for(int i=0;i<=n;++i) for(int j=0;j<=n;++j) for(int a=0;a<=n;++a) for(int b=0;b<=n;++b) if( (i!=a || j!=b ) && (i*b!=a*j)) { int u=a-i,v=b-j,cnt=0; for(int x=0;x<=n;++x) { int t=(x-i)*v+j*u; // y * u > t if(u==0) cnt+=t<0?n+1:0; else if(u>0) cnt+=t<0?n+1:n-min(t/u,n); else cnt+= t<0?min((t+1)/u,n) +1:0; } if(cnt!=0) ans+=(i*b-a*j)*pow(p,all-abs(gcd(u,v))-cnt-1)*(1-pow(p,cnt)); } return ans/(2*all); } int main() { cout << sol(1) << "\n" ; cout << sol(5) << "\n" ; cout << sol(10) << "\n" ; } In Mathematica like this (absolutely the same instructions)
f[n_] := Module[ {p = n/(n + 1), ans = 0., all = (n + 1) (n + 1), i, j, a, b, x, u, v, cnt, t}, For[i = 0, i <= n, ++i, For[j = 0, j <= n, ++j, For[a = 0, a <= n, ++a, For[b = 0, b <= n, ++b, If[(i != a || j != b) && (i b != a j), u = a - i; v = b - j; cnt = 0; For[x = 0, x <= n, ++x, t = (x - i) v + j u; Which[ u == 0, cnt += If[t < 0, n + 1, 0], u > 0, cnt += If[t < 0, n + 1, n - Min[t/u, n]], True, cnt += If[t < 0, Min[(t + 1)/u, n] + 1, 0] ]; ]; If[cnt != 0, ans += (i b - a j)* Power[p, all - Abs[GCD[u, v]] - cnt - 1]*(1 - Power[p, cnt])] ]]]]]; Return[ans/(2 all) // N] ] I obtain the same result with n = 1 but different results whith n != 1.
Somebody can explain where I am wrong ?
Thanks for help
