[seqfan] Re: Sequence related to A180919

[Joerg Arndt]

* Robert Israel <israel at math.ubc.ca> [Jan 18. 2011 19:00]:
> I'm not "heavy duty" by any means, but if n is even (say n = 2*m),
> x^2 = k^2 + 2*m*k + 1 becomes
> (k+m+x)(k+m-x) = m^2 - 1
> and the number of nonnegative integer solutions is the number of
> ways to write m^2 - 1 as the (unordered) product of two nonnegative
> integers that are either both odd (if m is even) or both even (if m
> is odd).  For example, with n=10, m=5, 24 = 2*12 = 4*6 corresponding
> to the solutions (k=2,x=5) and (k=0,x=1), while with n=16, m=8, 63 =
> 1*63 = 3*21 = 7*9
> corresponding to the solutions (k=24,x=31), (k=4,x=9) and (k=0,x=1).
> 
> If n is odd (say n=2*m+1), x^2=k^2+(2*m+1)*k+1 becomes
> (2*k+2*m+1-2*x)*(2*k+2*m+1+2*x) = (2*m+3)*(2*m-1)
> and solutions correspond to ways to write (2*m+3)*(2*m-1) as the
> product of two nonnegative odd integers.  For example, with
> n=11, m=5, 117 = 1*117 = 3*39 = 9*13 corresponding to (k=24,x=29),
> (k=5,x=9) and (k=0,x=1).
> 
> Robert Israel                                israel at math.ubc.ca
> Department of Mathematics        http://www.math.ubc.ca/~israel
> University of British Columbia            Vancouver, BC, Canada
> 
> On Tue, 18 Jan 2011, Joerg Arndt wrote:
> 

/* Failed attempt to program this: */

ev(n)=
{
    local(ct, m, pr);
    ct = 0; /* count squares */
    m = n / 2;
    pr = m^2 - 1;
    if ( bitand(m,1)==0, /* m even */
        fordiv (pr, d,  if (  bitand(d,1)==1 && bitand(pr\d,1)==1, ct += 1) );
        , /* else  (n even) */
        fordiv (pr, d,  if (  bitand(d,1)==0 && bitand(pr\d,1)==0, ct += 1) );
    );
    return( ct/2 );  /* unordered products */
} /* ----- */

od(n)=
{
    local(ct, m, pr);
    ct = 0; /* count squares */
    m = n \ 2;  /* div without remainder */
    pr = (2*m+3)*(2*m-1);
    fordiv (pr, d,  if (  bitand(d,1)==1 && bitand(pr\d,1)==1, ct += 1) );
    return( ct/2 );  /* unordered products */
} /* ----- */

N=100;  /* max N */
{ for (n=3, N,
    if ( n%2==0, ct=ev(n), ct=od(n) );
    print1(ct,", ");
); }

/* Gives: */
1, 1, 2, 1, 3, 2, 2, 2, 3, 2, 4, 3, 2, 3, 4, 3, 4, 3, 2, 4, 6, 2, 4,
4, 3, 4, 4, 4, 4, 4, 4, 6, 4, 2, 8, 6, 2, 4, 4, 4, 6, 4, 2, 6, 9, 3,
4, 6, 3, 6, 8, 4, 4, 4, 4, 8, 4, 2, 6, 10, 4, 4, 6, 5, 8, 8, 2, 6, 4,
4, 6, 6, 4, 4, 6, 6, 10, 4, 2, 12, 10, 2, 4, 8, 4, 6, 8, 4, 4, 6, 8,
8, 4, 3, 12, 12, 2, 6,