Факторизация числа при помощи квантового алгоритма Гровера (гостевой пост Романа Душкина)

14 января 2015

Сегодня я хотел бы показать один интересный трюк из области квантовых вычислений, который может иметь многочисленные последствия в применении к теории сложности вычислений. Все мы помним про алгоритм Гровера, который даёт хоть и не сверхполиномиальное ускорение для решения задачи неструктурированного поиска, но всё же является более эффективным по сравнению с классическим алгоритмом поиска грубой силы. Собственно, вдумчивый читатель уже должен был всё понять :)

Эта статья является продолжением в цикле статей по модели квантовых вычислений, поэтому если начать читать её без ознакомления с предыдущими статьями цикла, что-то может показаться непонятным. Так что читателя, который впервые набрёл на эти мои заметки, я отправляю к первым статьям:

В теории сложности вычислений выделен класс задач, которые сложно решить, но легко проверить результат решения (так называемый класс NP). Например, считается, что задача факторизации целого числа находится в этом классе, поскольку до сих пор неизвестно эффективного алгоритма для разложения целого числа на множители, а вот проверить результат довольно легко — достаточно просто перемножить список множителей и сравнить с исходным числом. И таких задач много. Собственно, почему бы не использовать алгоритм Гровера для повышения эффективности решения таких задач? Давайте попробуем…

Подготовка оракула

Как обычно прежде всего реализуем оракул. Но в отличие от алгоритма Шора, когда мы реализовали оракул для факторизации конкретного числа, и он не подходил больше ни для чего, сегодня мы попробуем написать оракул для произвольного числа, которое подходит под критерии Шора для факторизации. Критерии просты — это число должно быть произведением двух и только двух различающихся простых чисел, никакое из которых не равно 2.

Кратко напомню, что представляет собой оракул для алгоритма Гровера. Оракул — это квадратная матрица необходимого размера, на диагонали которой стоят единицы (1) для тех элементов, которые не являются целевыми, и стоит значение -1 для элемента, который необходимо найти. На недиагональных позициях стоят нули (0). Необходимый размер определяется размером индекса того элемента, который необходимо найти при помощи неструктурированного поиска. Количество необходимых кубитов рассчитывается как логарифм по основанию 2 этого индекса (округляется вверх), а размер матрицы равен степени числа 2 от количества кубитов. Другими словами, размер матрицы должен равняться следующей за индексом степенью двойки.

Что значит в контексте данной задачи индекс элемента при неструктурируемом поиске? Можно построить список пар, первым элементом которых будет пара простых чисел, а вторым — их произведение. Прямая задача проста, поскольку для нахождения произведения нам надо просто перемножить два числа в паре. А вот обратная сложна как с точки зрения факторизации (алгоритм Шора), так и с точки зрения неструктурированного поиска (алгоритм Гровера). В общем, сейчас мы будем строить оракул для второго алгоритма в применении к задачи первого.

Для начала составим список всех пар простых чисел, которые больше 2 и не равны друг другу. Чтобы пары не повторялись, первый элемент всегда будет меньше второго. Этот список (бесконечный) построим при помощи следующей функции:

pairsOfPrimes :: Integral a => [(a, a)]
pairsOfPrimes =
[(x,y) | y <- primesWO2, x <- takeWhile (< y) primesWO2]
where
primesWO2 = tail primes

Определение говорит само за себя. Теперь подготовим соответствующий этому список, в котором собраны произведения чисел в парах. В итоге у нас получится два списка, на одинаковых позициях в которых будут находиться пара простых множителей и их произведение соответственно. Объединять их в список пар нет никакой необходимости.

goodNumbers :: Integral a => [a]
goodNumbers = map (uncurry (*)) pairsOfPrimes

Теперь определим функцию для построения оракула. Она должна получать на вход количество кубитов и факторизуемое число, а на выходе возвращать матрицу комплексных чисел (её структура описана выше). Вот определение:

makeOracle :: Integral a => Int -> a -> Matrix (Complex Double)
makeOracle q m = matrixToComplex $ map makeLine [1..2^q]
where
zeroes = replicate (2^q) 0
index = 1 + fromJust (elemIndex m goodNumbers)

makeLine :: Int -> [Int]
makeLine i = changeElement zeroes i ((-1)^(if i == index then 1
else 0))

Здесь локальное определение makeLine меняет элемент в заданном списке по заданной позиции на заданное значение. С таким определением мы уже сталкивались при рассмотрении алгоритма Дойча-Йожи, поэтому определения функции changeElementприводить здесь не будем. Остальное всё в соответствии с приведённым ранее определением оракула.

Поскольку здесь используется функция fromJust, которая небезопасна, то необходимо гарантировать, что ей на вход будет подаваться число, которое заведомо есть в списке хороших чисел goodNumbers. Но именно в этой функции никогда не произойдёт аварийного останова программы из-за некорректного вызова функции fromJust, поскольку список хороших чисел бесконечен, и какое бы число для факторизации мы не передали в функцию makeOracle, она будет работать. Если это будет плохое число, произойдёт перебор бесконечного списка со всеми вытекающими последствиями.

