%I #31 Jan 02 2024 09:35:57
%S 3,7,8,13,17,18,21,30,31,32,41,43,46,47,50,55,57,68,72,73,75,76,83,91,
%T 93,98,99,100,105,111,112,117,119,122,123,128,129,132,133,142,144,155,
%U 157,162,172,173,174,177,182,183,185,187,189,192,193,200,203,211
%N Numbers n such that arctan(1/n) = arctan(1/x) - arctan(1/y) for some integers 0<x<y<n.
%C 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.
%H Robert Israel, <a href="/A252496/b252496.txt">Table of n, a(n) for n = 1..10000</a>
%e 8 is in the sequence since arctan(1/8) = arctan(1/3) - arctan(1/5)
%p N:= 1000: # to get all terms <= N
%p A:= {}:
%p for x from 1 to N/2 do
%p ds:= select(d -> (d <= x and d >= (x^2+1)/(N-x)), numtheory:-divisors(x^2+1));
%p A:= A union map(d -> x + (x^2+1)/d, ds);
%p od:
%p A;
%p # if using Maple 11 or earlier, uncomment the next line
%p # sort(convert(A,list));
%p # _Robert Israel_, Dec 19 2014
%o (SageMath)
%o S = []
%o bound = 50
%o for b in (1 bound-1):
%o bb = b*b+1
%o for d in divisors(bb):
%o if (2*b < d) & (d-b < 2*bound):
%o c = d-b
%o a = (b*c-1)/(b+c)
%o S.append((c,b,a))
%o S.sort()
%o print(S)
%K nonn
%O 1,1
%A _Matthijs Coster_, Dec 17 2014