2018-4-2

Project Euler Problem 17 - Number Letter Counts - solved with R

Solution to the 17th Project Euler Problem - Number Letter Counts - coded and visualised with R.

The 15th Project Euler problem - Number Letter Counts - is stated as follows. If the numbers 1 to 5 are written out in words: one, two, three, four, five, then there are 3 + 3 + 5 + 4 + 4 = 19 letters used in total. If all the numbers from 1 to 1000 were written out in words, how many letters would be used?

Do not count spaces or hyphens. For example, 342 (three hundred and forty-two) contains 23 letters and 115 (one hundred and fifteen) contains 20 letters. The use of "and" when writing out numbers is in compliance with British usage.

nchars.singles = nchar(c("one", "two", "three", "four", "five", "six", "seven", "eight", "nine"))
nchars.teens = nchar(c("ten", "eleven", "twelve", "thirteen", "fourteen", "fifteen", "sixteen", "seventeen", "eighteen", "nineteen"))
nchars.tens =  nchar(c("twenty", "thirty", "forty", "fifty", "sixty", "seventy", "eighty", "ninety"))

# One, two, three...nine.
nchars.001.009 <- sum(nchars.singles)

# Ten, eleven, twelve...ninteen.
nchars.010.019 <- sum(nchars.teens)

# Twentyone, twentytwo...nintynine.
nchars.020.099 <- 10*sum(nchars.tens) + 8 * nchars.001.009

# One, two, three...ninetynine.
nchars.001.099 <- nchars.001.009 + nchars.010.019 + nchars.020.099

# When counting from 1 - 999 we will say 1-9 a hundred times each
# while quantifying hundreds.
# ONE hundred, ONE hundred and one, ONE hundred and two.
a <- nchars.001.009 * 100

# When counting from 1 - 999 we will count from 1-99 ten times.
b <- nchars.001.099 * 10

# Nine of the numbers will include just "hundred"
# One hundred, two hundred, three hundred...nine hundred.
c <- length(nchars.singles) * nchar("hundred")

# 891 (99*9) of the numbers will include "hundred and".
# One hundred and one, one hundred and two... nine hundred and ninety nine.
d <- 891 * nchar("hundredand")

# ...and finnish with One Thousand.
e <- sum(nchar("onethousand"))

sum <- a+b+c+d+e;
sprintf("Number of letters: %s", sum)
2018-4-2

Project Euler Problem 16 - Power Digit Sum - solved with Javascript

Solution to the 16th Project Euler Problem - Power Digit Sum - in Javascript.

The 16th Project Euler problem - Power Digit Sum - is stated as follows. 215 = 32768 and the sum of its digits is 3 + 2 + 7 + 6 + 8 = 26. What is the sum of the digits of the number 21000?

function addIntegerArrays(a = [], b = [], memo = '', carry = 0) {

    if (a.length + b.length + carry === 0) {
        return memo
    }

    // Get current number from the end of both
    const currentNumberA = a[a.length - 1] || 0;
    const currentNumberB = b[b.length - 1] || 0;

    // Sum the carry and the last digits or default to 0
    const sum = carry + currentNumberA + currentNumberB;

    const nextNumber = sum % 10;
    const nextCarry = sum > 9 ? 1 : 0;
    const nextA = a.slice(0, a.length - 1);
    const nextB = b.slice(0, b.length - 1);

    const result = `${nextNumber}${memo}`;

    return addIntegerArrays(nextA, nextB, result, nextCarry);

}

function addStrings(a, b) {
    // Convert strings to arrays.
    const preparedA = a.split('')
        .map(v => parseInt(v, 10));

    const preparedB = b.split('')
        .map(v => parseInt(v, 10));

    return addIntegerArrays(preparedA, preparedB);
}

// Create a range of length N - 1.
const result = new Array(1000 - 1).fill(0)
    // Starting at 2, hence N-1, double the value for every step.
    .reduce(memo => addStrings(memo, memo), "2")
    // Split the result into single chars.
    .split('')
    // Convert the chars to integers.
    .map(str => parseInt(str))
    // Sum the integers.
    .reduce((memo, item) => memo + item);
2018-4-2

Project Euler Problem 15 - Lattice Paths - solved with R

Solution to the 15th Project Euler Problem coded and visualised with R.

The 15th Project Euler problem - Lattice Paths - is stated as follows. Starting in the top left corner of a 2×2 grid, and only being able to move to the right and down, there are exactly 6 routes to the bottom right corner. How many such routes are there through a 20×20 grid?

Number of Lattice Paths

countRoutes <- function(grid.size){
  # Get the binomial coefficient i.e. the number of ways
  # of picking k unordered outcomes from n possibilities.
  k <- grid.size
  n <- k*2
  result <- factorial(n) / (factorial(k)*factorial(n-k))
  return(result)
}

# Call the device driver.
png(file="graph.png",width=1080,height=1080,res=120)

# Set margins.
par(mar=c(5,3,3,3))

# Create a range from 1 to 20 for the x axis.
x <- c(1:20)

