Anonymous user
P-Adic numbers, basic: Difference between revisions
→{{header|Haskell}}: added Haskell solution
m (Changed comment.) |
(→{{header|Haskell}}: added Haskell solution) |
||
Line 1,093:
47099/10977
</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}}==
| |||