[seqfan] A055744

[David Corneth]

Hi all,
Brace yourself, it's a long one.It's about A055744, with the name "Numbers n such that n and phi(n) have same prime factors". I've submitted two sequences that may help to add elements to A055744. Those sequences are A256430 and A256431. A256430 is a subset of A256431. 
Some terms: Let HPF(n) be the highest prime factor of n. So HPF(12) = 3. Let P be a set of (distinct) primes, for example {2} or {3,5} and Let L(P) be the least element of A055744 divisible by that set of primes. A256431 has all and only those elements found by this function. Let LCE(m,n) be the least common element of m and n, m and n both elements of A256431 and one of them, say n, is an element of A256430. If m = L({3, 5, 7}) (Then, m = 7350 = 2 * 3 * 5^2 * 7^2), and n = L({13}) (Then, n = 1014 = 2 * 3 * 13^2). Then LCE(m,n) = 1242150. I described how to find these numbers; how LCE is defined here: https://oeis.org/history/view?seq=A256430&v=18  (it's still proposed while I write this). 
Sometimes, one gets duplicate values. For example, LCE(L({2,3,11}, L({5})) = LCE(18090, 50) = 18090. So 18090 is the least element of A055744 divisible by 2, 3 and 11, but also by 2, 3, 5 and 11. 
What I wonder about is whether the LCE is commutative. From the definition, I'd say yes. It would save half of the computation work; Since it would imply LCE(L({2,3,5}, L({11})) = LCE(L({11}), L({2,3,5})). Even more, is the function associative. I mean, is LCE(L({2,3,53}, L({11})) = LCE(L({2,3,11}, L({53})) without the need to compute these values to check if so. This would help since we could just compute new LCE's where HPF(n) > HPF(m). Again, I'd say yes though when I coded this up I noticed some elements where missed while assuming so.  
To actually compute some elements of A055744, I set an upperlimit. Say, I'd like to find all elements up to 10000. At some point, we have, to grasp back to my initial example, m = 7350 and n = 1014. We could just say that there is no need to compute their LCE since we can tell straight away it exceeds 10000. It's certainly >= 7350 and two factors 13 are added. Maybe some prime factor is divided away but that doesn't bring it under 13? Can we at some point before the highest prime <= 10000 say that we don't need to check anymore? Now, I stop after the highest prime beyond sqrt(floor(upperlimit / 6)). So in my latest edit in A055744 I checked n for all all primes from 2 to up to 129097. That's way to more than needed in general I think. As I couldn't find a proper stop, I just checked them all. 
Now the final bit, the multiples of elements of A256431. If we know 450 = 2^1 * 3^2 * 5^2 is in A256431 then 450 * k is in A055744 where k is a number of the form 2^x * 3^y * 5^z where x, y, z non-negative integers. To find these multiples upto say, 10000, I've constructed a numbersystem. Digits are exponents x, y and z. As we look for 450 * k < 10000, k integer, k < floor(10000 / 45) = 22. So we just need numbers upto 22 of the form 2^x * 3^y * 5^z. Start with the tuple (x, y, z) = (1, 0, 0) which represent k = 2. We increase z to 2 and now k = 4. That's okay. So we pick (3, 0, 0) and then (4, 0, 0). (5, 0, 0) is too high so we pick (0, 1, 0) then (1, 1, 0). At some point, the number would get 4 digits at which we stop searching and we have all desired numbers of k; (2, 4, 8, 16, 3, 6, 12, 9, 18, 5, 10, 20, 15) It finds them in that order. Starting with tuples (z, y, x) might be a more suitable choice; less exceeds of the upperlimit. I have some tight rules to choose a next tuple but I'll just give the code below.
\\------------------------------------------------------------------/*print all numbers <= n having primefactors just of g. in the Example I gave, n = 22 and g is for example, 2*3*5 = 30 (as long as it only has distinct primefactors 2, 3 and 5 for this example). */b(n,{g=n})={\\vector with distinct primefactors of gpfs = factor(g)[,1];t=1;c=0;
\\maximum exponents for primefactors. In the example, the vector would be [4, 2, 1]. x cannot exceed 4, y not 2 and z not \\1.mex = vector(#pfs,b,maxexp(n,pfs[b]));\\factors for numbersystem.fts=vector(#pfs,i,prod(j=i+1,#pfs,1+maxexp(n,pfs[j])));v=vector(#pfs);v[1]=1;

if(#pfs==1,for(i=1,log(u)\log(pfs[1]),print(pfs[1]^i)),while(v[#v]<=mex[#v],	   if(prod(i=1,#pfs,pfs[i]^v[i])<=n,t++;\\c++;       print(prod(i=1,#pfs,pfs[i]^v[i]));v=newLvfromleft(v),v=newSvfromleft(v);	   ));t)
}
\\next vector of exponents hopefully giving a lower numbernewSvfromleft(v)={i=1;while(v[i]==0,i++);if(i==#v,v[#v]++;v,v[i+1]++;v[i]=0);v}
newLvfromleft(v)={my(i=1);while(v[i]==mex[i],i++);for(j=1,i-1,v[j]=0);v[i]++;v}
\\highest exponent e s.t. b^e <= n.maxexp(n,b) = log(n)\log(b);\\----------------------------------------------------------------------------
I modified this from an algorithm I wrote for A244052 and the related A010846. In the latter, Charles Greathouse wrote an algorithm I couldn't modify to what I need here but this works good enough.
Enough issues as it seems and it was getting a bit over my head. Improving this method could, except the insight it gives, result in finding more elements to A256248. Any ideas?
Best,David