P108, P109: Added Python and Mathematica solutions.

P109, P119: Added Python solutions, tweaked Java solutions for clarity.
Readme: Updated solution count.

Author: nayuki

Readme.markdown

@@ -9,7 +9,7 @@ mathematical explanation/proof in the source code to justify the program logic.
 
 All problems from #1 to #100 have a Java program, and #1 to #50 have Python and Mathematica programs.
 The project codebase contains at least 185 solutions in Java, at least 135 in Python,
-at least 115 in Mathematica, and at least 65 in Haskell.
+at least 120 in Mathematica, and at least 65 in Haskell.
 
 Java solutions require JDK 5 or above. Python solutions are tested to work on CPython 2.7.10 and 3.4.3. Mathematica solutions are tested to work on Mathematica 5.1.
 

java/p109.java

@@ -19,21 +19,21 @@ public static void main(String[] args) {
 
 
 	public String run() {
-		// Initialization
-		points = new ArrayList<Integer>();  // Orderless, but duplicates are important
+		// Both lists are orderless but duplicates are important; they are sort of like multisets
+		points = new ArrayList<Integer>();
 		for (int i = 1; i <= 20; i++) {
-			points.add(i * 1);
-			points.add(i * 2);
-			points.add(i * 3);
+			for (int j = 1; j <= 3; j++)
+				points.add(i * j);
 		}
 		points.add(25);
 		points.add(50);
 
-		doublePoints = new ArrayList<Integer>();  // Orderless, but duplicates are important
+		List<Integer> doublePoints = new ArrayList<Integer>();  // Orderless
 		for (int i = 1; i <= 20; i++)
 			doublePoints.add(i * 2);
 		doublePoints.add(25 * 2);
 
+		// Memoization array
 		ways = new int[3][101][points.size()];
 		for (int[][] x : ways) {
 			for (int[] y : x)
@@ -54,7 +54,6 @@ public String run() {
 
 
 	private List<Integer> points;
-	private List<Integer> doublePoints;
 
 	private int[][][] ways;
 

java/p119.java

@@ -8,8 +8,8 @@
 
 import java.math.BigInteger;
 import java.util.ArrayList;
-import java.util.Collections;
-import java.util.List;
+import java.util.SortedSet;
+import java.util.TreeSet;
 
 
 public final class p119 implements EulerSolution {
@@ -33,23 +33,18 @@ public static void main(String[] args) {
 
 	public String run() {
 		for (BigInteger limit = BigInteger.ONE; ; limit = limit.shiftLeft(8)) {
-			List<BigInteger> candidates = new ArrayList<BigInteger>();
-			for (int k = 2; BigInteger.valueOf(2).pow(k).compareTo(limit) < 0; k++) {
-				BigInteger n = BigInteger.valueOf(2);
-				while (true) {
-					BigInteger pow = n.pow(k);
-					if (pow.compareTo(limit) >= 0 && BigInteger.valueOf(pow.toString().length() * 9).compareTo(n) < 0)
+			SortedSet<BigInteger> candidates = new TreeSet<BigInteger>();
+			for (int k = 2; BigInteger.valueOf(1).shiftLeft(k).compareTo(limit) < 0; k++) {
+				for (int n = 2; ; n++) {
+					BigInteger pow = BigInteger.valueOf(n).pow(k);
+					if (pow.compareTo(limit) >= 0 && pow.toString().length() * 9 < n)
 						break;
-					
-					if (pow.compareTo(BigInteger.TEN) >= 0 && !candidates.contains(pow) && isDigitSumPower(pow))
+					if (pow.compareTo(BigInteger.TEN) >= 0 && isDigitSumPower(pow))
 						candidates.add(pow);
-					n = n.add(BigInteger.ONE);
 				}
 			}
-			if (candidates.size() >= INDEX) {
-				Collections.sort(candidates);
-				return candidates.get(INDEX - 1).toString();
-			}
+			if (candidates.size() >= INDEX)
+				return new ArrayList<BigInteger>(candidates).get(INDEX - 1).toString();
 		}
 	}
 
@@ -64,7 +59,7 @@ private static boolean isDigitSumPower(BigInteger x) {
 		BigInteger pow = base;
 		while (pow.compareTo(x) < 0)
 			pow = pow.multiply(base);
-		return pow.compareTo(x) == 0;
+		return pow.equals(x);
 	}
 
 

mathematica/p108.mathematica

@@ -0,0 +1,27 @@
+(* 
+ * Solution to Project Euler problem 108
+ * by Project Nayuki
+ * 
+ * https://www.nayuki.io/page/project-euler-solutions
+ * https://github.com/nayuki/Project-Euler-solutions
+ *)
+
+
+(* 
+ * Rewrite the equation with x = n+i, y = n+j, and manipulate it:
+ *   1/n = 1/x + 1/y
+ *       = 1/(n+i) + 1/(n+j)
+ *       = (2n+i+j) / ((n+i)(n+j)).
+ *   n(2n+i+j) = (n+i)(n+j).
+ *   2n^2 + ni + nj = n^2 + ni + nj + ij.
+ *   n^2 = ij.
+ * Hence i and j are divisors of n^2. To ensure unique solutions,
+ * we impose that x <= y, so i <= j. Also, i > 0, otherwise no j exists.
+ * We have i <= j = n^2 / i, thus i^2 <= n^2. With i being positive, we get that i <= n.
+ * Therefore the number of solutions for i is the number of divisors of n^2 in the range [1, n].
+ * n^2 always has an odd number of divisors. One of them is n. As for the remainder of them, half of them are below n
+ * and half of them are above n. So if n^2 has m divisors, then we want (m+1)/2 of them as solutions for i.
+ *)
+
+For[n = 1, (DivisorSigma[0, n^2] + 1) / 2 <= 1000, n++]
+n

mathematica/p109.mathematica

@@ -0,0 +1,26 @@
+(* 
+ * Solution to Project Euler problem 109
+ * by Project Nayuki
+ * 
+ * https://www.nayuki.io/page/project-euler-solutions
+ * https://github.com/nayuki/Project-Euler-solutions
+ *)
+
+
+(* Both lists are orderless but duplicates are important; they are sort of like multisets *)
+points = Join[Flatten[Table[i * j, {i, 20}, {j, 3}]], {25, 50}];
+doublePoints = Join[Table[i * 2, {i, 20}], {25 * 2}];
+
+Ways[throws_, total_, maxIndex_] := Ways[throws, total, maxIndex] =  (* Memoization *)
+  If[throws == 0,
+    Boole[total == 0],
+    If[maxIndex > 1, Ways[throws, total, maxIndex - 1], 0] +
+      If[points[[maxIndex]] <= total, Ways[throws - 1, total - points[[maxIndex]], maxIndex], 0]]
+
+Sum[
+  If[doublePoints[[i]] <= remainingPoints,
+    Ways[throws, remainingPoints - doublePoints[[i]], Length[points]],
+    0],
+  {remainingPoints, 99},
+  {throws, 0, 2},
+  {i, Length[doublePoints]}]

python/eulertest.py

@@ -111,12 +111,15 @@ def main():
 	101: "37076114526",
 	102: "228",
 	104: "329468",
+	108: "180180",
+	109: "38182",
 	112: "1587000",
 	113: "51161058134250",
 	114: "16475640049",
 	115: "168",
 	116: "20492570929",
 	117: "100808458960497",
+	119: "248155780267521",
 	120: "333082500",
 	123: "21035",
 	124: "21417",

python/p108.py

@@ -0,0 +1,53 @@
+# 
+# Solution to Project Euler problem 108
+# by Project Nayuki
+# 
+# https://www.nayuki.io/page/project-euler-solutions
+# https://github.com/nayuki/Project-Euler-solutions
+# 
+
+import eulerlib, itertools
+
+
+# Rewrite the equation with x = n+i, y = n+j, and manipulate it:
+#   1/n = 1/x + 1/y
+#       = 1/(n+i) + 1/(n+j)
+#       = (2n+i+j) / ((n+i)(n+j)).
+#   n(2n+i+j) = (n+i)(n+j).
+#   2n^2 + ni + nj = n^2 + ni + nj + ij.
+#   n^2 = ij.
+# Hence i and j are divisors of n^2. To ensure unique solutions,
+# we impose that x <= y, so i <= j. Also, i > 0, otherwise no j exists.
+# We have i <= j = n^2 / i, thus i^2 <= n^2. With i being positive, we get that i <= n.
+# Therefore the number of solutions for i is the number of divisors of n^2 in the range [1, n].
+# n^2 always has an odd number of divisors. One of them is n. As for the remainder of them, half of them are below n
+# and half of them are above n. So if n^2 has m divisors, then we want (m+1)/2 of them as solutions for i.
+def compute():
+	for n in itertools.count(1):
+		if (count_divisors_squared(n) + 1) // 2 > 1000:
+			return str(n)
+
+
+# Returns the number of divisors of n^2
+def count_divisors_squared(n):
+	count = 1
+	end = eulerlib.sqrt(n)
+	for i in itertools.count(2):
+		if i > end:
+			break
+		if n % i == 0:
+			j = 0
+			while True:
+				n //= i
+				j += 1
+				if n % i != 0:
+					break
+			count *= j * 2 + 1
+			end = eulerlib.sqrt(n)
+	if n != 1:  # Remaining largest prime factor
+		count *= 3
+	return count
+
+
+if __name__ == "__main__":
+	print(compute())

python/p109.py

@@ -0,0 +1,43 @@
+# 
+# Solution to Project Euler problem 109
+# by Project Nayuki
+# 
+# https://www.nayuki.io/page/project-euler-solutions
+# https://github.com/nayuki/Project-Euler-solutions
+# 
+
+
+def compute():
+	# Both lists are orderless but duplicates are important; they are sort of like multisets
+	points = [i * j for i in range(1, 21) for j in range(1, 4)] + [25, 50]
+	doublepoints = [i * 2 for i in range(1, 21)] + [25 * 2]
+	
+	# Memoization array, with dimensions (3, 101, len(points))
+	ways = [[[None] * len(points) for j in range(101)] for i in range(3)]
+	
+	# Number of ways to get exactly 'total' points in exactly 'throwz' throws, using
+	# items (unordered) from the 'points' list with index less than or equal to 'maxIndex'.
+	def calc_ways(throws, total, maxindex):
+		if ways[throws][total][maxindex] is None:
+			if throws == 0:
+				result = 1 if total == 0 else 0
+			else:
+				result = 0
+				if maxindex > 0:
+					result += calc_ways(throws, total, maxindex - 1)
+				if points[maxindex] <= total:
+					result += calc_ways(throws - 1, total - points[maxindex], maxindex)
+			ways[throws][total][maxindex] = result
+		return ways[throws][total][maxindex]
+	
+	checkouts = 0
+	for remainingpoints in range(1, 100):
+		for throws in range(0, 3):
+			for p in doublepoints:
+				if p <= remainingpoints:
+					checkouts += calc_ways(throws, remainingpoints - p, len(points) - 1)
+	return str(checkouts)
+
+
+if __name__ == "__main__":
+	print(compute())

python/p119.py

@@ -0,0 +1,49 @@
+# 
+# Solution to Project Euler problem 119
+# by Project Nayuki
+# 
+# https://www.nayuki.io/page/project-euler-solutions
+# https://github.com/nayuki/Project-Euler-solutions
+# 
+
+import itertools
+
+
+# Candidates have the form n^k, where n >= 2, k >= 2, n^k >= 10, and isDigitSumPower(n^k) == true.
+# We also impose n^k < limit. If there are at least 30 candidates under 'limit',
+# then the 30th smallest candidate is the answer. Otherwise we raise the limit and search again.
+# 
+# We only need to try the exponents k until 2^k exceeds the limit.
+# We only need to try the bases n until the power of the digit sum is too small to match n^k.
+# The power of the digit sum is digitSum(n^k)^k, which is at most (9 * digitLength(n^k))^k.
+def compute():
+	INDEX = 30  # 1-based
+	limit = 1
+	while True:
+		candidates = set()
+		k = 2
+		while (1 << k) < limit:
+			for n in itertools.count(2):
+				pow = n**k
+				if pow >= limit and len(str(pow)) * 9 < n:
+					break
+				if pow >= 10 and is_digit_sum_power(pow):
+					candidates.add(pow)
+			k += 1
+		if len(candidates) >= INDEX:
+			return str(sorted(candidates)[INDEX - 1])
+		limit <<= 8
+
+
+def is_digit_sum_power(x):
+	digitsum = sum(int(c) for c in str(x))
+	if digitsum == 1:  # Powers of 10 are never a power of 1
+		return False
+	pow = digitsum
+	while pow < x:
+		pow *= digitsum
+	return pow == x
+
+
+if __name__ == "__main__":
+	print(compute())