Знакомство с p-адическими числами. Часть 2, практическая24.01.2022 21:47

Эта часть продолжает неформальный рассказ о p-адических числах и она посвящена практическим аспектам работы с этой числовой системой и, в частности, некоторым деталям реализации p-адической алгебры на языке Haskell. О том, что это за система и зачем она может понадобиться, читайте в предыдущей части.

Мы поговорим об эффективном внутреннем представлении p-адических чисел, о базовых алгоритмах и методах работы с ними, а также о двух классных инструментах в системе типов языка Haskell: о типах-литералах (type literals) и семействах типов (type families).

Как готовить p-адические числа. Рецепты и хитрости.

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

Мы уже достаточно подкованы теоретически, чтобы понять, что p-адические числа — это, всё-таки, не цепочка цифр, а числовая система. Эффективнее всего компьютер работает с целыми числами, и в роли системы, представляющей кольцо целых p-адических чисел с ограниченной точностью в $k$ разрядов, может выступать кольцо $\mathbb{Z}/p^k\mathbb{Z}$ , содержащее обычные целые числа с обычной модулярной арифметикой. Например, вместо «честных» 2-адических чисел можно использовать 64-разрядное число типа Word64, с которым можно обращаться, как с представителем кольца $\mathbb{Z}/2^{64}\mathbb{Z}$ . Также 7-адические числа можно представить как числа целые числа, вычисляемые по модулю $7^{22}$ , которые тоже все помещаются в Word64.

По существу, понижая $\mathbb{Z}_p$ до $\mathbb{Z}/p^k\mathbb{Z}$ , мы ограничиваемся лишь первой цифрой канонического разложения в кольце $\mathbb{Z}_{p^k}$ , в котором очень много элементов и они могут быть очень большими. Существование изоморфизмов между числовыми системами, упомянутых в первой части статьи, приводит к тому, что сложение, вычитание, умножение и даже деление для этой первой цифры будет в полной мере соответствовать p-адическим операциям. А поскольку это просто большое целое число, то арифметика с ним может быть реализована очень эффективно, например, с использованием GMP или подобных ей библиотек. Таким образом, на практике стоит рассматривать p-адические числа (и рациональные и трансцендентные) как большие целые, вся арифметика с которыми выполняется по модулю $p^k$ .

Ограничение разрядности требует осознанной работы с точностью чисел. И нам нужно научиться переходить от вещественных чисел, выраженных с конечной точностью, в $\mathbb{Z}/p^k\mathbb{Z}$ и возвращаться обратно так, чтобы по-возможности, не потерять необходимой информации при этом переходе.

В нашем контексте точность вещественных чисел определяется не числом значащих цифр десятичного разложения, а максимальными абсолютными значениями для числителя и знаменателя дроби, приближающей вещественные числа. Например, для каких-то вычислений мы можем ограничиться дробями, в которых числитель и знаменатель не превышают $1000$ , а для более точных расчетов повысить верхнюю границу используемых числителей и знаменателей, скажем, до максимального значения допускаемого типом Word64 — огромного числа $2^{64} \approx 10^{20}$ .

Приведу два базовых алгоритма для перевода чисел из $\mathbb{Q}$ в $\mathbb{Q}_p$ и наоборот.

Отображение рационального числа в p-адическое

Для перевода числа $x = \frac{r}{s}$ из множества $\mathbb{Q}$ в $\mathbb{Z}/p^k\mathbb{Z}$ достаточно найти обратный элемент для $s$ в $\mathbb{Z}/p^k\mathbb{Z}$ , то есть элемент $q = s^{-1}$ , такой, что $q\cdot s \equiv 1 \mod p^k$ . Это сравнение эффективно решается расширенным алгоритмом Евклида или через цепные дроби, еcли $\gcd(s, p) = 1$ .

Вот как выглядит алгоритм для нахождения округлённого p-адического образа числа $x \in Q$ , как элемента $\mathbb{Z}/p^k\mathbb{Z}$ :

