2

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

QuantumDot
  • 19,601
  • 7
  • 45
  • 121
nomis
  • 21
  • 1

1 Answers1

13

In C, a/b for integers a and b does something different than in Mathematica.

To translate the code faithfully, we could use IntegerPart[t/u] and IntegerPart[(t+1)/u], or as Szabolcs points out, Quotient[t, u] and Quotient[t+1, u].

With this modification, the results are

f[1]
f[5]
f[10]

(* 0.1875
   9.76895
   55.0301 *)

which is in line with what is printed by the C code.

Community
  • 1
ilian
  • 25,474
  • 4
  • 117
  • 186
  • Or also Quotient. – Szabolcs Sep 05 '15 at 07:38
  • What do you think about this? Bug or not? What bothers me is now how JSON behaves but the inconsistency between JSON and RawJSON. – Szabolcs Sep 05 '15 at 08:30
  • @Szabolcs Exactly, UTF-8 encoded output is not a bug, but I don't know why the newer & separate implementation (Developer\WriteJSONString`) chooses to not do it. – ilian Sep 05 '15 at 16:12
  • ilian, could I please have your attention on this issue to determine if it is intentional or (as I presume) a bug?: (92578) – Mr.Wizard Sep 05 '15 at 18:35
  • @Mr.Wizard It does seem like a bug to me, not that I'm any kind of expert on the frontend's mysterious ways. I also couldn't observe it on my Windows machine so maybe I have some lucky combination of video card, monitor and color profile... I'd recommend sending this issue to support with your SystemInformation[] so they can file a bug report. – ilian Sep 05 '15 at 20:39