Add more solutions

Added solutions for problems 61, 62, 63, 65 and 65 in C and python.
This commit is contained in:
daniele 2019-09-28 16:22:03 +02:00
parent 68b33e946c
commit 8e8ac9931a
Signed by: fuxino
GPG Key ID: 6FE25B4A3EE16FDA
13 changed files with 1067 additions and 0 deletions

195
C/p061.c Normal file
View File

@ -0,0 +1,195 @@
/* Triangle, square, pentagonal, hexagonal, heptagonal, and octagonal numbers are all figurate (polygonal) numbers and are generated
* by the following formulae:
*
* Triangle P3,n=n(n+1)/2 1, 3, 6, 10, 15, ...
* Square P4,n=n2 1, 4, 9, 16, 25, ...
* Pentagonal P5,n=n(3n1)/2 1, 5, 12, 22, 35, ...
* Hexagonal P6,n=n(2n1) 1, 6, 15, 28, 45, ...
* Heptagonal P7,n=n(5n3)/2 1, 7, 18, 34, 55, ...
* Octagonal P8,n=n(3n2) 1, 8, 21, 40, 65, ...
*
* The ordered set of three 4-digit numbers: 8128, 2882, 8281, has three interesting properties.
*
* The set is cyclic, in that the last two digits of each number is the first two digits of the next number (including the last number with the first).
* Each polygonal type: triangle (P3,127=8128), square (P4,91=8281), and pentagonal (P5,44=2882), is represented by a different number in the set.
* This is the only set of 4-digit numbers with this property.
*
* Find the sum of the only ordered set of six cyclic 4-digit numbers for which each polygonal type: triangle, square, pentagonal,
* hexagonal, heptagonal, and octagonal, is represented by a different number in the set.*/
#include <stdio.h>
#include <stdlib.h>
#include <time.h>
int find_set(int step, int *sum);
int polygonal[6][10000] = {0};
int chain[6], flags[6] = {0};
int main(int argc, char **argv)
{
int i, n, sum = 0;
double elapsed;
struct timespec start, end;
clock_gettime(CLOCK_MONOTONIC, &start);
i = 1;
n = 1;
/* Generate all triangle numbers smaller than 10000*/
do
{
polygonal[0][n] = 1;
i++;
n = i * (i + 1) / 2;
}while(n < 10000);
i = 1;
n = 1;
/* Generate all square numbers smaller than 10000.*/
do
{
polygonal[1][n] = 1;
i++;
n = i * i;
}while(n < 10000);
i = 1;
n = 1;
/* Generate all pentagonal numbers smaller than 10000.*/
do
{
polygonal[2][n] = 1;
i++;
n = i * (3 * i - 1) / 2;
}while(n < 10000);
i = 1;
n = 1;
/* Generate all hexagonal numbers smaller than 10000.*/
do
{
polygonal[3][n] = 1;
i++;
n = i * (2 * i - 1);
}while(n < 10000);
i = 1;
n = 1;
/* Generate all heptagonal numbers smaller than 10000.*/
do
{
polygonal[4][n] = 1;
i++;
n = i * (5 * i - 3) / 2;
}while(n < 10000);
i = 1;
n = 1;
/* Generate all octagonal numbers smaller than 10000.*/
do
{
polygonal[5][n] = 1;
i++;
n = i * (3 * i - 2);
}while(n < 10000);
/* Find the requested set of numbers.*/
if(find_set(0, &sum) == 0)
{
printf("Set not found\n");
}
clock_gettime(CLOCK_MONOTONIC, &end);
elapsed = (end.tv_sec - start.tv_sec) + (double)(end.tv_nsec - start.tv_nsec) / 1000000000;
printf("Project Euler, Problem 61\n");
printf("Answer: %d\n", sum);
printf("Elapsed time: %.9lf seconds\n", elapsed);
return 0;
}
/* Recursive function to find the required set. It finds a polygonal number,
* check if it can be part of the chain, then use recursion to find the next
* number. If a solution can't be found with the current numbers, it uses
* backtracking and tries the next polygonal number.*/
int find_set(int step, int *sum)
{
int i, j;
/* Use one polygonal number per type, starting from triangular.*/
for(i = 0; i < 6; i++)
{
/* If the current type has not been used yet, try it.*/
if(!flags[i])
{
/* Set a flag to record that the current polygonal type has been used.*/
flags[i] = 1;
/* Start from 1010 because numbers finishing with 00, 01, ..., 09 can't
* be part of the chain.*/
for(j = 1010; j < 10000; j++)
{
/* If the number doesn't finish with 00, 01, ..., 09 and is poligonal,
* try adding it to the chain and add its value to the total sum.*/
if(j % 100 > 9 && polygonal[i][j])
{
/* If it's the first number, just add it as first step in the chain.*/
if(step == 0)
{
chain[step] = j;
*sum += j;
/* Recursively try to add other numbers to the chain. If a solution
* is found, return 1.*/
if(find_set(step+1, sum))
{
return 1;
}
/* If a solution was not found, backtrack, subtracting the value of
* the number from the total.*/
*sum -= j;
}
/* If this is the last step and the current number can be added to the chain,
* add it, update the sum and return 1. A solution has been found.*/
else if(step == 5 && j % 100 == chain[0] / 100 && j / 100 == chain[step-1] % 100)
{
chain[step] = j;
*sum += j;
return 1;
}
/* For every other step, add the number to the chain if possible, then recursively
* try to add other numbers.*/
else if(step < 5 && j / 100 == chain[step-1] % 100)
{
chain[step] = j;
*sum += j;
if(find_set(step+1, sum))
{
return 1;
}
/* If a solution was not found, backtrack.*/
*sum -= j;
}
}
}
/* Remove the flag for the current polygonal type.*/
flags[i] = 0;
}
}
return 0;
}