Вроде бы всё. Теперь можно перейти к реализации квантового алгоритма.

Реализация алгоритма

Определим функцию, которая реализует квантовую схему алгоритма Гровера. Поскольку теперь у нас производится поиск в неструктурированном списке произвольной длины, в эту функцию необходимо передавать как количество используемых кубитов, так и количество итераций (всё это должно быть вычислено вне функции на основании числа для факторизации и, возможно. иных факторов). Ну и, само собой разумеется, необходимо передать оракул.

Вот определение:

grover :: Matrix (Complex Double) -> Int -> Int -> IO String
grover f q n = initial |> gateHn q
|> iteration
>>> (measure . fromVector q)
where
initial = toVector $ foldr1 entangle $ replicate q qubitZero
iteration = foldr1 (<**>) $ replicate n (f <**> diffusion q)

Оно практически тождественно тому, что было рассмотрено в соответствующей статье про алгоритм Гровера. При помощи количества кубитов строится начальное квантовое состояние соответствующего размера, равно как и подготавливается преобразование Адамара. Затем производится необходимое количество итераций Гровера (локальное определение iteration). Здесь есть одна тонкость.

Если ранее мы использовали оператор компоновки квантовой схемы (|>), то в этом случае он не подходит. Но на выручку приходит математика — поскольку квантовые вычисления по своей сути представляют собой перемножение матриц, а сама операция умножения матриц ассоциативна, то общую матрицу для всех итераций Гровера можно подготовить при помощи перемножения необходимого количества матриц. Это и делается в замыкании iteration.

Ну а оператор диффузии определён практически так же, как и раньше, за исключением передачи ему количества кубитов для формирования матрицы правильной размерности:

diffusion :: Int -> Matrix (Complex Double)
diffusion q = 2 <*> (qubitPlusQ |><| qubitPlusQ) <-> gateIn q
where
qubitPlusQ = toVector $ foldl1 entangle $ replicate q qubitPlus

В общем-то, всё. Сам алгоритм готов. Необходимо теперь сделать некую обвязку для получения от пользователя числа и разложения его на множители.

Обвязка для взаимодействия с пользователем

Поскольку мы сейчас хотим реализовать что-то большее, чем просто квантовую схему, надо определить функции для организации цикла взаимодействия с пользователем. Что необходимо делать? В целом, всё просто. Во-первых, надо получить от пользователя число, которое он хотел бы факторизовать. Затем надо проверить, число ли вообще ввёл пользователь, а если число, то натуральное ли, а если натуральное, то подходит ли оно под критерии Шора. Далее, если всё нормально, надо запустить алгоритм Гровера, получить результат и вывести его пользователю. А ещё было бы прикольно, если бы программа считала и сообщала сравнение по эффективности (хотя бы теоретической) классического алгоритма поиска и алгоритма Гровера.

Попробуем реализовать такую функцию main:

main :: IO ()
main = do
putStr "Введите число для факторизации: "
mStr <- getLine
unless (map toLower mStr `isPrefixOf` "quit") $
if not (all isDigit mStr) || null mStr
then putStrLn "Надо вводить число, а не что-то иное." >> main
else do
let m = read mStr
if m `mod` 2 == 0
then putStrLn "Введённое число чётно." >> main
else
if isPrime m
then putStrLn "Введённое число простое." >> main
else
if length (primeFactors m) > 2
then do
putStrLn "У введённого числа больше 2-х простых делителей."
main
else
if isSquare m
then do
putStrLn "Введённое число является квадратом."
main
else do
let (cc, cq) = calculateComplexity m
((x, y), count) <- repeatUntil 0
(calculatePrimeFactors m)
(\(x, y) -> x * y == m)
putStrLn ("Решение найдено: " ++ show m ++ " = " ++
show x ++ " * " ++ show y)
putStrLn ("Алгоритм Гровера был вызван " ++ show count ++
" раз.")
showComplexity (cc, cq * count)
putStrLn ""
main

Как видно, здесь функция main бесконечно зациклена, но пользователю оставляется возможность выйти из цикла при помощи слова quit. Далее осуществляется проверка на целочисленность ввода. Если число чётное, то его раскладывать на множители смысла нет. Затем проверяем число на простоту (зачем факторизовать простое число?). Затем проверяется количество простых делителей у числа (должно быть 2 делителя), а после этого и то, не является ли оно квадратом. Если все эти проверки успешно пройдены, то делается следующее.

