contour_plot issue with large numbers
I was trying to whip up a variant on Tupper's "self-referential" formula graph (nice exposition here if you're not familiar with it) in Sage, but I ran into a graphing issue with largish numbers. In an attempt to nail down the problem, I've tried to reduce this to a minimal case:
x,y = var('x,y')
def mod2(x): return floor(x)-2*floor(x/2)
def mod4(x): return floor(x)-4*floor(x/4)
N=2^53+1
contour_plot(mod2((N+floor(y/4)) / 2^(4*floor(x)+mod4(y))),(x,0,14),(y,0,4))
(I wrote the mod2/4 functions because I kept on running into issues with the % operator and type coercion.)
If you understand the idea behind the Tupper formula, then you know that the binary representation of N is basically a readout of the "pixels" of the eventual graph generated. Now, 2^53+1 is 1000...001 in binary, so there should be TWO "pixels" lit.
But, when I actually execute the above code on my home computer or in the Sage Cloud, something funny happens. I get a graph where only ONE "pixel" is lit -- the one corresponding to the least significant bit of 2^53+1, which should be in the bottom left of the graph (i.e., near the origin) is mysteriously absent.
However, if you check for...
N=2^52+1
...the graph correctly shows two "pixels". So I first presumed this was some sort of floating point or rounding error. But if you actually check the values returned by the function being graphed, the returned values are always correct. For example, setting (x,y) to a value squarely within the "missing pixel" and manually evaluating...
x,y=0.5,0.5
mod2((N+floor(y/4)) / 2^(4*floor(x)+mod4(y)))
...yields the correct value of 1, even for N=2^53+1! Similarly, setting (x,y) to (0.5,1.5), i.e., outside said "pixel", appropriately yields 0. So I don't get at all why the graphing isn't working.
Any insight would be appreciated.
