mata
function totient(n_) {
	n = n_
	if (n<4) {
		if (n<1) return(.)
		else if (n>1) return(n-1)
		else return(1)
	}
	else {
		r = 1
		if (mod(n,2)==0) {
			n = floor(n/2)
			while (mod(n,2)==0) {
				n = floor(n/2)
				r = r*2
			}
		}
		for (k=3; k*k<=n; k=k+2) {
			if (mod(n,k)==0) {
				r = r*(k-1)
				n = floor(n/k)
				while (mod(n,k)==0) {
					n = floor(n/k)
					r = r*k
				}
			}
		}
		if (n>1) r = r*(n-1)
		return(r)
	}
}

function map(f,a) {
	n = rows(a)
	p = cols(a)
	b = J(n,p,.)
	for (i=1; i<=n; i++) {
		for (j=1; j<=p; j++) {
			b[i,j] = (*f)(a[i,j])
		}
	}
	return(b)
}

function farey_length(n) {
	return(1+sum(map(&totient(),1::n)))
}

function farey(n) {
	m = 1+sum(map(&totient(),1::n))
	r = J(m,2,.)
	r[1,.] = 0,1
	a = 0
	b = 1
	c = 1
	d = n
	i = 1
	while (c<=n) {
		k = floor((n+b)/d)
		a = k*c-a
		b = k*d-b
		swap(a,c)
		swap(b,d)
		r[++i,.] = a,b
	}
	return(r)
}

for (n=1; n<=11; n++) {
	a = farey(n)
	m = rows(a)
	for (i=1; i<=m; i++) printf("%f/%f ",a[i,1],a[i,2])
	printf("\n")
}

map(&farey_length(),100*(1..10))
end