commit 8afba98c4434f978ad66674c02d2f7a91dc7d61e Author: filifa Date: Mon Mar 24 00:06:46 2025 -0400 write up the first few problems, some wip diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..5dc97b5 --- /dev/null +++ b/.gitignore @@ -0,0 +1,2 @@ +.ipynb_checkpoints/ +html/ diff --git a/Makefile b/Makefile new file mode 100644 index 0000000..9ba22d9 --- /dev/null +++ b/Makefile @@ -0,0 +1,14 @@ +.PHONY: all clean + +HTMLS = $(addprefix html/,problem0001.html problem0002.html problem0003.html problem0004.html problem0005.html problem0006.html problem0007.html problem0008.html problem0009.html problem0010.html problem0012.html problem0014.html) + +all: $(HTMLS) + +clean: + rm -r html + +html/%.html: %.ipynb | html + jupyter nbconvert --to html $< --output $@ --execute + +html: + mkdir -p html diff --git a/problem0001.ipynb b/problem0001.ipynb new file mode 100644 index 0000000..11fbec1 --- /dev/null +++ b/problem0001.ipynb @@ -0,0 +1,59 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "fdd3f4eb", + "metadata": {}, + "source": [ + "# Multiples of 3 or 5\n", + "> If we list all the natural numbers below 10 that are multiples of 3 or 5, we get 3, 5, 6, and 9. The sum of these multiples is 23.\n", + ">\n", + "> Find the sum of all the multiples of 3 or 5 below 1000.\n", + "\n", + "Lots of different approaches to this problem! Writing a program may be the easiest, but I always find it interesting when these problems can be solved without a computer doing the heavy lifting, and with a little bit of [summation](https://en.wikipedia.org/wiki/Summation) knowledge and the [inclusion-exclusion principle](https://en.wikipedia.org/wiki/Inclusion%E2%80%93exclusion_principle), we can calculate this answer by hand.\n", + "\n", + "First off, the sum of the multiples of 3 below 1000 is\n", + "$$3 + 6 + 9 + \\cdots + 999 = \\sum_{k=1}^{333} 3k$$\n", + "Likewise, the sum of the multiples of 5 below 1000 is\n", + "$$5 + 10 + 15 + \\cdots + 995 = \\sum_{k=1}^{199} 5k$$\n", + "Although these sums may seem tedious to calculate by hand, we can apply the distributive property and factor out the 3 and 5, respectively, from each summation.\n", + "$$\\sum_{k=1}^{333} 3k = 3\\sum_{k=1}^{333} k$$\n", + "$$\\sum_{k=1}^{199} 5k = 5\\sum_{k=1}^{199} k$$\n", + "We are then left with a summation of the form\n", + "$$\\sum_{k=1}^{n} k$$\n", + "in both instances. If we successively evaluate this summation at $n = 1,2,3,\\ldots$ we get the [triangular numbers](https://en.wikipedia.org/wiki/Triangular_number), which have a well-known closed formula:\n", + "$$\\sum_{k=1}^{n} k = \\frac{n(n+1)}{2}$$\n", + "(This formula can be proven valid with [mathematical induction](https://en.wikipedia.org/wiki/Mathematical_induction), but you don't need to formally prove it to gain an intuition for it. Click the link on triangular numbers to find other explanations.)\n", + "\n", + "Therefore, the first summation equals $3(333(334)/2) = 166833$, and the second equals $5(199(200)/2) = 99500$.\n", + "\n", + "However, we can't just sum these two values. That's because any number that's a multiple of both 3 and 5 is included in *both* of these summations, so if we add the results up, we'll count those values twice! To adjust for this, we'll need to subtract the sum of those numbers from our total so they are only counted once. Luckily, a number is a multiple of both 3 and 5 if and only if it is a multiple of 15, so we can apply our triangular number formula from above.\n", + "$$15 + 30 + 45 + \\cdots + 990 = \\sum_{k=1}^{66} 15k = 15(66(67)/2) = 33165$$\n", + "\n", + "Our answer is then\n", + "$$\\sum_{k=1}^{333} 3k + \\sum_{k=1}^{199} 5k - \\sum_{k=1}^{66} 15k = 233168$$" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "SageMath 9.5", + "language": "sage", + "name": "sagemath" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.11.2" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/problem0002.ipynb b/problem0002.ipynb new file mode 100644 index 0000000..79f0b64 --- /dev/null +++ b/problem0002.ipynb @@ -0,0 +1,84 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "1899e933", + "metadata": {}, + "source": [ + "# Even Fibonacci Numbers\n", + "> Each new term in the Fibonacci sequence is generated by adding the previous two terms. By starting with 1 and 2, the first terms will be:\n", + "> $$1, 2, 3, 5, 8, 13, 21, 34, 55, 89, \\ldots$$\n", + "> By considering the terms in the Fibonacci sequence whose values do not exceed four million, find the sum of the even-valued terms.\n", + "\n", + "There's a lot to learn about the [Fibonacci numbers](https://en.wikipedia.org/wiki/Fibonacci_sequence), and this will not be the last time we see them in a Project Euler problem, but for now this problem keeps it relatively simple.\n", + "\n", + "Since this is a SageMath notebook, we'll use its functionality so we can keep our code simple like the problem." + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "id": "57b68465", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "4613732" + ] + }, + "execution_count": 1, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "sum(n for n in fibonacci_xrange(1, 4000000) if n % 2 == 0)" + ] + }, + { + "cell_type": "markdown", + "id": "eef788da", + "metadata": {}, + "source": [ + "Even without SageMath, it's easy to implement our own generator of the sequence up to $n$." + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "id": "24a7c257", + "metadata": {}, + "outputs": [], + "source": [ + "def fib(n):\n", + " a = 0\n", + " b = 1\n", + " while a <= n:\n", + " yield a\n", + " a, b = b, a + b" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "SageMath 9.5", + "language": "sage", + "name": "sagemath" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.11.2" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/problem0003.ipynb b/problem0003.ipynb new file mode 100644 index 0000000..9cdc5fa --- /dev/null +++ b/problem0003.ipynb @@ -0,0 +1,115 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "1a8d8609", + "metadata": {}, + "source": [ + "# Largest Prime Factor\n", + "> The prime factors of 13195 are 5, 7, 13, and 29.\n", + "> \n", + "> What is the largest prime factor of the number 600851475143?\n", + "\n", + "The [fundamental theorem of arithmetic](https://en.wikipedia.org/wiki/Fundamental_theorem_of_arithmetic) says that every integer greater than 1 has a unique prime factorization.\n", + "\n", + "There are numerous factorization algorithms. Some of the easiest to implement are [trial division](https://en.wikipedia.org/wiki/Trial_division), [wheel factorization](https://en.wikipedia.org/wiki/Wheel_factorization) and [Pollard's rho algorithm](https://en.wikipedia.org/wiki/Pollard%27s_rho_algorithm).\n", + "\n", + "In general, as an integer gets larger, [factoring it gets harder](https://en.wikipedia.org/wiki/Integer_factorization). We don't have any efficient algorithm for factoring in general, but somewhat strangely, we don't know for certain that no such algorithm exists.\n", + "\n", + "Fortunately for us, this number is not too large for any of these algorithms, or for SageMath (you can even factor it by hand if you have time to kill)." + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "id": "d2f8884c", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "[71, 839, 1471, 6857]" + ] + }, + "execution_count": 1, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "prime_factors(600851475143)" + ] + }, + { + "cell_type": "markdown", + "id": "d4343a42", + "metadata": {}, + "source": [ + "If we wanted to implement our own factorization algorithm, though, Pollard's rho is relatively easy to implement and more efficient (on average) than trial division. It's also got some interesting math at its foundation, so let's take a look at it.\n", + "\n", + "## Pollard's rho algorithm\n", + "Let $n=pq$ be a positive integer where $p$ and $q$ are non-trivial factors of $n$ (i.e. not 1 or $n$), and let $\\{x_i\\}$ be a sequence such that $x_0$ is some random initial value and\n", + "$$x_{i+1} = (x_i^2 - 1) \\bmod{n}$$\n", + "Since the number of values this sequence can take on is finite, and each term is determined entirely by the previous term, it must ultimately become periodic, although there may be some initial terms before we reach a term that will repeat (refer to [this image](https://en.wikipedia.org/wiki/File:Pollard_rho_cycle.svg) for a visual aid - the fact this diagram looks like the Greek letter $\\rho$ is where the algorithm gets its name).\n", + "\n", + "How many terms will we see before the sequence repeats? If we assume the sequence is random, the [birthday paradox](https://en.wikipedia.org/wiki/Birthday_problem) implies that it will take an average of $O(\\sqrt{n})$ iterations before it repeats. Of course, the sequence is not truly random, but it *is* [pseudorandom](https://en.wikipedia.org/wiki/Pseudorandomness), and experimental results seem to corroborate this analysis despite the somewhat faulty assumption.\n", + "\n", + "Now consider the related sequence $\\{x_i \\bmod p\\}$. Obviously, we can't calculate this sequence directly since we don't know the value of $p$, but similarly to $\\{x_i\\}$, this sequence is eventually periodic, with an average of $O(\\sqrt{p})$ iterations before a repetition. In other words, this sequence is expected to repeat earlier than $\\{x_i\\}$.\n", + "\n", + "When $\\{x_i \\bmod p\\}$ cycles, it implies that $x_j \\equiv x_k \\pmod{p}$ for some $j \\neq k$. This implies that $p | (x_j - x_k)$. Since $p$ also divides $n$, this means that $\\gcd(x_j - x_k, n)$ will be some multiple of $p$.\n", + "\n", + "This all means that we simply need to iterate through terms in $\\{x_i\\}$ until we find $x_j$ and $x_k$ such that $\\gcd(x_j - x_k, n)$ equals something other than 1. When it does, by definition this will be a factor of $n$. However, this factor may be $n$ itself, which is trivial. In this case the algorithm fails, but it can be attempted again with a different value for $x_0$.\n", + "\n", + "To facilitate finding $x_j$ and $x_k$, [Floyd's cycle-finding algorithm](https://en.wikipedia.org/wiki/Cycle_detection) is used. Briefly, two variables are used to iterate through $\\{x_i\\}$ at different rates. Because of this, the two variables will never have the same index, but eventually $x_i \\equiv x_j \\pmod{p}$, which we will detect with the GCD." + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "id": "d12993bf", + "metadata": {}, + "outputs": [], + "source": [ + "from itertools import count\n", + "\n", + "def factorization(n):\n", + " x = randint(0, n-1)\n", + " y = x\n", + " k = 2\n", + " for i in count(1):\n", + " x = (x ** 2 - 1) % n\n", + " d = gcd(y - x, n)\n", + " \n", + " assert d != n\n", + " \n", + " if d != 1:\n", + " return d\n", + " \n", + " if i == k:\n", + " y = x\n", + " k *= 2" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "SageMath 9.5", + "language": "sage", + "name": "sagemath" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.11.2" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/problem0004.ipynb b/problem0004.ipynb new file mode 100644 index 0000000..4e02b8f --- /dev/null +++ b/problem0004.ipynb @@ -0,0 +1,63 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "9c7a281d", + "metadata": {}, + "source": [ + "# Largest Palindrome Product\n", + "> A palindromic number reads the same both ways. The largest palindrome made from the product of two 2-digit numbers is $9009 = 91 \\times 99$.\n", + ">\n", + "> Find the largest palindrome made from the product of two 3-digit numbers.\n", + "\n", + "Our search space is only $\\binom{900}{2} = 404550$ 3-digit pairs - a lot to check by hand but peanuts to a modern computer. " + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "id": "e5fd66f6", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "906609" + ] + }, + "execution_count": 1, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "from itertools import combinations\n", + "\n", + "is_palindrome = lambda x: str(x) == str(x)[::-1]\n", + "three_digit_pairs = combinations(range(100, 1000), 2)\n", + "max(x*y for (x,y) in three_digit_pairs if is_palindrome(x*y))" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "SageMath 9.5", + "language": "sage", + "name": "sagemath" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.11.2" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/problem0005.ipynb b/problem0005.ipynb new file mode 100644 index 0000000..112d047 --- /dev/null +++ b/problem0005.ipynb @@ -0,0 +1,77 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "088fe896", + "metadata": {}, + "source": [ + "# Smallest Multiple\n", + "> 2520 is the smallest number that can be divided by each of the numbers from 1 to 10 without any remainder.\n", + "> \n", + "> What is the smallest positive number that is evenly divisible by all of the numbers from 1 to 20?\n", + "\n", + "Another problem we can solve by hand! We're looking for the [least common multiple](https://en.wikipedia.org/wiki/Least_common_multiple), or LCM, of all the numbers from 1 to 20. Obviously,\n", + "$$20! = 20 \\times 19 \\times 18 \\times \\cdots \\times 1 = 2432902008176640000$$\n", + "is *a* common multiple, but the LCM is indeed smaller than this.\n", + "\n", + "One way to find the LCM is to first determine the prime factorizations of these numbers.\n", + "* $1 = 1$\n", + "* $2 = 2$\n", + "* $3 = 3$\n", + "* $4 = 2^2$\n", + "* $5 = 5$\n", + "* $6 = 2 \\times 3$\n", + "* $7 = 7$\n", + "* $8 = 2^3$\n", + "* $9 = 3^2$\n", + "* $10 = 2 \\times 5$\n", + "* $11 = 11$\n", + "* $12 = 2^2 \\times 3$\n", + "* $13 = 13$\n", + "* $14 = 2 \\times 7$\n", + "* $15 = 3 \\times 5$\n", + "* $16 = 2^4$\n", + "* $17 = 17$\n", + "* $18 = 2 \\times 3^2$\n", + "* $19 = 19$\n", + "* $20 = 2^2 \\times 5$\n", + "\n", + "Now, for each prime factor that appears in these factorizations, we'll take note of the largest power they're raised to. For instance, $2^4$ (the prime factorization of 16) is the largest power of 2 appearing.\n", + "\n", + "The LCM is the product of all these \"largest prime powers.\"\n", + "$$2^4 \\times 3^2 \\times 5 \\times 7 \\times 11 \\times 13 \\times 17 \\times 19 = 232792560$$\n", + "How do we know this is the LCM? Well, it's trivial to see it is *a* common multiple - just divide by each of the above numbers to confirm. Furthermore, we know it's the *least* common multiple because if we divided by any of these prime factors to get a smaller result, it would no longer be divisible by one or more of our numbers. For instance, if we divide by 3, our result is no longer divisible by $9 = 3^2$ or $18 = 2 \\times 3^2$.\n", + "\n", + "## Alternate method\n", + "I personally find the prime factorization method above the easiest way to solve this specific problem by hand. However, if there were a lot more numbers or the numbers were larger, I wouldn't want to do all those factorizations. Fortunately, there is another method that doesn't involve factoring at all, and can be implemented efficiently on a computer.\n", + "\n", + "For two numbers $m$ and $n$,\n", + "$$\\mathrm{lcm}(m,n) = \\frac{mn}{\\gcd(m,n)}$$\n", + "where $\\gcd$ is the [greatest common divisor](https://en.wikipedia.org/wiki/Greatest_common_divisor). The [Euclidean algorithm](https://en.wikipedia.org/wiki/Euclidean_algorithm) gives an efficient method for calculating the GCD.\n", + "\n", + "This formula is only defined for two numbers, but if we want to find an LCM of a set of three or more numbers, we can simply find the LCM of any two numbers in the set, then proceed to find the LCM of that value and another number from the set, repeating until we have used each number. For example, to find the LCM of 5, 8, and 14, you can first find the LCM of 5 and 8 using the above formula (40), then find the LCM of 40 and 14 (280)." + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "SageMath 9.5", + "language": "sage", + "name": "sagemath" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.11.2" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/problem0006.ipynb b/problem0006.ipynb new file mode 100644 index 0000000..8713b28 --- /dev/null +++ b/problem0006.ipynb @@ -0,0 +1,52 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "a39ba505", + "metadata": {}, + "source": [ + "# Sum Square Difference\n", + "> The sum of the squares of the first ten natural numbers is,\n", + "> $$1^2 + 2^2 + \\cdots + 10^2 = 385$$\n", + "> The square of the sum of the first ten natural numbers is,\n", + "> $$(1 + 2 + \\cdots + 10)^2 = 55^2 = 3025$$\n", + "> Hence the difference between the sum of the squares of the first ten natural numbers and the square of the sum is $3025 - 385 = 2640$.\n", + "> \n", + "> Find the difference between the sum of the squares of the first one hundred natural numbers and the square of the sum.\n", + "\n", + "In [problem 1](https://projecteuler.net/problem=1), we applied the following formula for [triangular numbers](https://en.wikipedia.org/wiki/Triangular_number):\n", + "$$\\sum_{k=1}^n k = \\frac{n(n+1)}{2}$$\n", + "We can apply it again here and determine that\n", + "$$(1 + 2 + 3 + \\cdots + 100)^2 = \\left(\\frac{100(101)}{2}\\right)^2 = 25502500$$\n", + "\n", + "A similar formula exists for computing [sums of squares](https://en.wikipedia.org/wiki/Square_pyramidal_number):\n", + "$$\\sum_{k=1}^n k^2 = \\frac{n(n+1)(2n+1)}{6}$$\n", + "(In fact, [Faulhaber's formula](https://en.wikipedia.org/wiki/Faulhaber%27s_formula) gives a formula for the sum of $k$th powers, but we obviously only need the cases $k=1$ and $k=2$ for this problem.) Consequently,\n", + "$$1^2 + 2^2 + 3^2 + \\cdots + 100^2 = \\frac{100(101)(201)}{6} = 338350$$\n", + "\n", + "Therefore, the difference is $25502500 - 338350 = 25164150$." + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "SageMath 9.5", + "language": "sage", + "name": "sagemath" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.11.2" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/problem0007.ipynb b/problem0007.ipynb new file mode 100644 index 0000000..5bc6f0f --- /dev/null +++ b/problem0007.ipynb @@ -0,0 +1,181 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "87fa1305", + "metadata": {}, + "source": [ + "# 100001st Prime\n", + "> By listing the first six prime numbers: 2, 3, 5, 7, 11, and 13, we can see that the 6th prime is 13.\n", + ">\n", + "> What is the 10001st prime number?\n", + "\n", + "Unlike integer factorization, we *can* efficiently determine whether a number is prime. One of the simplest and fastest methods is the [Miller-Rabin primality test](https://en.wikipedia.org/wiki/Miller%E2%80%93Rabin_primality_test), which is a probabilistic test.\n", + "\n", + "Of course, SageMath has tools to let us compute the 10001st prime without implementing a primality test ourselves! " + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "id": "3a27f454", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "104743" + ] + }, + "execution_count": 1, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "# Primes.unrank is 0-indexed, so the 10001st prime is at index 10000\n", + "Primes().unrank(10000)" + ] + }, + { + "cell_type": "markdown", + "id": "6f644292", + "metadata": {}, + "source": [ + "Let's look into how the Miller-Rabin test works anyway, though.\n", + "\n", + "## Fermat's little theorem\n", + "First, a little bit of elementary number theory. [Fermat's little theorem](https://en.wikipedia.org/wiki/Fermat%27s_little_theorem) says that if $p$ is a prime number, then for any integer $a$,\n", + "$$a^p \\equiv a \\pmod{p}$$\n", + "It follows from the [contrapositive](https://en.wikipedia.org/wiki/Contraposition) that if there exists an integer $a$ such that $a^p \\not\\equiv a \\pmod{p}$, then $p$ is not prime. We then call $a$ a *witness* that proves $p$ is composite.\n", + "\n", + "Consider the fact that $3^{6} \\equiv 3 \\pmod{6}$, but 6 is composite. In cases like this, we say that $p$ is a *base $a$ Fermat pseudoprime*, e.g. 6 is a base 3 Fermat pseudoprime. 6 is also a base 4 Fermat pseudoprime, but 2 and 5 both serve as witnesses to 6 being composite (we need only consider bases modulo $p$).\n", + "\n", + "This may lead us to believe that we can just check every base $a$, and if there are no witnesses, then $p$ is prime. However, (along with this method being inefficient for large $p$), there are numbers that are Fermat pseudoprimes for *every* base, called [Carmichael numbers](https://en.wikipedia.org/wiki/Carmichael_number). Therefore, the [converse](https://en.wikipedia.org/wiki/Converse_(logic)) of Fermat's little theorem does not hold.\n", + "\n", + "Nevertheless, Carmichael numbers are relatively rare, so if you pick a random number to test and a few random bases, and find no witnesses, you can be pretty sure the number is in fact prime.\n", + "\n", + "## Miller-Rabin test\n", + "\n", + "The Miller-Rabin test employs another useful fact about prime numbers: if $p$ is prime, then the only values of $x$ such that \n", + "$$x^2 \\equiv 1 \\pmod{p}$$\n", + "are $x \\equiv \\pm 1 \\pmod{p}$. (*Proof: if $x^2 \\equiv 1 \\pmod{p}$, then $p \\vert (x+1)(x-1)$. Then by [Euclid's lemma](https://en.wikipedia.org/wiki/Euclid%27s_lemma), $p$ must divide one of $x+1$ or $x-1$, implying $x \\equiv \\pm 1 \\pmod{p}$.*)\n", + "\n", + "Similarly to Fermat's little theorem, the converse of this statement doesn't hold. For example, $x^2 \\equiv 1 \\pmod{9}$ only has $x=1$ and $x=-1$ as solutions, but 9 is composite.\n", + "\n", + "We now have two questions to lead an investigation into the primality of a number $n$:\n", + "* Is there an integer $a$ such that $a^n \\not\\equiv a \\pmod{n}$?\n", + "* Does 1 have a non-trivial square root modulo $n$?\n", + "\n", + "If the answer to either of these is yes, then $n$ is *definitely* composite. If both answers are no, this is strong evidence that $n$ is prime, *but not proof*, since we have already demonstrated that both of these conditions can have false positives. A composite number that is erroneously suggested to be prime by these conditions is called a *base $a$ strong pseudoprime*. Fortunately, no number (including Carmichael numbers) is a base $a$ strong pseudoprime for all bases $a$. In fact, Rabin (one of the namesakes of the algorithm) proved that if $n$ is composite, at least 3/4 of the numbers 1 to $n-1$ are witnesses.\n", + "\n", + "Let's finally look at the actual test. Suppose we want to test if an odd number $n$ is prime (testing whether an even number is prime is trivial). There must exist $t \\geq 1$ and odd $u$ such that $n - 1 = 2^t u$." + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "id": "720f7e8c", + "metadata": {}, + "outputs": [], + "source": [ + "def factor_twos(n):\n", + " t = 0\n", + " u = n - 1\n", + " while u % 2 == 0:\n", + " u //= 2\n", + " t += 1\n", + " \n", + " return t, u" + ] + }, + { + "cell_type": "markdown", + "id": "980bd2cf", + "metadata": {}, + "source": [ + "It follows that\n", + "$$a^{n-1} \\equiv a^{2^t u} \\pmod{n}$$\n", + "for any base $a$, so we'll randomly select a base $1 \\leq a \\leq n-1$.\n", + "\n", + "The test proceeds by repeatedly squaring $a^u$ modulo $n$ until we reach $a^{n-1}$. If at any point in this squaring, we find that $a^{2^{i-1} u}$ does not equal 1 or -1, but the following squaring gives $a^{2^i u} = 1$, then we have found a non-trivial square root of 1 modulo $n$. Therefore, $n$ must be composite.\n", + "\n", + "After $t$ squarings, we will reach $a^{n-1}$. If this equals anything other than 1, then by the contrapositive of Fermat's last theorem, this also tells us that $n$ is composite.\n", + "\n", + "If, after $t$ squarings, we have not found a non-trivial square root of 1 modulo $n$ and $a^{n-1} \\equiv 1 \\pmod{n}$, then either $n$ is prime, or is a base $a$ strong pseudoprime." + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "id": "123d1068", + "metadata": {}, + "outputs": [], + "source": [ + "def witness(a, n):\n", + " t, u = factor_twos(n)\n", + " x = pow(a, u, n)\n", + " for _ in range(0, t):\n", + " x, y = pow(x, 2, n), x\n", + " if x == 1 and y != 1 and y != n - 1:\n", + " return True\n", + " \n", + " return x != 1" + ] + }, + { + "cell_type": "markdown", + "id": "f8454603", + "metadata": {}, + "source": [ + "You can reduce the chance of a false report of primality by repeating the procedure with different bases. If we test several bases and do not find any witnesses, it is exceedingly unlikely that $n$ is composite." + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "id": "e71ba6ec", + "metadata": {}, + "outputs": [], + "source": [ + "from random import randint\n", + "\n", + "def is_prime(n, k=10):\n", + " if n < 2:\n", + " return False\n", + " elif n == 2:\n", + " return True\n", + " elif n % 2 == 0:\n", + " return False\n", + " \n", + " for _ in range(0, k):\n", + " a = randint(1, n - 1)\n", + " if witness(a, n):\n", + " return False\n", + " \n", + " return True" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "SageMath 9.5", + "language": "sage", + "name": "sagemath" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.11.2" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/problem0008.ipynb b/problem0008.ipynb new file mode 100644 index 0000000..a4ae8eb --- /dev/null +++ b/problem0008.ipynb @@ -0,0 +1,107 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "8744af43", + "metadata": {}, + "source": [ + "# Largest Product in a Series\n", + "> The four adjacent digits in the 1000-digit number that have the greatest product are $9 \\times 9 \\times 8 \\times 9 = 5832$.\n", + "> ```\n", + "73167176531330624919225119674426574742355349194934\n", + "96983520312774506326239578318016984801869478851843\n", + "85861560789112949495459501737958331952853208805511\n", + "12540698747158523863050715693290963295227443043557\n", + "66896648950445244523161731856403098711121722383113\n", + "62229893423380308135336276614282806444486645238749\n", + "30358907296290491560440772390713810515859307960866\n", + "70172427121883998797908792274921901699720888093776\n", + "65727333001053367881220235421809751254540594752243\n", + "52584907711670556013604839586446706324415722155397\n", + "53697817977846174064955149290862569321978468622482\n", + "83972241375657056057490261407972968652414535100474\n", + "82166370484403199890008895243450658541227588666881\n", + "16427171479924442928230863465674813919123162824586\n", + "17866458359124566529476545682848912883142607690042\n", + "24219022671055626321111109370544217506941658960408\n", + "07198403850962455444362981230987879927244284909188\n", + "84580156166097919133875499200524063689912560717606\n", + "05886116467109405077541002256983155200055935729725\n", + "71636269561882670428252483600823257530420752963450\n", + "> ```\n", + "> Find the thirteen adjacent digits in the 1000-digit number that have the greatest product. What is the value of this product?\n", + "\n", + "First, let's load the number as a tuple." + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "id": "699aa447", + "metadata": {}, + "outputs": [], + "source": [ + "num = tuple(int(d) for d in \"7316717653133062491922511967442657474235534919493496983520312774506326239578318016984801869478851843858615607891129494954595017379583319528532088055111254069874715852386305071569329096329522744304355766896648950445244523161731856403098711121722383113622298934233803081353362766142828064444866452387493035890729629049156044077239071381051585930796086670172427121883998797908792274921901699720888093776657273330010533678812202354218097512545405947522435258490771167055601360483958644670632441572215539753697817977846174064955149290862569321978468622482839722413756570560574902614079729686524145351004748216637048440319989000889524345065854122758866688116427171479924442928230863465674813919123162824586178664583591245665294765456828489128831426076900422421902267105562632111110937054421750694165896040807198403850962455444362981230987879927244284909188845801561660979191338754992005240636899125607176060588611646710940507754100225698315520005593572972571636269561882670428252483600823257530420752963450\")" + ] + }, + { + "cell_type": "markdown", + "id": "8b53af56", + "metadata": {}, + "source": [ + "Now we'll just do a sliding window over the number to find the thirteen adjacent digits with the greatest product." + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "id": "6b760673", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "(5, 5, 7, 6, 6, 8, 9, 6, 6, 4, 8, 9, 5)" + ] + }, + "execution_count": 2, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "subnums = (num[i:i+13] for i in range(0, len(num)-13))\n", + "max(subnums, key=prod)" + ] + }, + { + "cell_type": "markdown", + "id": "4148d160", + "metadata": {}, + "source": [ + "Their product is 23514624000." + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "SageMath 9.5", + "language": "sage", + "name": "sagemath" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.11.2" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/problem0009.ipynb b/problem0009.ipynb new file mode 100644 index 0000000..5aeb6bc --- /dev/null +++ b/problem0009.ipynb @@ -0,0 +1,119 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "df906ad5", + "metadata": {}, + "source": [ + "# Special Pythagorean Triplet\n", + "> A Pythagorean triplet is a set of three natural numbers, $a < b < c$, for which, \n", + "> $$a^2 + b^2 = c^2$$\n", + "> For example, $3^2 + 4^2 = 9 + 16 = 25 = 5^2$.\n", + ">\n", + "> There exists exactly one Pythagorean triplet for which $a + b + c = 1000$.\n", + ">\n", + "> Find the product $abc$.\n", + "\n", + "The key mathematical insight for this problem is that Euclid gave a parameterization for generating [Pythagorean triplets](https://en.wikipedia.org/wiki/Pythagorean_triple): given integers $m, n$ such that $m > n > 0$, a triplet is given by\n", + "$$\n", + "\\begin{align}\n", + "a &= m^2 - n^2 \\\\\n", + "b &= 2mn \\\\\n", + "c &= m^2 + n^2\n", + "\\end{align}\n", + "$$\n", + "Furthermore, if $m$ and $n$ are [coprime](https://en.wikipedia.org/wiki/Coprime_integers) (i.e. $\\gcd(m,n)=1$) and exactly one of $m$ and $n$ is even, they will generate a primitive triplet, i.e. a triplet such that $\\gcd(a,b,c)=1$." + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "id": "644fb73e", + "metadata": {}, + "outputs": [], + "source": [ + "from itertools import count\n", + "\n", + "def primitive_pythagorean_triplets():\n", + " for m in count(2):\n", + " for n in range(1, m):\n", + " if not ((m % 2) ^^ (n % 2)) or gcd(m, n) != 1:\n", + " continue\n", + "\n", + " a = m ** 2 - n ** 2\n", + " b = 2 * m * n\n", + " c = m ** 2 + n ** 2\n", + "\n", + " yield (a, b, c)" + ] + }, + { + "cell_type": "markdown", + "id": "db3c3419", + "metadata": {}, + "source": [ + "We can simply iterate through these primitive triplets, and for each one, scale them until we either find our special triplet or their sum is greater than 1000, at which point we'll try the next primitive triplet." + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "id": "7b2c3efa", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "(375, 200, 425)" + ] + }, + "execution_count": 2, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "def find_special_triplet():\n", + " for (a, b, c) in primitive_pythagorean_triplets():\n", + " for k in count(1):\n", + " if k * (a + b + c) == 1000:\n", + " return (k * a, k * b, k * c)\n", + " \n", + " if k * (a + b + c) > 1000:\n", + " break\n", + " \n", + " \n", + "find_special_triplet()" + ] + }, + { + "cell_type": "markdown", + "id": "8809ed3b", + "metadata": {}, + "source": [ + "So our product is $abc=31875000$." + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "SageMath 9.5", + "language": "sage", + "name": "sagemath" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.11.2" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/problem0010.ipynb b/problem0010.ipynb new file mode 100644 index 0000000..bb9dde6 --- /dev/null +++ b/problem0010.ipynb @@ -0,0 +1,88 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "ec6f67ee", + "metadata": {}, + "source": [ + "# Summation of Primes\n", + "> The sum of the primes below 10 is 2 + 3 + 5 + 7 = 17.\n", + ">\n", + "> Find the sum of all the primes below two million.\n", + "\n", + "Once again, SageMath comes with functions that do the hard part for us." + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "id": "127873d5", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "142913828922" + ] + }, + "execution_count": 1, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "sum(prime_range(2000000))" + ] + }, + { + "cell_type": "markdown", + "id": "46f100db", + "metadata": {}, + "source": [ + "If `prime_range` didn't exist, we may have attempted using a primality test (like in [problem 7](https://projecteuler.net/problem=7)) to test each number less than 2,000,000, but for this problem, there's an even simpler way than implementing a test: the [sieve of Eratosthenes](https://en.wikipedia.org/wiki/Sieve_of_Eratosthenes), which can even be carried out on pen and paper (although maybe not for primes up to 2,000,000).\n", + "\n", + "Given an upper limit $n$, the sieve determines all the primes less than or equal to $n$. It does this by first assuming every number from 2 to $n$ is prime. Then, since 2 is assumed to be prime, every multiple of 2 up to $n$ is marked as composite. After this, we skip to the next number still assumed prime (3 for the second iteration) and mark all its multiples as composite. Repeat for each successive \"assumed\" prime. At the end, any number not marked composite is prime. We can then add those primes up to get our answer." + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "id": "35ec90ad", + "metadata": {}, + "outputs": [], + "source": [ + "def sieve(n):\n", + " # leave out multiples of 2 to save memory\n", + " P = set(range(3, n + 1, 2))\n", + " \n", + " for i in range(3, isqrt(n) + 1, 2):\n", + " if i in P:\n", + " P -= set(range(i ** 2, n + 1, 2 * i))\n", + " \n", + " P.add(2)\n", + " return P" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "SageMath 9.5", + "language": "sage", + "name": "sagemath" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.11.2" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/problem0012.ipynb b/problem0012.ipynb new file mode 100644 index 0000000..30fbb19 --- /dev/null +++ b/problem0012.ipynb @@ -0,0 +1,95 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "802ccaa1", + "metadata": {}, + "source": [ + "# Highly Divisible Triangular Number\n", + "> The sequence of triangle numbers is generated by adding the natural numbers. So the 7th triangle number would be $1+2+3+4+5+6+7=28$. The first ten terms would be:\n", + "> $$1,3,6,10,15,21,28,36,45,55,\\ldots$$\n", + "> Let us list the factors of the first seven triangle numbers:\n", + "> * 1: 1\n", + "> * 3: 1,3\n", + "> * 6: 1,2,3,6\n", + "> * 10: 1,2,5,10\n", + "> * 15: 1,3,5,15\n", + "> * 21: 1,3,7,21\n", + "> * 28: 1,2,4,7,14,28\n", + "> \n", + "> We can see that 28 is the first triangle number to have over five divisors.\n", + ">\n", + "> What is the value of the first triangle number to have over five hundred divisors?\n", + "\n", + "Yet another problem with a straightforward solution thanks to SageMath's preinstalled functions. We'll also once again apply the closed formula for triangle numbers." + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "id": "c6249c71", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "76576500\n" + ] + } + ], + "source": [ + "for n in PositiveIntegers():\n", + " t = n * (n+1) / 2\n", + " if number_of_divisors(t) > 500:\n", + " print(t)\n", + " break" + ] + }, + { + "cell_type": "markdown", + "id": "b504d3c6", + "metadata": {}, + "source": [ + "If we wanted to implement our own [divisor counting function](https://en.wikipedia.org/wiki/Divisor_function) (denoted $\\sigma_0$), we could apply a factorization algorithm (as discussed for [problem 3](https://projecteuler.net/problem=3)) to compute each triangle number's prime factorization\n", + "$$n = 2^{r_1} 3^{r_2} 5^{r_3} 7^{r_4} \\cdots$$\n", + "(Note that every number can be written in this form, but some exponents will be 0 if the corresponding prime number is not a divisor of $n$.)\n", + "From there, the number of divisors is simply\n", + "$$\\sigma(2^{r_1} 3^{r_2} 5^{r_3} 7^{r_4} \\cdots) = (r_1 + 1)(r_2 + 1)(r_3 + 1)\\cdots$$\n", + "For example, $12 = 2^2 \\times 3$, so $\\sigma_0(12) = (2 + 1)(1 + 1) = 6$ (the factors are 1, 2, 3, 4, 6, and 12).\n", + "\n", + "If this formula doesn't make sense to you, try thinking about it like this: a divisor of $n$ will have the same construction $2^{s_1} 3^{s_2} 5^{s_3} 7^{s_4} \\cdots$, but each exponent must be less than or equal to its corresponding exponent in the factorization of $n$ (if it's larger, it won't divide $n$). This means $s_1$ can be any value from $0,1,2,\\ldots,r_1$. There are therefore $r_1 + 1$ options for the value of $s_1$, $r_2 + 1$ options for the value of $s_2$, and so on. Each choice is independent of the others, so we can multiply the number of options for each exponent together to get the total number of options, which gets us the total number of divisors.\n", + "\n", + "Continuing with our previous example, to construct a divisor of $12=2^2 \\times 3$, our options for the power of 2 are 0, 1, or 2, and our options for the power of 3 are 0 or 1. This gives us 6 options altogether for the two exponents. Sure enough, here are the 6 divisors of 12, each composed of one of $2^0$, $2^1$, or $2^2$, and one of $3^0$ or $3^1$:\n", + "\n", + "* $1=1$\n", + "* $2=2$\n", + "* $3=3$\n", + "* $4=2^2$\n", + "* $6=2 \\times 3$\n", + "* $12=2^2 \\times 3$\n" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "SageMath 9.5", + "language": "sage", + "name": "sagemath" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.11.2" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/problem0014.ipynb b/problem0014.ipynb new file mode 100644 index 0000000..f97a962 --- /dev/null +++ b/problem0014.ipynb @@ -0,0 +1,98 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "43f43a42", + "metadata": {}, + "source": [ + "# Longest Collatz Sequence\n", + "> The following iterative sequence is defined for the set of positive integers:\n", + "> * $n \\to n/2$ ($n$ is even)\n", + "> * $n \\to 3n + 1$ ($n$ is odd)\n", + "> \n", + "> Using the rule above and starting with 13, we generate the following sequence: \n", + "> $$13 \\to 40 \\to 20 \\to 10 \\to 5 \\to 16 \\to 8 \\to 4 \\to 2 \\to 1$$\n", + "> It can be seen that this sequence (starting at 13 and finishing at 1) contains 10 terms. Although it has not been proved yet (Collatz Problem), it is thought that all starting numbers finish at 1.\n", + "> \n", + "> Which starting number, under one million, produces the longest chain?\n", + "> \n", + "> **NOTE:** Once the chain starts the terms are allowed to go above one million.\n", + "\n", + "The [Collatz conjecture](https://en.wikipedia.org/wiki/Collatz_conjecture) is a famous unsolved problem, and a classic example of a seemingly simple question that has proven very difficult, if not impossible, to answer.\n", + "\n", + "It's easy enough to *define* a [recursive](https://en.wikipedia.org/wiki/Recursion_(computer_science%29) function to calculate the chain length for a starting number $n$.\n", + "$$\n", + "f(n) = \\begin{cases}\n", + "1 & n = 1 \\\\\n", + "1+f(n/2) & n \\equiv 0 \\pmod{2} \\\\\n", + "1+f(3n+1) & n \\neq 1\\ \\text{and}\\ n \\equiv 1 \\pmod{2}\n", + "\\end{cases}\n", + "$$\n", + "However, we want its *implementation* to be efficient. We can optimize greatly if we cache the outputs we compute (the computer science term for this is [memoization](https://en.wikipedia.org/wiki/Memoization)). For instance, if we store the fact that $f(4) = 3$ after we've computed it, when we later compute $f(8) = 1 + f(4)$, the program can use the stored value of 3 rather than recomputing $f(4)$. For large inputs, this will save us (or really the computer, I guess) from redoing work.\n", + "\n", + "Python has a nice decorator called [cache](https://docs.python.org/3/library/functools.html) that will automagically memoize our function." + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "id": "d1017296", + "metadata": {}, + "outputs": [], + "source": [ + "from functools import cache\n", + "\n", + "@cache\n", + "def collatz_length(n):\n", + " if n == 1:\n", + " return 1\n", + " elif n % 2 == 0:\n", + " return 1 + collatz_length(n // 2)\n", + " else:\n", + " return 1 + collatz_length(3 * n + 1)" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "id": "fd145a5b", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "837799" + ] + }, + "execution_count": 2, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "max(range(1, 1000000), key=collatz_length)" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "SageMath 9.5", + "language": "sage", + "name": "sagemath" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.11.2" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +}