An integral with a fractional part in 3 dimensions

The analytic answer is $$ %\sum_{n=1}^\infty \frac{1}{2n(n+1)^2}\sum_{n=1}^\infty \frac{1+3n}{6n^2(n+1)^3}+\sum_{n=1}^\infty \frac{3n^2-1}{6n^2(n+1)^3}=\\ 1+\pi ^2\frac{2 \zeta (3)-9}{72} \approx 0.0958502 $$

Therefore, Mathematica is correct.

Proof

The 3D integral

NIntegrate[FractionalPart[x/y] FractionalPart[y/z] FractionalPart[z/x], 
   {x, 0, 1}, {y, 0, 1}, {z, 0, 1}]

0.0958741

depends only on ratios between x, y, and z. Therefore, it can be reduced to the 2D integral

NIntegrate[FractionalPart[x/y] y FractionalPart[1/x], {x, 0, 1}, {y, 0, 1}]

0.0958383

It has the following structure

ArrayPlot@Table[FractionalPart[x/y] y FractionalPart[1/x], 
    {x, 0.001, 1, 0.001}, {y, 0.001, 1, 0.001}]

We can split the integral to two parts

NIntegrate[FractionalPart[x/y] y FractionalPart[1/x], {x, 0, 1}, {y, 0, x}]
NIntegrate[FractionalPart[x/y] y FractionalPart[1/x], {x, 0, 1}, {y, x, 1}]

0.0176263

0.0782212

After substitution $y/x = ξ$ the first one becomes

NIntegrate[x FractionalPart[1/x], {x, 0, 1}] NIntegrate[x^2 FractionalPart[1/x], {x, 0, 1}]

0.0176233

The second one is equivalent to

NIntegrate[x (1 - x) FractionalPart[1/x], {x, 0, 1}]

0.0781823

These integrals can be written as a sum of

Integrate[x (1/x - n), {x, 1/(n + 1), 1/n}] // Simplify
Integrate[x^2 (1/x - n), {x, 1/(n + 1), 1/n}] // Simplify
Integrate[x (1 - x) (1/x - n), {x, 1/(n + 1), 1/n}] // Simplify

1/(2 n (1 + n)^2)
(1 + 3 n)/(6 n^2 (1 + n)^3)
(-1 + 3 n^2)/(6 n^2 (1 + n)^3)

Sum[1/(2 n (1 + n)^2), {n, 1, ∞}]
Sum[(1 + 3 n)/(6 n^2 (1 + n)^3), {n, 1, ∞}]
Sum[(3 n^2 - 1)/(6 n^2 (1 + n)^3), {n, 1, ∞}]

1/2 (2 - π^2/6)
1/6 (3 - 2 Zeta[3])
1/12 (6 - π^2 + 4 Zeta[3])

Finally, we obtain

1/2 (2 - π^2/6) 1/6 (3 - 2 Zeta[3]) + 1/12 (6 - π^2 + 4 Zeta[3]) // Simplify
N[%]

1 + 1/72 π^2 (-9 + 2 Zeta[3])
0.0958502

The Maple image you showed in your comment above was done in Maple 2D math (document mode). little hard to replicate it using classical interface (worksheet mode) which I use since one does not see the actual Maple commands that way.

The Maple commands you ave also might be using another package that was loaded before and not shown, such as Maple's Student. It is best that you post complete self contained Maple code to be able to compare.

But I have Maple 17, and not able to get the same result. Maple is not able to do this integral either.

restart:
f:=x->x-floor(x);
Int(f(x/y)*f(y/z)*f(z/x),[x=0..1,y=0..1,z=0..1]);
value(%);

Error, (in floor) too many levels of recursion

There is one thing I noticed though, is that Maple has a rule to integrate floor function, and Mathematica does not do it: Maple's definition is same as M:

floor - greatest integer less than or equal to a number

hence

restart:
int(floor(x),x);

    (*  floor(x)*x  *)

and in Mathematica, it remains unevaluated

Integrate[Floor[x], x]

But do you want to see something bizarre? I gave Integrate a rule integrate floor(x/y) to x*floor(x/y), which is also how it is defined in Maple:

 int(floor(x/y),x);
 (*  floor(x/y)*x *)

Unprotect Integrate and add the rule:

 (Unprotect[Integrate]; Integrate[Floor[x_/y_],x_] := x Floor[x/y];Protect[Integrate];)
 Integrate[x-Floor[x/y] , {x, 0, 1}]

and this is the result

Where did the 99999999999999999999/(100000000000000000000 y) come from? This should be 1

 N[%]

Very strange result. Here is Maple's result:

 int(x-floor(x/y),x=0..1);
 (* -floor(1/y)+1/2 *)

So same final result, but M exact answer has this 99999999999999999999/(100000000000000000000 y) term which only becomes 1.0 with using N on it. This might be related to the WorkingPrecision messages you are seeing.

The bottom line is: Maple does not give 1 as you showed in your screen image.