Выделяем из числителя или знаменателя множитель , то есть, представляем входное число как произведение

$x = \frac{r}{s}\cdot p^v.$

Если $\gcd(s, p) = c \neq 1$ , то нужно избавиться от делителя $p$ в знаменателе:

$x = \frac{p^{v'}}{p^{v'}} \frac{c^{v'}}{c^{v'}}\frac{r}{s} \cdot p^{v} = \frac{r\cdot (p/c)^{v'}}{s / c^{v'}}\cdot p^{(v-v')}.$

Полученная дробь будет обратимым элементом в $\mathbb{Z}_p$ . Для знаменателя обратимой части находим обратный элемент $q = s^{-1} \mod p^k$ расширенным алгоритмом Евклида или библиотечной функцией, наподобие mpz_invert из библиотеки GMP.
Результатом будет p-адическое число с обратимым элементом $u = r\cdot q$ и порядком $v$ .

Пример для простого основания

Пример для составного основания

Восстановление рационального числа из p-адического

Как мы уже говорили, из-за того, что множество p-адических чисел существенно больше множества рациональных, корректная реконструкция рационального числа возможна не всегда. При этом особенности p-адической топологии приводят к тому, что постепенное увеличение точности для нерационального p-адического числа не порождает последовательности приближений, сходящейся к какому-то рациональному числу, а работает скорее, как генератор псевдослучайных рациональных чисел. Однако, как в 1981 году показал Пол Ванг, если рациональная реконструкция существует, то она может быть однозначно получена алгоритмом, основанным на расширенном алгоритме Евклида, более известном, как метод Лагранжа для выделения базиса целочисленной решётки.

Вот как выглядит алгоритм восстановления рационального числа $x = \frac{r}{s}, r < N, w < D$ из p-адического числа, округлённого до $k$ разрядов, представленного числом $n \in \mathbb{Z}/p^k\mathbb{Z}$ . Необходимое (но не достаточное) условие корректной работы алгоритма: $p^k > 2ND$» />.</p>

<p><em>Инициализируем два вектора <img src= . Пока $w_1 \geq N$ выполняем итерации

$w \to v,\quad v - \left\lfloor\frac{v_1}{w_1}\right\rfloor w \to w.$

Если по окончании цикла получаем $w_2 < D$ и $\gcd(w_1, w_2) = 1$ , то результатом будет число $x = w_1/w_2$ . В противном случае восстановление невозможно, и необходимо увеличивать число разрядов $k$ .

Пример восстановления

Лемма Гензеля

Один из основных инструментов при работе в p-адическом мире был описан самим основателем этой теории Куртом Гензелем. Процесс построения проективного предела из начального кольца вычетов, о котором шла речь в первой части статьи, можно применять ко всей числовой системе сразу, а можно сосредоточиться только на некоторой её части: на множествах корней алгебраических уравнений. Он позволяет для уравнения $f(x) = 0$ , которое очень легко решить в самом маленьком кольце вычетов $\mathbb{Z}/p\mathbb{Z}$ , повысить эти элементарные решения до любого кольца $\mathbb{Z}/p^k\mathbb{Z}$ , в котором такие решения существуют.

Алгоритм поднятия решений алгебраического уравнения $f(x) = 0$ с целыми (или p-адическими целыми) коэффициентами с помощью леммы Гензеля состоит в следующем:

Находим решение в кольце $\mathbb{Z}/p\mathbb{Z}$ , например, перебором, если кольцо имеет небольшой порядок. На этом этапе надо убедиться в том, что в этом кольце $f'(x_1) \neq 0$ . Для $0 < m \leq k$ , имея решение в $\mathbb{Z}/p^k\mathbb{Z}$ находим решение $x_{k+m}$ в $\mathbb{Z}/p^{k+m}\mathbb{Z}$ , следующим образом: $x_{k+m} = x_k - a\cdot f(x_k)$ , где $a = f'(x_k) \mod p^m$ .

Пример поднятия решения уравнения

По существу, это метод Ньютона, хорошо известный в вещественном мире. На практике, если p-адические числа реализованы, как кольцо $\mathbb{Z}/p^k\mathbb{Z}$ , проще вычислять корни уравнений напрямую методом Ньютона в этом большом кольце, без возни с промежуточными кольцами. Главное при этом, начинать итерации с корней, существующих в базовом кольце $\mathbb{Z}/p\mathbb{Z}$ .

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

Вычисление квадратных корней

Коль скоро в нашем арсенале появился метод Ньютона, мы можем вычислять p-адические квадратные корни (если они существуют), при этом шаг итераций не будет отличаться от традиционного $x \to (x + a/x)/2$ . Однако, он точно не сможет работать в $\mathbb{Z}_2$ , поскольку производная целевой функции $f(x) = x ^2 - a$ тождественно равна нулю на всём $\mathbb{Z}/2\mathbb{Z}$ . Для этого случая можно использовать оригинальный рекурсивный метод:

$\sqrt{a} = g\left(\frac{a-1}{8}\right),\\ g(x) = (1 + 4x) \cdot g\left[-2\left(\frac{x}{1 + 4x}\right)^2\right].$

Каждый новый множитель добавляет две новых верные цифры в каноническом разложении результата.

Инструменты для работы с p-адическими числами есть во многих математических пакетах, таких как Sage, Maple, Magma и им подобных, на GitHub выложены библиотеки для многих языков, а теперь такая библиотека есть и в Hackage.

Пара слов о библиотеке

Для Haskell Андреем Лелеченко написан отличный пакет модулярной арифметики mod, с очень удобным интерфейсом и эффективной реализацией. Он и послужил основой для моей библиотеки. Каноническое разложение в список цифр используется только для текстового представления чисел, а вся арифметика производится в целочисленной модулярной арифметике.

Простейшая работа с библиотекой выглядит так:

> :set -XDataKinds
> 45 :: Z 5 -- целое число просто превращается в пятиричную запись
140
> -45 :: Z 5 -- отрицательное целое число, в скобках -- периодическая последовательность
(4)310
> toInteger it -- восстанавливаем целое из предыдущего результата
-45
> 13 / 7 :: Q 10 -- рациональное число
(714285)9.0
> toRational it 
13 % 7
> sqrt 11 :: Q 5 -- квадратный корень
…231012244200433234102330200211.0
> sqrt 11 :: Z' 5 40 -- увеличиваем разрядность
…00104441102231221020231012244200433234102330200211.0
> toRational (it ^ 2) -- проверяем, что корень найден верно
11 % 1
> sin 49 :: Q 7 -- пример трансцендентной функции
…013021253020521111000100.0
> fromDigits [1,2,3,0,3,2,1] :: Z 4 -- конструирование числа из цифр
1230321

Всё достаточно удобно и интуитивно.

Продемонстрируем накопление ошибки округления при работе с IEEE и её отсутствие при работе с 10-адическими числами:

> iterate (1 / 3 + ) 0 !! 300000 :: Double
99999.99999968921
> iterate (1 `div` 3 + ) 0 !! 300000 :: Z 10
100000

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

> henselLifting (\x -> x^3 - 2*x +3) (\x -> 3*x^2 -2) :: [Z 7]
[…106254154414566525205522]

Производную надо задавать явно, но можно воспользоваться автоматическим дифференцированием, которое я не захотел тащить в библиотеку, чтобы не утяжелять её зависимости.

Можно найти 10-адические автоморфные числа и убедиться в том, что в какую бы степень мы их не возвели, они останутся неизменными:

> aut = henselLifting (\x -> x^2 - x) (\x -> 2*x - 1) :: [Z 10]
> aut
[0,1,…392256259918212890625,…607743740081787109376]
> and [ x == x^i | x <- aut, i <- [2..200] ]
True

А вот ещё числа, которых не встретишь в вещественном мире: $n$ различных корней $n$ -ной степени из единицы, как в комплексных числах (правда, они существуют не для всех $n$ ).

> unityRoots 3 :: [Q 7]
[1.0,…053116412125443426203642.0,…613550254541223240463024.0]
> (^3) <$> it
[1.0,1.0,1.0]
> unityRoots 6 :: [Q 7]
[1.0,…053116412125443426203642.0,…053116412125443426203643.0,
…613550254541223240463024.0,…613550254541223240463025.0,(6).0]
> (^6) <$> it
[1.0,1.0,1.0,1.0,1.0,1.0]

Методом Гензеля можно решать и трансцендентные уравнения тоже. Например, вот 7-адическое число, синус которого равен $7^2$ :

> x = head $ henselLifting (\x -> sin x - 7^2) cos :: Q 7
> x
…313125366542105556000100.0
> sin x
100.0 -- это 49 в 7-ричном представлении.

Из-за особенностей округления в р-адических числах в них можно работать с невероятно большими числами, правда, только с младшими разрядами, наименее значащими в вещественном мире, но в модулярных арифметиках и они способны нести полезную информацию. Например, нам не составит особого труда вычислить младшие разряды знаменитого числа Грэма. Оно определяется через т.н. гипероперации — обобщения арифметических операций. Определим одну из них — тетрацию, то есть башню возведений в степень a ^^^ b = a ^ (a ^ (a ... a ^ a)) состоящую из b этажей.

(^^^) :: Radix p prec => Z' p prec -> Z' p prec -> Z' p prec 
a ^^^ 0 = 1
a ^^^ 1 = a
a ^^^ b = a `zPow` (a ^^^ (b - 1))

Здесь функция zPow — это p-адическое возведение в степень для кольца целых p-адических чисел. Значение этой функции быстро становятся такими большими, что в обычных вещественных числах с ними уже не поработаешь, тогда как в p-адической арифметике вычисления не занимают сколь-нибудь заметного времени.

Число Грэма строится на основе тетрации тройки. Нетрудно убедиться, что последовательность $a _i = 3 \text{^^^} i$ сходится в $\mathbb{Z}_{10}$ к некоторому предельному числу $\mathscr{G}$ , которое решает уравнение $3^\mathscr{G} = \mathscr{G}$ :

> λ> 3 ^^^ 2 :: Z 10
27
λ> 3 ^^^ 3 :: Z 10
7625597484987
λ> 3 ^^^ 4 :: Z 10
…206738945776100739387
λ> 3 ^^^ 5 :: Z 10
…315006939489660355387
λ> 3 ^^^ 6 :: Z 10
…913333445425126595387
λ> 3 ^^^ 10 :: Z 10
…254100942846464195387
λ> 3 ^^^ 100 :: Z 10
…104575627262464195387
λ> g = 3 ^^^ 1000 :: Z 10
…104575627262464195387

Существование этого предела означает, что какую бы огромную башню степеней тройки мы бы ни построили (а для числа Грэма она больше чем число планковских расстояний, укладывающихся в диаметр видимой Вселенной), мы получим всего лишь конечное округление числа $\mathscr{G}$ ! В русскоязычной статье из Википедии приводится 50 младших разрядов числа Грэма. Работая в p-адической арифметике, мы легко можем получить хоть 50, хоть 100 последних цифр этого монстра:

> g = 3 ^^^ 1000 :: Z' 10 100
> g
…9404248265018193851562535796399618993967905496638003222348723967018485186439059104575627262464195387
> 3 `zPow` g == g
True

Некоторые детали реализации

Работа над библиотекой padic доставила мне массу удовольствия!

Особенно приятно, что при реализации этой библиотеки удалось задействовать элементы зависимых типов, обеспечив типобезопасность и непредставимость некорректных данных, а также переведя часть вычислений на этап компиляции. На уровне типов определяются как модуль в котором мы работаем, так и задаваемая точность. Если точность не указана, она подбирается такой, чтобы можно было корректно реконструировать рациональные числа с числителями и знаменателями в пределах типа Int32 из p-адического представления. На уровне типов так же происходит и согласование модуля p-адического числа $p$ и модуля его внутреннего представления $p^k$ c требуемой точностью. Для этого используются зависимые типы, предоставляемые типами-литералами (type literals) и семействами типов (type families).

Типы-литералы, представленные библиотекой GHC.TypeLits позволяют оперировать с натуральными числами и со строками, как с фантомными типами, параметризующими обычные типы данных, и к которым есть доступ в рантайме (лишь на чтение, изменять фантомы можно только на уровне типов). В качестве простого примера их использования создадим минимальный работающий прототип для целочисленной p-адической арифметики.

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

{-# LANGUAGE DataKinds #-}
{-# LANGUAGE KindSignatures #-}

import GHC.TypeLits hiding (Mod)

newtype Mod (p :: Nat) = Mod Integer

Тип-фантом p присутствует только в левой части определения типа-обёртки над типом Integer. Его роль только помечать на уровне типов порядок кольца вычетов. Доступ к типу-параметру можно получить с помощью функции natVal:

> :t natVal
> natVal (Mod 5 :: Mod 10)
10
> :t it
Natural

В результате мы получаем «нормальное» натуральное число, с которым можно проводить все возможные вычисления. Например, мы можем определить толковый способ вывода чисел в кольце вычетов на печать:

instance KnownNat p => Show (Mod p) where
  show n@(Mod x) = show (x `mod` p) ++ " mod " ++ show p
    where p = natVal n

Класс KnownNat содержит определение для функции natVal, поэтому его необходимо указывать в качестве условия (constraint).

> Mod 5 :: Mod 10
5 mod 10
> Mod 17 :: Mod  10
7 mod 10

Получить доступ к фантомному типу-параметру можно и у ещё не созданного объекта. Вот, например, как можно реализовать базовый арифметический функционал:

instance KnownNat p => Num (Mod p) where
  fromInteger n = let res = Mod $ n `mod` natVal res in res
  Mod a + Mod b = let res = Mod $ (a + b) `mod` natVal res in res
  Mod a - Mod b = let res = Mod $ (a - b) `mod` natVal res in res
  Mod a * Mod b = let res = Mod $ (a * b) `mod` natVal res in res
  negate (Mod a) = let res = Mod (natVal res - a) in res
  abs = id
  signum _ = 1

Этот нехитрый приём позволяет нам создавать параметризованные типы из чисел, которые «узнают» о том, в какую модулярную арифметику их погрузили от вызывающей стороны.

> 5 :: Mod 3
2 mod 3
> 5 + (14 :: Mod 10)
9 mod 10

Корректное представление p-адических чисел подразумевает ряд условий. Во-первых, основание $p$ должно быть числом больше единицы. Во-вторых, внутреннее представление в виде кольца вычетов должно иметь модуль $p^k$ , где $k$ должно вычисляться исходя и заданной точности округления p-адических чисел. Всё это можно сделать на уровне типов, используя семейства типов.

-- добавляем к преамбуле
{-# LANGUAGE TypeFamilies #-}
{-# LANGUAGE TypeOperators #-}
{-# LANGUAGE NoStarIsType #-}
{-# LANGUAGE UndecidableInstances #-}

-- добавляем к списку импортов
import Data.Constraint (Constraint)
import Data.Word

-- Ограничение на допустимые значения основания
type family ValidRadix (m :: Nat) :: Constraint where
  ValidRadix 0 = TypeError ('Text "Zero radix!")
  ValidRadix 1 = TypeError ('Text "Radix should be more then 1!")
  ValidRadix m = ()

-- Основание для внутреннего представления p-адического числа
type family LiftedRadix p prec where
  LiftedRadix p prec = p ^ (2*prec + 1)

-- Ограничение для основания p-адического числа с точностью prec
type family Radix p prec :: Constraint where
  Radix p prec =
    ( KnownNat prec
    , KnownNat p, ValidRadix p
    , KnownNat (LiftedRadix p prec), ValidRadix (LiftedRadix p prec) )

-- Точность, необходимая для восстановления рациональных чисел 
-- с числителем и знаменателем, принадлежащим  типу num
type family SufficientPrecision num (p :: Nat) :: Nat where
  SufficientPrecision Word32 2 = 64
  SufficientPrecision Word32 3 = 43
  SufficientPrecision Word32 5 = 30
  SufficientPrecision Word32 6 = 26
  SufficientPrecision Word32 7 = 24
  SufficientPrecision t p = Div (SufficientPrecision t 2) (Log2 p)

Эти условия и ограничения позволяют нам определить тип для целых p-адических чисел

-- добавляем в преамбулу
{-# LANGUAGE DerivingVia #-}
{-# LANGUAGE StandaloneDeriving #-}

newtype Z' (p :: Nat) (prec :: Nat) = Z' (R prec p)
newtype R (prec :: Nat ) (p :: Nat) = R (Mod (LiftedRadix p prec))
type Z p = Z' p (SufficientPrecision Word32 p)

deriving via (Mod p) instance Radix p prec => Num (R p prec)
deriving via (Mod (LiftedRadix p prec)) instance Radix p prec => Num (R p prec)

Чехарда с двумя вложенными типами Z и R понадобилась для того, чтобы обеспечить доступ к обоим фантомным типам, при том, что функция natVal способна «считать» только внешний тип-фантом.

precision :: (Radix p prec, Z' p prec, Integral i) => Z' p prec -> i
precision z = fromIntegral (natVal z)

radix :: (Radix p prec, Z' p prec, Integral i) => Z' p prec -> i
radix (Z' r) = fromIntegral (natVal r)

На этапе компиляции все эти типы-обертки исчезнут так что в результирующем коде останутся только операции над числами типа Integer заточенные для работы с p-адическими числами.
Осталось только сделать так, чтобы p-адические числа выводились в своей системе счисления. Для этого напишем примитивную реализацию метода show:

instance Radix p prec => Show (Z' p prec) where
  show n@(Z' (R (Mod x))) = 
  foldMap show $ reverse $ take (precision n) $ toRadix (radix n) x

toRadix :: Integer -> Integer -> [Int]
toRadix _ 0 = [0]
toRadix p n = unfoldr go n
  where
    go 0 = Nothing
    go x = let (q, r) = quotRem x p in Just (fromIntegral r, q)

> 42 :: Z 5
132

-- посмотрим на внутреннее представление этого числа
> Z' (R x) = 42 :: Z' 5 10
> x
42 mod 476837158203125
-- поменяем точность
> Z' (R x) = 42 :: Z' 5 5
> x
42 mod 48828125

-- работа с отрицательными числами:
> -42 :: Z 5
444444444444444444444444444313
> -42 :: Z' 5 5
44313
> (-42 :: Z 5) + 52
20

-- получим последовательность, сходящуюся к автоморфному числу
> mapM_ print $ take 10 $ iterate (^2) (5 :: Z' 10 10)
5
25
625
390625
2587890625
6962890625
5712890625
3212890625
8212890625
8212890625

По большому счёту, того, что мы сделали, уже достаточно для реализации простейшей p-адической арифметики. Для полноценной работы нужно ещё реализовать корректные переходы от рациональных чисел к p-адическим и обратно, а также многие другие функции. Те, кому интересны, эти детали, могут посмотреть код библиотеки на GitHub.

Habrahabr.ru прочитано 5870 раз