114
C/p062.c Normal file
View File

@ -0,0 +1,114 @@
/* The cube, 41063625 (345^3), can be permuted to produce two other cubes: 56623104 (384^3) and 66430125 (405^3).
* In fact, 41063625 is the smallest cube which has exactly three permutations of its digits which are also cube.
*
* Find the smallest cube for which exactly five permutations of its digits are cube.*/
#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#define N 10000
int check_digits(long int a, long int b);
int count_digits(long int a);
int main(int argc, char **argv)
{
int count;
long int i, j, cubes[N];
double elapsed;
struct timespec start, end;
clock_gettime(CLOCK_MONOTONIC, &start);
/* Calculate i^3 for all i smaller than 10000.*/
for(i = 0; i < N; i++)
{
cubes[i] = i * i * i;
}
/* For each cube, check if there are four other cubes which are also
* a permutation of the first cube.*/
for(i = 0; i < N - 5; i++)
{
count = 1;
/* Stop when the limit has been reached, when 5 values have been found or
* when j^3 has more digits than i^3 (if they don't have the same digits,
* they can't be permutations).*/
for(j = i + 1; j < N && count_digits(cubes[j]) == count_digits(cubes[i]); j++)
{
if(check_digits(cubes[i], cubes[j]))
{
count++;
}
if(count == 5)
{
break;
}
}
if(count == 5)
{
break;
}
}
clock_gettime(CLOCK_MONOTONIC, &end);
elapsed = (end.tv_sec - start.tv_sec) + (double)(end.tv_nsec - start.tv_nsec) / 1000000000;
printf("Project Euler, Problem 62\n");
printf("Answer: %ld\n", cubes[i]);
printf("Elapsed time: %.9lf seconds\n", elapsed);
return 0;
}
int count_digits(long int a)
{
int count = 0;
if(a == 0)
{
return 1;
}
while(a > 0)
{
count++;
a /= 10;
}
return count;
}
int check_digits(long int a, long int b)
{
int i;
int digits1[10] = {0}, digits2[10] = {0};
while(a > 0)
{
digits1[a%10]++;
a /= 10;
}
while(b > 0)
{
digits2[b%10]++;
b /= 10;
}
for(i = 0; i < 10; i++)
{
if(digits1[i] != digits2[i])
{
return 0;
}
}
return 1;
}

81
C/p063.c Normal file
View File

