Euler Project 64

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	= 4 +	1
	1 23—4		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?


My Solution
This time, I use the Mathematica

Length[Select[Table[ContinuedFraction[Sqrt[i]], {i, 1, 10000}], Length[#] == 2 && OddQ [Length[#[[2]]]] &]]

Others' Solution using LISP

(defun list-as-for-sqrt (n)
(iter (for i first (isqrt n) then (truncate x))
        (for di first i then (- (- di (* i d))))
        (for d first (- n (* di di))
                then (/ (- n (* di di)) d))
        (for dl previous d initially -1)
        (for x = (/ (+ di (sqrt n)) d))
              (collecting i)
              (until (= dl 1))))
(defun 64-count-odd-periods (max)
    (iter (for n from 2 to max)
           (when (integer-root-p n) (next-iteration))
;; This checks evenness, because it counts the ;; initial digit, which isn't in the repeating ;; section.
           (when (evenp (length (list-as-for-sqrt n))) (sum 1))))

First Solution

Second Solution

(defun cont-frac-odd-p (n) 

    (do* ((b0 (isqrt n))
          (l b0 (- (* b p) l))
          (bflag t (not bflag))
          (p (- n (expt b0 2)) (/ (- n (expt l 2)) p)) 

          (b (if (zerop p) 

                 (return nil) 

                 (floor (+ b0 l) p)) 

            (floor (+ b0 l) p))) ((= 1 p) bflag))) 

(defun main (toplim) 

    (loop :for n :from 1 :to toplim :count (cont-frac-odd-p n))) 

(princ (main 10000)) (terpri)