Calculating a large sum of positive and negative integers gives wrong answer

General Usage
,
1

Hello! I am trying to solve a Project Euler problem (#193) to get the squarefree integers less than the large number 2^50. I normally wouldn’t just post the code here directly because of spoilers for the problem except I keep getting the wrong answer in base Julia and in Pluto.jl

using Primes
limit193 = Int64(2)^50-1
sqrt193 = Int64(2)^25-1
function mobius(max_n)
	res = ones(Int8, max_n)
	for p in primes(floor(Int32, sqrt(Int64(max_n))))
		for i = p*p:p*p:max_n
			res[i] = 0
		end
		for i = p:p:max_n
			res[i] *= -1
		end
	end
	res
end
m25 = mobius(Int64(2)^25)
sum(m25[k]*limit193÷k^2 for k = 1:sqrt193)

This last line returns 684489637771331 which is not the correct answer (except in the first few digits, which may or may not be a factor). I am summing 33 million positive and negative integers here so maybe I’m hitting some sort of limitation, but I can’t figure out what. Do you guys have any ideas?

2

this is probably integer overflow. Try doing this with Int128

2 Likes
3

Not sure I understand the task here but mobius is calculating the Mobius function, right? Which is only ever -1, 0, or 1?

The sum is doing an integer division (÷) which will rarely be exact, is this intended?

4

Yes, that is the Mobius function and the integer division is intended.

1 Like
5

Your Möbius function does not give the same result as

term 5801 is 1 in your version and -1 with the other (and from there on a lot of other terms are different)

replacing the Möbius function, the final result is 684465067343069 on my machine (no idea if that is actually the right answer)

2 Likes
6

I hate figuring out integer overflows, so try sum(m25[k]*limit193÷k^2 for k = 1:sqrt193; init=BigInt(0) and see if you’re lucky enough to get a verified answer. I would try this after you figure out the mobius implementation concern JM_Beckers raised.

7

Aw heck. That’s probably it. Is it okay etiquette if I delete this thread since it’s just a bug on my own part but also it’s a spoiler for a Project Euler problem?
EDIT: Yeah it was the mobius function. Thanks guys! :smiley:
EDIT 2: Gonna leave this thread up because there’s already spoilers elsewhere on the Internet if someone wants to find them. Again, thanks

4 Likes
8

Just a few comments:

Are you on a 32-bit system, because on mine, there’s no primes method for Int32:

julia> primes(Int32(11))
ERROR: MethodError: no method matching primes(::Int32)

Closest candidates are:
  primes(::Int64)
   @ Primes C:\Users\DNF\.julia\packages\Primes\Yo1YT\src\Primes.jl:130
  primes(::Int64, ::Int64)
   @ Primes C:\Users\DNF\.julia\packages\Primes\Yo1YT\src\Primes.jl:114

I was wondering why you wrapped your literals like this: Int64(2), since that’s not normally needed.

Secondly, it’s safer to use isqrt(n) rather than floor(Int, sqrt(n)). I struggled a bit to find an example, but here’s one from an old thread using an Int128:

1 Like