# Map the range using the countRoutes function.
y <- countRoutes(x)

# Plot the graph without axes and annotations.
plot(y, type="o", col="blue", axes=FALSE, ann=FALSE, pch=16)

# Draw a box around the plot.
box()

# Add gridd lines.
abline(h=tail(y,3), v=x, col="gray", lty=3)

# Add labels for the last three markers.
text(tail(x,3),tail(y,3),labels=tail(y,3), pos=2)

# Add the X axis.
axis(1, at=x)

# Add a title for the x axis.
title(xlab= "Lattice Grid Size")

# End call to device driver.
dev.off()
2018-4-1

Finding the 10001st prime with R

Solution to the seventh Project Euler problem – 10001st prime - in R.

The seventh Project Euler problem - 10001st prime - is stated as follows. By listing the first six prime numbers: 2, 3, 5, 7, 11, and 13, we can see that the 6th prime is 13. What is the 10001st prime number?

isPrime <- function(n) {
  if (n == 2 || n == 3) {
    # Handle the edge cases 2 and 3.
    return(TRUE)                    
  } else if (n < 2 || n %% 2 == 0) {
    # Check if n is less than 2 or equal to a multiple of 2.
    return(FALSE)
  } else {
    # Factors come in pairs and one of the factors has
    # to be smaller than the square root of N.
    upper.limit <- floor(sqrt(n))
    range <- seq(3, upper.limit)

    # Check if N is divisible 
    d <- n %% (range)

    # PRIME = Nothing in range is divisble by N.
    # NOT PRIME = Something in range is divisible by N.
    return(! 0%in%d)
  }
}

# Init target prime, a counter and an empty vector.
target <- 10001
current <- 1
primes <- c()

# Find primes untill we have TARGET primes.
while (length(primes) < target) {
  primes = if (isPrime(current)) c(primes, current) else primes
  current = current + 1
}

# Print the prime we are looking for.
sprintf("The %sth prime is %s", target, primes[target])
2018-3-24

Find longest Collatz sequence using Javascript

Solution to the 14th Project Euler problem – Longest Collatz Sequence - in both Kotlin and Javascript.

The 14th Project Euler problem - Longest Collatz Sequence - is stated as follows. Which starting number, under one million, produces the longest Collatz conjecture sequence?

const memoize = fn => {
    let cache = {};
    return n => cache[n] = cache[n] ? cache[n] : fn(n);
};

const getColatzLength = memoize(n => n <= 1 ? 1 : 1 + getColatzLength(n % 2 === 0 ? n / 2 : n * 3 + 1));

const [first] = new Array(1000000).fill()
    .map((_, index) => index)
    .map(n => getColatzLength(n))
    .map((len, n) => ({len, n}))
    .sort((a, b) => b.len - a.len);
2018-3-20

Sum large numbers with Kotlin versus Javascript

Solution to the 13th Project Euler problem – Large sum - in both Kotlin and Javascript.

The 13th Project Euler problem - Large sum - is stated as follows. Work out the first ten digits of the sum of the following one-hundred 50-digit numbers.

37107287533902102798797998220837590246510135740250
46376937677490009712648124896970078050417018260538
[...]
53503534226472524250874054075591789781264330331690

This problem is trivial in Kotlin due to the built-in class BigInteger which lets you do arithmetic on arbitrary big integers.

import java.math.BigInteger

fun main(args: Array<String>) {

    var numbers = """
    37107287533902102798797998220837590246510135740250
    53503534226472524250874054075591789781264330331690
        """

    var numbersAsArray = numbers.split("\n")    
        .filter { !it.isBlank() }
        .map { it.trim() }
        .map { BigInteger(it) }
        .fold(BigInteger("0"), { sum, n -> sum.add(n) })


    println(numbersAsArray);

}

This problem is less trivial in languages where you don't have support for 1arge enough integers, e.g., Javascript. In those languages, you have to implement your own add-function.

const numbers = `
    37107287533902102798797998220837590246510135740250
    53503534226472524250874054075591789781264330331690`;


function addIntegerArrays(a = [], b = [], memo = '', carry = 0) {

    if (a.length + b.length + carry === 0) {
        return memo
    }

    // Get current number from the end of both
    const currentNumberA = a[a.length - 1] || 0;
    const currentNumberB = b[b.length - 1] || 0;

    // Sum the carry and the last digits or default to 0
    const sum = carry + currentNumberA + currentNumberB;

    const nextNumber = sum % 10;
    const nextCarry = sum > 9 ? 1 : 0;
    const nextA = a.slice(0, a.length - 1);
    const nextB = b.slice(0, b.length - 1);

    const result = `${nextNumber}${memo}`;

    return addIntegerArrays(nextA, nextB, result, nextCarry);

}

function addStrings(a, b) {
    // Convert strings to arrays.
    const preparedA = a.split('')
        .map(v => parseInt(v, 10));

    const preparedB = b.split('')
        .map(v => parseInt(v, 10));

    return addIntegerArrays(preparedA, preparedB);
}


