# Matlab sum is wrong: Double symsum gives incorrect result

I'm trying to calculate the double sum

$$\frac{1}{10} \sum_{x=1}^{10} \left( \frac{1}{x} \sum_{n=0}^{floor(log_{10}x)} 10^n \right).$$

In MATLAB, my result is

>> syms n x
>> vpa(symsum(symsum(10^n, n, 0, floor(log(x)/log(10)))/x, x, 1, 10)/10)
ans =
0.29289682539682539682539682539683


This is incorrect. Am I using wrong syntax? Is it a MATLAB bug? Wolfram Alpha gives

sum(sum(10^n, n, 0, floor(log(x)/log(10)))/x, x, 1, 10)/10)
ans = 0.392897


WxMaxima gives

float((sum(sum(10^n, n, 0, floor(log(x)/log(10)))/x, x, 1, 10)/10));
0.3928968253968254


The results from Wolfram Alpha and WxMaxima are correct, or at least they match my hand calculations.

NOTE: I'm using 10 here for the upper limit on the first sum since it's the smallest value for which this discrepancy appears. I'm guessing it has something to do with the upper limit on the inner sum, but when I test $floor(\log_{10}x) = floor(\log(x)/\log(10))$ for different values of $x$ in MATLAB, I get the expected results.

• Are you sure you got the same limits in the sum in Matlab? – mathreadler May 11 '16 at 21:15
• Try using multiple lines so that you can isolate the problem. It is most likely in your code, not in Matlab. (The only issue I can see that would have that size of an impact would be if somehow $\lfloor \log_{10}(x) \rfloor$ was being numerically computed slightly wrong.) – Ian May 11 '16 at 21:17
• Numerically you can't be sure that a floating point value of functions like log behave perfectly because of rounding effects. Using log10 could be better but still not robust. Mathematically log10(10) is 1.0000..., and floor of that is 1. But floor of 0.9999 is 0. – mathreadler May 11 '16 at 21:28
• @mathreadler: I'm not sure I understand your first question. I literally copy and paste the first two lines of code given above into Matlab, and I get the result shown. I too, was wondering, if perhaps it's because of it being computed numerically. – L. Hostetler May 11 '16 at 21:36

You're using symbolic math, but the log(10) in the middle of your code is being evaluating numerically before any symbolic calculation, which leads to imprecision. If you evaluate just your upper bound, you'll see that

floor(log(x)/log(10))


returns floor((1125899906842624*log(x))/2592480341699211) rather than floor(log(x)/log(10)). Order of operations matter when coding for symbolic math and the log function has no way of knowing that the surrounding context is symbolic so it defaults to numeric evaluation.

syms n x
vpa(symsum(symsum(10^n, n, 0, floor(log(x)/log(sym(10))))/x, x, 1, 10)/10)


or just:

syms n x
vpa(symsum(symsum(10^n, n, 0, floor(log10(x)))/x, x, 1, 10)/10)


This can also be solved numerically via:

s = 0;
for x = 1:10
s = s+sum(10.^(0:floor(log10(x))))/x;
end
s = s/10


Finally, WxMaxima and Wolfram Alpha are both computer algebra systems that work symbolically by default. Matlab is numeric by default.

• Okay, that makes sense. It seemed like something that might happen with numerical evaluation, but I mistakenly believed symsum would cause everything to be evaluated symbolically. – L. Hostetler May 11 '16 at 22:46
• Great answer! Impressive you took the time to write all that. +1 – mathreadler May 13 '16 at 13:41

The funny thing about your sum is that it has an exact closed form that is trivial to compute:

Sum[10^n/x, {x, 1, 10}, {n, 0, Floor[Log[10, x]]}]/10


gives the result $9901/25200$ in Mathematica, because for all $1 \le x < 10$, $\lfloor \log_{10} x \rfloor = 0$, thus the inner sum is simply $10^0 = 1$ for all but the last value of $x$, in which case the sum is $10^0 + 10^1 = 11$. Therefore, the given double sum is simply $$\frac{1 + H_{10}}{10},$$ where $H_n = \sum_{k=1}^n 1/k$ is the $n^{\rm th}$ harmonic number.

This is likely a numerical error coming from the fact that $\log_{10}(10) = 1$ and the approximation $$\log_{10}(x) \approx \frac{\log(x)}{\log(10)}$$ Feeding this function $x=10$ would make the slightest numerical floating point rounding error (downwards) would give us a value just slightly lower than 1 ( even if by as little as $10^{-15}$ ) which then the floor function would round down to 0 instead of 1.

• Except that if you actually run floor(log(10)/log(10)) within Matlab you get 1. – Ian May 11 '16 at 21:41
• I don't know if that is actually run internally in this symsum function. Anyway you should definitely not rely on it. Also the difference looks very close to 0.1 which is what the difference would be if that error occured : $\frac{1}{10}\frac{1}{10} 10^1$. – mathreadler May 11 '16 at 21:44