Anonymous user
P-Adic numbers, basic: Difference between revisions
→{{header|Haskell}}: changed to more complete version
(→{{header|Haskell}}: added Haskell solution) |
(→{{header|Haskell}}: changed to more complete version) |
||
Line 1,095:
=={{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 with infinite part going to the left.
p-Adic arithmetics and conversion between rationals is implemented as instances of <code>Eq</code>, <code>Num</code>, <code>Fractional</code> and <code>Real</code> classes, so, they could be treated as usual real numbers (up to existence of some rationals for non-prime bases).
<lang haskell>{-# LANGUAGE KindSignatures, DataKinds #-}
module Padic where
import Data.Ratio
import Data.
import GHC.TypeLits
data Padic
-- valuation of the base
modulo :: (KnownNat p, Integral i) => Padic p -> i
modulo = fromIntegral . natVal
-- Constructor for zero value
pZero :: KnownNat p => Padic p
pZero = Padic (repeat 0) 0
-- Smart constructor, adjusts trailing zeros with the order.
mkPadic :: (KnownNat p,
mkPadic u k = go 0 (fromIntegral <$> u)
where
go 17 _ = pZero
go i (0:u) = go (i+1) u
go i u = Padic
-- Constructor for
-- Zero test (up to 1/p^17)
isZero :: KnownNat p => Padic p -> Bool
isZero (Padic
isZero _ = False
-- p-adic norm
pNorm :: KnownNat p => Padic p -> Ratio Int
pNorm
pNorm p = fromIntegral (modulo p) ^^ (- order p)
-- test for an integerness up to p^-17
isInteger :: KnownNat p => Padic p -> Bool
isInteger Null = False
isInteger (Padic s k) = case splitAt k s of
([],i) -> length (takeWhile (==0) $ reverse (take 20 i)) > 3
_ -> False
-- p-adics are shown with 1/p^17 precision
instance KnownNat p => Show (Padic p) where
show
show (modulo x) ++ "-adic: " ++
(case si of {[] -> "0"; _ -> si})
++ "." ++
Line 1,142 ⟶ 1,166:
showD n = [(['0'..'9']++['a'..'z']) !! n]
instance KnownNat p => Eq (Padic p) where
a == b = isZero (a - b)
instance
compare = error "Ordering is undefined fo p-adics."
instance KnownNat p => Num (Padic p) where
fromInteger
fromInteger n = pAdic (fromInteger n)
x@(Padic
where
k = ka `max` kb
s = addMod
(replicate (k-ka) 0 ++ a)
(replicate (k-kb) 0 ++ b)
_ + _ = Null
x@(Padic
mkPadic (mulMod (modulo x) a b) (ka + kb)
_ * _ = Null
negate x@(Padic u k) =
case map (\y -> modulo x - 1 - y) u of
n:ns -> Padic ((n+1):ns) k
[] -> pZero
negate _ = Null
abs p = pAdic (pNorm p)
signum = undefined
------------------------------------------------------------
-- conversion from rationals to p-adics
instance KnownNat p => Fractional (Padic p) where
fromRational = pAdic
recip Null = Null
recip x@(Padic (u:us) k)
| isZero x = Null
| gcd p u /= 1 = Null
| otherwise = mkPadic res (-k)
where
p = modulo x
res = longDivMod p (1:repeat 0) (u:us)
pAdic :: (Show i, Integral i, KnownNat p)
=> Ratio i -> Padic p
pAdic 0 = pZero
pAdic x = res
where
p = modulo res
(k, q) = getUnit p x
(n, d) = (numerator q, denominator q)
res = maybe Null process $ recipMod p d
process r = mkPadic (series n) k
where
series n
| n == 0 = repeat 0
| n `mod` p == 0 = 0 : series (n `div` p)
| otherwise =
in m : series ((n - m * d) `div` p)
------------------------------------------------------------
-- conversion from p-adics to rationals
-- works for relatively small denominators
instance KnownNat p => Real (Padic p) where
toRational Null = error "no rational representation!"
toRational x@(Padic s k) = res
where
p = modulo x
res = case break isInteger $ take 10000 $ iterate (x +) x of
(_,[]) -> - toRational (- x)
(d, i:_) -> (fromBase p (unit i) * (p^(- order i))) % (genericLength d + 1)
fromBase p = foldr (\x r -> r*p + x) 0 .
take 20 . map fromIntegral
--------------------------------------------------------------------------------
-- helper functions
-- extracts p-adic unit from a rational number
getUnit ::
getUnit p x = (
where
(k1,b:_) = span (\n -> denominator n `mod` p == 0) $
Line 1,204 ⟶ 1,260:
iterate (/ fromIntegral p) b
-- Reciprocal of a number modulo p (extended Euclidean
-- For non-prime p returns Nothing non-
recipMod ::
recipMod p 1 = Just 1
recipMod p a | gcd p a == 1 = Just $ go 0 1 p a
Line 1,236 ⟶ 1,284:
in r : go q xs ys
--
subMod p a (b:bs) = addMod p a $ (p-b) : ((p - 1 -) <$> bs)
-- Multiplication of two sequences modulo p
mulMod p as [b] =
mulMod p as bs = case as of
[1] ->
[a] -> go 0 bs
where
go s [] = [s]
go s (b:bs) =
let (q, r) = (a * b + s) `divMod` p
in r : go q bs
as -> go bs
where
go [] = []
go (b:bs) =
let c:cs = mulMod p [b] as
in c : addMod p (go bs) cs
-- Division of two sequences modulo p
longDivMod p a (b:bs) = case recipMod p b of
Nothing -> error $
show b ++ " is not invertible modulo " ++ show p
Just r -> go a
where
go [] = []
go (0:xs) = 0 : go xs
go (x:xs) =
let m = (x*r) `mod` p
_:zs = subMod p (x:xs) (mulMod p [m] (b:bs))
in m : go zs</lang>
Convertation between rationals and p-adic numbers
<pre>λ>
5-adic: 0.0
λ>
5-adic: 3.0
λ> 1/3 :: Padic 5
5-adic: 13131313131313132.0
λ>
5-adic: 0.02
λ>
2-adic: 1011100001010010.0
λ>
7-adic: 36303630363036304.0
λ>
2 % 25
λ>
10-adic: 25.0
λ>
Null
λ> -12/23 :: Padic 3
3-adic: 21001011200210010.0
λ> toRational it
(-12) % 23</pre>
Arithmetic:
<pre>λ> pAdic (12/25) + pAdic (23/56) :: Padic 7
7-adic: 32523252325232524.2
λ> pAdic (12/25 + 23/56) :: Padic 7
7-adic: 32523252325232524.2
λ> toRational it
1247 % 1400
λ> 12/25 + 23/56 :: Rational
1247 % 1400
λ> let x = 2/7 :: Padic 13
λ>
> toRational it
8 % 49
λ>
=={{header|Julia}}==
| |||