@ -0,0 +1,81 @@
/* The 5-digit number, 16807=7^5, is also a fifth power. Similarly, the 9-digit number, 134217728=8^9, is a ninth power.
*
* How many n-digit positive integers exist which are also an nth power?*/
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>
#include <gmp.h>
int count_digits(mpz_t n);
int main(int argc, char **argv)
{
int i, j, count = 0, finished = 0;
double elapsed;
struct timespec start, end;
mpz_t p;
mpz_init(p);
clock_gettime(CLOCK_MONOTONIC, &start);
i = 1;
while(!finished)
{
/* When j=10, j^i will have i+1 digits (e.g if i=3, 10^3=1000).*/
for(j = 1; j < 10; j++)
{
mpz_ui_pow_ui(p, j, i);
if(count_digits(p) == i)
{
count++;
}
}
/* When 9^i has less than i digits, all the numbers have been found.*/
if(count_digits(p) < i)
{
finished = 1;
}
i++;
}
mpz_clear(p);
clock_gettime(CLOCK_MONOTONIC, &end);
elapsed = (end.tv_sec - start.tv_sec) + (double)(end.tv_nsec - start.tv_nsec) / 1000000000;
printf("Project Euler, Problem 63\n");
printf("Answer: %d\n", count);
printf("Elapsed time: %.9lf seconds\n", elapsed);
return 0;
}
int count_digits(mpz_t n)
{
int count = 0;
mpz_t tmp;
if(mpz_cmp_ui(n, 0) == 0)
{
return 1;
}
mpz_init_set(tmp, n);
while(mpz_cmp_ui(tmp, 0))
{
mpz_tdiv_q_ui(tmp, tmp, 10);
count++;
}
return count;
}

97
C/p064.c Normal file
View File

@ -0,0 +1,97 @@
/* All square roots are periodic when written as continued fractions and can be written in the form:
*
* N=a0+1/(a1+1/(a2+1/(a3+
*
* For example, let us consider 23:
*
* 23=4+234=4+1/(1/(234))=4+1/(1+(233)/7)
*
* If we continue we would get the following expansion:
*
* 23=4+1/(1+1/(3+1/(1+1/(8+
*
* The process can be summarised as follows:
* a0=4,1/(234)=(23+4)/7=1+(233)/7
* a1=1,7/(233)=7(23+3)/14=3+(233)/2
* a2=3,2/(233)=2(23+3)/14=1+(234)/7
* a3=1,7/(234)=7(23+4)/7=8+234
* a4=8,1/(234)=(23+4)/7=1+(233)/7
* a5=1,7/(233)=7(23+3)/14=3+(233)/2
* a6=3,2/(233)=2(23+3)/14=1+(234)/7
* a7=1,7/(234)=7(23+4)/7=8+234
*
* It can be seen that the sequence is repeating. For conciseness, we use the notation 23=[4;(1,3,1,8)], to indicate that the block (1,3,1,8)
* repeats indefinitely.
*
* The first ten continued fraction representations of (irrational) square roots are:
*
* 2=[1;(2)], period=1
* 3=[1;(1,2)], period=2
* 5=[2;(4)], period=1
* 6=[2;(2,4)], period=2
* 7=[2;(1,1,1,4)], period=4
* 8=[2;(1,4)], period=2
* 10=[3;(6)], period=1
* 11=[3;(3,6)], period=2
* 12=[3;(2,6)], period=2
* 13=[3;(1,1,1,1,6)], period=5
*
* Exactly four continued fractions, for N13, have an odd period.
*
* How many continued fractions for N10000 have an odd period?*/
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>
#include "projecteuler.h"
int is_square(int n);
int main(int argc, char **argv)
{
int i, count = 0;
double elapsed;
struct timespec start, end;
clock_gettime(CLOCK_MONOTONIC, &start);
for(i = 2; i <= 10000; i++)
{
/* Perfect squares are obviously not represented as continued fractions.
* For all other numbers, calculate their period and check if it's odd.*/
if(!is_square(i) && period_cf(i) % 2 != 0)
{
count++;
}
}
clock_gettime(CLOCK_MONOTONIC, &end);
elapsed = (end.tv_sec - start.tv_sec) + (double)(end.tv_nsec - start.tv_nsec) / 1000000000;
printf("Project Euler, Problem 64\n");
printf("Answer: %d\n", count);
printf("Elapsed time: %.9lf seconds\n", elapsed);
return 0;
}
int is_square(int n)
{
int m;
double p;
p = sqrt(n);
m = p;
if(p == m)
{
return 1;
}
else
{
return 0;
}
}

88
C/p065.c Normal file
View File

