|
%I #28 Mar 13 2023 15:40:43
%S 1,9,36,100,225,441,784,900,1225,1296,2025,3025,4356,6084,7056,8281,
%T 11025,14400,18496,23409,29241,32400,36100,41616,44100,53361,64009,
%U 76176,88209,90000,105625,108900,123201,142884,164836,189225,216225,246016,278784,298116
%N Perfect squares that are a product of two triangular numbers.
%C Includes (except for 0) A000537 and 3/2*x*(x+1) for x in A132596. - _Robert Israel_, Jan 16 2015
%H Alois P. Heinz, <a href="/A169835/b169835.txt">Table of n, a(n) for n = 1..1000</a>
%H Erich Friedman, <a href="https://erich-friedman.github.io/numbers.html">What's Special About This Number?</a> (See entry 7056.)
%p N:= 10^6: # to get all terms <= N
%p A:= select(issqr,{seq(seq(a*(a+1)*b*(b+1)/4,
%p b = a .. floor(sqrt(1/4+4*N/a/(a+1))-1/2)),a=1..floor(sqrt(4*N)))});
%p # if using Maple 11 or earlier, uncomment the next line
%p # sort(convert(A, list)); # _Robert Israel_, Jan 16 2015
%t M = 10^6; (* to get all terms <= M *)
%t A = Union[Select[Flatten[Table[Table[(1/4) a (a+1) b (b+1), {b, a, Floor[ Sqrt[1/4 + 4M/(a (a+1))] - 1/2]}], {a, 1, Floor[Sqrt[4M]]}]], IntegerQ[ Sqrt[#]]&]] (* _Jean-François Alcover_, Mar 09 2019, after _Robert Israel_ *)
%o (PARI) istriangular(n)=issquare(8*n+1) \\ now one can use ispolygonal(n, 3)
%o isok(n) = {if (issquare(n), fordiv(n, d, if (d > sqrtint(n), break); if (istriangular(d) && istriangular(n/d), return (1)););); return (0);} \\ _Michel Marcus_, Jul 24 2013
%o (Haskell)
%o a169835 n = a169835_list !! (n-1)
%o a169835_list = f [] (tail a000217_list) (tail a000290_list) where
%o f ts us'@(u:us) vs'@(v:vs)
%o | u <= v = f (u : ts) us vs'
%o | any p $ map (divMod v) ts = v : f ts us' vs
%o | otherwise = f ts us' vs
%o where p (q, r) = r == 0 && a010054 q == 1
%o -- _Reinhard Zumkeller_, Mar 03 2015
%o (Python)
%o from itertools import count, islice, takewhile
%o from sympy import divisors
%o from sympy.ntheory.primetest import is_square
%o def A169835_gen(): # generator of terms
%o return filter(lambda k:any(map(lambda d: is_square((d<<3)+1) and is_square((k//d<<3)+1), takewhile(lambda d:d**2<=k,divisors(k)))),(m**2 for m in count(0)))
%o A169835_list = list(islice(A169835_gen(),20)) # _Chai Wah Wu_, Mar 13 2023
%Y Superset of A000537. Cf. A000217, A132596, A169836.
%Y Cf. A000290, A010054.
%K nonn
%O 1,2
%A _R. J. Mathar_, May 30 2010
%E Corrected (missing terms inserted) by _R. J. Mathar_, Jun 04 2010
| 