RosettaCodeData/Task/Farey-sequence/Stata/farey-sequence.stata
2023-07-01 13:44:08 -04:00

77 lines
1 KiB
Text

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