@ -0,0 +1,88 @@
/* The square root of 2 can be written as an infinite continued fraction.
*
* 2=1+1/(2+1/(2+1/(2+1/(2+...
*
* The infinite continued fraction can be written, 2=[1;(2)], (2) indicates that 2 repeats ad infinitum. In a similar way, 23=[4;(1,3,1,8)].
* It turns out that the sequence of partial values of continued fractions for square roots provide the best rational approximations.
* Let us consider the convergents for 22.
*
* 1+1/2=3/2
* 1+1/(2+1/2)=7/5
* 1+1/(2+1/(2+1/2))=17/12
* 1+1/(2+1/(2+1/(2+1/2)))=41/29
*
* Hence the sequence of the first ten convergents for 2 are:
*
* 1,3/2,7/5,17/12,41/29,99/70,239/169,577/408,1393/985,3363/2378,...
*
* What is most surprising is that the important mathematical constant,
*
* e=[2;1,2,1,1,4,1,1,6,1,...,1,2k,1,...].
*
* The first ten terms in the sequence of convergents for e are:
*
* 2,3,8/3,11/4,19/7,87/32,106/39,193/71,1264/465,1457/536,...
*
* The sum of digits in the numerator of the 10th convergent is 1+4+5+7=17.
*
* Find the sum of digits in the numerator of the 100th convergent of the continued fraction for e.*/
#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <gmp.h>
int main(int argc, char **argv)
{
int i, ai[3] = {1, 2, 1}, count = 4, sum = 0;
mpz_t n0, n1, n2;
double elapsed;
struct timespec start, end;
clock_gettime(CLOCK_MONOTONIC, &start);
mpz_init_set_ui(n0, 3);
mpz_init_set_ui(n1, 8);
mpz_init_set_ui(n2, 11);
/* For a continued fractions [a_0; a_1, a_2, ...], the numerator of the
* next convergent N_n=a_n*N_(n-1)+N_(n-2). The first three values for e are
* 3, 8 and 11, the next ones are easily calculated, considering that a_n
* follows a simple pattern:
* a_1=1, a_2=2, a_3=1
* a_4=1, a_5=4, a_6=1
* a_7=1, a_8=6, a_9=1
* and so on.*/
while(count < 100)
{
ai[1] += 2;
for(i = 0; i < 3 && count < 100; i++)
{
mpz_set(n0, n1);
mpz_set(n1, n2);
mpz_mul_ui(n2, n1, ai[i]);
mpz_add(n2, n2, n0);
count++;
}
}
/* Sum the digits.*/
while(mpz_cmp_ui(n2, 0))
{
sum += mpz_tdiv_q_ui(n2, n2, 10);
}
mpz_clears(n0, n1, n2, NULL);
clock_gettime(CLOCK_MONOTONIC, &end);
elapsed=(end.tv_sec-start.tv_sec)+(double)(end.tv_nsec-start.tv_nsec)/1000000000;
printf("Project Euler, Problem 65\n");
printf("Answer: %d\n", sum);
printf("Elapsed time: %.9lf seconds\n", elapsed);
return 0;
}

View File

@ -466,3 +466,35 @@ int is_pentagonal(long int n)
return 0;
}
}
/* Function implementing the iterative algorithm taken from Wikipedia
* to find the continued fraction for sqrt(S). The algorithm is as
* follows:
*
* m_0=0
* d_0=1
* a_0=floor(sqrt(n))
* m_(n+1)=d_n*a_n-m_n
* d_(n+1)=(S-m_(n+1)^2)/d_n
* a_(n+1)=floor((sqrt(S)+m_(n+1))/d_(n+1))=floor((a_0+m_(n+1))/d_(n+1))
* if a_i=2*a_0, the algorithm ends.*/
int period_cf(int n)
{
int mn = 0, mn1, dn = 1, dn1, a0, an, an1, count = 0;
a0 = floor(sqrt(n));
an = a0;
do
{
mn1 = dn * an - mn;
dn1 = (n - mn1 * mn1) / dn;
an1 = floor((a0+mn1)/dn1);
mn = mn1;
dn = dn1;
an = an1;
count++;
}while(an != 2 * a0);
return count;
}

View File

@ -19,5 +19,6 @@ void quick_sort(void **array, int l, int r, int (*cmp)(void *lv, void *rv));
int next_permutation(void **perm, int n, int (*cmp)(void *a, void *b));
int is_pandigital(int value, int n);
int is_pentagonal(long int n);
int period_cf(int n);
#endif

180
Python/p061.py Normal file
View File

