Последовательный и параллельный (CUDA) поиски диофант. ур-ий 2 степени

Диофантово уравнение 2-го порядка с двумя неизвестными имеет вид a*x^2 + b*y^2 = c, 3 коэффициента и 2 переменные - все они неизвестны, но лежат в диапазоне [-1000000;1000000 ] и принадлежат множеству целых чисел. Нужно написать 2 скрипта - последовательный и параллельный (CUDA) вариант нахождения всех связок {a,b,c,x,y}, при которых уравнение существует.
По идеям есть расширенный алгоритм Евклида и поиск в лоб.
В лоб на псевдокоде представляется так:
for a in [-10^6; 10^6]:
for b in [-10^6; 10^6]:
for c in [-10^6; 10^6]:
for x in [-10^6; 10^6]:
y = |sqrt(c-a*x^2)/b|
if y - целое, выводим {a, b, c, x, y} и сохраняем в файл.
Перебор всех вариантов составит (2*10^6)^4. Как распараллелить - вопрос открытый.

Цена договорная.