a propos divisors...

[Peter Pein]

Peter Pein schrieb:
> Dear seqfans,
> 
> I#ve been surprised not to find sequences of the following form in the OEIS:
> 
> a(n)=min(k in N: sigma(r,n)=sigma(r,k)) with sigma(r,n)=sum of the r-th
> power of the divisors of n:
...
> Are these of interest? And if so, up to which exponent r?
> 
> Peter
> 

Although the interst seems to be bounded by a fairly small n (0 until
now), I would like to ask you if you could please help me.

I decided to publish the sequences in two forms:

a) n such there is at least one x<n such that sigma(r,x)=sigma(r,n),
because there exists such a series for r=1 (A069822)

and

b) a(n)=min({k>0: sigma(r,k)=sigma(r,n)})


The issue with these is that they behave in a surprising way for r=3 (I
think, sequences for higher exponents are too hard to calculate without
 tricks)

for r = 0 the type-a-seq starts:
3, 5, 7, 8, 9, 10, 11, 13, 14, 15, 17, 18, 19, 20
for r = 1:
11, 15, 17, 23, 25, 26, 31, 35, 38, 39, 41, 46, 47, 51
for r = 2
7, 26, 35, 47, 77, 91, 119, 130, 133, 141, 157, 161, 175

but for r = 3 the air becomes very thin. If I did not make a mistake,
this seq starts(!):

194315, 295301, 590602, 1181204, 1476505, 1886920, 2067107, 2362408,
2526095, 2953010, 3248311, 3691985, 3838913, 4134214, 4469245, 4724816,
5020117, 5610719, 5635135, 5906020

If any of you with some Mathematica knowledge could please confirm, that
the following lines calculate the sequence described above?

Timing[Block[{spa, nmax = 6*10^6, expo = 3},
   Reap[For[n = 1, n <= nmax, n++,
      (If[Head[#1] === spa, #1 = n, Sow[{n, #1}]] & )[
       spa[DivisorSigma[expo, n]]]]][[2,1]]]]

the name spa is an artefact; I tried this with SparseArrays, but the
allowed range of indices has not been sufficient. I use it as an
initially undefined function (Block[{spa..}]). For each n to test I look
wether spa[sigma(r,n)] has been defined. If not, the Head is still spa
and I set spa[sigma(r,n)] to n; else the remembered value together with
n will go to the result via the Sow-Reap mechanism. This way I get
type-a and a hint for type-b sequences in one run of the proggie.

And if you want to make me really happy (I have had birthday on July
24th ;-) (http://www.stevesbeatles.com/songs/when_im_sixty_four.asp this
age will be reached in 20 years, but the symptoms... )):

 If you've got Mathematica and more RAM (4GB or so) than I do (1.5 GB),
could you please run this code  with, say nmax=10 or 20 million? On my
machine it swapped heavily with nmax=6 million and I had to kill
MathKernel as I tried nmax=10^7. The lines above took ~181 seconds to
evaluate (nmax=10^7 has been stopped by me after 15 minutes). I do not
expect any runtimes of more than 7 minutes. Would this be possible, please?

 Alternatively any hints how to calculate these sequences more efficient
would be highly appreciated (AFAIK there exists no kind of "inverse
function" to sigma(r,n) w.r.t. n which could be calculated without this
brute-force method).

Thank you for your attention and in advance for CPU-time,

Peter



* Peter Pein <petsie at dordos.net> [Jul 26. 2007 08:51]:
> The link to a Mma-script to easy put sequences in the OEIS-format on the
> page http://www.research.att.com/~njas/sequences/Submit.html which
> should link to a Mma-file http://www.seqfan.net/EISFormat.m leads me to
> a strange (hijacking?) site full of advertising which states that
> "seqfan.net expired on 07/12/2007 and is pending renewal or deletion."

This is not a domain grabber but "network solutions"

> 
> If I remember correctly, there's a package at mathworld.wolfram.com with
> some utilities needed for Eric's project which include some EISFormat
> functionality.

http://www.aboutus.org/Seqfan.net
ogerard at ext.jussieu.fr (CC)

The domain should simply be renewed.
Else at some point in time a domain grabber will get it.