diff --git a/CHANGELOG.md b/CHANGELOG.md
new file mode 100644
--- /dev/null
+++ b/CHANGELOG.md
@@ -0,0 +1,23 @@
+# Changelog
+
+All notable changes to this project will be documented in this file.
+
+## 0.1.0.0 — Initial release
+
+* Initial public release.
+* Added discrete probability distributions:
+
+  * Binomial
+  * Negative Binomial
+  * Geometric
+  * Multinomial
+* Added continuous probability distributions:
+
+  * Normal
+  * Exponential
+  * Gamma
+  * Uniform
+* Added probability mass, density, and cumulative distribution functions where applicable.
+* Added explicit error handling through `Either`.
+* Added numerical support for logarithmic factorial, combination, Gamma, Beta, and error functions.
+* Added Cabal test suite covering numerical results, boundary conditions, and invalid inputs.
diff --git a/LICENSE b/LICENSE
new file mode 100644
--- /dev/null
+++ b/LICENSE
@@ -0,0 +1,29 @@
+BSD 3-Clause License
+
+Copyright (c) 2026, Barış Barış
+
+Redistribution and use in source and binary forms, with or without
+modification, are permitted provided that the following conditions are met:
+
+1. Redistributions of source code must retain the above copyright notice,
+   this list of conditions and the following disclaimer.
+
+2. Redistributions in binary form must reproduce the above copyright notice,
+   this list of conditions and the following disclaimer in the documentation
+   and/or other materials provided with the distribution.
+
+3. Neither the name of the copyright holder nor the names of its
+   contributors may be used to endorse or promote products derived from
+   this software without specific prior written permission.
+
+THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
+AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
+LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
+CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
+SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
+INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
+CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
+ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
+POSSIBILITY OF SUCH DAMAGE.
diff --git a/README.md b/README.md
new file mode 100644
--- /dev/null
+++ b/README.md
@@ -0,0 +1,351 @@
+# probability-dist
+
+A Haskell library providing probability distributions and related probability functions.
+
+The library is designed around a small, explicit public API. Discrete and continuous probability distributions are exposed through separate modules, allowing applications to depend only on the functionality they need.
+
+## Features
+
+* Discrete probability distributions
+* Continuous probability distributions
+* Probability mass functions (PMF)
+* Probability density functions (PDF)
+* Cumulative distribution functions (CDF)
+* Distribution-specific moments where implemented
+* Numerically stable logarithmic calculations for combinatorial and special functions
+* Explicit error handling through `Either`
+
+## Requirements
+
+* GHC 9.10.3 or compatible GHC version
+* Cabal 3.16.1.0 or compatible Cabal version
+* `base >= 4.20 && < 5`
+
+## Installation
+
+Clone the repository and build it with Cabal:
+
+```bash
+git clone <repository-url>
+cd probability-dist
+cabal build
+```
+
+Run the test suite with:
+
+```bash
+cabal test
+```
+
+The package can also be built as a source distribution:
+
+```bash
+cabal sdist
+```
+
+## Public API
+
+The library exposes two public modules.
+
+### Discrete distributions
+
+```haskell
+import Probability.Discrete
+```
+
+This module contains the discrete probability distributions implemented by the package.
+
+Current distributions include:
+
+* Binomial
+* Negative Binomial
+* Geometric
+* Multinomial
+
+### Continuous distributions
+
+```haskell
+import Probability.Continuous
+```
+
+This module contains the continuous probability distributions implemented by the package.
+
+Current distributions include:
+
+* Normal
+* Exponential
+* Gamma
+* Uniform
+
+The modules are intentionally separated so that an application requiring only discrete or only continuous distributions does not need to import both APIs.
+
+## Basic Usage
+
+### Binomial distribution
+
+The binomial PMF is exposed through `binomialPMF`.
+
+```haskell
+import Probability.Discrete
+
+main :: IO ()
+main = do
+    print (binomialPMF 3 10 0.5)
+```
+
+The parameters are:
+
+```text
+k  n  p
+```
+
+where:
+
+* `k` is the number of successes
+* `n` is the number of trials
+* `p` is the probability of success
+
+The result is returned as an `Either` value, allowing invalid parameters to be handled explicitly.
+
+For example:
+
+```haskell
+binomialPMF 3 10 0.5
+```
+
+returns a `Right` value containing the probability.
+
+An invalid success count produces an error:
+
+```haskell
+binomialPMF 11 10 0.5
+```
+
+which returns a `Left` value.
+
+## Negative Binomial and Geometric Distributions
+
+The geometric distribution is implemented as the special case of the negative binomial distribution where:
+
+```text
+r = 1
+```
+
+For example:
+
+```haskell
+import Probability.Discrete
+
+main :: IO ()
+main = do
+    print (negativeBinomialPMF 2 1 0.5)
+    print (geometricPMF 2 0.5)
+```
+
+The geometric distribution is therefore provided as a convenience function rather than requiring users to manually express it as a negative binomial distribution with `r = 1`.
+
+## Multinomial Distribution
+
+The multinomial PMF accepts the total number of trials, a vector of category counts, and a corresponding probability vector.
+
+```haskell
+import Probability.Discrete
+
+main :: IO ()
+main = do
+    print $
+        multinomialPMF
+            10
+            [4, 3, 3]
+            [0.4, 0.3, 0.3]
+```
+
+The count vector and probability vector must be dimensionally compatible.
+
+Invalid input is represented through the package's error type rather than silently producing an invalid probability.
+
+## Normal Distribution
+
+The normal distribution provides PDF and CDF functionality.
+
+```haskell
+import Probability.Continuous
+
+main :: IO ()
+main = do
+    print (normalPDF 0.0 0.0 1.0)
+    print (normalCDF 0.0 0.0 1.0)
+```
+
+For the standard normal distribution:
+
+```text
+μ = 0
+σ = 1
+```
+
+the CDF at zero is:
+
+```text
+0.5
+```
+
+The standard deviation is validated explicitly; non-positive values result in an error.
+
+## Exponential Distribution
+
+The exponential distribution is parameterized by its rate `λ`.
+
+```haskell
+import Probability.Continuous
+
+main :: IO ()
+main = do
+    print (exponentialPDF 1.0 2.0)
+    print (exponentialCDF 1.0 2.0)
+```
+
+The implementation uses numerically appropriate calculations for expressions such as:
+
+```text
+1 - exp(-λx)
+```
+
+in order to improve numerical behavior for small values.
+
+## Gamma Distribution
+
+The Gamma distribution is parameterized by shape and rate parameters.
+
+```haskell
+import Probability.Continuous
+
+main :: IO ()
+main = do
+    print (gammaPDF 1.0 2.0 1.0)
+    print (gammaMean 2.0 1.0)
+    print (gammaVar 2.0 1.0)
+```
+
+The package also provides the corresponding mean and variance functions.
+
+The exponential distribution can be viewed as a special case of the Gamma distribution with shape parameter equal to one.
+
+## Uniform Distribution
+
+The continuous uniform distribution is parameterized by its lower and upper bounds.
+
+```haskell
+import Probability.Continuous
+
+main :: IO ()
+main = do
+    print (uniformPDF 0.5 0.0 1.0)
+    print (uniformCDF 0.5 0.0 1.0)
+    print (uniformMean 0.0 1.0)
+    print (uniformVar 0.0 1.0)
+```
+
+Invalid bounds are reported through the package's error handling mechanism.
+
+## Error Handling
+
+Public distribution functions return results using `Either`.
+
+This makes invalid statistical parameters explicit instead of relying on exceptions or silently returning meaningless numerical values.
+
+For example:
+
+```haskell
+case binomialPMF 11 10 0.5 of
+    Right probability ->
+        print probability
+
+    Left err ->
+        print err
+```
+
+The package defines distribution-related errors in its internal error module and exposes the resulting error values through the public functions.
+
+## Numerical Design
+
+Several calculations involved in probability distributions can become numerically unstable when performed directly.
+
+The package therefore uses logarithmic forms where appropriate, particularly for:
+
+* factorial-related calculations
+* combinations
+* Gamma functions
+* Beta functions
+* probability expressions involving products of many terms
+
+Internally, the package provides mathematical support functions such as:
+
+* `logFactorial`
+* `logCombination`
+* `logGamma`
+* `logBeta`
+* `erf`
+
+These functions are implementation details and are kept outside the public API.
+
+This separation allows the public distribution modules to remain focused on probability distributions while the numerical machinery remains internal to the package.
+
+## Testing
+
+The project includes a Cabal test suite covering:
+
+* ordinary distribution values
+* boundary conditions
+* invalid parameters
+* probability constraints
+* dimensional consistency
+* numerical results
+* degenerate cases
+* error handling
+
+Run all tests with:
+
+```bash
+cabal test
+```
+
+The package is also checked with:
+
+```bash
+cabal check
+```
+
+and can be packaged with:
+
+```bash
+cabal sdist
+```
+
+## Project Structure
+
+```text
+probability-dist/
+├── src/
+│   └── Probability/
+│       ├── Continuous.hs
+│       ├── Discrete.hs
+│       ├── Error.hs
+│       └── Math.hs
+├── test/
+│   └── Main.hs
+├── LICENSE
+├── README.md
+├── probability-dist.cabal
+└── CHANGELOG.md
+```
+
+`Probability.Discrete` and `Probability.Continuous` form the public API.
+
+`Probability.Math` and `Probability.Error` are internal implementation modules.
+
+## License
+
+This project is licensed under the BSD 3-Clause License.
+
+See the `LICENSE` file for the complete license text.
diff --git a/probability-dist.cabal b/probability-dist.cabal
new file mode 100644
--- /dev/null
+++ b/probability-dist.cabal
@@ -0,0 +1,47 @@
+cabal-version: 3.0
+
+name:                probability-dist
+version:             0.1.0.0
+synopsis:            Probability distributions in Haskell
+description:
+    A Haskell library providing discrete and continuous probability
+    distributions.
+category:            Mathematics
+license:             BSD-3-Clause
+license-file:        LICENSE
+author:              BARIŞ BARIŞ
+maintainer:          barisbaris2005@gmail.com
+build-type:          Simple
+extra-doc-files:
+        README.md
+        CHANGELOG.md
+
+library
+    hs-source-dirs:      src
+
+    exposed-modules:
+        Probability.Discrete
+        Probability.Continuous
+
+    other-modules:
+        Probability.Error
+        Probability.Math
+
+    build-depends:
+        base >= 4.18 && < 5
+    default-language:    Haskell2010
+
+test-suite probability-dist-test
+    type:         exitcode-stdio-1.0
+    hs-source-dirs:      test
+    main-is:             Main.hs
+
+    build-depends:
+        base,
+        probability-dist
+
+    default-language:    Haskell2010
+
+source-repository head
+    type:     git
+    location: https://github.com/barisbarisgithub/probability-dist
diff --git a/src/Probability/Continuous.hs b/src/Probability/Continuous.hs
new file mode 100644
--- /dev/null
+++ b/src/Probability/Continuous.hs
@@ -0,0 +1,340 @@
+module Probability.Continuous
+(
+    normalPDF
+    ,normalLogPDF
+    ,normalCDF
+    ,normalMean
+    ,normalVar
+    ,exponentialLogPDF
+    ,exponentialPDF
+    ,exponentialCDF
+    ,exponentialMean
+    ,exponentialVar
+    ,gammaLogPDF
+    ,gammaPDF
+    ,gammaMean
+    ,gammaVar
+    ,uniformLogPDF
+    ,uniformPDF
+    ,uniformCDF
+    ,uniformMean
+    ,uniformVar
+)where
+
+import Probability.Error (ProbabilityError(..))
+import Probability.Math
+    (logFactorial
+    , logCombination
+    , logGamma
+    , logBeta
+    ,erf
+    ,expm1)
+normalLogPDF
+    :: Double
+    -> Double
+    -> Double
+    -> Either ProbabilityError Double
+normalLogPDF x mu sigma
+    | sigma <= 0 =
+        Left (InvalidStandardDeviation sigma)
+
+    | otherwise =
+        let z = (x - mu) / sigma
+
+            logP =
+                - log sigma
+                - 0.5 * log (2.0 * pi)
+                - 0.5 * z * z
+
+        in Right logP
+
+
+normalPDF
+    :: Double
+    -> Double
+    -> Double
+    -> Either ProbabilityError Double
+normalPDF x mu sigma
+    | sigma <= 0 =
+        Left (InvalidStandardDeviation sigma)
+
+    | otherwise = do
+        logP <- normalLogPDF x mu sigma
+        pure (exp logP)
+
+
+normalCDF
+    :: Double
+    -> Double
+    -> Double
+    -> Either ProbabilityError Double
+normalCDF x mu sigma
+    | sigma <= 0 =
+        Left (InvalidStandardDeviation sigma)
+
+    | otherwise =
+        let z = (x - mu) / (sigma * sqrt 2.0)
+
+        in Right
+            (0.5 * (1.0 + erf z))
+
+
+normalMean
+    :: Double
+    -> Double
+    -> Either ProbabilityError Double
+normalMean mu sigma
+    | sigma <= 0 =
+        Left (InvalidStandardDeviation sigma)
+
+    | otherwise =
+        Right mu
+
+
+normalVar
+    :: Double
+    -> Double
+    -> Either ProbabilityError Double
+normalVar mu sigma
+    | sigma <= 0 =
+        Left (InvalidStandardDeviation sigma)
+
+    | otherwise =
+        Right (sigma * sigma)
+
+
+exponentialLogPDF
+    :: Double
+    -> Double
+    -> Either ProbabilityError Double
+exponentialLogPDF x lambda
+    | lambda <= 0 =
+        Left (InvalidRate lambda)
+
+    | x < 0 =
+        Right (-1.0 / 0.0)
+
+    | otherwise =
+        Right (log lambda - lambda * x)
+
+
+exponentialPDF
+    :: Double
+    -> Double
+    -> Either ProbabilityError Double
+exponentialPDF x lambda
+    | lambda <= 0 =
+        Left (InvalidRate lambda)
+
+    | x < 0 =
+        Right 0.0
+
+    | otherwise = do
+        logP <- exponentialLogPDF x lambda
+        pure (exp logP)
+
+
+exponentialCDF
+    :: Double
+    -> Double
+    -> Either ProbabilityError Double
+exponentialCDF x lambda
+    | lambda <= 0 =
+        Left (InvalidRate lambda)
+
+    | x < 0 =
+        Right 0.0
+
+    | otherwise =
+        Right (- expm1 (-lambda * x))
+
+
+exponentialMean
+    :: Double
+    -> Either ProbabilityError Double
+exponentialMean lambda
+    | lambda <= 0 =
+        Left (InvalidRate lambda)
+
+    | otherwise =
+        Right (1.0 / lambda)
+
+
+exponentialVar
+    :: Double
+    -> Either ProbabilityError Double
+exponentialVar lambda
+    | lambda <= 0 =
+        Left (InvalidRate lambda)
+
+    | otherwise =
+        Right (1.0 / (lambda * lambda))
+
+gammaLogPDF
+    :: Double
+    -> Double
+    -> Double
+    -> Either ProbabilityError Double
+gammaLogPDF x alpha lambda
+    | alpha <= 0 =
+        Left (InvalidShape alpha)
+
+    | lambda <= 0 =
+        Left (InvalidRate lambda)
+
+    | x < 0 =
+        Right (-1.0 / 0.0)
+
+    | x == 0 && alpha < 1 =
+        Right (1.0 / 0.0)
+
+    | x == 0 && alpha == 1 =
+        Right (log lambda)
+
+    | x == 0 =
+        Right (-1.0 / 0.0)
+
+    | otherwise = do
+        logG <- logGamma alpha
+
+        let logP =
+                alpha * log lambda
+                - logG
+                + (alpha - 1.0) * log x
+                - lambda * x
+
+        pure logP
+
+
+gammaPDF
+    :: Double
+    -> Double
+    -> Double
+    -> Either ProbabilityError Double
+gammaPDF x alpha lambda
+    | alpha <= 0 =
+        Left (InvalidShape alpha)
+
+    | lambda <= 0 =
+        Left (InvalidRate lambda)
+
+    | x < 0 =
+        Right 0.0
+
+    | x == 0 && alpha < 1 =
+        Right (1.0 / 0.0)
+
+    | x == 0 && alpha == 1 =
+        Right lambda
+
+    | x == 0 =
+        Right 0.0
+
+    | otherwise = do
+        logP <- gammaLogPDF x alpha lambda
+        pure (exp logP)
+
+
+gammaMean
+    :: Double
+    -> Double
+    -> Either ProbabilityError Double
+gammaMean alpha lambda
+    | alpha <= 0 =
+        Left (InvalidShape alpha)
+
+    | lambda <= 0 =
+        Left (InvalidRate lambda)
+
+    | otherwise =
+        Right (alpha / lambda)
+
+
+gammaVar
+    :: Double
+    -> Double
+    -> Either ProbabilityError Double
+gammaVar alpha lambda
+    | alpha <= 0 =
+        Left (InvalidShape alpha)
+
+    | lambda <= 0 =
+        Left (InvalidRate lambda)
+
+    | otherwise =
+        Right (alpha / (lambda * lambda))
+
+uniformLogPDF
+    :: Double
+    -> Double
+    -> Double
+    -> Either ProbabilityError Double
+uniformLogPDF x a b
+    | b <= a =
+        Left (InvalidBounds a b)
+
+    | x < a || x > b =
+        Right (-1.0 / 0.0)
+
+    | otherwise =
+        Right (- log (b - a))
+
+
+uniformPDF
+    :: Double
+    -> Double
+    -> Double
+    -> Either ProbabilityError Double
+uniformPDF x a b
+    | b <= a =
+        Left (InvalidBounds a b)
+
+    | x < a || x > b =
+        Right 0.0
+
+    | otherwise =
+        Right (1.0 / (b - a))
+
+
+uniformCDF
+    :: Double
+    -> Double
+    -> Double
+    -> Either ProbabilityError Double
+uniformCDF x a b
+    | b <= a =
+        Left (InvalidBounds a b)
+
+    | x < a =
+        Right 0.0
+
+    | x > b =
+        Right 1.0
+
+    | otherwise =
+        Right ((x - a) / (b - a))
+
+
+uniformMean
+    :: Double
+    -> Double
+    -> Either ProbabilityError Double
+uniformMean a b
+    | b <= a =
+        Left (InvalidBounds a b)
+
+    | otherwise =
+        Right ((a + b) / 2.0)
+
+
+uniformVar
+    :: Double
+    -> Double
+    -> Either ProbabilityError Double
+uniformVar a b
+    | b <= a =
+        Left (InvalidBounds a b)
+
+    | otherwise =
+        let width = b - a
+        in Right (width * width / 12.0)
diff --git a/src/Probability/Discrete.hs b/src/Probability/Discrete.hs
new file mode 100644
--- /dev/null
+++ b/src/Probability/Discrete.hs
@@ -0,0 +1,420 @@
+module Probability.Discrete
+    ( bernoulliPMF
+    , bernoulliMean
+    , bernoulliVar
+    ,binomialPMF
+    ,binomialLogPMF
+    ,binomialMean
+    ,binomialVar
+    ,poissonPMF
+    ,poissonLogPMF
+    ,poissonMean
+    ,poissonVar
+    ,negativeBinomialPMF
+    ,negativeBinomialLogPMF
+    ,negativeBinomialMean
+    ,negativeBinomialVar
+    ,geometricLogPMF
+    ,geometricPMF
+    ,geometricMean
+    ,geometricVar
+    ,hypergeometricPMF
+    ,hypergeometricLogPMF
+    ,hypergeometricMean
+    ,hypergeometricVar
+    ,multinomialLogPMF
+    ,multinomialPMF
+    ) where
+
+import Probability.Error (ProbabilityError(..))
+import Probability.Math
+    (logCombination
+    ,logFactorial)
+
+bernoulliPMF :: Int -> Double -> Either ProbabilityError Double
+bernoulliPMF x p
+    | p < 0 || p > 1 = Left (ProbabilityOutOfRange p)
+    | x == 0         = Right (1 - p)
+    | x == 1         = Right p
+    | otherwise      = Right 0.0
+
+
+bernoulliMean :: Double -> Either ProbabilityError Double
+bernoulliMean p
+    | p < 0 || p > 1 = Left (ProbabilityOutOfRange p)
+    | otherwise      = Right p
+
+
+bernoulliVar :: Double -> Either ProbabilityError Double
+bernoulliVar p
+    | p < 0 || p > 1 = Left (ProbabilityOutOfRange p)
+    | otherwise      = Right (p * (1 - p))
+
+
+
+binomialLogPMF
+    :: Int
+    -> Int
+    -> Double
+    -> Either ProbabilityError Double
+binomialLogPMF k n p
+    | n < 0 = Left (InvalidSampleSize n)
+    | p < 0 || p > 1 = Left (ProbabilityOutOfRange p)
+    | k < 0 || k > n = Right 0.0
+    | otherwise = do
+        logC <- logCombination n k
+
+        let termK =
+                if k == 0
+                    then 0.0
+                    else fromIntegral k * log p
+
+            termNMinusK =
+                if k == n
+                    then 0.0
+                    else fromIntegral (n - k) * log (1 - p)
+
+        return (logC + termK + termNMinusK)
+
+
+binomialPMF
+    :: Int
+    -> Int
+    -> Double
+    -> Either ProbabilityError Double
+binomialPMF k n p
+    | n < 0 = Left (InvalidSampleSize n)
+    | p < 0 || p > 1 = Left (ProbabilityOutOfRange p)
+    | k < 0 || k > n = Left(InvalidSuccessCount k)
+    | otherwise = do
+        logP <- binomialLogPMF k n p
+        pure (exp logP)
+
+binomialMean
+    :: Int
+    -> Double
+    -> Either ProbabilityError Double
+binomialMean n p
+    | n < 0 = Left (InvalidSampleSize n)
+    | p < 0 || p > 1 = Left (ProbabilityOutOfRange p)
+    | otherwise = Right (fromIntegral n * p)
+
+binomialVar
+    :: Int
+    -> Double
+    -> Either ProbabilityError Double
+binomialVar n p
+    | n < 0 = Left (InvalidSampleSize n)
+    | p < 0 || p > 1 = Left (ProbabilityOutOfRange p)
+    | otherwise = Right (fromIntegral n * p * (1 - p))
+
+
+poissonLogPMF
+    :: Int
+    -> Double
+    -> Either ProbabilityError Double
+poissonLogPMF k lambda
+    | lambda <= 0 = Left (InvalidRate lambda)
+    | k < 0       = Right 0.0
+    | otherwise   = do
+        logFact <- logFactorial k
+
+        let logP =
+                -lambda
+                + if k == 0
+                    then 0.0
+                    else fromIntegral k * log lambda
+                - logFact
+
+        pure logP
+
+
+poissonPMF
+    :: Int
+    -> Double
+    -> Either ProbabilityError Double
+poissonPMF k lambda
+    | lambda <= 0 = Left (InvalidRate lambda)
+    | k < 0       = Right 0.0
+    | otherwise   = do
+        logP <- poissonLogPMF k lambda
+        pure (exp logP)
+
+
+poissonMean
+    :: Double
+    -> Either ProbabilityError Double
+poissonMean lambda
+    | lambda <= 0 = Left (InvalidRate lambda)
+    | otherwise   = Right lambda
+
+
+poissonVar
+    :: Double
+    -> Either ProbabilityError Double
+poissonVar lambda
+    | lambda <= 0 = Left (InvalidRate lambda)
+    | otherwise   = Right lambda
+
+
+negativeBinomialLogPMF
+    :: Int
+    -> Int
+    -> Double
+    -> Either ProbabilityError Double
+negativeBinomialLogPMF k r p
+    | r <= 0 = Left (InvalidSuccessCount r)
+    | p <= 0 || p > 1 = Left (ProbabilityOutOfRange p)
+    | k < 0 = Left (InvalidFailureCount k)
+    | otherwise = do
+        logC <- logCombination (r + k - 1) k
+
+        let logP =
+                logC
+                + fromIntegral r * log p
+                + if k == 0
+                    then 0.0
+                    else fromIntegral k * log (1 - p)
+
+        pure logP
+
+
+negativeBinomialPMF
+    :: Int
+    -> Int
+    -> Double
+    -> Either ProbabilityError Double
+negativeBinomialPMF k r p
+    | r <= 0 = Left (InvalidSuccessCount r)
+    | p <= 0 || p > 1 = Left (ProbabilityOutOfRange p)
+    | k < 0 = Left (InvalidFailureCount k)
+    | otherwise = do
+        logP <- negativeBinomialLogPMF k r p
+        pure (exp logP)
+
+
+negativeBinomialMean
+    :: Int
+    -> Double
+    -> Either ProbabilityError Double
+negativeBinomialMean r p
+    | r <= 0 = Left (InvalidSuccessCount r)
+    | p <= 0 || p > 1 = Left (ProbabilityOutOfRange p)
+    | otherwise =
+        Right (fromIntegral r * (1 - p) / p)
+
+
+negativeBinomialVar
+    :: Int
+    -> Double
+    -> Either ProbabilityError Double
+negativeBinomialVar r p
+    | r <= 0 = Left (InvalidSuccessCount r)
+    | p <= 0 || p > 1 = Left (ProbabilityOutOfRange p)
+    | otherwise =
+        Right (fromIntegral r * (1 - p) / (p * p))
+
+
+geometricLogPMF
+    :: Int
+    -> Double
+    -> Either ProbabilityError Double
+geometricLogPMF k p =
+    negativeBinomialLogPMF k 1 p
+
+
+geometricPMF
+    :: Int
+    -> Double
+    -> Either ProbabilityError Double
+geometricPMF k p =
+    negativeBinomialPMF k 1 p
+
+
+geometricMean
+    :: Double
+    -> Either ProbabilityError Double
+geometricMean p =
+    negativeBinomialMean 1 p
+
+
+geometricVar
+    :: Double
+    -> Either ProbabilityError Double
+geometricVar p =
+    negativeBinomialVar 1 p
+
+
+hypergeometricLogPMF
+    :: Int
+    -> Int
+    -> Int
+    -> Int
+    -> Either ProbabilityError Double
+hypergeometricLogPMF k sampleSize successes populationSize
+    | populationSize <= 0 =
+        Left (InvalidPopulation populationSize) --size
+
+    | successes < 0 || successes > populationSize =
+        Left (InvalidSuccessCount successes)
+
+    | sampleSize < 0 || sampleSize > populationSize =
+        Left (InvalidSampleSize sampleSize)
+
+    | k < 0
+        || k > sampleSize
+        || k > successes
+        || sampleSize - k > populationSize - successes =
+        Right 0.0
+
+    | otherwise = do
+        logC1 <- logCombination successes k
+        logC2 <- logCombination
+                    (populationSize - successes)
+                    (sampleSize - k)
+        logC3 <- logCombination populationSize sampleSize
+
+        pure (logC1 + logC2 - logC3)
+
+
+hypergeometricPMF
+    :: Int
+    -> Int
+    -> Int
+    -> Int
+    -> Either ProbabilityError Double
+hypergeometricPMF k sampleSize successes populationSize
+    | populationSize <= 0 =
+        Left (InvalidPopulation populationSize) --size
+
+    | successes < 0 || successes > populationSize =
+        Left (InvalidSuccessCount successes)
+
+    | sampleSize < 0 || sampleSize > populationSize =
+        Left (InvalidSampleSize sampleSize)
+
+    | k < 0
+        || k > sampleSize
+        || k > successes
+        || sampleSize - k > populationSize - successes =
+        Right 0.0
+
+    | otherwise = do
+        logP <- hypergeometricLogPMF
+                    k
+                    sampleSize
+                    successes
+                    populationSize
+
+        pure (exp logP)
+
+
+hypergeometricMean
+    :: Int
+    -> Int
+    -> Int
+    -> Either ProbabilityError Double
+hypergeometricMean sampleSize successes populationSize
+    | populationSize <= 0 =
+        Left (InvalidPopulation populationSize) -- size kısımını düzelttik
+
+    | successes < 0 || successes > populationSize =
+        Left (InvalidSuccessCount successes)
+
+    | sampleSize < 0 || sampleSize > populationSize =
+        Left (InvalidSampleSize sampleSize)
+
+    | otherwise =
+        Right
+            ( fromIntegral sampleSize
+            * fromIntegral successes
+            / fromIntegral populationSize
+            )
+
+
+hypergeometricVar
+    :: Int
+    -> Int
+    -> Int
+    -> Either ProbabilityError Double
+hypergeometricVar sampleSize successes populationSize
+    | populationSize <= 1 =
+        Left (InvalidPopulation populationSize)   -- düzeltme yaptık
+
+    | successes < 0 || successes > populationSize =
+        Left (InvalidSuccessCount successes)
+
+    | sampleSize < 0 || sampleSize > populationSize =
+        Left (InvalidSampleSize sampleSize)
+
+    | otherwise =
+        let n = fromIntegral sampleSize
+            k = fromIntegral successes
+            nPop = fromIntegral populationSize
+
+            p = k / nPop
+
+        in Right
+            ( n
+            * p
+            * (1 - p)
+            * ((nPop - n) / (nPop - 1))
+            )
+
+multinomialLogPMF
+    :: Int
+    -> [Int]
+    -> [Double]
+    -> Either ProbabilityError Double
+multinomialLogPMF n counts probabilities
+    | n < 0 =
+        Left (InvalidSampleSize n)
+
+    | null probabilities =
+        Left (InvalidProbabilityVector probabilities)
+
+    | length counts /= length probabilities =
+        Left (DimensionMismatch
+                (length counts)
+                (length probabilities))
+
+    | any (< 0) counts =
+        Left (InvalidCountVector counts)
+
+    | sum counts /= n =
+        Left (InvalidSampleSize n)
+
+    | any (< 0) probabilities =
+        Left (InvalidProbabilityVector probabilities)
+
+    | abs (sum probabilities - 1.0) > 1.0e-12 =
+        Left (InvalidProbabilityVector probabilities)
+
+    | otherwise = do
+        logNFact <- logFactorial n
+
+        logCountsFact <- mapM logFactorial counts
+
+        let logCountTerm =
+                logNFact - sum logCountsFact
+
+            logProbabilityTerm =
+                sum
+                    [ if k == 0
+                        then 0.0
+                        else fromIntegral k * log p
+                    | (k, p) <- zip counts probabilities
+                    ]
+
+        pure (logCountTerm + logProbabilityTerm)
+
+
+multinomialPMF
+    :: Int
+    -> [Int]
+    -> [Double]
+    -> Either ProbabilityError Double
+multinomialPMF n counts probabilities
+    | otherwise = do
+        logP <- multinomialLogPMF n counts probabilities
+        pure (exp logP)
diff --git a/src/Probability/Error.hs b/src/Probability/Error.hs
new file mode 100644
--- /dev/null
+++ b/src/Probability/Error.hs
@@ -0,0 +1,19 @@
+module Probability.Error
+(ProbabilityError(..)) where
+
+data ProbabilityError
+    = ProbabilityOutOfRange Double
+    | NegativeValue Double
+    | InvalidBounds Double Double
+    | InvalidSampleSize Int
+    | InvalidSuccessCount Int
+    | InvalidFailureCount Int
+    | InvalidRate Double
+    | InvalidScale Double
+    | InvalidShape Double
+    | InvalidPopulation Int
+    |InvalidProbabilityVector [Double]
+    |InvalidCountVector [Int]
+    |DimensionMismatch Int Int
+    |InvalidStandardDeviation Double
+    deriving(Show, Eq)
diff --git a/src/Probability/Math.hs b/src/Probability/Math.hs
new file mode 100644
--- /dev/null
+++ b/src/Probability/Math.hs
@@ -0,0 +1,128 @@
+module Probability.Math
+    ( logFactorial
+    , logCombination
+    , logGamma
+    , logBeta
+    ,erf
+    ,expm1
+    ) where
+
+import Probability.Error (ProbabilityError(..))
+
+
+logFactorial :: Int -> Either ProbabilityError Double
+logFactorial n
+    | n < 0  = Left (NegativeValue (fromIntegral n))
+    | n < 35 = Right (sum [log (fromIntegral i) | i <- [1 .. n]])
+    | otherwise = logGamma (fromIntegral n + 1)
+
+
+logCombination :: Int -> Int -> Either ProbabilityError Double
+logCombination n r
+    | n < 0 = Left (NegativeValue (fromIntegral n))
+    | r < 0 = Left (NegativeValue (fromIntegral r))
+    | r > n = Left (InvalidSuccessCount r)
+    | otherwise = do
+        lnN  <- logFactorial n
+        lnR  <- logFactorial r
+        lnNR <- logFactorial (n - r)
+
+        pure (lnN - lnR - lnNR)
+
+
+logGamma :: Double -> Either ProbabilityError Double
+logGamma z
+    | z <= 0 = Left (NegativeValue z)
+    | otherwise = Right result
+  where
+    g :: Double
+    g = 7.0
+
+    coefficients :: [Double]
+    coefficients =
+        [ 0.99999999999980993
+        , 676.5203681218851
+        , -1259.1392167224028
+        , 771.32342877765313
+        , -176.61502916214059
+        , 12.507343278686905
+        , -0.13857109526572012
+        , 9.9843695780195716e-6
+        , 1.5056327351493116e-7
+        ]
+
+    x :: Double
+    x = z - 1.0
+
+    a :: Double
+    a =
+        foldl
+            (\acc (i, c) ->
+                acc + c / (x + fromIntegral i))
+            (head coefficients)
+            (zip [1 ..] (tail coefficients))
+
+    t :: Double
+    t = x + g + 0.5
+
+    result :: Double
+    result =
+        0.5 * log (2.0 * pi)
+            + (x + 0.5) * log t
+            - t
+            + log a
+
+
+logBeta :: Double -> Double -> Either ProbabilityError Double
+logBeta a b
+    | a <= 0 = Left (InvalidShape a)
+    | b <= 0 = Left (InvalidShape b)
+    | otherwise = do
+        logA   <- logGamma a
+        logB   <- logGamma b
+        logAB  <- logGamma (a + b)
+
+        pure (logA + logB - logAB)
+
+erf
+    :: Double
+    -> Double
+erf x
+    | x == 0.0 = 0.0
+    | otherwise =
+        let sign = if x < 0.0 then -1.0 else 1.0
+            ax   = abs x
+
+            p    = 0.3275911
+            t    = 1.0 / (1.0 + p * ax)
+
+            poly =
+                t * exp
+                    ( -ax * ax
+                        -1.26551223
+                        + t * ( 1.00002368
+                        + t * ( 0.37409196
+                        + t * ( 0.09678418
+                        + t * (-0.18628806
+                        + t * (0.27886807
+                        + t * (-1.13520398
+                        + t * (1.48851587
+                        + t * (-0.82215223
+                        + t * 0.17087277)))))))))
+
+        in sign * (1.0 - poly)
+
+expm1
+    :: Double
+    -> Double
+expm1 x
+    | abs x < 1.0e-5 =
+        x
+        + x^2 / 2.0
+        + x^3 / 6.0
+        + x^4 / 24.0
+        + x^5 / 120.0
+        + x^6 / 720.0
+
+    | otherwise =
+        exp x - 1.0
diff --git a/test/Main.hs b/test/Main.hs
new file mode 100644
--- /dev/null
+++ b/test/Main.hs
@@ -0,0 +1,393 @@
+module Main where
+
+import Control.Monad (unless)
+import Data.Either (isLeft)
+import System.Exit (exitFailure)
+
+--import Probability.Math
+import Probability.Discrete
+import Probability.Continuous
+
+
+-- ============================================================
+-- Test helpers
+-- ============================================================
+
+epsilon :: Double
+epsilon = 1.0e-12
+
+
+approxEqual :: Double -> Double -> Bool
+approxEqual x y =
+    abs (x - y) < epsilon
+
+
+assert :: Bool -> String -> IO ()
+assert condition message =
+    unless condition $ do
+        putStrLn ("FAIL: " ++ message)
+        exitFailure
+
+
+assertApprox :: Double -> Double -> String -> IO ()
+assertApprox actual expected message =
+    assert (approxEqual actual expected) message
+
+
+assertRightApprox
+    :: Either e Double
+    -> Double
+    -> String
+    -> IO ()
+assertRightApprox result expected message =
+    case result of
+        Right value ->
+            assertApprox value expected message
+        Left _ ->
+            assert False (message ++ " returned Left")
+
+
+assertLeft
+    :: Either e a
+    -> String
+    -> IO ()
+assertLeft result message =
+    assert (isLeft result) message
+
+{-
+-- ============================================================
+-- Math
+-- ============================================================
+
+testMath :: IO ()
+testMath = do
+    putStrLn "Testing Math..."
+
+    -- logGamma
+    assertApprox
+        (logGamma 4)
+        1.791759469228055
+        "logGamma 4"
+
+    assertApprox
+        (logGamma 5)
+        3.178053830347944
+        "logGamma 5"
+
+    assertApprox
+        (logGamma 10)
+        12.801827480081474
+        "logGamma 10"
+
+    -- logFactorial
+    assertApprox
+        (logFactorial 0)
+        0.0
+        "logFactorial 0"
+
+    assertApprox
+        (logFactorial 10)
+        15.104412573075516
+        "logFactorial 10"
+
+    assertApprox
+        (logFactorial 40)
+        110.32063971475739
+        "logFactorial 40"
+
+    -- logCombination
+    assertRightApprox
+        (logCombination 10 0)
+        0.0
+        "logCombination 10 0"
+
+    assertRightApprox
+        (logCombination 10 10)
+        0.0
+        "logCombination 10 10"
+
+    assertLeft
+        (logCombination 10 11)
+        "logCombination 10 11"
+
+    assertLeft
+        (logCombination (-1) 2)
+        "logCombination -1 2"
+
+    -- logBeta
+    assertApprox
+        (logBeta 2 2)
+        (-1.791759469228055)
+        "logBeta 2 2"
+
+    assertApprox
+        (logBeta 1 2)
+        (-0.6931471805599472)
+        "logBeta 1 2"
+
+    -- erf
+    assertApprox
+        (erf 0.0)
+        0.0
+        "erf 0"
+
+    -- expm1
+    assertApprox
+        (expm1 0.0)
+        0.0
+        "expm1 0"
+
+    assertApprox
+        (expm1 1.0)
+        (exp 1.0 - 1.0)
+        "expm1 1"
+
+    putStrLn "Math: OK"
+
+-}
+-- ============================================================
+-- Discrete
+-- ============================================================
+
+testDiscrete :: IO ()
+testDiscrete = do
+    putStrLn "Testing Discrete..."
+
+    -- --------------------------------------------------------
+    -- Binomial
+    -- --------------------------------------------------------
+
+    assertRightApprox
+        (binomialPMF 0 10 0.5)
+        0.0009765625
+        "binomialPMF 10 0 0.5"
+
+    assertRightApprox
+        (binomialPMF 10 10 0.5)
+        0.0009765625
+        "binomialPMF 10 10 0.5"
+
+    assertLeft
+        (binomialPMF 11 10 0.5)
+        "binomialPMF invalid success count"
+
+    assertLeft
+        (binomialPMF (-1) 2 0.5)
+        "binomialPMF negative population"
+
+    -- --------------------------------------------------------
+    -- Geometric
+    -- --------------------------------------------------------
+
+    assertRightApprox
+        (geometricPMF 1 0.5)
+        0.25
+        "geometricPMF 1 0.5"
+
+    -- --------------------------------------------------------
+    -- Negative Binomial
+    -- --------------------------------------------------------
+
+    assertRightApprox
+        (negativeBinomialPMF 0 1 0.5)
+        0.5
+        "negativeBinomialPMF r=1"
+
+    -- --------------------------------------------------------
+    -- Multinomial
+    -- --------------------------------------------------------
+
+    assertRightApprox
+        (multinomialPMF
+            10
+            [10, 0]
+            [1.0, 0.0])
+        1.0
+        "multinomialPMF degenerate case"
+
+    -- --------------------------------------------------------
+    -- Poisson
+    -- --------------------------------------------------------
+
+    assertRightApprox
+        (poissonPMF 0 2.0)
+        (exp (-2.0))
+        "poissonPMF 0 2"
+
+    -- --------------------------------------------------------
+    -- Bernoulli
+    -- --------------------------------------------------------
+
+    assertRightApprox
+        (bernoulliPMF 1 0.7)
+        0.7
+        "bernoulliPMF success"
+
+    assertRightApprox
+        (bernoulliPMF 0 0.7)
+        0.3
+        "bernoulliPMF failure"
+
+    putStrLn "Discrete: OK"
+
+
+-- ============================================================
+-- Continuous
+-- ============================================================
+
+testContinuous :: IO ()
+testContinuous = do
+    putStrLn "Testing Continuous..."
+
+    -- --------------------------------------------------------
+    -- Normal
+    -- --------------------------------------------------------
+
+    assertRightApprox
+        (normalPDF 0.0 0.0 1.0)
+        0.3989422804014327
+        "normalPDF standard normal at 0"
+
+    assertRightApprox
+        (normalCDF 0.0 0.0 1.0)
+        0.5
+        "normalCDF standard normal at 0"
+
+    -- Symmetry
+    let leftCDF  = normalCDF (-1.0) 0.0 1.0
+        rightCDF = normalCDF 1.0 0.0 1.0
+
+    case (leftCDF, rightCDF) of
+        (Right l, Right r) ->
+            assertApprox
+                (l + r)
+                1.0
+                "normalCDF symmetry"
+        _ ->
+            assert False "normalCDF symmetry returned Left"
+
+    -- Invalid standard deviation
+    assertLeft
+        (normalPDF 0.0 0.0 (-1.0))
+        "normalPDF invalid standard deviation"
+
+    -- --------------------------------------------------------
+    -- Exponential
+    -- --------------------------------------------------------
+
+    assertRightApprox
+        (exponentialPDF 0.0 2.0)
+        2.0
+        "exponentialPDF x=0"
+
+    assertRightApprox
+        (exponentialCDF 0.0 2.0)
+        0.0
+        "exponentialCDF x=0"
+
+    assertRightApprox
+        (exponentialCDF 1.0 2.0)
+        (1.0 - exp (-2.0))
+        "exponentialCDF x=1 lambda=2"
+
+    assertLeft
+        (exponentialPDF 1.0 (-1.0))
+        "exponentialPDF invalid lambda"
+
+    -- --------------------------------------------------------
+    -- Gamma
+    -- --------------------------------------------------------
+
+    -- Gamma(1, lambda) = Exponential(lambda)
+    assertRightApprox
+        (gammaPDF 0.0 1.0 2.0)
+        2.0
+        "gammaPDF alpha=1 x=0"
+
+    assertRightApprox
+        (gammaPDF 1.0 2.0 1.0)
+        0.36787944117144233
+        "gammaPDF alpha=2 lambda=1 x=1"
+
+    assertRightApprox
+        (gammaMean 2.0 1.0)
+        2.0
+        "gammaMean"
+
+    assertRightApprox
+        (gammaVar 2.0 1.0)
+        2.0
+        "gammaVar"
+
+    assertLeft
+        (gammaPDF 1.0 (-1.0) 1.0)
+        "gammaPDF invalid shape"
+
+    assertLeft
+        (gammaPDF 1.0 1.0 (-1.0))
+        "gammaPDF invalid rate"
+
+    -- --------------------------------------------------------
+    -- Uniform
+    -- --------------------------------------------------------
+
+    assertRightApprox
+        (uniformPDF 0.5 0.0 1.0)
+        1.0
+        "uniformPDF U(0,1)"
+
+    assertRightApprox
+        (uniformCDF 0.5 0.0 1.0)
+        0.5
+        "uniformCDF U(0,1)"
+
+    assertRightApprox
+        (uniformMean 0.0 1.0)
+        0.5
+        "uniformMean U(0,1)"
+
+    assertRightApprox
+        (uniformVar 0.0 1.0)
+        (1.0 / 12.0)
+        "uniformVar U(0,1)"
+
+    -- CDF boundaries
+    assertRightApprox
+        (uniformCDF (-1.0) 0.0 1.0)
+        0.0
+        "uniformCDF below lower bound"
+
+    assertRightApprox
+        (uniformCDF 2.0 0.0 1.0)
+        1.0
+        "uniformCDF above upper bound"
+
+    -- Invalid bounds
+    assertLeft
+        (uniformPDF 0.0 1.0 1.0)
+        "uniformPDF invalid bounds"
+
+    assertLeft
+        (uniformCDF 0.0 2.0 1.0)
+        "uniformCDF invalid bounds"
+
+    putStrLn "Continuous: OK"
+
+
+-- ============================================================
+-- Main
+-- ============================================================
+
+main :: IO ()
+main = do
+    putStrLn "========================================"
+    putStrLn " probability-dist test suite"
+    putStrLn "========================================"
+
+    --testMath
+    testDiscrete
+    testContinuous
+
+    putStrLn "========================================"
+    putStrLn " All tests passed."
+    putStrLn "========================================"