OFFSET
1,1
COMMENTS
arctan(1/a(n)) = arctan(1/x) - arctan(1/y) for some integers x and y where 0 < x < y < a(n). We use the formula tan(a+b) = (tan a + tan b)/(1 - tan a.tan b) which implies that 1/a(n) = (1/x - 1/y)/(1+1/(xy)) or a(n) = (xy+1)/(y-x) = x + (x^2+1)/(y-x). So we look for divisors of x^2+1.
LINKS
Robert Israel, Table of n, a(n) for n = 1..10000
EXAMPLE
8 is in the sequence since arctan(1/8) = arctan(1/3) - arctan(1/5)
MAPLE
N:= 1000: # to get all terms <= N
A:= {}:
for x from 1 to N/2 do
ds:= select(d -> (d <= x and d >= (x^2+1)/(N-x)), numtheory:-divisors(x^2+1));
A:= A union map(d -> x + (x^2+1)/d, ds);
od:
A;
# if using Maple 11 or earlier, uncomment the next line
# sort(convert(A, list));
# Robert Israel, Dec 19 2014
PROG
(SageMath)
S = []
bound = 50
for b in (1 bound-1):
bb = b*b+1
for d in divisors(bb):
if (2*b < d) & (d-b < 2*bound):
c = d-b
a = (b*c-1)/(b+c)
S.append((c, b, a))
S.sort()
print(S)
CROSSREFS
KEYWORD
nonn
AUTHOR
Matthijs Coster, Dec 17 2014
STATUS
approved