Ask Your Question
0

contour_plot issue with large numbers

asked 2013-09-13 08:25:34 -0500

sappidus gravatar image

updated 2013-09-13 10:17:01 -0500

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.

edit retag flag offensive close merge delete

1 answer

Sort by ยป oldest newest most voted
1

answered 2013-09-13 14:09:24 -0500

sappidus gravatar image

Ah, I did not realize plotting is by default done with floats. By setting:

f(x,y)=mod2((N+floor(y/4)) / 2^(4*floor(x)+mod4(y)))
g=fast_callable(f, vars=[x,y], domain=RealField(100))

and contour_plot'ing g instead, which forces the use of higher-precision arithmetic, I get what I want.

edit flag offensive delete link more

Comments

That's a very good observation. Your solution is also a potential solution to ticket 13559. There is also ticket 15030 that aims to switch the plot functions to use `fast_callable` instead of `fast_float`.

ppurka gravatar imageppurka ( 2013-09-14 05:01:39 -0500 )edit

Your Answer

Please start posting anonymously - your entry will be published after you log in or create a new account.

Add Answer

Question Tools

Stats

Asked: 2013-09-13 08:25:34 -0500

Seen: 184 times

Last updated: Sep 13 '13