A169835 - OEIS

%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