const res = numbers.split('\n')
    .map(v => v.trim())
    .filter(v => /^[0-9]+/.test(v))
    .reduce((mem, n) => addStrings(mem, n), '0');
2018-3-18

What is a binary tree and how to invert it using Kotlin

Explanation of what a binary tree is, what it means to invert a binary tree and how to do it in Kotlin.

A binary tree is a data structure in which each node has at most two children. The children are referred to as the left child and the right child.

                    1
                  /   \    
   (Left Child)  2     3  (Right Child)
                / \
 (Left Child)  4   5  (Right Child)             

To invert, sometimes also called to mirror or reverse, a binary tree is to interchange all left and right nodes.

   (Original)     (Inverted)

       1              1
     /   \          /   \      
    2     3        3     2
   / \   / \      / \   / \
  4   5 6   7    7   6 5   4
 /                          \
8                            8

Given that a node has this shape.

class TreeNode(var value: Int = 0) {
    var left: TreeNode? = null
    var right: TreeNode? = null
}

Inverting can be done by recursivly walking through the tree and on every level change left for right and right for left.

fun invert(root: TreeNode?): TreeNode? {
    // Check if we have a node.
    if (root == null) {
        return null
    }

    // Invert the sub structures and save the result
    // to temporary variables.
    val right: TreeNode? = invert(root.right)
    val left: TreeNode? = invert(root.left)

    // Switch places.
    root.left = right
    root.right = left

    // Return the root.
    return root
}
2018-3-14

Find Highly Divisible Triangular Numbers with Kotlin

A solution for the Project Euler problem number 12. Written in the Java-like language Kotlin. Performance optimizations have primarily been made in the function for finding divisors.

Project Euler problem number 12 - Highly Divisible Triangular Number - is stated as follows.

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 number 28 has six divisors - 1, 2, 4, 7, 14, 28.

What is the value of the first triangle number to have over five hundred divisors?

fun main(args: Array<String>) {
    // Set the target to 500 divisors.
    val target = 500

    // Create a sequence that steps 1 at the time.
    val triangleNumber = generateSequence(1, { it + 1 })
            // Get the triangle number for every number.
            .map { getTriangleNumber(it) }
            // Stop when we find a triangle number with more divisors
            // than our target.
            .first { getDivisors(it).size > target }

    // Print the answer.
    println("$triangleNumber is the first triangle number with more than $target divisors.")
}

fun getTriangleNumber(n: Int): Int {
    // Triangle numbers are given by the following formula.
    // n(n+1)/2
    return (n*(n+1)/2.0).toInt()
}

fun getDivisors(n: Int): List<Int> {
    // Factors comes in pairs, you only need to iterate to the
    // square root of n and get the paired factor using n / it.
    val limit = kotlin.math.sqrt(n.toDouble()).toInt()
    return (1..limit).filter { n % it == 0 }.flatMap {
        val squaredIsN = it * it == n
        if (squaredIsN) listOf(it) else listOf(it, n / it)
    }
}

You can run and modify the code at https://try.kotlinlang.org.

2018-2-27

Find the Largest Product in a Grid with Rust

Solution for the Project Euler Problem 11 - Largest product in a grid. Written in last years most popular language, according to Github statistics, Rust.

Project Euler problem number 11 - Largest product in a grid - is stated as follows. In a defined 20x20 grid what is the greatest product of four adjacent numbers in the same direction (up, down, left, right, or diagonally)?

fn main() {
    let lens_size: usize = 4;
    let string: &str = "\
        08 02 22 97 38 15 00 40 00 75 04 05 07 78 52 12 50 77 91 08
        49 49 99 40 17 81 18 57 60 87 17 40 98 43 69 48 04 56 62 00
        81 49 31 73 55 79 14 29 93 71 40 67 53 88 30 03 49 13 36 65
        52 70 95 23 04 60 11 42 69 24 68 56 01 32 56 71 37 02 36 91
        22 31 16 71 51 67 63 89 41 92 36 54 22 40 40 28 66 33 13 80
        24 47 32 60 99 03 45 02 44 75 33 53 78 36 84 20 35 17 12 50
        32 98 81 28 64 23 67 10 26 38 40 67 59 54 70 66 18 38 64 70
        67 26 20 68 02 62 12 20 95 63 94 39 63 08 40 91 66 49 94 21
        24 55 58 05 66 73 99 26 97 17 78 78 96 83 14 88 34 89 63 72
        21 36 23 09 75 00 76 44 20 45 35 14 00 61 33 97 34 31 33 95
        78 17 53 28 22 75 31 67 15 94 03 80 04 62 16 14 09 53 56 92
        16 39 05 42 96 35 31 47 55 58 88 24 00 17 54 24 36 29 85 57
        86 56 00 48 35 71 89 07 05 44 44 37 44 60 21 58 51 54 17 58
        19 80 81 68 05 94 47 69 28 73 92 13 86 52 17 77 04 89 55 40
        04 52 08 83 97 35 99 16 07 97 57 32 16 26 26 79 33 27 98 66
        88 36 68 87 57 62 20 72 03 46 33 67 46 55 12 32 63 93 53 69
        04 42 16 73 38 25 39 11 24 94 72 18 08 46 29 32 40 62 76 36
        20 69 36 41 72 30 23 88 34 62 99 69 82 67 59 85 74 04 36 16
        20 73 35 29 78 31 90 01 74 31 49 71 48 86 81 16 23 57 05 54
        01 70 54 71 83 51 54 69 16 92 33 48 61 43 52 01 89 19 67 48";

    // Convert the multiline string to a two dimensional array.
    //  "01 02
    //  "03 04" --> [[1, 2], [3, 4]]
    let matrix = get_matrix_from_string(string);

    // Get the greatest product of N (lens_size) adjacent numbers
    // in the same direction (up, down, left, right, or diagonally) in the matrix.
    let result = (0..matrix.len()-lens_size)
        .map(|x| (0..matrix.len()-lens_size).map(|y| get_greatest_product_in_lens(&matrix, x,y, lens_size)).max())
        .max()
        .unwrap();

    println!("{:?}", result);
}

