How to handle float error for plots near discontinuities in Sage

Last week I read this blog post by +Patrick Honner. In the post +Patrick Honner plots a graph of a function with a removable discontinuity on Desmos and when zooming in enough he got some errors.

I was waiting around to start this (ridiculously fun) Hangout on Air with a bunch of mathematicians hosted by +Amy Robinson of +Science on Google+:

While waiting I rushed to write this blog post claiming that if you did the same thing with +Sage Mathematical Software System you did not get any errors. It was quickly pointed out to me on twitter and in the comments that I just had not zoomed in enough.

I edited the blog post to first of all change the title (it was originally 'When Sage doesn't fail' but now reads 'When Sage also fails') and also to include some code that shows that the exact same error appears.

On G+, +Robert Jacobson (who's the owner of the Mathematics community which you should check out if you haven't already) pointed out that you could surely use Sage's exact number fields to avoid this error.

He put together some code and shared it with me on +The Sagemath Cloud that does exactly this. Here's a slight tweak of the code Robert wrote (hopefully you haven't changed your mind and still don't mind if I blog this Robert!):

f(x) = (x + 2) / (x ^ 2 + 3 * x + 2) # Define the function
discontinuity = -1  # The above function has two discontinuities, this one I don't want to plot
hole = -2  # The hole described by Patrick Honner

def make_list_for_plot(f, use_floats=False, zoom_level=10^7, points=1001):
    count = 0  # Adding this to count how many tries fail
    z = zoom_level
    xmin = hole - 10/z # Setting lower bound for plot
    xmax = min(hole + 10/z, discontinuity - 1/10) # Setting upper bound for plot only up until the second (messy) discontinuity
    x_vals = srange(start=xmin, end=xmax, step=(xmax-xmin)/(points-1), universe=QQ, check=True, include_endpoint=True)

    # If we are using floating point arithmetic, cast all QQ numbers to floating point numbers using the n() function.
    if use_floats:
        x_vals = map(n, x_vals)

    lst = []
    for x in x_vals:
        if x != hole and x != discontinuity:  # Robert originally had a try/except statement here to pick up ANY discontinuities. This is not as good but I thought was a bit fairer...
            y = f(x)
            lst.append((x, y))

    return lst

The code above makes sure we stay away from the discontinuity but also allows us to swap over to floating point arithmetic to see the effect. The following plots the functions using exact arithmetic:

exact_arithmetic = make_list_for_plot(f)

p = list_plot(exact_arithmetic, plotjoined=True)  # Plot f
p += point([hole, -1], color='red', size=30)  # Add a point
show(p)

We see the plot here (with no errors):

To call the plots with floating point arithmetic:

float_arithmetic = make_list_for_plot(f, use_floats=True)

p = list_plot(float_arithmetic, plotjoined=True)  # Plot f
p += point([hole, -1], color='red', size=30)  # Add a point
show(p)

We see that we now get the numerical error:

Just to confirm here is the same two plots with an even higher zoom:

To change the zoom, try out the code in the sage cell linked here: simply change the zoom_level which was set to $10^12$ for the last two plots.

(Going any higher than $10^14$ seems to bring in another error that does not get picked up by my if statement in my function definition: Robert originally had a try except method but I thought that in a way this was a 'fairer' way of doing things. Ultimately though it's very possible and easy to get an error-less plot.)