Re: sigma

[Karim Belabas <Karim.Belabas@math.u-psud.fr>]

* Bill Allombert [2004-05-30 12:58]:
> On Sun, May 30, 2004 at 12:18:57PM +0200, Karim Belabas wrote:
> > * Jon Perry [2004-03-22 17:10]:
> > > ? sigma(4,1.5)
> > >   ***   this should be an integer: sigma(4,1.5)
> > >                                            ^----
> > You may use something like
> > 
> > SIGMA(n, k) = 
> > { local(R, S, f, pk); 
> >   R = 1; f = factor(n);
> >   for (i = 1, matsize(f)[1],
> >     pk = f[i, 1]^k; S = 1 + pk; 
> >     for (j = 2, f[i, 2], S = 1 + pk*S);
> >     R *= S
> >   ); R 
> > }
> > 
> > it's about 2 to 10 times slower than the built-in sigma [ for smooth n ! ].
> > The latter makes heavy use of immediate small integers which is tough to
> > emulate in GP.
> 
> Much simpler, just use sumdiv
> SIGMA(n, k)=sumdiv(n,d,d^k)

I considered this at first, but it's much slower and less robust:

(13:05) gp > n = round( prodeuler(p = 2, 70, p) )
%1 = 7858321551080267055879090

(13:05) gp > for (i=1,1000, SIGMA(n,2))  \\ my routine
time = 910 ms.

(13:05) gp > for (i=1,1000, sigma(n,2))
time = 40 ms.

(13:05) gp > sumdiv(n, d, d^2);   \\ default stack
  ***   the PARI stack overflows !  

(13:07) gp > sumdiv(n, d, d^2);   \\ 80MB stack
time = 10,960 ms.

    Karim.
-- 
Karim Belabas                     Tel: (+33) (0)1 69 15 57 48
Dep. de Mathematiques, Bat. 425   Fax: (+33) (0)1 69 15 60 19
Universite Paris-Sud              http://www.math.u-psud.fr/~belabas/
F-91405 Orsay (France)            http://pari.math.u-bordeaux.fr/  [PARI/GP]