fn get_greatest_product_in_lens(matrix: &[Vec<u32>], x:usize, y:usize, lens_size:usize) -> u32 {
    // Get a sub set of the matrix based on lens size.
    // Matrix   +   Lens    =   Result
    // [            [           [
    //   [1,2,3]     [0,0]        [1,2]
    //   [4,5,6]     [1,1]        [3,4]
    //   [7,8,9]    ]           ]
    // ]
    let submatrix = get_matrix_subset(&matrix, [[x, y], [x+lens_size-1, y+lens_size-1]]);

    let array: [u32; 4] = [
        get_largest_row_product(&submatrix),
        get_largest_column_product(&submatrix),
        get_lr_diagonal_product(&submatrix),
        get_rl_diagonal_product(&submatrix)
    ];

    *array.iter().max().unwrap_or(&0)
}

// Split a string into an int array.
// "01 02" --> [1,2]
fn string_to_intarray(string: &str) -> Vec<u32> {
    // Split the string into a string array.
    // "01 02" --> ["01", "02"]
    string.split(' ')
        // Remove leading zeroes
        // ["01", "02"] --> ["1", "2"]
        .map(|s: &str| s.trim_left_matches('0'))
        // Convert the strings into integers.
        // ["1", "2"] --> [1,2]
        .map(|s: &str| s.parse().unwrap_or(0))
        // Transform the iterator into a collection.
        .collect()
}

fn get_matrix_subset(matrix: &[Vec<u32>], lens: [[usize; 2]; 2]) -> Vec<Vec<u32>> {
    // Remove rows.
    matrix[lens[0][0]..lens[1][0] + 1]
        .iter()
        // Remove columns.
        .map(|row| row[lens[0][1]..lens[1][1] + 1].to_owned())
        // Transform the iterator into a collection.
        .collect()
}

fn get_largest_row_product(matrix: &[Vec<u32>]) -> u32 {
    (0..matrix.len())
        .map(|x| (0..matrix[x].len()).fold(1, |sum, n| sum * matrix[x][n]))
        .max()
        .unwrap_or(0)
}

fn get_largest_column_product(matrix: &[Vec<u32>]) -> u32 {
    (0..matrix.len())
        .map(|x| (0..matrix[x].len()).fold(1, |sum, n| sum * matrix[n][x]))
        .max()
        .unwrap_or(0)
}

fn get_lr_diagonal_product(matrix: &[Vec<u32>]) -> u32 {
    (0..matrix.len()).fold(1, |sum, n| sum * matrix[n][n])
}

fn get_rl_diagonal_product(matrix: &[Vec<u32>]) -> u32 {
    (0..matrix.len()).fold(1, |sum, n| sum * matrix[n][matrix.len() - n - 1])
}

fn get_matrix_from_string(s: &str) -> Vec<Vec<u32>> {
    // Convert the multiline string to a two dimensional array.
    //  "01 02
    //  "03 04" --> [[1, 2], [3, 4]]
    s.split('\n')
        // Creates rows by splitting the string on newline.
        // ["01 02", "        03 04"]
        .map(|x: &str| x.trim())
        // Removes leading whitespaces.
        // ["01 02", "03 04"]
        .map(string_to_intarray)
        // Splits the row strings into integer arrays.
        // [[1, 2],[3, 4]]
        .collect()
}
2018-2-18

Summation of the First Two Million Primes with Rust

This is a one-liner if you have a function that can identify primes.

Project Euler problem number 10 - Summation of primes - is stated as follows. The sum of the primes below 10 is 2 + 3 + 5 + 7 = 17. Find the sum of all the primes below two million.

#![feature(iterator_step_by)]

fn is_prime(n: u64) -> bool {
    if n == 2 || n == 3 {
        return true;
    } else if (n < 2) || (n % 2 == 0) {
        return false;
    }
    let upper_limit = (n as f64).sqrt() as u64;
    (3..upper_limit + 1)
        .step_by(2)
        .find(|i| n % i == 0)
        .is_none()
}

