This is an alternative site for discovering Elm packages. You may be looking for the official Elm package site instead.

# Arithmetic

A library that provides useful number-theoretical functions for dealing with integers, primes, divisibility, et cetera.

# Primes

isPrime : Int -> Bool

Test whether an integer is a positive prime.

``````isPrime 2357 == True
isPrime 500 == False
``````
primesBelow : Int -> List Int

Get all primes in the given range `[0..n-1]`, using the Sieve of Eratosthenes.

``````primesBelow 4 == [2, 3]
primesBelow 17 == [2, 3, 5, 7, 11, 13]
``````
primeFactors : Int -> List Int

Return a list of all prime factors for a given positive integer, in ascending order. If the input is less than 2, the empty list is returned.

``````primeFactors 24 == [2, 2, 2, 3]
primeFactors 767 == [13, 59]
primeFactors 1 == []
``````
primeExponents : Int -> List ( Int, Int )

Return a list of all prime-exponent pairs for a given positive integer's prime decomposition, with the primes in ascending order. If the input is less than 2, the empty list is returned.

``````primeExponents 24 == [(2, 3), (5, 2)]                -- 2^3 * 5^2
primeExponents 531764 == [(2, 1), (11, 2), (13, 3)]  -- 2^1 * 11^2 * 13^3
primeExponents 1 == []                               -- empty product
``````

# Parity

isEven : Int -> Bool

Test whether an integer is even.

``````isEven 2 == True
isEven 3 == False
``````
isOdd : Int -> Bool

Test whether an integer is odd.

``````isOdd 2 == False
isOdd 3 == True
``````

# Divisors

divides : Int -> Int -> Bool

Test whether one number divides another.

``````10 `divides` 120 == True
10 `divides` 125 == False
``````
isMultipleOf : Int -> Int -> Bool

Test whether one number is a multiple of another.

``````120 `isMultipleOf` 10 == True
125 `isMultipleOf` 10 == False
``````
divisors : Int -> List Int

Get all divisors of a number, in ascending order.

``````divisors 20 == [1, 2, 4, 5, 10, 20]
``````
properDivisors : Int -> List Int

Get all proper divisors (i.e., divisors less than the input) of a number, in ascending order.

``````properDivisors 20 == [1, 2, 4, 5, 10]
``````
divisorCount : Int -> Int

Get the number of divisors of a number (counting itself).

# GCD and LCM

gcd : Int -> Int -> Int

Calculate the greatest common divisor of two integers. `gcd x 0` and `gcd 0 x` both return `x`. Negative arguments are made positive first.

``````gcd 56 80 == 8
``````
lcm : Int -> Int -> Int

Calculate the least common multiple of two integers. `lcm x 0` and `lcm 0 x` both return `0`. Negative arguments are made positive first.

``````lcm 56 80 == 560
``````
isCoprimeTo : Int -> Int -> Bool

Test whether two integers are coprime.

``````56 `isCoprimeTo` 80 == False
5 `isCoprimeTo` 8
``````
totient : Int -> Int

Compute Euler's totient function `φ(n)`: the number of positive integers `1 <= k <= n` for which `gcd(n, k) == 1`. The input is made positive first.

``````totient 99 == 60
totient 1450 == 560
``````
extendedGcd : Int -> Int -> ( Int, Int, Int )

