Re^2: Error-finding program

[Richard Mathar]

> From seqfan-owner at ext.jussieu.fr  Wed Oct 25 21:07:37 2006
> Return-Path: <seqfan-owner at ext.jussieu.fr>
> To: seqfan at ext.jussieu.fr
> Subject: Re: Error-finding program
> Date: Wed, 25 Oct 2006 15:05:48 -0400
> ...
>     ... are suspicious:
> A007769 It looks like terms were added here that belong with an
> entirely different sequence.  I don't know if this is a computational
> error in the original sequence, the addition, or possibly the terms
> were added to the wrong sequence.


The first error seems to be in the final digit of A007769(14). A007769(15)
looks alright, then A007769(16...) do not.
Just plugging in the formula given in the OEIS, we generate

n A007769(n)
0 1
1 1
2 2
3 5
4 18
5 105
6 902
7 9749
8 127072
9 1915951
10 32743182
11 624999093
12 13176573910
13 304072048265
14 7623505722158
15 206342800616597
16 5996837126024824
17 186254702826289089
18 6156752656678674792
19 215810382466145354405
20 7995774669504366055054

The PARI program is

doublefactorial(n)={
	local(resul) ;
	resul=1 ;
	forstep(i=n,2,-2,
		resul *= i ;
	) ;
	return(resul) ;
}
alpha(n,q)={
	if(q %2,
		return( q^(p/2)*doublefactorial(p-1)) ,
		return( sum(k=0,p/2,binomial(p,2*k)*q^k*doublefactorial(2*k-1)) ) ;
	) ;
}

A007769(n)={
	local(resul,q) ;
	if(n==0,
		return(1),
		resul=0 ;
		fordiv(2*n,p,
			q=2*n/p ;
			resul += alpha(p,q)*eulerphi(q) ;
		);
		return(resul/(2*n)) ;
	) ;
}

{
	for(n=0,20,
		print(n," ",A007769(n)) ;
	) ;
}