How to exactly solve quadratic equations with large integer coefficients (over integers)?

Question

I read a problem about bullseyes in Google Code Jam. (The contest is over now, so it's okay to talk about it)

Maria starts with t millilitres of black paint, which she will use to draw rings of thickness 1cm (one centimetre). A ring of thickness 1cm is the space between two concentric circles whose radii differ by 1cm.

Maria draws the first black ring around a white circle of radius r cm.

The area of a disk with radius 1cm is π cm2. One millilitre of paint is required to cover area π cm2. What is the maximum number of black rings that Maria can draw?

By my calculations on paper, the area of paint to draw a bullseye with n rings, inner radius r, as a multiple of pi is 2*n**2 + n*(2*r-1)

So given t*pi millitres of paint the problem is to find the greatest n such that f(n,r) <= t.

This morning I solved that with binary search https://github.com/hickford/codejam/blob/master/2013/1A/bullseye/bullseye.py

I chose binary search over the quadratic equation because I'm very wary of floating point imprecision — in this problem t and r are integers as big as 10**18). Arithmetic imprecision bit me in a previous Code Jam.

But I'm curious. Can you shore up the quadratic equation to give the correct answer for equations with large integer coefficients? Do maths libraries like Sympy or Numpy have anything to offer me?

---

Demonstration that quadratic equation gives wrong answer for large inputs. For example, with r=308436464205151562 and t=1850618785230909388 . The quadratic equation to solve is

2*n**2 + 616872928410303123*n -1850618785230909388 <= 0

ie. the coefficients are

a = 2
b = 616872928410303123
c = -1850618785230909388

Computing in Python

    > int((-b + math.sqrt(b**2 - 4*a*c)) / (2*a))
    0

This is the wrong answer! The right answer (found by binary search) is 3

>>> n = 3
>>> 2*n**2 + 616872928410303123*n -1850618785230909388 <= 0
True

Answer 1

For symbolic exact manipulations, there's sympy.

If you paste the following:

a, b, c = 2, 616872928410303123, -1850618785230909388
x = Symbol('x')
int(max(solve(a*x**2 + b*x + c, x)))

here, you get 3.

[edited following OP comment].

Answer 2

The roundoff precision killed me on this problem... but you could keep everything at 64 bit integer precision, and do a binary search on the resulting quadratic equation. I outlined my approach here.

Answer 3

If (t / r^2) > 10000 calculate via sqrt.
If (t / r^2) < 10000 try every n starting from 0 increasing by 1.

Answer 4

Solution using integer square roots as suggested by @Vaughn

def solve2(r,t):
    """Maximum number of black rings that Maria can draw, with inner radius r, given pi*t of paint. Solve using quadratic equation"""
    import gmpy
    from gmpy import mpz

    a = 2
    b = 2*r - 1
    c = -t

    x = (-b + mpz(b**2 - 4*a*c).sqrt()) // (2*a) 
    return int(x)