fn sum_primes_below(n: u64) -> u64 {
    (0..n).filter(|x| is_prime(*x)).sum()
}

fn main() {
    let n = 2_000_000;
    let sum = sum_primes_below(n);
    println!("The Sum of the first {} primes is {}.", n, sum);
}
2018-2-18

Finding the largest product in a series with Rust

Map, map, fold, max and unwrap.

The eight Project Euler problem - Largest Product in a Series - is stated as follows. The four adjacent digits in the 1000-digit number below that have the greatest product are 9 × 9 × 8 × 9 = 5832.

73167176531330624919225119674426574742355349194934
96983520312774506326239578318016984801869478851843
85861560789112949495459501737958331952853208805511
12540698747158523863050715693290963295227443043557
66896648950445244523161731856403098711121722383113
62229893423380308135336276614282806444486645238749
30358907296290491560440772390713810515859307960866
70172427121883998797908792274921901699720888093776
65727333001053367881220235421809751254540594752243
52584907711670556013604839586446706324415722155397
53697817977846174064955149290862569321978468622482
83972241375657056057490261407972968652414535100474
82166370484403199890008895243450658541227588666881
16427171479924442928230863465674813919123162824586
17866458359124566529476545682848912883142607690042
24219022671055626321111109370544217506941658960408
07198403850962455444362981230987879927244284909188
84580156166097919133875499200524063689912560717606
05886116467109405077541002256983155200055935729725
71636269561882670428252483600823257530420752963450

Find the thirteen adjacent digits in the 1000-digit number that have the greatest product. What is the value of this product?

fn char_to_u64(c: char) -> u64 {
    // Convert the char to a digit in base 10.
    c.to_digit(10) 
        // Return the value or default to 0.
        .unwrap_or(0) as u64
}

fn product_of_string(s: &str) -> u64 {
    // Returns an iterator over the chars.
    s.chars()
        // Convert chars to u64 and get their product.
        .fold(1, |product, n| char_to_u64(n) * product)
}

fn largest_product_in_a_series(series: &str, len: usize) -> u64 {
    // Iterate from 0 to the series length minus the length
    // of the number of adjacent digits we are looking for.
    (0..series.len() - (len - 1))
        // Get the digits we are looking for as a string.
        .map(|x| &series[x..x + len])
        // Convert the chars to u64 and get their product.
        .map(|x| product_of_string(x))
        // Get the largest product.
        .max()
        // Return the value or default to 0.
        .unwrap_or(0)
}

fn main() {
    let l = 13;
    let s = "73167176531330624919225119674426574742355349194934\
             96983520312774506326239578318016984801869478851843\
             85861560789112949495459501737958331952853208805511\
             12540698747158523863050715693290963295227443043557\
             66896648950445244523161731856403098711121722383113\
             62229893423380308135336276614282806444486645238749\
             30358907296290491560440772390713810515859307960866\
             70172427121883998797908792274921901699720888093776\
             65727333001053367881220235421809751254540594752243\
             52584907711670556013604839586446706324415722155397\
             53697817977846174064955149290862569321978468622482\
             83972241375657056057490261407972968652414535100474\
             82166370484403199890008895243450658541227588666881\
             16427171479924442928230863465674813919123162824586\
             17866458359124566529476545682848912883142607690042\
             24219022671055626321111109370544217506941658960408\
             07198403850962455444362981230987879927244284909188\
             84580156166097919133875499200524063689912560717606\
             05886116467109405077541002256983155200055935729725\
             71636269561882670428252483600823257530420752963450";

    println!(
        "{} adjacent digits = {}",
        l,
        largest_product_in_a_series(s, l)
    );
}
2018-2-18

Finding the 10001st prime with Rust

If a number N has a prime factor larger than the square root of N, then it also has a prime factor smaller than square root of N.

The seventh Project Euler problem - 10001st prime - is stated as follows. By listing the first six prime numbers: 2, 3, 5, 7, 11, and 13, we can see that the 6th prime is 13. What is the 10001st prime number?

#![feature(iterator_step_by)]

fn is_prime(n: u64) -> bool {
    if n == 2 || n == 3 {
        return true;
    } else if (n < 2) || (n % 2 == 0) {
        return false;
    }
    let upper_limit = (n as f64).sqrt() as u64;
    (3..upper_limit + 1)
        .step_by(2)
        .find(|i| n % i == 0)
        .is_none()
}

fn main() {
    let n = 10_001;
    let mut primes = 0u64;
    let mut i = 0u64;

    while primes < n {
        i += 1;
        if is_prime(i) {
            primes += 1;
        }
    }

    println!("{} is prime number {}.", i, n);
}
2018-2-15

The complete list of rational numbers with Stern-Brocot and Javascript

We can get infinity minus infinity to equal any real number. Therefore, infinity subtracted from infinity is undefined and a tl:dr; not possible.

The Stern-Brocot sequence is a number sequence not too different from the Fibonacci series but created using a slightly different rule.

Just like with the Fibonacci sequence you sum the last two numbers to get the next, but then you also copy the very last digit and add that as well.

