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

In a previous code jam, I got bitten by a floating point precision error http://stackoverflow.com/q/15978781/284795 — Colonel Panic, Apr 27 '13 at 13:35
You shouldn't have any problem using math.sqrt on an integer that small. — Vaughn Cato, Apr 27 '13 at 13:43
My mistake, for some reason I was reading the maximum value as 2**18 instead of 10**18. — Vaughn Cato, Apr 27 '13 at 14:07
Integer square roots aren't helpful for solving the quadratic equation, because we need to add `-b/2a`, a fraction, to the square root *before* any rounding. — Colonel Panic, Apr 27 '13 at 14:24
Once you have the integer sqrt, you subtract the integer b to get the integer x. Then you just need to determine `x/4` rounded down to an integer. `int(x/4) = int(int(x)/4)`. — Vaughn Cato, Apr 27 '13 at 14:37

Haile · Accepted Answer · 2013-04-27T14:21:01.310

2

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].

edited Apr 27 '13 at 14:21

answered Apr 27 '13 at 14:15

Haile

3,120
3
24
40

More like 3.0? `int(max(solve(a*x**2 + b*x + c, x)))` returns 3 which is the right answer. Thanks! Neat REPL – Colonel Panic Apr 27 '13 at 14:18

score 1 · Answer 2 · answered Apr 27 '13 at 19:28

1

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.

answered Apr 27 '13 at 19:28

tbischel

6,337
11
51
73

score 0 · Answer 3 · answered Apr 27 '13 at 14:22

0

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

answered Apr 27 '13 at 14:22

Egor Skriptunoff

23,359
2
34
64

score 0 · Answer 4 · answered Apr 27 '13 at 17:15

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)

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

4 Answers4