5 files changed, 146 insertions, 4 deletions
diff --git a/challenge-167/paulo-custodio/perl/ch-1.pl b/challenge-167/paulo-custodio/perl/ch-1.pl
index 6ee40c97f8..0b61d7bcf4 100644
--- a/challenge-167/paulo-custodio/perl/ch-1.pl
+++ b/challenge-167/paulo-custodio/perl/ch-1.pl
@@ -1,4 +1,4 @@
-#!/usr/bin/perl
+#!/usr/bin/env perl
 
 # Challenge 167
 #
diff --git a/challenge-167/paulo-custodio/perl/ch-2.pl b/challenge-167/paulo-custodio/perl/ch-2.pl
index 65e303772e..2a804b6757 100644
--- a/challenge-167/paulo-custodio/perl/ch-2.pl
+++ b/challenge-167/paulo-custodio/perl/ch-2.pl
@@ -1,4 +1,4 @@
-#!/usr/bin/perl
+#!/usr/bin/env perl
 
 # Challenge 167
 #
@@ -65,4 +65,4 @@ sub gamma {
     return drop_imag($y);
 }
 
-say gamma(Math::Complex->new(shift));
+say sprintf("%.0f", gamma(Math::Complex->new(shift)));
diff --git a/challenge-167/paulo-custodio/python/ch-1.py b/challenge-167/paulo-custodio/python/ch-1.py
new file mode 100644
index 0000000000..d32a0d33ba
--- /dev/null
+++ b/challenge-167/paulo-custodio/python/ch-1.py
@@ -0,0 +1,81 @@
+#!/usr/bin/env python3
+
+# Challenge 167
+#
+# Task 1: Circular Prime
+# Submitted by: Mohammad S Anwar
+#
+# Write a script to find out first 10 circular primes having at least 3 digits
+# (base 10).
+#
+# Please checkout wikipedia for more information.
+#
+#     A circular prime is a prime number with the property that the number
+#     generated at each intermediate step when cyclically permuting its (base 10)
+#     digits will also be prime.
+#
+# Output
+#
+# 113, 197, 199, 337, 1193, 3779, 11939, 19937, 193939, 199933
+
+def is_prime(n):
+    if n <= 1:
+        return False
+    if n <= 3:
+        return True
+    if n % 2 == 0 or n % 3 == 0:
+        return False
+    i = 5
+    while i * i <= n:
+        if n % i == 0 or n % (i + 2) == 0:
+            return False
+        i += 6
+    return True
+
+def next_prime(n):
+    if n <= 1:
+        return 2
+    p = n
+    while not is_prime(p := p + 1):
+        pass
+    return p
+
+seen = {}
+def is_circular_prime(n):
+    global seen
+
+    if n in seen:
+        return False
+    seen[n] = True
+    if n < 10:
+        return False
+    if any(digit in '024568' for digit in str(n)):
+        return False
+    rotations = []
+    for _ in range(len(str(n))):
+        rotations.append(n)
+        if not is_prime(n):
+            return False
+        n = int(str(n)[1:] + str(n)[0])
+    for rotation in rotations:
+        seen[rotation] = True
+    return True
+
+def next_circular_prime(n):
+    n += 1
+    while not is_circular_prime(n):
+        n = next_prime(n)
+    return n
+
+import sys
+
+if len(sys.argv) != 2:
+    raise ValueError("usage: ch-1.py n")
+n = int(sys.argv[1])
+cprimes = []
+cp = 1
+while len(cprimes) < n:
+    cp = next_circular_prime(cp)
+    if cp >= 100:
+        cprimes.append(cp)
+print(", ".join(map(str, cprimes)))
diff --git a/challenge-167/paulo-custodio/python/ch-2.py b/challenge-167/paulo-custodio/python/ch-2.py
new file mode 100644
index 0000000000..1a8fc16a90
--- /dev/null
+++ b/challenge-167/paulo-custodio/python/ch-2.py
@@ -0,0 +1,61 @@
+#!/usr/bin/env perl
+
+# Challenge 167
+#
+# Task 2: Gamma Function
+# Submitted by: Mohammad S Anwar
+#
+# Implement subroutine gamma() using the Lanczos approximation method.
+#
+# [2022-05-31]
+#
+# Ryan Thompson wrote an interesting blog explaining the subject in details.
+# Highly recommended if you are looking for more information. BLOG.
+#
+# Example
+#
+# print gamma(3); # 1.99
+# print gamma(5); # 24
+# print gamma(7); # 719.99
+
+import sys
+import cmath
+import math
+
+# Constants
+g = 7
+p = [
+    0.99999999999980993227684700473478,
+    676.520368121885098567009190444019,
+    -1259.13921672240287047156078755283,
+    771.3234287776530788486528258894,
+    -176.61502916214059906584551354,
+    12.507343278686904814458936853,
+    -0.13857109526572011689554707,
+    9.984369578019570859563e-6,
+    1.50563273514931155834e-7
+]
+
+EPSILON = 1e-07
+
+def drop_imag(z):
+    if abs(z.imag) <= EPSILON:
+        z = z.real
+    return z
+
+def gamma(z):
+    if z.real < 0.5:
+        y = math.pi / (math.sin(math.pi * z) * gamma(1 - z))  # Reflection formula
+    else:
+        z -= 1
+        x = p[0]
+        for i in range(1, len(p)):
+            x += p[i] / (z + i)
+        t = z + g + 0.5
+        t = drop_imag(t)
+        y = math.sqrt(2 * math.pi) * t ** (z + 0.5) * math.exp(-t) * x
+    return drop_imag(y)
+
+# Example usage
+z = complex(float(sys.argv[1]))
+print(f"{gamma(z):.0f}")
diff --git a/challenge-167/paulo-custodio/t/test-2.yaml b/challenge-167/paulo-custodio/t/test-2.yaml
index e575042550..0a1a933c29 100644
--- a/challenge-167/paulo-custodio/t/test-2.yaml
+++ b/challenge-167/paulo-custodio/t/test-2.yaml
@@ -12,4 +12,4 @@
   cleanup:
   args:     7
   input:
-  output:   720.000000000001
+  output:   720