cp's OEIS Frontend

A054994 Numbers of the form q1^b1 * q2^b2 * q3^b3 * q4^b4 * q5^b5 * ... where q1=5, q2=13, q3=17, q4=29, q5=37, ... (A002144) and b1 >= b2 >= b3 >= b4 >= b5 >= ....

%I A054994 #53 Feb 16 2025 08:32:42
%S A054994 1,5,25,65,125,325,625,1105,1625,3125,4225,5525,8125,15625,21125,
%T A054994 27625,32045,40625,71825,78125,105625,138125,160225,203125,274625,
%U A054994 359125,390625,528125,690625,801125,1015625,1185665,1221025,1373125,1795625
%N A054994 Numbers of the form q1^b1 * q2^b2 * q3^b3 * q4^b4 * q5^b5 * ... where q1=5, q2=13, q3=17, q4=29, q5=37, ... (A002144) and b1 >= b2 >= b3 >= b4 >= b5 >= ....
%C A054994 This sequence is related to Pythagorean triples regarding the number of hypotenuses which are in a particular number of total Pythagorean triples and a particular number of primitive Pythagorean triples.
%C A054994 Least integer "mod 4 prime signature" values that are the hypotenuse of at least one primitive Pythagorean triple. - _Ray Chandler_, Aug 26 2004
%C A054994 See A097751 for definition of "mod 4 prime signature"; terms of A097752 with all prime factors of form 4*k+1.
%C A054994 Sequence A006339 (Least hypotenuse of n distinct Pythagorean triangles) is a subset of this sequence. - _Ruediger Jehn_, Jan 13 2022
%H A054994 Ray Chandler, <a href="/A054994/b054994.txt">Table of n, a(n) for n = 1..10000</a>
%H A054994 Eric Weisstein's World of Mathematics, <a href="https://mathworld.wolfram.com/PythagoreanTriple.html">Pythagorean Triple</a>
%F A054994 Sum_{n>=1} 1/a(n) = Product_{n>=1} 1/(1 - 1/A006278(n)) = 1.2707219403... - _Amiram Eldar_, Oct 20 2020
%e A054994 1=5^0, 5=5^1, 25=5^2, 65=5*13, 125=5^3, 325=5^2*13, 625=5^4, etc.
%t A054994 maxTerm = 10^15;(* this limit gives ~ 500 terms *) maxNumberOfExponents = 9;(* this limit has to be increased until the number of reaped terms no longer changes *) bmax = Ceiling[ Log[ maxTerm]/Log[q]]; q = Reap[For[k = 0 ; cnt = 0, cnt <= maxNumberOfExponents, k++, If[PrimeQ[4*k + 1], Sow[4*k + 1]; cnt++]]][[2, 1]]; Clear[b]; b[maxNumberOfExponents + 1] = 0; iter = Sequence @@ Table[{b[k], b[k + 1], bmax[[k]]}, {k, maxNumberOfExponents, 1, -1}]; Reap[ Do[an = Product[q[[k]]^b[k], {k, 1, maxNumberOfExponents}]; If[an <= maxTerm, Print[an]; Sow[an]], Evaluate[iter]]][[2, 1]] // Flatten // Union (* _Jean-François Alcover_, Jan 18 2013 *)
%o A054994 (PARI) list(lim)=
%o A054994 {
%o A054994   my(u=[1], v=List(), w=v, pr, t=1);
%o A054994   forprime(p=5,,
%o A054994     if(p%4>1, next);
%o A054994     t*=p;
%o A054994     if(t>lim, break);
%o A054994     listput(w,t)
%o A054994   );
%o A054994   for(i=1,#w,
%o A054994     pr=1;
%o A054994     for(e=1,logint(lim\=1,w[i]),
%o A054994       pr*=w[i];
%o A054994       for(j=1,#u,
%o A054994         t=pr*u[j];
%o A054994         if(t>lim, break);
%o A054994         listput(v,t)
%o A054994       )
%o A054994     );
%o A054994     if(w[i]^2<lim, u=Set(concat(Vec(v),u)); v=List());
%o A054994   );
%o A054994   Set(concat(Vec(v),u));
%o A054994 } \\ _Charles R Greathouse IV_, Dec 11 2016
%o A054994 (Python)
%o A054994 def generate_A054994():
%o A054994     """generate arbitrarily many elements of the sequence.
%o A054994     TO_DO is a list of pairs (radius, exponents) where
%o A054994     "exponents" is a weakly decreasing sequence, and
%o A054994     radius == prod(prime_4k_plus_1(i)**j for i,j in enumerate(exponents))
%o A054994     An example entry is (5525, (2, 1, 1)) because 5525 = 5**2 * 13 * 17.
%o A054994     """
%o A054994     TO_DO = {(1,())}
%o A054994     while True:
%o A054994         radius, exponents = min(TO_DO)
%o A054994         yield radius #, exponents
%o A054994         TO_DO.remove((radius, exponents))
%o A054994         TO_DO.update(successors(radius,exponents))
%o A054994 def successors(radius,exponents):
%o A054994     # try to increase each exponent by 1 if possible
%o A054994     for i,e in enumerate(exponents):
%o A054994         if i==0 or exponents[i-1]>e:
%o A054994             # can add 1 in position i without violating monotonicity
%o A054994             yield (radius*prime_4k_plus_1(i), exponents[:i]+(e+1,)+exponents[i+1:])
%o A054994     if exponents==() or exponents[-1]>0: # add new exponent 1 at the end:
%o A054994         yield (radius*prime_4k_plus_1(len(exponents)), exponents+(1,))
%o A054994 from sympy import isprime
%o A054994 primes_congruent_1_mod_4 = [5] # will be filled with 5,13,17,29,37,...
%o A054994 def prime_4k_plus_1(i): # the i-th prime that is congruent to 1 mod 4
%o A054994     while i>=len(primes_congruent_1_mod_4): # generate primes on demand
%o A054994         n = primes_congruent_1_mod_4[-1]+4
%o A054994         while not isprime(n): n += 4
%o A054994         primes_congruent_1_mod_4.append(n)
%o A054994     return primes_congruent_1_mod_4[i]
%o A054994 for n,radius in enumerate(generate_A054994()):
%o A054994     if n==34:
%o A054994         print(radius)
%o A054994         break # print the first 35 elements
%o A054994     print(radius, end=", ")
%o A054994 # _Günter Rote_, Sep 12 2023
%Y A054994 Subsequence of A097752.
%Y A054994 Cf. A002144, A006278, A006339, A071383, A097751, A097752, A097753, A097754, A097755, A097756.
%K A054994 easy,nonn
%O A054994 1,2
%A A054994 Bernard Altschuler (Altschuler_B(AT)bls.gov), May 30 2000
%E A054994 More terms from _Henry Bottomley_, Mar 14 2001