Сначала вычисляется теоретическая сложность поиска пары простых делителей в списке классическим и квантовым алгоритмами. Это делается при помощи функции calculateComplexity. Затем посредством функции repeatUntil осуществляется поиска простых делителей заданного числа. Эта функция не только возвращает результат, но и указывает, сколько раз пришлось вызвать функцию поиска простых делителей. А критерием успешности выступает задаваемый последним аргументом предикат. После этого на экран выводится результат, пишется количество вызовов алгоритма Гровера, а также указывается сравнительная эффективность классического и квантового алгоритмов. Для квантового алгоритма теоретическая сложность умножается на количество вызовов.

Теперь рассмотрим определения всех функций, задействованных в этом процессе. Вот функция calculateComplexity:

calculateComplexity :: Integral a => a -> (Int, Int)
calculateComplexity m = (2^q `div` 2, nofIterations (2^q))
where
q' = log2 $ toInteger (1 + length (takeWhile (/= m) goodNumbers))
q = if q' == 0 then 1 else q'

Сложность классического алгоритма вычисляется как N/2, где N — это 2q, а q — количество требуемых для вычисления кубитов. Квантовая сложность вычисляется при помощи следующей функции:

nofIterations :: Integer -> Int
nofIterations i = ceiling (pi/4 * sqrt (fromInteger i))

Это просто реализация формулы для оценки количества итераций Гровера, которая дана в описании алгоритма. Вынесена в отдельное определение только потому, что она нам ещё потребуется.

Вот определение монадической функции repeatUntil. Она не относится как таковая к квантовым вычислениям, определена только для организации циклических монадических вычислений, а потому сделана полиморфной для произвольной монады:

repeatUntil :: Monad m => Int -> m a -> (a -> Bool) -> m (a, Int)
repeatUntil counter action condition = do
let counter' = counter + 1
result <- action
if condition result
then return (result, counter')
else repeatUntil counter' action condition

При помощи этой функции мы запускаем процесс поиска пары простых делителей. Вот определение функции calculatePrimeFactors:

calculatePrimeFactors :: Integral a => a -> IO (a, a)
calculatePrimeFactors m = do
let q' = log2 $ toInteger (1 + length (takeWhile (/= m) goodNumbers))
q = if q' == 0 then 1 else q'
oracle = makeOracle q m
result <- grover oracle q $ nofIterations (2^q)
return $ pairsOfPrimes !! binToNumber result

Вот здесь второй раз используется количество итераций, поскольку мы явно указываем, сколько раз надо произвести итерацию Гровера. Функция binToNumber преобразует строку битов в целое число (не будем приводить определение, оно тривиально; равно как и не будет далее приведено определений прочих мелких служебных функций).

Теперь осталось определить функцию для вывода на экран соотношения эффективностей двух алгоритмов. Она выделена отдельно исключительно для того, чтобы не загромождать дальше и так уже загромождённую по самое некуда функцию main. Вот определение:

showComplexity :: (Int, Int) -> IO ()
showComplexity (c, q)
| r < 1.0 =
putStrLn ("Квантовый алгоритм отработал менее эффективно " ++
"классического в " ++
show (ceiling $ recip r) ++ " раз.")
| r == 1.0 =
putStrLn "Эффективность квантового и классического алгоритмов " ++
"одинакова."
| otherwise =
putStrLn ("Превышение эффективности квантового алгоритма над " ++
"классическим: " ++
show (ceiling r) ++ " раз.")
where
r = ((/) `on` (fromInteger . toEnum)) c q

На этом, видимо, всё. Посмотрим, как это дело работает…

GF> main
Введите число для факторизации: 51
Решение найдено: 51 = 3 * 17
Алгоритм Гровера был вызван 1 раз.
Превышение эффективности квантового алгоритма над классическим: 2 раз.

Ну вот как-то так…

Заключение

Предложенный подход в каком-то смысле является универсальным. Алгоритм Гровера позволяет решать обратную задачу с квадратичным ускорением, и в этой статье было показано, как это сделать. Основная проблема заключается в правильном и эффективном построении оракула, а также точном вычислении количества итераций.

Последнее очень важно для алгоритма Гровера, поскольку, как было показано в статье про этот алгоритм, если сделать больше итераций, чем требуется, то результаты встанут с ног на голову. И этот эффект можно наблюдать и при помощи этой программы. Если ввести какое-нибудь большое число (но не слишком большое, а то размеры матриц будут запредельными, и ваш компьютер начнёт поглощать материю из параллельных вселенных), то количество потребных итераций может посчитаться неправильно (из-за округлений), а потому сама квантовая схема будет вызываться большое количество раз, что приведёт в итоге категорическое падение эффективности перед классическим вариантом.

Всякий желающий, как обычно, может получить полный исходный код модуля здесь. Все заинтересованные приглашаются к обсуждению, равно как и к скачиванию моей книги про квантовые вычисления.

Дополнение: Квантовый зоопарк: карта отношений квантовых алгоритмов

Метки: Haskell, Алгоритмы, Функциональное программирование.

Вы можете прислать свой комментарий мне на почту, или воспользоваться комментариями в Telegram-группе.

Записки программиста