diff --git a/C/p061.c b/C/p061.c
new file mode 100644
index 0000000..0185e3d
--- /dev/null
+++ b/C/p061.c
@@ -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(3n−1)/2	 	1, 5, 12, 22, 35, ...
+ * Hexagonal	 	P6,n=n(2n−1)	 	1, 6, 15, 28, 45, ...
+ * Heptagonal	 	P7,n=n(5n−3)/2	 	1, 7, 18, 34, 55, ...
+ * Octagonal	 	P8,n=n(3n−2)	 	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;
+}
diff --git a/C/p062.c b/C/p062.c
new file mode 100644
index 0000000..e368524
--- /dev/null
+++ b/C/p062.c
@@ -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;
+}
diff --git a/C/p063.c b/C/p063.c
new file mode 100644
index 0000000..9e1840f
--- /dev/null
+++ b/C/p063.c
@@ -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;
+}
diff --git a/C/p064.c b/C/p064.c
new file mode 100644
index 0000000..337a3fd
--- /dev/null
+++ b/C/p064.c
@@ -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+√23−4=4+1/(1/(√23−4))=4+1/(1+(√23−3)/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/(√23−4)=(√23+4)/7=1+(√23−3)/7
+ * a1=1,7/(√23−3)=7(√23+3)/14=3+(√23−3)/2
+ * a2=3,2/(√23−3)=2(√23+3)/14=1+(√23−4)/7
+ * a3=1,7/(√23−4)=7(√23+4)/7=8+√23−4
+ * a4=8,1/(√23−4)=(√23+4)/7=1+(√23−3)/7
+ * a5=1,7/(√23−3)=7(√23+3)/14=3+(√23−3)/2
+ * a6=3,2/(√23−3)=2(√23+3)/14=1+(√23−4)/7
+ * a7=1,7/(√23−4)=7(√23+4)/7=8+√23−4
+ *
+ * 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?*/
+
+#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;
+   }
+}
diff --git a/C/p065.c b/C/p065.c
new file mode 100644
index 0000000..1102ef6
--- /dev/null
+++ b/C/p065.c
@@ -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;
+}
diff --git a/C/projecteuler.c b/C/projecteuler.c
index 2467c96..8801c54 100644
--- a/C/projecteuler.c
+++ b/C/projecteuler.c
@@ -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;
+}
diff --git a/C/projecteuler.h b/C/projecteuler.h
index 55260cb..e87cf0c 100644
--- a/C/projecteuler.h
+++ b/C/projecteuler.h
@@ -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
diff --git a/Python/p061.py b/Python/p061.py
new file mode 100644
index 0000000..4525763
--- /dev/null
+++ b/Python/p061.py
@@ -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(3n−1)/2	    1, 5, 12, 22, 35, ...
+# Hexagonal	 	P6,n=n(2n−1)	    1, 6, 15, 28, 45, ...
+# Heptagonal	 	P7,n=n(5n−3)/2	    1, 7, 18, 34, 55, ...
+# Octagonal	 	P8,n=n(3n−2)	    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()
diff --git a/Python/p062.py b/Python/p062.py
new file mode 100644
index 0000000..9118fc9
--- /dev/null
+++ b/Python/p062.py
@@ -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()
diff --git a/Python/p063.py b/Python/p063.py
new file mode 100644
index 0000000..e376a9b
--- /dev/null
+++ b/Python/p063.py
@@ -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()
diff --git a/Python/p064.py b/Python/p064.py
new file mode 100644
index 0000000..e02fbda
--- /dev/null
+++ b/Python/p064.py
@@ -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+√23−4=4+1/(1/(√23−4))=4+1/(1+(√23−3)/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/(√23−4)=(√23+4)/7=1+(√23−3)/7
+# a1=1,7/(√23−3)=7(√23+3)/14=3+(√23−3)/2
+# a2=3,2/(√23−3)=2(√23+3)/14=1+(√23−4)/7
+# a3=1,7/(√23−4)=7(√23+4)/7=8+√23−4
+# a4=8,1/(√23−4)=(√23+4)/7=1+(√23−3)/7
+# a5=1,7/(√23−3)=7(√23+3)/14=3+(√23−3)/2
+# a6=3,2/(√23−3)=2(√23+3)/14=1+(√23−4)/7
+# a7=1,7/(√23−4)=7(√23+4)/7=8+√23−4
+#
+# 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()
diff --git a/Python/p065.py b/Python/p065.py
new file mode 100644
index 0000000..650b8e2
--- /dev/null
+++ b/Python/p065.py
@@ -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()
diff --git a/Python/projecteuler.py b/Python/projecteuler.py
index 67c83ee..56e8449 100644
--- a/Python/projecteuler.py
+++ b/Python/projecteuler.py
@@ -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