P-Adic numbers, basic: Difference between revisions
Content added Content deleted
m (Changed comment.) |
(→{{header|Haskell}}: added Haskell solution) |
||
Line 1,093: | Line 1,093: | ||
47099/10977 |
47099/10977 |
||
</pre> |
</pre> |
||
=={{header|Haskell}}== |
|||
p-Adic numbers in this implementation are represented in floating point manner, with p-adic unit as mantissa and p-adic norm as an exponent. Textual presentation is given in traditional form whith infinite part going to the left. |
|||
<lang haskell>module Padic where |
|||
import Data.Ratio |
|||
import Data.Maybe |
|||
data Padic = Padic { pBase :: Int |
|||
, pUnit :: [Int] |
|||
, pOrder :: Int } |
|||
-- Smart constructor, adjusts trailing zeros with the order. |
|||
mkPadic p u k = go 0 u |
|||
where |
|||
go 17 _ = pZero p |
|||
go i (0:u) = go (i+1) u |
|||
go i u = Padic p u (k-i) |
|||
-- Constructor for zero value |
|||
pZero :: Int -> Padic |
|||
pZero p = Padic p [] 0 |
|||
-- Zero test (up to 1/p^17) |
|||
isZero :: Padic -> Bool |
|||
isZero (Padic _ u _) = all (== 0) (take 17 u) |
|||
-- p-adic norm |
|||
pNorm :: Padic -> Ratio Int |
|||
pNorm (Padic p _ k) = fromIntegral p ^^ (-k) |
|||
-- p-adics are shown with 1/p^17 precision |
|||
instance Show Padic where |
|||
show (Padic p u k) = |
|||
show p ++ "-adic: " ++ |
|||
(case si of {[] -> "0"; _ -> si}) |
|||
++ "." ++ |
|||
(case f of {[] -> "0"; _ -> sf}) |
|||
where |
|||
(f,i) = case compare k 0 of |
|||
LT -> ([], replicate (-k) 0 ++ u) |
|||
EQ -> ([], u) |
|||
GT -> splitAt k (u ++ repeat 0) |
|||
sf = foldMap showD $ reverse $ take 17 f |
|||
si = foldMap showD $ dropWhile (== 0) $ reverse $ take 17 i |
|||
el s = if length s > 16 then "…" else "" |
|||
showD n = [(['0'..'9']++['a'..'z']) !! n] |
|||
-------------------------------------------------------------------------------- |
|||
-- basic arythmetic |
|||
instance Eq Padic where |
|||
a == b = isZero (a - b) |
|||
instance Num Padic where |
|||
fromInteger _ = undefined |
|||
Padic p a ka + Padic _ b kb = mkPadic p s k |
|||
where |
|||
k = ka `max` kb |
|||
s = addMod p (replicate (k-ka) 0 ++ a) (replicate (k-kb) 0 ++ b) |
|||
Padic p a ka * Padic _ b kb = mkPadic p (mulMod p a b) (ka + kb) |
|||
negate (Padic p u k) = |
|||
case map (\x -> p - 1 - x) u of |
|||
n:ns -> Padic p ((n+1):ns) k |
|||
[] -> pZero p |
|||
abs p = fromRat (pBase p) (pNorm p) |
|||
signum = undefined |
|||
-------------------------------------------------------------------------------- |
|||
-- conversion between rationals and p-adics |
|||
fromRat :: Int -> Ratio Int -> Padic |
|||
fromRat p 0 = pZero p |
|||
fromRat p x = mkPadic p (series n) k |
|||
where |
|||
(k, q) = getUnit p x |
|||
(n, d) = (numerator q, denominator q) |
|||
r = fromMaybe |
|||
(error $ concat [show x, " can't be represented as " |
|||
, show p, "-adic."]) |
|||
(recipMod p d) |
|||
series n |
|||
| n == 0 = repeat 0 |
|||
| n `mod` p == 0 = 0 : series (n `div` p) |
|||
| otherwise = |
|||
let m = (n * r) `mod` p |
|||
in m : series ((n - m * d) `div` p) |
|||
-- rational reconstruction |
|||
toRat :: Padic -> Ratio Int |
|||
toRat x@(Padic p s k) = |
|||
(fromBase p (pUnit i) * (p^(- pOrder i))) % length d |
|||
where |
|||
(d, i:_) = break isInteger $ iterate (x +) $ pZero p |
|||
fromBase p = foldr (\x r -> r*p + x) 0 . take 20 |
|||
-- extracts p-adic unit from a rational number |
|||
getUnit :: Int -> Ratio Int -> (Int, Ratio Int) |
|||
getUnit p x = (length k1 - length k2, c) |
|||
where |
|||
(k1,b:_) = span (\n -> denominator n `mod` p == 0) $ |
|||
iterate (* fromIntegral p) x |
|||
(k2,c:_) = span (\n -> numerator n `mod` p == 0) $ |
|||
iterate (/ fromIntegral p) b |
|||
-- naive test for an integerness up to p^-17 |
|||
isInteger :: Padic -> Bool |
|||
isInteger (Padic p s k) = case splitAt k s of |
|||
([],i) -> length (takeWhile (==0) $ reverse (take 20 i)) > 3 |
|||
_ -> False |
|||
-------------------------------------------------------------------------------- |
|||
-- helper functions |
|||
-- Reciprocal of a number modulo p (extended Euclidean algorythm). |
|||
-- For non-prime p returns Nothing non-inverible element of the ring. |
|||
recipMod :: Int -> Int -> Maybe Int |
|||
recipMod p 1 = Just 1 |
|||
recipMod p a | gcd p a == 1 = Just $ go 0 1 p a |
|||
| otherwise = Nothing |
|||
where |
|||
go t _ _ 0 = t `mod` p |
|||
go t nt r nr = |
|||
let q = r `div` nr |
|||
in go nt (t - q*nt) nr (r - q*nr) |
|||
-- Addition of two sequences modulo p |
|||
addMod p = go 0 |
|||
where |
|||
go 0 [] ys = ys |
|||
go 0 xs [] = xs |
|||
go s [] ys = go 0 [s] ys |
|||
go s xs [] = go 0 xs [s] |
|||
go s (x:xs) (y:ys) = |
|||
let (q, r) = (x + y + s) `divMod` p |
|||
in r : go q xs ys |
|||
-- Multiplication of sequence by a number modulo p |
|||
scaleMod p a = go 0 |
|||
where |
|||
go s [] = [s] |
|||
go s (b:bs) = |
|||
let (q, r) = (a * b + s) `divMod` p |
|||
in r : go q bs |
|||
-- Multiplication of two sequences modulo p |
|||
mulMod p as = go |
|||
where |
|||
go [] = [] |
|||
go (b:bs) = |
|||
let c:cs = scaleMod p b as |
|||
in c : addMod p (go bs) cs</lang> |
|||
Convertation between rationals and p-adic numbers |
|||
<pre>λ> pZero 5 |
|||
5-adic: 0.0 |
|||
λ> fromRat 5 (1/3) |
|||
5-adic: 13131313131313132.0 |
|||
λ> fromRat 5 (25/3) |
|||
5-adic: 13131313131313200.0 |
|||
λ> fromRat 5 (2/75) |
|||
5-adic: 31313131313131313.14 |
|||
λ> fromRat 5 (2/25) |
|||
5-adic: 0.02 |
|||
λ> fromRat 2 (2/25) |
|||
2-adic: 1011100001010010.0 |
|||
λ> fromRat 7 (2/25) |
|||
7-adic: 36303630363036304.0 |
|||
λ> toRat it |
|||
2 % 25 |
|||
λ> fromRat 10 (25) |
|||
10-adic: 25.0 |
|||
λ> fromRat 10 (2/25) |
|||
*** Exception: 2 % 25 can't be represented as 10-adic.</pre> |
|||
Addition and subtraction. |
|||
<pre>λ> fromRat 7 (12/25) + fromRat 7 (23/567) |
|||
7-adic: 51251554136620421.4 |
|||
λ> fromRat 7 (12/25 + 23/567) |
|||
7-adic: 51251554136620421.4 |
|||
λ> toRat it |
|||
7379 % 14175 |
|||
λ> 12/25 + 23/567 :: Rational |
|||
7379 % 14175 |
|||
λ> fromRat 3 (12/25) - fromRat 3 (23/567) |
|||
3-adic: 20210021212110012.1121 |
|||
λ> fromRat 3 (23/567) - fromRat 3 (12/25) |
|||
3-adic: 2012201010112210.1102 |
|||
λ> fromRat 17 1 == fromRat 17 (1/2) + fromRat 17 (1/2) |
|||
True</pre> |
|||
=={{header|Julia}}== |
=={{header|Julia}}== |