@ -0,0 +1,180 @@
#!/usr/bin/python3
# Triangle, square, pentagonal, hexagonal, heptagonal, and octagonal numbers are all figurate (polygonal) numbers and are generated
# by the following formulae:
#
# Triangle P3,n=n(n+1)/2 1, 3, 6, 10, 15, ...
# Square P4,n=n2 1, 4, 9, 16, 25, ...
# Pentagonal P5,n=n(3n1)/2 1, 5, 12, 22, 35, ...
# Hexagonal P6,n=n(2n1) 1, 6, 15, 28, 45, ...
# Heptagonal P7,n=n(5n3)/2 1, 7, 18, 34, 55, ...
# Octagonal P8,n=n(3n2) 1, 8, 21, 40, 65, ...
#
# The ordered set of three 4-digit numbers: 8128, 2882, 8281, has three interesting properties.
#
# The set is cyclic, in that the last two digits of each number is the first two digits of the next number (including the last number with the first).
# Each polygonal type: triangle (P3,127=8128), square (P4,91=8281), and pentagonal (P5,44=2882), is represented by a different number in the set.
# This is the only set of 4-digit numbers with this property.
#
# Find the sum of the only ordered set of six cyclic 4-digit numbers for which each polygonal type: triangle, square, pentagonal,
# hexagonal, heptagonal, and octagonal, is represented by a different number in the set.
from numpy import zeros
from timeit import default_timer
# Recursive function to find the required set. It finds a polygonal number,
# check if it can be part of the chain, then use recursion to find the next
# number. If a solution can't be found with the current numbers, it uses
# backtracking and tries the next polygonal number.
def find_set(step):
global polygonal
global chain
global flags
global sum_
# Use one polygonal number per type, starting from triangular.
for i in range(6):
# If the current type has not been used yet, try it.
if flags[i] == 0:
# Set a flag to record that the current polygonal type has been used.
flags[i] = 1
# Start from 1010 because numbers finishing with 00, 01, ..., 09 can't
# be part of the chain.
for j in range(1010, 10000):
# If the number doesn't finish with 00, 01, ..., 09 and is poligonal,
# try adding it to the chain and add its value to the total sum.
if j % 100 > 9 and polygonal[i, j] == 1:
# If it's the first number, just add it as first step in the chain.
if step == 0:
chain[step] = j
sum_ = sum_ + j
# Recursively try to add other numbers to the chain. If a solution
# is found, return 1.
if find_set(step+1):
return 1
# If a solution was not found, backtrack, subtracting the value of
# the number from the total.
sum_ = sum_ - j
# If this is the last step and the current number can be added to the chain,
# add it, update the sum and return 1. A solution has been found.
elif step == 5 and j % 100 == chain[0] // 100 and j // 100 == chain[step-1] % 100:
chain[step] = j
sum_ = sum_ + j
return 1
# For every other step, add the number to the chain if possible, then recursively
# try to add other numbers.
elif step < 5 and j // 100 == chain[step-1] % 100:
chain[step] = j
sum_ = sum_ + j
if find_set(step+1):
return 1
# If a solution was not found, backtrack.
sum_ = sum_ - j
# Remove the flag for the current polygonal type.
flags[i] = 0
return 0
def main():
start = default_timer()
global polygonal
global chain
global flags
global sum_
polygonal = zeros((6, 10000), int)
chain = [0] * 6
flags = [0] * 6
sum_ = 0
i = 1
n = 1
# Generate all triangle numbers smaller than 10000
while True:
polygonal[0, n] = 1
i = i + 1
n = i * (i + 1) // 2
if n >= 10000:
break
i = 1
n = 1
# Generate all square numbers smaller than 10000
while True:
polygonal[1, n] = 1
i = i + 1
n = i * i
if n >= 10000:
break
i = 1
n = 1
# Generate all pentagonal numbers smaller than 10000
while True:
polygonal[2, n] = 1
i = i + 1
n = i * (3 * i - 1) // 2
if n >= 10000:
break
i = 1
n = 1
# Generate all hexagonal numbers smaller than 10000
while True:
polygonal[3, n] = 1
i = i + 1
n = i * (2 * i - 1)
if n >= 10000:
break
i = 1
n = 1
# Generate all heptagonal numbers smaller than 10000
while True:
polygonal[4, n] = 1
i = i + 1
n = i * (5 * i - 3) // 2
if n >= 10000:
break
i = 1
n = 1
# Generate all octagonal numbers smaller than 10000
while True:
polygonal[5, n] = 1
i = i + 1
n = i * (3 * i - 2)
if n >= 10000:
break
# Find the requested set of numbers
if find_set(0) == 0:
print('Set not found')
end = default_timer()
print('Project Euler, Problem 61')
print('Answer: {}'.format(sum_))
print('Elapsed time: {:.9f} seconds'.format(end - start))
if __name__ == '__main__':
main()