____
1, 1
// Sum the current pair 1+1=2 and copy the second number 1.
   ____
1, 1, 2, 1
// Sum the current pair 1+2=3 and copy the second number 2.
      ____
1, 1, 2, 1, 3, 2
// Sum the current pair 2+1=3 and copy the second number 1.
         ____
1, 1, 2, 1, 3, 2, 3, 1
// Sum the current pair 1+3=4 and copy the second number 3.

1, 1, 2, 1, 3, 2, 3, 1, 4, 3

The exciting thing about the Stern-Brocot sequence is that it can be used to generate every possible fraction in their most simplified form. No misses, no duplicates.

____
1, 1 = 1/1
// 1 as the numerator and 1 as the denominator.
   ____
1, 1, 2 = 1/1, 1/2
// 1 as the numerator and 2 as the denominator.
      ____
1, 1, 2, 1 = 1/1, 1/2, 2/1
// 2 as the numerator and 1 as the denominator.
         ____
1, 1, 2, 1, 3 = 1/1, 1/2, 2/1, 1/3 
// 1 as the numerator and 3 as the denominator.

Below is a Stern-Brocot generator written in Javascript as a recursive function.

const addSternBrocot = (seq, n) => {
  const n1 = seq[seq.length/2-1];
  const n2 = seq[seq.length/2];
  const next = seq.concat(n1+n2, n2);
  return --n ? addSternBrocot(next, n) : next;
};

addSternBrocot([1,1], 25); // ==> First 25 numbers.
2018-2-14

Project Euler number six solved with Rust

Solved by taking the 100th triangular number minus a 2000-year-old formula by Archimedes.

The sixth Project Euler problem - Sum Square Difference - is stated as follows.

The sum of the squares of the first ten natural numbers is 385. The square of the sum of the first ten natural numbers is 3025. 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.

Find the difference between the sum of the squares of the first one hundred natural numbers and the square of the sum.

fn sum_square_difference(n: u32) -> u32 {
    (n*(n+1)/2).pow(2) - ((n*(n+1)*(2*n+1))/6)
}

fn main() {
    let n = 100;
    println!(
        "The difference is {} with n = {}.",
        sum_square_difference(n),
        n,
    );
}
2018-2-14

Smallest positive number that is evenly divisible by all of the numbers from 1 to 20 with Rust

The smallest common multiple of 1, 2 and 3 is the smallest common multiple of 1 and the smallest common multiple of 2 and 3.

The fifth Project Euler problem - Smallest Multiple - is stated as follows.

2520 is the smallest number that can be divided by each of the numbers from 1 to 10 without any remainder. What is the smallest positive number that is evenly divisible by all of the numbers from 1 to 20?

fn lcm(a: u64, b: u64) -> u64 {
    if a <= 0 || b <= 0 {
        0
    } else {
        a * b / gcd(a, b)
    }
}

fn gcd(a: u64, b: u64) -> u64 {
    if b <= 0 {
        a
    } else {
        gcd(b, a % b)
    }
}

fn main() {
    let n = 20;
    println!(
        "{} is evenly divisible by all of the numbers from 1 to {}.",
        (1..n).fold(n, |a, b| lcm(a, b)),
        n
    );
}
2018-1-31

Finding the largest palindrome product with Rust

An easy way to check whether a number is a palindrome is to convert the number to a string, reverse it and check for equality. Neither string conversions or pure brute-force will lead to an optimal solution.

The fourth Project Euler problem - Largest Palindrome Product - is stated as follows.

A palindromic number reads the same both ways. The largest palindrome made from the product of two 2-digit numbers is 9009 = 91 × 99. Find the largest palindrome made from the product of two 3-digit numbers.

Pure brute-force with string conversions

fn is_palindrome(product: u32) -> bool {
    let str = product.to_string();
    let rev_str = str.chars().rev().collect::<String>();
    str == rev_str
}

fn get_largest_palindrome(n: u32) -> u32 {
    let start = 1 * u32::pow(10, n - 1);
    let end = start * 10;

    (start..end)
        .flat_map(|x| (x..end).map(move |y| x * y))
        .filter(|&p| is_palindrome(p))
        .max()
        .unwrap()
}

fn main() {
    println!(
        "Largest palindrome n=3: {}",
        get_largest_palindrome(3)
    );
}
N Execution Time Improvement
2 9.8 ms Baseline
3 933.5 ms Baseline
4 103286.6 ms Baseline

Pure brute-force with modulo operation

fn is_palindrome(n: u32, num: u32) -> bool {
    !(0..n).any(|a| (num / 10u32.pow(2 * n - 1 - a)) % 10 != (num / 10u32.pow(a)) % 10)
}

fn get_largest_palindrome(n: u32) -> u32 {
    let start = 1 * u32::pow(10, n - 1);
    let end = start * 10;

    (start..end)
        .flat_map(|x| (x..end).map(move |y| x * y))
        .filter(|&p| is_palindrome(n, p))
        .max()
        .unwrap()
}

