diff --git a/integration_tests/run_tests.py b/integration_tests/run_tests.py
index aa70e4bd27..d144ea8d5a 100755
--- a/integration_tests/run_tests.py
+++ b/integration_tests/run_tests.py
@@ -16,6 +16,7 @@
     "modules_01.py",
     #"modules_02.py",
     "test_math.py",
+    "test_random.py",
     "test_builtin.py",
     "test_builtin_abs.py",
     "test_builtin_bool.py",
diff --git a/integration_tests/test_random.py b/integration_tests/test_random.py
new file mode 100644
index 0000000000..2ad5e15a54
--- /dev/null
+++ b/integration_tests/test_random.py
@@ -0,0 +1,43 @@
+from ltypes import i32, f64
+import random
+
+
+def test_random():
+    r: f64
+    r = random.random()
+    print(r)
+    assert r >= 0.0 and r < 1.0
+    r = random.random()
+    print(r)
+    assert r >= 0.0 and r < 1.0
+
+def test_randrange():
+    r: i32
+    r = random.randrange(0, 10) # [0, 10)
+    print(r)
+    assert r >= 0 and r < 10
+    r = random.randrange(-50, 76) # [-50, 76)
+    print(r)
+    assert r >= -50 and r < 76
+
+def test_randint():
+    ri1: i32
+    ri2: i32
+    ri1 = random.randint(0, 10) # [0, 10]
+    print(ri1)
+    assert ri1 >= 0 and ri1 <= 10
+    ri2 = random.randint(-50, 76) # [-50, 76]
+    print(ri2)
+    assert ri2 >= -50 and ri2 <= 76
+
+def test_paretovariate():
+    r: f64
+    r = random.paretovariate(2.0)
+    print(r)
+    r = random.paretovariate(-5.6)
+    print(r)
+
+test_random()
+test_randrange()
+test_randint()
+test_paretovariate()
diff --git a/src/runtime/impure/lfortran_intrinsics.c b/src/runtime/impure/lfortran_intrinsics.c
index 8ab63d6961..c228ebde27 100644
--- a/src/runtime/impure/lfortran_intrinsics.c
+++ b/src/runtime/impure/lfortran_intrinsics.c
@@ -30,6 +30,23 @@ LFORTRAN_API void _lfortran_random_number(int n, double *v)
     }
 }
 
+LFORTRAN_API float _lfortran_random_float()
+{
+    return ((float) rand() / (float) RAND_MAX);
+}
+
+LFORTRAN_API int _lfortran_randrange(int lower, int upper)
+{
+    int rr = lower + (rand() % (upper - lower));
+    return rr;
+}
+
+LFORTRAN_API int _lfortran_random_int(int lower, int upper)
+{
+    int randint = lower + (rand() % (upper - lower + 1));
+    return randint;
+}
+
 LFORTRAN_API void _lfortran_printf(const char* format, ...)
 {
     va_list args;
diff --git a/src/runtime/impure/lfortran_intrinsics.h b/src/runtime/impure/lfortran_intrinsics.h
index a72f7eb1a4..aff7f50664 100644
--- a/src/runtime/impure/lfortran_intrinsics.h
+++ b/src/runtime/impure/lfortran_intrinsics.h
@@ -31,6 +31,9 @@ typedef double _Complex double_complex_t;
 
 LFORTRAN_API double _lfortran_sum(int n, double *v);
 LFORTRAN_API void _lfortran_random_number(int n, double *v);
+LFORTRAN_API float _lfortran_random_float();
+LFORTRAN_API int _lfortran_randrange(int lower, int upper);
+LFORTRAN_API int _lfortran_random_int(int lower, int upper);
 LFORTRAN_API void _lfortran_printf(const char* format, ...);
 LFORTRAN_API void _lfortran_complex_add(struct _lfortran_complex* a,
         struct _lfortran_complex* b, struct _lfortran_complex *result);
diff --git a/src/runtime/random.py b/src/runtime/random.py
new file mode 100644
index 0000000000..2b8da20ed4
--- /dev/null
+++ b/src/runtime/random.py
@@ -0,0 +1,40 @@
+from ltypes import i32, f64, ccall
+
+
+def random() -> f64:
+    """
+    Returns a random floating point number in the range [0.0, 1.0)
+    """
+    return _lfortran_random_float()
+
+@ccall
+def _lfortran_random_float() -> f64:
+    pass
+
+def randrange(lower: i32, upper: i32) -> i32:
+    """
+    Return a random integer N such that `lower <= N < upper`.
+    """
+    return _lfortran_randrange(lower, upper)
+
+@ccall
+def _lfortran_randrange(lower: i32, upper: i32) -> i32:
+    pass
+
+def randint(lower: i32, upper: i32) -> i32:
+    """
+    Return a random integer N such that `lower <= N <= upper`.
+    """
+    return _lfortran_random_int(lower, upper)
+
+@ccall
+def _lfortran_random_int(lower: i32, upper: i32) -> i32:
+    pass
+
+def paretovariate(alpha: f64) -> f64:
+    """
+    Return a random number from a Pareto distribution with parameter `alpha`.
+    """
+    u: f64
+    u = 1.0 - random()
+    return u ** (-1.0 / alpha)