53
Python/p062.py Normal file
View File

@ -0,0 +1,53 @@
#!/usr/bin/python
# The cube, 41063625 (345^3), can be permuted to produce two other cubes: 56623104 (384^3) and 66430125 (405^3).
# In fact, 41063625 is the smallest cube which has exactly three permutations of its digits which are also cube.
#
# Find the smallest cube for which exactly five permutations of its digits are cube.
from numpy import zeros
from timeit import default_timer
def main():
start = default_timer()
N = 10000
cubes = zeros(N, int)
# Calculate i^3 for all i smaller than 10000
for i in range(N):
cubes[i] = i * i * i
# For each cube, check if there are four other cubes which are also
# a permutation of the first cube.
for i in range(N-5):
count = 1
# Stop when the limit has been reached, when 5 values have been found or
# when j^3 has more digits than i^3 (if they don't have the same digits,
# they can't be permutations).
j = i + 1
while j < N and len(str(cubes[j])) == len(str(cubes[i])):
if ''.join(sorted(str(cubes[i]))) == ''.join(sorted(str(cubes[j]))):
count = count + 1
if count == 5:
break
j = j + 1
if count == 5:
break
end = default_timer()
print('Project Euler, Problem 62')
print('Answer: {}'.format(cubes[i]))
print('Elapsed time: {:.9f} seconds'.format(end - start))
if __name__ == '__main__':
main()

38
Python/p063.py Normal file
View File

@ -0,0 +1,38 @@
#!/usr/bin/python
# The 5-digit number, 16807=7^5, is also a fifth power. Similarly, the 9-digit number, 134217728=8^9, is a ninth power.
#
# How many n-digit positive integers exist which are also an nth power?
from timeit import default_timer
def main():
start = default_timer()
i = 1
count = 0
finished = 0
while not finished:
# When j=10, j^i will have i+1 digits (e.g. if i=3, 10^3=1000).
for j in range(1, 10):
p = j ** i
if len(str(p)) == i:
count = count + 1
# When 9^i has less than i digits, all the numbers have been found.
if len(str(p)) < i:
finished = 1
i = i + 1
end = default_timer()
print('Project Euler, Problem 63')
print('Answer: {}'.format(count))
print('Elapsed time: {:.9f} seconds'.format(end - start))
if __name__ == '__main__':
main()

78
Python/p064.py Normal file
View File

@ -0,0 +1,78 @@
#!/usr/bin/python
# All square roots are periodic when written as continued fractions and can be written in the form:
#
# √N=a0+1/(a1+1/(a2+1/(a3+…
#
# For example, let us consider √23:
#
# √23=4+√234=4+1/(1/(√234))=4+1/(1+(√233)/7)
#
# If we continue we would get the following expansion:
#
# √23=4+1/(1+1/(3+1/(1+1/(8+…
#
# The process can be summarised as follows:
# a0=4,1/(√234)=(√23+4)/7=1+(√233)/7
# a1=1,7/(√233)=7(√23+3)/14=3+(√233)/2
# a2=3,2/(√233)=2(√23+3)/14=1+(√234)/7
# a3=1,7/(√234)=7(√23+4)/7=8+√234
# a4=8,1/(√234)=(√23+4)/7=1+(√233)/7
# a5=1,7/(√233)=7(√23+3)/14=3+(√233)/2
# a6=3,2/(√233)=2(√23+3)/14=1+(√234)/7
# a7=1,7/(√234)=7(√23+4)/7=8+√234
#
# It can be seen that the sequence is repeating. For conciseness, we use the notation √23=[4;(1,3,1,8)], to indicate that the block (1,3,1,8)
# repeats indefinitely.
#
# The first ten continued fraction representations of (irrational) square roots are:
#
# √2=[1;(2)], period=1
# √3=[1;(1,2)], period=2
# √5=[2;(4)], period=1
# √6=[2;(2,4)], period=2
# √7=[2;(1,1,1,4)], period=4
# √8=[2;(1,4)], period=2
# √10=[3;(6)], period=1
# √11=[3;(3,6)], period=2
# √12=[3;(2,6)], period=2
# √13=[3;(1,1,1,1,6)], period=5
#
# Exactly four continued fractions, for N≤13, have an odd period.
#
# How many continued fractions for N≤10000 have an odd period?
from math import floor, sqrt
from timeit import default_timer
from projecteuler import period_cf
def is_square(n):
p = sqrt(n)
m = int(p)
if p == m:
return 1
else:
return 0
def main():
start = default_timer()
count = 0
for i in range(2, 10000):
# Perfect squares are obviously not represented as continued fractions.
# For all other numbers, calculate their period and check if it's odd.
if not is_square(i) and period_cf(i) % 2 != 0:
count = count + 1
end = default_timer()
print('Project Euler, Problem 64')
print('Answer: {}'.format(count))
print('Elapsed time: {:.9f} seconds'.format(end - start))
if __name__ == '__main__':
main()