fn main() {
    println!(
        "Largest palindrome n=3: {}",
        get_largest_palindrome(3)
    );
}
N Execution Time Improvement
2 1.4 ms 7X
3 104.8 ms 9X
4 10568.0 ms 10X

More efficient algorithm with string conversions

fn is_palindrome(product: u32) -> bool {
    let str = product.to_string();
    let rev_str = str.chars().rev().collect::<String>();
    str == rev_str
}

fn get_largest_palindrome(n: u32) -> u32 {
    let a_start = 10u32.pow(n - 1) / 11 + 1;
    let a_end = 10u32.pow(n) / 11 + 1;
    let b_start = 10u32.pow(n - 1);
    let b_end = 10u32.pow(n);

    (a_start..a_end)
        .flat_map(|x| (b_start..b_end).map(move |y| 11 * x * y))
        .filter(|&prod| prod > 10u32.pow(2 * n - 1) && is_palindrome(prod))
        .max()
        .unwrap()
}

fn main() {
    println!(
        "Largest palindrome n=3: {}",
        get_largest_palindrome(3)
    );
}
N Execution Time Improvement
2 1.3 ms 8X
3 125.8 ms 7X
4 14348.6 ms 7X

More efficient algorithm with modulo operation

fn is_palindrome(n: u32, num: u32) -> bool {
    !(0..n).any(|a| (num / 10u32.pow(2 * n - 1 - a)) % 10 != (num / 10u32.pow(a)) % 10)
}

fn get_largest_palindrome(n: u32) -> u32 {
    let a_start = 10u32.pow(n - 1) / 11 + 1;
    let a_end = 10u32.pow(n) / 11 + 1;
    let b_start = 10u32.pow(n - 1);
    let b_end = 10u32.pow(n);

    (a_start..a_end)
        .flat_map(|x| (b_start..b_end).map(move |y| 11 * x * y))
        .filter(|&prod| prod > 10u32.pow(2 * n - 1) && is_palindrome(n, prod))
        .max()
        .unwrap()
}

fn main() {
    println!(
        "Largest palindrome n=3: {}",
        get_largest_palindrome(3)
    );
}
N Execution Time Improvement
2 0.3 ms 30X
3 20.0 ms 47X
4 2149.2 ms 48X

The optimized versions of my pure brute-force algorithm is more or less a product of the friendly community over at the Rust Discrord. If you need help with the programming language Rust, I couldn't recommend the Discord channel enough.

2018-1-21

Get the largest prime factor with Reason

The simplest solution to this problem is to list all prime factors, sort, reverse and print the first value in the list.

The third Project Euler problem - Largest Prime Factor - is stated as follows. What is the largest prime factor of the number 600851475143?

let getPrimeFactors = (num) => {
  let rec next = (num, i, res) =>
    switch (num > 1, num mod i === 0) {
    | (true, true) => next(num / i, 2, [i, ...res])
    | (true, false) => next(num, i + 1, res)
    | _ => res
    };
  next(num, 2, [])
};

getPrimeFactors(600851475143)
    |> List.sort(compare)
    |> List.rev
    |> List.hd
    |> print_int;
2018-1-18

Find the sum of all even Fibonacci numbers below four million with OCaml

The sum of all even Fibonacci numbers below four million is easily found with a recursive function and a switch.

The second Project Euler problem - Even Fibonacci Numbers - is stated as follows - By considering the terms in the Fibonacci sequence whose values do not exceed four million, find the sum of the even-valued terms.

let fib target =
  let rec next x y sum =
    match (y mod 2 == 0, y < target) with
    | (true, true) -> next y (x + y) (sum + y)
    | (false,true) -> next y (x + y) sum
    | _ -> sum in
  next 1 1 0

let _ = print_int (fib 4000000)
2018-1-18

Finding the sum of all multiples of 3 or 5 below 1000 with OCaml

The easiest way to find all the multiples of 3 or 5 below 1000 with Ocaml is to use a recursive function with a switch.

The first Project Euler problem - Multiples of 3 and 5 - is stated as follows. Find the sum of all the multiples of 3 or 5 below 1000.

(* Recursive switch *)
let multiples target =
  let rec next x sum =
    match (x mod 3 = 0 || x mod 5 = 0, x > 0) with
    | (false, true) -> next (x - 1) sum
    | (true, true) -> next (x - 1) (sum + x)
    | _ -> sum in
  next (target - 1) 0
(* Recursive switch + one function argument  *)
let multiples target =
  let rec next sum =
    function
    | 0 -> sum
    | n when n mod 3 = 0 || n mod 5 = 0 -> next (sum + n) (n - 1)
    | n -> next sum (n - 1) in
  next 0 target
2018-1-17

Special Pythagorean Triplet solved with Reason

Brute force is a good enough approach when finding small Pythagorean triplets, but unsuitable with greater target sums. When brute force takes minutes, a smarter approach will take microseconds.

The ninth Project Euler problem - Special Pythagorean Triplet - is stated as follows. A Pythagorean triplet is a set of three natural numbers, a < b < c, for which a2 + b2 = c2. There exists exactly one Pythagorean triplet for which a + b + c = 1000. Find the product abc.