Given `a` and `b`, compute integers `(d, u, v)` so that ```a * u + b * v == d``` where `d == gcd a b`. (These are known as [Bézout coefficients]( https://en.wikipedia.org/wiki/B%C3%A9zout%27s_identity). If the inputs are both positive, the solution returned satisfies `abs u < b // gcd a b` and `abs v < a // gcd a b`.)

``````extendedGcd 1215 465 == (15, -13, 34)
-- because gcd 1215 465 == 15 == -13 * 1215 + 34 * 465
``````

# Base conversion

toBase : Int -> Int -> List Int

Convert a number to a list of digits in the given base. The input number is made positive first.

``````toBase 2 42 = [1, 0, 1, 0, 1, 0]  -- 42 in binary
``````
fromBase : Int -> List Int -> Int

Interpret a list of digits as a number in the given base. The input is expected to consist of integers `d` for which `0 <= d < base`.

``````fromBase 2 [1, 0, 1, 0, 1, 0] = 42
``````

# Squares

squareRoot : Float -> Float

Take the square root of a number. Return `NaN` (not a number) for negative arguments.

``````squareRoot 5.76 == 2.4
squareRoot (-1) |> isNaN
``````
safeSquareRoot : Float -> Maybe Float

Safely take the square root of a number: return `Just (squareRoot n)` if the input `n` is nonnegative; otherwise, return `Nothing`.

``````squareRoot 5.76 == Just 2.4
squareRoot (-1) == Nothing
``````
intSquareRoot : Int -> Int

Take the square root, rounding down. Return `NaN` (not a number) for negative arguments.

``````intSquareRoot 20 == 4
intSquareRoot 25 == 5
``````
exactIntSquareRoot : Int -> Maybe Int

Return `Just s` if the given integer is a square, and `s` is its square root; otherwise, return `Nothing`.

``````exactIntSquareRoot 20 == Nothing
exactIntSquareRoot 25 == Just 5
``````
isSquare : Int -> Bool

Test whether a number is a square.

``````isSquare 20 == False
isSquare 25 == True
``````

# Cubes

cubeRoot : Float -> Float

Take the cube root of a number.

``````cubeRoot 15.625 == 2.5
``````
intCubeRoot : Int -> Int

Integer cube root, rounding down.

``````intCubeRoot 800 == 9
intCubeRoot 1000 == 10
``````
exactIntCubeRoot : Int -> Maybe Int

Return `Just s` if the given integer is a cube, and `s` is its cube root; otherwise, return `Nothing`.

``````exactIntCubeRoot 800 == Nothing
exactIntCubeRoot 1000 == Just 10
``````
isCube : Int -> Bool

Test whether a number is a cube.

``````isCube 800 == False
isCube 1000 == True
``````

# Modular arithmetic

powerMod : Int -> Int -> Int -> Int

`powerMod b e m` efficiently calculates `b ^ e` (modulo `m`). It assumes `b >= 0`, `e >= 0` and `m >= 1`.

For example, to compute `4147 ^ 8671 % 1000` efficiently:

``````powerMod 4147 8671 1000 == 803
``````
modularInverse : Int -> Int -> Maybe Int

Given a number `a` and a modulus `n`, return the multiplicative inverse of `a` modulo `n`, if it exists. That is: try to return `Just b`, with `0 <= b < n`, so that `a * b == 1` modulo `n`, but return `Nothing` if no such `b` exists. (`b` exists precisely when `a` and the modulus `n` are coprime.)

``````modularInverse 3 11 == Just 4    -- 3 * 4 == 12 == 1 (mod 11)
modularInverse 3 15 == Nothing   -- 3 and 15 aren't coprime
``````
chineseRemainder : List ( Int, Int ) -> Maybe Int

Given a list of residue-modulus pairs `[(r1, m1), (r2, m2), ...]`, solve the system of linear congruences:

``````x = r1 (mod m1)
x = r2 (mod m2)
...
``````

Let `M` be the product of all moduli in the list. The Chinese remainder theorem tells us that

• if all of the moduli are pairwise coprime (their `gcd` is 1), there are infinitely many solutions, all congruent `mod M`;
• if there is a pair of non-coprime moduli in the list, there is no solution.

The result is a solution `Just x` with `0 <= x < M` in the first case; or `Nothing` if the system is unsolvable.

``````chineseRemainder [(10, 11), (4, 12), (12, 13)] == Just 1000
-- Solution to x = 10 (mod 11), x = 4 (mod 12), x = 12 (mod 13).

chineseRemainder [(2, 3), (4, 6)] == Nothing
-- 3 and 6 are not coprime, so there is no solution.

chineseRemainder [] == Just 0
-- The trivial solution, modulo M = 1.
``````
``````module Arithmetic
exposing
( isEven
, isOdd
, toBase
, fromBase
, squareRoot
, safeSquareRoot
, intSquareRoot
, exactIntSquareRoot
, isSquare
, cubeRoot
, intCubeRoot
, exactIntCubeRoot
, isCube
, divides
, isMultipleOf
, divisors
, properDivisors
, divisorCount
, gcd
, lcm
, isCoprimeTo
, totient
, extendedGcd
, powerMod
, modularInverse
, chineseRemainder
, isPrime
, primesBelow
, primeFactors
, primeExponents
)

{-| A library that provides useful number-theoretical functions for dealing
with integers, primes, divisibility, et cetera.

# Primes
@docs isPrime, primesBelow, primeFactors, primeExponents

# Parity
@docs isEven, isOdd

# Divisors
@docs divides, isMultipleOf, divisors, properDivisors, divisorCount

# GCD and LCM
@docs gcd, lcm, isCoprimeTo, totient, extendedGcd

# Base conversion
@docs toBase, fromBase

# Squares
@docs squareRoot, safeSquareRoot, intSquareRoot, exactIntSquareRoot, isSquare

# Cubes
@docs cubeRoot, intCubeRoot, exactIntCubeRoot, isCube

# Modular arithmetic
@docs powerMod, modularInverse, chineseRemainder

-}

import Array

{- Parity -}

{-| Test whether an integer is even.

isEven 2 == True
isEven 3 == False
-}
isEven : Int -> Bool
isEven n =
n % 2 == 0

{-| Test whether an integer is odd.

isOdd 2 == False
isOdd 3 == True
-}
isOdd : Int -> Bool
isOdd n =
n % 2 /= 0

{- Base conversion -}

{-| Convert a number to a list of digits in the given base. The input number is

toBase 2 42 = [1, 0, 1, 0, 1, 0]  -- 42 in binary
-}
toBase : Int -> Int -> List Int
toBase base n =
let
n_ =
abs n

go x acc =
if x <= 0 then
acc
else
go (x // base) ((x % base) :: acc)
in
go n_ []

{-| Interpret a list of digits as a number in the given base. The input is
expected to consist of integers `d` for which `0 <= d < base`.

fromBase 2 [1, 0, 1, 0, 1, 0] = 42
-}
fromBase : Int -> List Int -> Int
fromBase base =
List.foldl (\x acc -> acc * base + x) 0

{- Squares -}

{-| Take the square root of a number. Return `NaN` (not a number) for negative
arguments.

squareRoot 5.76 == 2.4
squareRoot (-1) |> isNaN
-}
squareRoot : Float -> Float
squareRoot =
sqrt

{-| Safely take the square root of a number: return `Just (squareRoot n)` if
the input `n` is nonnegative; otherwise, return `Nothing`.

squareRoot 5.76 == Just 2.4
squareRoot (-1) == Nothing
-}
safeSquareRoot : Float -> Maybe Float
safeSquareRoot n =
if n < 0 then
Nothing
else
Just (sqrt n)

{-| Take the square root, rounding down. Return `NaN` (not a number) for
negative arguments.

intSquareRoot 20 == 4
intSquareRoot 25 == 5
-}
intSquareRoot : Int -> Int
intSquareRoot =
toFloat >> sqrt >> round

{-| Return `Just s` if the given integer is a square, and `s` is its square
root; otherwise, return `Nothing`.

exactIntSquareRoot 20 == Nothing
exactIntSquareRoot 25 == Just 5
-}
exactIntSquareRoot : Int -> Maybe Int
exactIntSquareRoot n =
let
s =
intSquareRoot n
in
if s * s == n then
Just s
else
Nothing

{-| Test whether a number is a square.

isSquare 20 == False
isSquare 25 == True
-}
isSquare : Int -> Bool
isSquare n =
let
r =
n % 48
in
(r == 0 || r == 1 || r == 4 || r == 9 || r == 16 ||
r == 25 || r == 33 || r == 36) && intSquareRoot n ^ 2 == n

{- Cubes -}

{-| Take the cube root of a number.

cubeRoot 15.625 == 2.5
-}
cubeRoot : Float -> Float
cubeRoot n =
n ^ (1 / 3)

{-| Integer cube root, rounding down.

intCubeRoot 800 == 9
intCubeRoot 1000 == 10
-}
intCubeRoot : Int -> Int
intCubeRoot =
toFloat >> cubeRoot >> round

{-| Return `Just s` if the given integer is a cube, and `s` is its cube root;
otherwise, return `Nothing`.

exactIntCubeRoot 800 == Nothing
exactIntCubeRoot 1000 == Just 10
-}
exactIntCubeRoot : Int -> Maybe Int
exactIntCubeRoot n =
let
s =
intCubeRoot n
in
if s ^ 3 == n then
Just s
else
Nothing

{-| Test whether a number is a cube.

isCube 800 == False
isCube 1000 == True
-}
isCube : Int -> Bool
isCube n =
let
r =
n % 63
in
(r == 0 || r == 1 || r == 8 || r == 27 || r == 28 || r == 35 ||
r == 36 || r == 55 || r == 62) && intCubeRoot n ^ 3 == n

{- Divisors -}

{-| Test whether one number divides another.

10 `divides` 120 == True
10 `divides` 125 == False
-}
divides : Int -> Int -> Bool
divides a b =
b % a == 0

{-| Test whether one number is a multiple of another.

120 `isMultipleOf` 10 == True
125 `isMultipleOf` 10 == False
-}
isMultipleOf : Int -> Int -> Bool
isMultipleOf a b =
a % b == 0

{-| Get all divisors of a number, in ascending order.

divisors 20 == [1, 2, 4, 5, 10, 20]
-}
divisors : Int -> List Int
divisors =
let
f ( p, e ) =
List.concatMap (\a -> List.map (\x -> p ^ x * a) (List.range 0 e))
in
primeExponents >> List.foldr f [ 1 ] >> List.sort

{-| Get all proper divisors (i.e., divisors less than the input) of a number,
in ascending order.

properDivisors 20 == [1, 2, 4, 5, 10]
-}
properDivisors : Int -> List Int
properDivisors n =
divisors n
|> List.filter ((/=) n)

{-| Get the number of divisors of a number (counting itself).
-}
divisorCount : Int -> Int
divisorCount =
primeExponents >> List.map (\( _, e ) -> e + 1) >> List.product

{- GCD and LCM -}

{-| Calculate the greatest common divisor of two integers. `gcd x 0` and
`gcd 0 x` both return `x`. Negative arguments are made positive first.

gcd 56 80 == 8
-}
gcd : Int -> Int -> Int
gcd a b =
let
gcd_ a b =
if b == 0 then
a
else
gcd_ b (a % b)
in
gcd_ (abs a) (abs b)

{-| Calculate the least common multiple of two integers. `lcm x 0` and
`lcm 0 x` both return `0`. Negative arguments are made positive first.

lcm 56 80 == 560
-}
lcm : Int -> Int -> Int
lcm a b =
abs ((a // gcd a b) * b)

{-| Test whether two integers are coprime.

56 `isCoprimeTo` 80 == False
5 `isCoprimeTo` 8
-}
isCoprimeTo : Int -> Int -> Bool
isCoprimeTo a b =
gcd a b == 1

{-| Compute Euler's totient function `φ(n)`: the number of positive integers
`1 <= k <= n` for which `gcd(n, k) == 1`. The input is made positive first.

totient 99 == 60
totient 1450 == 560
-}
totient : Int -> Int
totient n =
let
n_ =
abs n

f p n =
n * (p - 1) // p
in
List.foldr f n_ (List.map Tuple.first (primeExponents n_))

{-| Given `a` and `b`, compute integers `(d, u, v)` so that `a * u + b * v ==
d` where `d == gcd a b`. (These are known as [Bézout coefficients](
https://en.wikipedia.org/wiki/B%C3%A9zout%27s_identity). If the inputs are both
positive, the solution returned satisfies `abs u < b // gcd a b` and
`abs v < a // gcd a b`.)

extendedGcd 1215 465 == (15, -13, 34)
-- because gcd 1215 465 == 15 == -13 * 1215 + 34 * 465
-}
extendedGcd : Int -> Int -> ( Int, Int, Int )
extendedGcd a b =
let
egcd n1 o1 n2 o2 r s =
if s == 0 then
( r, o1, o2 )
else
let
q =
r // s
in
egcd (o1 - q * n1) n1 (o2 - q * n2) n2 s (rem r s)

( d, x, y ) =
egcd 0 1 1 0 (abs a) (abs b)

u =
if a < 0 then
negate x
else
x

v =
if b < 0 then
negate y
else
y
in
( d, u, v )

{- Modular arithmetic -}

{-| `powerMod b e m` efficiently calculates `b ^ e` (modulo `m`). It assumes
`b >= 0`, `e >= 0` and `m >= 1`.

For example, to compute `4147 ^ 8671 % 1000` efficiently:

powerMod 4147 8671 1000 == 803
-}
powerMod : Int -> Int -> Int -> Int
powerMod base exponent modulus =
let
go b e r =
if e == 0 then
r
else
let
r_ =
if isOdd e then
(r * b) % modulus
else
r
in
go (b * b % modulus) (e // 2) r_
in
if modulus == 1 then
0
else
go (base % modulus) exponent 1

{- Halve this number until it is odd. Then, return a tuple `(k, m)`, where
`k` is the number of times the input was halved, and `m` is the resulting odd
number. In other words, the original number equals `(2 ^ k) * m`.

shiftToOdd 999 == (0, 999)
shiftToOdd 1000 == (3, 125)
-}

shiftToOdd : Int -> ( Int, Int )
shiftToOdd n =
let
f k m =
if m % 2 == 1 then
( k, m )
else
f (k + 1) (m // 2)
in
f 0 n

{-| Given a number `a` and a modulus `n`, return the multiplicative inverse of
`a` modulo `n`, if it exists. That is: try to return `Just b`, with
`0 <= b < n`, so that `a * b == 1` modulo `n`, but return `Nothing` if no such
`b` exists. (`b` exists precisely when `a` and the modulus `n` are coprime.)

modularInverse 3 11 == Just 4    -- 3 * 4 == 12 == 1 (mod 11)
modularInverse 3 15 == Nothing   -- 3 and 15 aren't coprime
-}
modularInverse : Int -> Int -> Maybe Int
modularInverse a modulus =
let
( d, x, _ ) =
extendedGcd a modulus
in
if d == 1 then
Just (x % modulus)
else
Nothing

{-| Given a list of residue-modulus pairs `[(r1, m1), (r2, m2), ...]`, solve
the system of linear congruences:

x = r1 (mod m1)
x = r2 (mod m2)
...

Let `M` be the product of all moduli in the list. The [Chinese remainder
theorem](https://en.wikipedia.org/wiki/Chinese_remainder_theorem) tells us that

* if all of the moduli are pairwise coprime (their `gcd` is 1), there are
infinitely many solutions, all congruent `mod M`;
* if there is a pair of non-coprime moduli in the list, there is no solution.

The result is a solution `Just x` with `0 <= x < M` in the first case; or
`Nothing` if the system is unsolvable.

chineseRemainder [(10, 11), (4, 12), (12, 13)] == Just 1000
-- Solution to x = 10 (mod 11), x = 4 (mod 12), x = 12 (mod 13).

chineseRemainder [(2, 3), (4, 6)] == Nothing
-- 3 and 6 are not coprime, so there is no solution.

chineseRemainder [] == Just 0
-- The trivial solution, modulo M = 1.
-}
chineseRemainder : List ( Int, Int ) -> Maybe Int
chineseRemainder equations =
let
( residues, moduli ) =
List.unzip equations

m =
List.product moduli

v =
List.map (\x -> m // x) moduli

fromJustCons x acc =
case x of
Just y ->
Maybe.map ((::) y) acc

Nothing ->
Nothing

fromJustList =
List.foldr fromJustCons (Just [])

inverses =
fromJustList (List.map2 modularInverse v moduli)
in
fromJustList (List.map2 modularInverse v moduli)
|> Maybe.map
(List.map2 (*) residues
>> List.map2 (*) v
>> List.sum
>> flip (%) m
)

{- Primes -}

{-| Test whether an integer is a positive prime.

isPrime 2357 == True
isPrime 500 == False
-}
isPrime : Int -> Bool
isPrime n =
if n < 13 then
n == 2 || n == 3 || n == 5 || n == 7 || n == 11
else if n % 2 == 0 || n % 3 == 0 || n % 5 == 0 then
False
else if n < 1373653 then
millerRabin n [ 2, 3 ]
else
-- n < 2152302898747
millerRabin n [ 2, 3, 5, 7, 11 ]

{-| Perform a Miller-Rabit pseudoprimality test on `n` with the given
witnesses, which should all be in the range `[2..n-2]`. Return `True` if the
integer is a probable prime, and `False` if it is definitely composite.

millerRabin 2357 [2, 3] == True
millerRabin 500 [2, 3] == False
-}
millerRabin : Int -> List Int -> Bool
millerRabin n witnesses =
let
( s, d ) =
shiftToOdd (n - 1)

check l x =
if l <= 0 then
True
else
let
y =
powerMod x 2 n
in
y == 1 || (y /= n - 1 && check (l - 1) y)

go witnesses =
case witnesses of
[] ->
True

a :: rest ->
let
x =
powerMod a d n
in
if x == 1 || x == n - 1 then
go rest
else if check (s - 1) x then
False
else
go rest
in
go witnesses

{-| Get all primes in the given range `[0..n-1]`, using the Sieve of
Eratosthenes.

primesBelow 4 == [2, 3]
primesBelow 17 == [2, 3, 5, 7, 11, 13]
-}
primesBelow : Int -> List Int
primesBelow n =
let
ps =
2 :: List.map (\x -> 2 * x + 1) (List.range 1 (intSquareRoot n // 2))

initial =
Array.repeat n True
|> Array.set 0 False
|> Array.set 1 False

sieve p arr =
let
mark i p n arr =
if i * p >= n then
arr
else
mark (i + 1) p n (Array.set (i * p) False arr)
in
if Array.get p arr == Just True then
mark 2 p n arr
else
arr

trueIndices =
let
f i pred =
if pred then
Just i
else
Nothing

g x acc =
case x of
Just x ->
x :: acc

Nothing ->
acc
in
Array.indexedMap f >> Array.foldr g []
in
trueIndices (List.foldr sieve initial ps)

{-| Return a list of all prime factors for a given positive integer, in
ascending order. If the input is less than 2, the empty list is returned.

primeFactors 24 == [2, 2, 2, 3]
primeFactors 767 == [13, 59]
primeFactors 1 == []
-}
primeFactors : Int -> List Int
primeFactors n =
let
go p n factors =
if p * p > n then
List.reverse factors ++ [ n ]
else if n % p == 0 then
go p (n // p) (p :: factors)
else
go (p + 1 + p % 2) n factors
in
if n <= 1 then
[]
else
go 2 n []

{-| Return a list of all prime-exponent pairs for a given positive integer's
prime decomposition, with the primes in ascending order. If the input is less
than 2, the empty list is returned.

primeExponents 24 == [(2, 3), (5, 2)]                -- 2^3 * 5^2
primeExponents 531764 == [(2, 1), (11, 2), (13, 3)]  -- 2^1 * 11^2 * 13^3
primeExponents 1 == []                               -- empty product
-}
primeExponents : Int -> List ( Int, Int )
primeExponents =
let
runLengthCons x acc =
case acc of
[] ->
[ ( x, 1 ) ]

( y, n ) :: rest ->
if x == y then
(( y, n + 1 ) :: rest)
else
(( x, 1 ) :: ( y, n ) :: rest)
in
primeFactors >> List.foldr runLengthCons []
```
```