78
Python/p065.py Normal file
View File

@ -0,0 +1,78 @@
#!/usr/bin/python
# The square root of 2 can be written as an infinite continued fraction.
#
# √2=1+1/(2+1/(2+1/(2+1/(2+...
#
# The infinite continued fraction can be written, √2=[1;(2)], (2) indicates that 2 repeats ad infinitum. In a similar way, √23=[4;(1,3,1,8)].
# It turns out that the sequence of partial values of continued fractions for square roots provide the best rational approximations.
# Let us consider the convergents for √22.
#
# 1+1/2=3/2
# 1+1/(2+1/2)=7/5
# 1+1/(2+1/(2+1/2))=17/12
# 1+1/(2+1/(2+1/(2+1/2)))=41/29
#
# Hence the sequence of the first ten convergents for √2 are:
#
# 1,3/2,7/5,17/12,41/29,99/70,239/169,577/408,1393/985,3363/2378,...
#
# What is most surprising is that the important mathematical constant,
#
# e=[2;1,2,1,1,4,1,1,6,1,...,1,2k,1,...].
#
# The first ten terms in the sequence of convergents for e are:
#
# 2,3,8/3,11/4,19/7,87/32,106/39,193/71,1264/465,1457/536,...
#
# The sum of digits in the numerator of the 10th convergent is 1+4+5+7=17.
#
# Find the sum of digits in the numerator of the 100th convergent of the continued fraction for e.
from timeit import default_timer
def main():
start = default_timer()
ai = [1, 2, 1]
count = 4
n0 = 3
n1 = 8
n2 = 11
# For a continued fractions [a_0; a_1, a_2, ...], the numerator of the
# next convergent N_n=a_n*N_(n-1)+N_(n-2). The first three values for e are
# 3, 8 and 11, the next ones are easily calculated, considering that a_n
# follows a simple pattern:
# a_1=1, a_2=2, a_3=1
# a_4=1, a_5=4, a_6=1
# a_7=1, a_8=6, a_9=1
# and so on.
while count < 100:
ai[1] = ai[1] + 2
for i in range(3):
n0 = n1
n1 = n2
n2 = n1 * ai[i] + n0
count = count + 1
if count >= 100:
break
sum_ = 0
while n2 != 0:
sum_ = sum_ + n2 % 10
n2 = n2 // 10
end = default_timer()
print('Project Euler, Problem 65')
print('Answer: {}'.format(sum_))
print('Elapsed time: {:.9f} seconds'.format(end - start))
if __name__ == '__main__':
main()

View File

@ -179,3 +179,35 @@ def is_pentagonal(n):
return i.is_integer()
# Function implementing the iterative algorithm taken from Wikipedia
# to find the continued fraction for sqrt(S). The algorithm is as
# follows:
#
# m_0=0
# d_0=1
# a_0=floor(sqrt(n))
# m_(n+1)=d_n*a_n-m_n
# d_(n+1)=(S-m_(n+1)^2)/d_n
# a_(n+1)=floor((sqrt(S)+m_(n+1))/d_(n+1))=floor((a_0+m_(n+1))/d_(n+1))
# if a_i=2*a_0, the algorithm ends.
def period_cf(n):
mn = 0
dn = 1
count = 0
a0 = floor(sqrt(n))
an = a0
while True:
mn1 = dn * an - mn
dn1 = (n - mn1 * mn1) // dn
an1 = floor((a0+mn1)/dn1)
mn = mn1
dn = dn1
an = an1
count = count + 1
if an == 2 * a0:
break
return count