Solving with brute force

The easiest way to solve this problem is by brute force and iterate through all possible values of a, b and c where their sum equals 1000.

let findTriplet = (sum) =>
  for (a in 1 to sum) {
    for (b in a to sum) {
      let c = sum - a - b;
      if (a * a + b * b === c * c) {
        print_int(a * b * c)
      }
    }
  };

findTriplet(1000); /* ==> 31875000 */

One of the first noticeable problems with the brute force approach is that it does not scale very well. When finding a triplet with the sum 1000 takes milliseconds it takes minutes to find a triplet with the sum 100000.

Sum Execution Time (s)
1000 0.009845
10000 0.964073
100000 92.694244

Brute force with smaller bounds

By applying some wetware to this problem, we can make the bounds for the loops smaller.

We can deduce that a < 1000/3 since if a is not less than 1000/3, then b and c can't be greater than a.

We can also deduce that b < 1000/2 since if b is greater than 1000/2, c can't be greater than b.

let findTriplet = (sum) => {
  let rec next = (a, b) => {
    let c = sum - a - b;
    switch (a * a + b * b === c * c) {
    | true => a * b * c
    | false when b < sum / 2 => next(a, b + 1)
    | false when a < sum / 3 => next(a + 1, a + 2)
    | _ => 0
    }
  };
  next(1, 2)
};

findTriplet(1000); /* ==> 31875000 */

The shorter loops significantly speed up the solution, but unfortunately, it still has the same scaling behavior.

Sum Execution Time (s)
1000 0.002457
10000 0.256682
100000 26.623702

Multiple possible triplets

Another flaw with the code above, when we start looking for a more general solution, is that we cannot be sure that there exists only one Pythagorean triplet where a + b + c = 1000 since it only returns the first triplet it encounters.

let findTriplets = (sum) => {
  let rec next = (a, b, memo) => {
    let c = sum - a - b;
    let nextMemo = a * a + b * b === c * c ? [a * b * c, ...memo] : memo;
    switch (b < sum / 2, a < sum / 3) {
    | (true, true) => next(a, b + 1, nextMemo)
    | (false, true) => next(a + 1, a + 2, nextMemo)
    | _ => memo
    }
  };
  next(1, 2, [])
};

findTriplets(990); /* ==> 34178760, 29452500, 28030860, 19645560 */

Adding that functionality takes a toll on the execution time and makes it even more apparent that this approach is flawed for greater sums.

Sum Execution Time (s)
1000 0.005372
10000 0.581307
100000 54.624186

The holy grail

After you turn in your answer to Project Euler you get access to a solution using parametrisation of Pythagorean triplets. I did this Reason implementation below just to see how it compare to my solitions - minutes became microseconds.

let getTriplet = (sum) => {
  let s2 = sum / 2;
  let mlimit = ceil(sqrt(float_of_int(s2))) -. 1.0;
  let rec gcd = (a, b) => b === 0 ? a : gcd(b, a mod b);
  let rec rss = (sm) => sm mod 2 === 0 ? rss(sm / 2) : sm;
  let search = (m, s2, memo) => {
    let sm = rss(s2 / m);
    let rec helper = (k) =>
      if (k < 2 * m && k <= sm) {
        if (sm mod k === 0 && gcd(k, m) === 1) {
          let d = s2 / (k * m);
          let n = k - m;
          let a = d * (m * m - n * n);
          let b = 2 * d * m * n;
          let c = d * (m * m + n * n);
          [a * b * c, ...memo]
        } else {
          helper(k + 2)
        }
      } else {
        memo
      };
    helper(m + m mod 2 + 1)
  };
  let rec helper = (m, memo) =>
    switch (s2 mod m === 0, m < int_of_float(mlimit)) {
    | (true, true) => helper(m + 1, search(m, s2, memo))
    | (true, false) => search(m, s2, memo)
    | _ => helper(m + 1, memo)
    };
  helper(2, [])
};

getTriplet(990) |> Array.of_list |> Array.map(string_of_int) |> Array.map(print_endline);
Sum Execution Time (s)
1000 0.000006
10000 0.000010
100000 0.000018
2017-12-29

New Year's Resolution 2018

My new year's resolution for 2018 is to become a professional artist, tattoo myself and dabble in a new programming language every month.

This year is the year I will become a professional artist. My style will be computer-generated and generative, and my tools will be Javascript, Gcode, and an XY plotter. I will consider the resolution kept if I sell at least one piece of art to a stranger.

During the year, I will also equip my plotter with a tattoo needle and create and get my first tattoo. It will probably be something minimalistic and practical like a ruler.

Last but not least, I want to broaden my knowledge of programming languages. I plan to dabble in a new programming language every month. Currently, the list is Ocaml/Reason, Rust, Python, Ruby, C, Go, R, Clojure, Wasm, Haskell, Kotlin, and C#.

,
Jan Järfalk — User experienced esthete, technician, geek and web worker. Aspiring artist and recreational mathematician. I indulge and travel plenty. Constraints are good.