User:Peter Luschny/OddsAndEnds

OddsAndEnds 

On this page I will gather some small observations on 
sequences and on OEIS as they appear in the course of time.
You are invited to comment, see the `discussion page´.

Synopsis

• What is a prime number?
  I was not satisfied how the term 'prime number' is defined in A000040. To say that a prime is a natural number with exactly two divisors is a correct but dumb definition. A prime number is a number which is conserved as a factor in all factorizations of all its multiples.
• How to compute the number of partitions of 1,000,000?
  Inspired by Fredrik Johansson's paper on the "Efficient implementation of the Hardy–Ramanujan–Rademacher formula" I gave a short implementation of A000041 and found an extension of Euler's pentagonal number theorem yet to be analyzed.
• How to compute sqrt(2+sqrt(2))+1)^n+(-sqrt(2)-1)^n+ ... ?
  Linear recurrences of higher order are not always easy to identify. In this example I didn't, but Colin Barker found the solution; see A215503.
• Gauss factorials.
  The term Gauss factorial was introduced by Cosgrave and Dilcher, a terminology suggested by a theorem of Gauss generalizing Wilsons's theorem; see A216919.
• The Scambler statistic on Dyck words.
  A new statistic on Dyck paths introduced by David Scambler has a bunch of connections with some well known sequences but unexpectedly also with two sequences I had submitted; see A217540.
• How to compute the 1000-th derivative of Lambert W(exp(x))?
  Admittedly this is a somewhat esoteric question. But it entails some nice relations to the second-order Eulerian numbers and polynomials A201637.
• What does 'frac' mean?
  When is the attribute 'frac' given to an integer sequence? I could not infer the meaning from its use on OEIS and tried to give a definition for this keyword.

What is a prime number?

Oct 09, 2012

We all know what prime numbers are, right? Well, I am not so sure. At least I got confused last week at some point. So I looked up A000040. Two definitions were given there as the first comments.

• (D1) A number n is prime if it is greater than 1 and has no positive divisors except 1 and n.
• (D2) A number n is prime if and only if it has exactly two positive divisors.

I will assume that in (D1) `if´ is to be read as `if and only if´ (symbol ⇔).

• (D1*) A number n is prime ⇔ n is greater than 1 and has no positive divisors except 1 and n.

First I had to know what a `number´ is. Numbers can be `natural numbers´, `integers´, `rational numbers´ or other types of numbers. From the name `OEIS´ I inferred that `number´ has on OEIS the default meaning `integer´ (ganze Zahlen, symbol: Z), and is defined as on Wikipedia. The natural numbers (natürliche Zahlen, symbol: N) are a subset of Z and have a different algebraic structure. The integers Z form a ring and the natural numbers N a semiring.

Often N is identified with the positive integers and this is a possible source of confusion. The expression `n is prime´ is meaningless. The universe of discourse in which the term `prime´ is applied has to be made explicit. The meaningful way is to say: `n is prime in R´, where R is a ring (or a semiring (or even a commutative monoid)).

So I rephrased the definitions:

• (D1z) n is prime in Z ⇔ n is greater than 1 and has no positive divisors except 1 and n.
• (D2z) n is prime in Z ⇔ n has exactly two positive divisors.

Then I made some tests.

If n = 2 then n is in Z, n>1, and only 1 and 2 divide 2. 
   (D1z passed).

If n = 2 then n is in Z, and {1, 2} are the only positive 
   divisors of 2. (D2z passed).

If n = -2 then n is in Z, n>1 is false and therefore -2 is not
   a prime. But this is false! -2 is prime in Z! (D1z failed).

If n = -2 then n is in Z, and {1, 2} are the only positive 
   divisors of 2. So -2 is prime (D2z passed). But this is 
   bad news because -2 is not in A000040.

Thus under the interpretation of number = element of Z both definitions (D1) and (D2) are wrong (or the sequence A000040 is incomplete)! I am sure that this fact has plunged more than one reader of A000040 into confusion.

So perhaps number = element of N is meant. So let's try this.

• (D1n) n is prime in N ⇔ n is greater than 1 and has no positive divisors except 1 and n.
• (D2n) n is prime in N ⇔ n has exactly two positive divisors.

Here my immediate critique is the use of the term `positive´. If n is in N then all ingredients are `positive´. The notion `negative element´ is not even defined in N! Thus the formulation is confusing right away (or better: suggests that using R=Z was indeed the intention of the author). Now, since a definition has, by definition of a definition, to focus on the basic essentials, I removed the inessential and thereby confusing attribute `positive´.

• (D1na) n is prime in N ⇔ n is greater than 1 and has no divisors except 1 and n.
• (D2na) n is prime in N ⇔ n has exactly two divisors.

[Edit start] As M. F. Hasler in his comment below rightly says, the use of an order (here: n>1) in the definition is not adequate; an order is not needed. Therefore (D1na) should be replaced by (D1nb):

• (D1nb) n is prime in N ⇔ n is not equal to 1 and has no divisors except 1 and n. [Edit end]

Since the two variants (D1) and (D2) are not so different I ultimately I agreed on the following two formulations:

• (Dnp) n is prime in N ⇔ n has exactly two divisors, which are {1, n}.
• (Dzp) n is prime in Z ⇔ n has exactly four divisors, which are {−n, −1, 1, n}.

However, these formulations are poor definitions (although they are correct). One is immediately willing to ask: "Why does this notion of a prime as a number with exactly two divisors only hold for natural numbers and not for integers?" "What is the common idea behind the different definitions in the cases N and Z?" "What do I have to expect next, an algebraic structure where a prime has exactly six divisors?"

These definitions are unconvincing, misleading and not capable of a generalization. They hide more than convey the meaning by rising a rather incidental phenomenon to a definition.

SSo is there a way to get rid of all these shortcomings? Yes. We have to switch to a slightly more abstract definition which applies equally for the integers and the natural numbers (and additionally also covers the general case of a commutative ring). It does not make any assertion on the number of divisors of n and reads:

• A number n is prime in R if and only if n is different from zero and different from a unit and each multiple of n decomposes into factors such that n divides at least one of the factors. /li>

NNo more confusion is possible. And the meaning of the term `prime number´ becomes much more transparent: a prime number is conserved in all its multiples. For example 12 is a multiple of 3 and 4, and in 12 you can get rid of the 4 in by factoring 12=2*6. 4 does not divide 2 and does not divide 6, thus 4 is not prime. But there is no way to get rid of the 3 by factoring 12 .

Comments:/p> <p>M. F. Hasler: "IMO, "greater than 1" is even worse ....

How to compute the number of partitions of 1,000,000?

Oct 12, 2012

I looked up A000041 to see what proposals are made in the code sections. For example the Mathematica section shows the command `PartitionsP[n]´ and indeed I got almost instantly from Wolfram Alpha the right answer. [1]

This is quite amazing. It isn't that trivial to compute p(n) for large n. But perhaps Mathematica just looked it up in a database? At any rate, I wanted to know the formula behind the magic and get at least an vague idea how it could have been implemented. So I looked for other options.

Joerg Arndt suggests in A000041 - despite the fact that Pari has for quite some time the command `numbpart´ which computes p(n) almost as fast as Mathematica - the formula

a(n) = if(n<0, 0, 
       polcoeff(exp(sum(k=1,n,x^k/(1-x^k)/k,x*O(x^n))),n)).

BBut I had no chance to verify this as Pari quits with the error message: "Overflow in valp()." The same happens when I try the also advertised formula:

a(n) = if(n<0, 0, polcoeff(1 / eta(x + x * O(x^n)), n)).

Ralf Stephan shared "the Hardy-Ramanujan-Rademacher exact formula" in Pari. Trying this, Pari quits again with an error message, this time saying "precision too low in truncr."

Now this is a very serious problem. In the lucky case you get an aborted computation; much more dangerous is to give a precision which is enough for `truncr´ but not enough for computing the right value.

This for example happens if you set default(realprecision, 1100) in Pari. The value returned is wrong. But how do I know what precision I have to use? Without an answer to this question the proposed code is useless and potential dangerous.

After several tries I managed ari to find the correct value with default(realprecision, 1200) in 20.8 seconds with Stephan's HRR implementation.

AAn editor-in-chief once wrote that a code snippet in OEIS should be able to compute the values of a sequence without any further addition by the user. Measured against this criterion this implementation clearly fails.

Omar ol's GWbasic program reminded me of Euler's pentagonal

number theorem. With Sage this can be implemented as:

@CachedFunction
def epnt(n):
    if n == 0 : return 1
    if n < 0 : return 0
    return add((-1)^(k+1)*(epnt(n-k*(3*k-1)/2)
               +epnt(n-k*(3*k+1)/2)) for k in (1..n))

This works nicely for n up to about 200. Then the maximum recursion depth of Python is exceeded (although it is a cached function!). And I do not want to tinker with sys.setrecursionlimit().

But Omar Pol's code is more sophisticated than that. He uses the generalized pentagonal numbers in a way, which is outlined by R. K. Guy in A001318. (Guy refers there to a relation of Conway which remains unexplained.) With a little more tuning I arrived at:

@CachedFunction
def A000041(n):
    if n == 0: return 1
    S = 0; J = n-1; k = 2
    while 0 <= J:
        T = A000041(J)
        S = S+T if is_odd(k//2) else S-T
        J -= k if is_odd(k) else k//2
        k += 1
     return S   

This function is well suited for computing the list of the first k partition numbers. On my laptop [A000041(n) for n in (0..1000)] executes in 0.6 seconds. However, it would exceed my patience to wait for (0..1000000).

Now it's time to make a serious effort. I went back to the HRR (Hardy-amanujan-Rademacher) formula, this time with Sage.

def A(n, k):
    if k <= 1: return k 
    if k == 2: return (-1)^n
    (s, r, m) = (0, 2, n%k)
    for l in range(2*k):
        if m == 0: s += (-1)^l*cos(pi*(6*l+1)/(6*k)) 
        m += r 
        if k <= m: m -= k
        r += 3 
        if k <= r: r -= k
    return (k/3)^(1/2)*s

def p(n):
    if n == 0: return 1
    U = lambda x: cosh(x)-sinh(x)/x
    C = lambda n: (pi/6)*sqrt(24*n-1)
    R = add(sqrt(3/k)*4/(24*n-1)*A(n,k)*U(C(n)/k) 
            for k in (1..isqrt(n))) 
    return R.round()

With the default precision of Sage all goes well up to n=235, it will fail thereafter because of the precision problem. So now I had to solve this nasty problem. After some trial and error and reading Fredrik Johansson's excellent account on the implementation of the Hardy-Ramanujan-Rademacher formula and from which I took the above algorithm for evaluating A(n, k) I came up with a simple formula:

    precision(n) = max(53, 3.7*sqrt(n)))

Precision is the number of bits used to represent the mantissa of a floating-point number, 53 is the default of Sage. I checked the validity numerically up to n=10000.

Thus all I had to do was to add to the above program the line R = RealField(max(53, 3.7*sqrt(n))) and make sure, that the functions are evaluated in this field.

And voilà, finally I could compute p(1000000) on my computer in 6.7 seconds, three times faster than the Pari implementation of HRR in A000041.

---

Fredrik Johansson, Efficient implementation of the Hardy–Ramanujan–Rademacher formula, LMS Journal of Computation and Mathematics (15), Oct. 2012, pp 341-359.

Comments:

Peter Luschny: Extending the implementation of Euler's pentagonal number theorem.

Susanne Wienand: Concerning the extension of Euler's pentagonal number theorem.

How to compute sqrt(2+sqrt(2))+ 1)^n + (-sqrt(2)-1)^n + ...?

Oct 19, 2012

$\left({\sqrt {2+{\sqrt {2}}}}+1\right)^{n}+\left(-{\sqrt {2}}-1\right)^{n}+\left({\sqrt {2-{\sqrt {2}}}}+1\right)^{n}+(-1)^{n}+\left(-{\sqrt {2-{\sqrt {2}}}}+1\right)^{n}+\left({\sqrt {2}}-1\right)^{n}+\left(-{\sqrt {2+{\sqrt {2}}}}+1\right)^{n}.$

Bad solution:

def a(n):
    return (sqrt(2+sqrt(2))+1)^n + (-sqrt(2)-1)^n
        + (sqrt(2-sqrt(2))+1)^n
        + (-1)^n + (-sqrt(2-sqrt(2))+1)^n
        + (sqrt(2)-1)^n + (-sqrt(2+sqrt(2))+1)^n

[a(i).round() for i in (0..28)]

CCan you find a better way to compute the sequence? (A solution can be found on the discussion page.)

Gauss factorials

Oct 21, 2012

For positive integers N and n let N_n! denote the product of all integers up to N that are relatively prime to n. The term Gauss factorial was introduced by Cosgrave and Dilcher, a terminology suggested by the following theorem of Gauss, a generalization of Wilsons's theorem.

For n in $\{2,\ 4,\ p^{a},\ 2p^{a}\}$, p odd prime, a positive integer $(n-1)_{n}!\ \equiv -1{\pmod {n}};$ 1 (mod n) otherwise.

To understand this theorem let us first consider this function:

def A(n):
    r = 1
    for k in (1..n):
        if gcd(n, k) == 1: r = r*k
    return mod(r, n) 

[A(n) for n in (1..20)]0, 1, 2, 3, 4, 5, 6, 1, 8, 9, 10, 1, 12, 13, 1, 1, 16, ...

ThThis is not an efficient implementation. First it builds a potentially huge number (which can be (n-1)! in the worst case) and then reduces via the modulus operation to something which is smaller than n. Much better is to reduce r in each step; this will keep r always < n.

def A160377(n):    r = 1
    for k in (1..n):
        if gcd(n, k) == 1: r = mod(r*k, n)
    return r

0, 1, 2, 3, 4, 5, 6, 1, 8, 9, 10, 1, 12, 13, 1, 1, ...

Now we can give A160377 a name which is easy to understand:

"The product of the residues in [1, n] relatively prime to n, taken modulo n."

(The actual name in the database is unattractive OEIS lingo (Phi-torial of n (A001783) modulo n) ) and the first Maple implementation is unfortunately the inefficient one.)

But wait, must not be something wrong here as the Gauss-Wilson theorem speaks only about -1 (mod n) and 1 (mod n)? The reason is this:p> <p>The value of the modulo function depends on which representation for integers modulo n is used. Here `mod´ means that the positive representation is used, which means that the values are reduced to integers in the range [0, abs(n)-1]. This `mod´ is also denoted by `modp´ or `pmod´. Alternatively the values can also be reduced to integers in the range [-(abs(n)-1)//2), abs(n)//2]. This `mod´ is denoted by `mods´ or `smod´ and called symmetric (or signed) modulo.

Let's see what happens if we use this mod operation.

def A103131(n):
    r = 1
    for k in (1..n):
        if gcd(n, k) == 1: r = smod(r*k, n)
    return r

0, -1, -1, -1, -1, -1, -1, 1, -1, -1, -1, 1, -1, -1, 1, ...

This is the way the Gauss-Wilson theorem is formulated. So now we can say:

A103131(n) = Gauss_factorial(n, n) signed_modulo n.

There remains a problem: As far as I can see the operation `smod´ is not implemented in Sage. Here is our version:

def smod(a,n): return a-n*ceil(a/n-1/2) if n <> 0 else a    

Now looking at -1, -1, -1, -1, -1, -1, 1, -1, ... it is natural to ask: For which n the Gauss_factorial(n, n) smod n is 1 and for which -1? The answer is given by the Gauss-Wilson theorem above.

More about the Gauss factorial can be found in A216919 and the

references given there.

The Scambler statistic on Dyck words

Oct 22, 2012

Last week I received an email from David Scambler.

<Scambler> -----------------------------

I hit a couple of your sequences when playing (as usual) 
with Dyck path statistics.
 
Consider the statistic S = NumReturns + NumHills – NumPeaks.
Count n-paths with S=-1, S=0, S=1 and S>=0.
 
   S   Seq
 =-1   A014531(n-2) [0,0,0],1,3,10,30,90,…    - Sloane
 = 0   A113682(n-2) [1,0],1,1,4,9,26,70,197,… - Barry
 = 1   A194588(n-1) [0],1,0,2,2,8,17,49,128,… - Luschny
>= 0   A189912(n-1) [1],1,2,4,10,25,66,177,…  - Luschny

</Scambler> -----------------------------

A Dyck path is a path in R² from (0, 0) to (2n, 0) consisting of a sequence of steps of length sqrt(2) and slope +/-1. These two types of steps are denoted by (1,1) and (1,-1), and called up steps and down steps, respectively.

• Number of hills is the number of up steps starting on the x-axis followed immediately by a down step.
• Number of peaks is the number of up steps that are followed by a down step.
• Number of returns is the number of points where D touches the x-axis y=0 and is not the origin.

WeWe will work here with Dyck words by letting 1’s represent steps in the direction (1, 1) and 0’s represent steps in the direction (1, -1).

def mpeaks(D): 
    return len(D.peaks())

def numret(D):   
    y = count = 0
    for d in D :
        y += 1 if d == 1 else -1
        if y == 0 : count += 1
    return count
    
def numhills(D):
    count = 0
    h = D.heights()
    for i in (1..len(D)-1):
        if h[i-1]==0 and h[i]==1 and h[i+1]==0: 
            count += 1
    return count 
    
def scambler_stat(n, S):
    count = 0
    D = DyckWords(n)
    for d in D:
        s = numret(d) + numhills(d) - numpeaks(d)
        if s == S: count += 1
    return count
    
def scambler_statSum(n):
    count_m = count_p = 0
    D = DyckWords(n)
    for d in D:
        s = numret(d) + numhills(d) - numpeaks(d)
        if s >= 0: count_p += 1
        else: count_m += 1
    return [count_m, count_p, count_p+count_m]   

print [scambler_statSum(n)[1] for n in (1..11)]
for S in (0..9): [scambler_stat(n, S) for n in (1..11)] 

S>=0 [1, 2, 4, 10, 25, 66, 177, 484, 1339, 3742, 10538]   A189912
S =0 [0, 1, 1,  4,  9, 26,  70, 197,  553, 1570,  4476]   A113682
S =1 [1, 0, 2,  2,  8, 17,  49, 128,  356,  983,  2759]   A194588
S =2 [0, 1, 0,  3,  3, 13,  27,  80,  207,  579,  1593]
S =3 [0, 0, 1,  0,  4,  4,  19,  39,  120,  310,   879]
S =4 [0, 0, 0,  1,  0,  5,   5,  26,   53,  170,   440]
S =5 [0, 0, 0,  0,  1,  0,   6,   6,   34,   69,   231]

print [scambler_statSum(n)[0] for n in (1..11)]
for S in (0..9): [scambler_stat(n, -S) for n in (1..11)] 

S<0  [0, 0, 1, 4, 17, 66, 252, 946, 3523, 13054, 48248]   A217539
S=-1 [0, 0, 1, 3, 10, 30,  90, 266,  784,  2304,  6765]   A014531
S=-2 [0, 0, 0, 1,  6, 25,  90, 301,  966,  3024,  9318]
S=-3 [0, 0, 0, 0,  1, 10,  56, 245,  938,  3318, 11166]
S=-4 [0, 0, 0, 0,  0,  1,  15, 112,  602,  2682, 10626]
S=-5 [0, 0, 0, 0,  0,  0,   1,  21,  204,  1344,  7017]

The Scambler triangle is the above array read by columns, S runs through (-n,..,n). The row sums are the Catalan numbers.

def scambler_row(n):
    return [scambler_stat(n, S) for S in (-n..n)]
fofor n in (0..8): scambler_row(n) 

Scambler's triangle

[n\k]-8,-7,-6,-5, -4,  -3,  -2,  -1,   0,   1,  2,  3,  4, 5, 6, 7, 8
---------------------------------------------------------------------
[0]                                    1,
[1]                               0,   0,   1,
[2]                          0,   0,   1,   0,  1,
[3]                     0,   0,   1,   1,   2,  0,  1,
[4]                0,   0,   1,   3,   4,   2,  3,  0,  1,
[5]           0,   0,   1,   6,  10,   9,   8,  3,  4,  0, 1,
[6]       0,  0,   1,  10,  25,  30,  26,  17, 13,  4,  5, 0, 1,
[7]    0, 0,  1,  15,  56,  90,  90,  70,  49, 27, 19,  5, 6, 0, 1,
[8] 0, 0, 1, 21, 112, 245, 301, 266, 197, 128, 80, 39, 26, 6, 7, 0, 1

Illustrating the row for n = 4:

[ScamblerWords(4, S) for S in (-4..4)]

[11010100] (()()())----
[11011000] (()(())) [11100100] ((())()) [11101000] ((()()))
----
[10110100] ()(()()) [11001100] (())(()) [11010010] (()())()
[11110000] (((())))
----
[10111000] ()((())) [11100010] ((()))()
----
[10101100] ()()(()) [10110010] ()(())() [11001010] (())()()
----
[10101010] ()()()()

I submitted this triangle because I think the connection with A014531, A113682, A194588, and A189912, which are defined in quite different ways, proves it as valuable. And it defines probably a new interpretation of the Catalan numbers. See A217540.

ToTo sum up, we can subsume the Scambler condition on Dyck paths in the following characterization:

def amblerWords(n, S):
    def characteristic(d):
        count = 1
        h = d.heights()
        for i in (1..len(d)-1):
            if d[i-1]==1 and d[i]==0: count -= 1
            if h[i]==0: count +=1
            else: 
                if h[i-1]==0 and h[i+1]==0: count += 1
        return count 
    
    count = 0
    D = DyckWords(n)
    for d in D:
        if S == characteristic(d): 
            count += 1 
    return count   

I also submitted a second sequence as a companion to A189912. In the Scambler interpretation A189912 is the "number of Dyck paths of semilength n which satisfy the condition: number of hills + number of returns ≥ number of peaks". The complement of this set, the paths which satisfy the condition: "number of hills + number of returns < number of peaks" now is A217539.

FoFor A217539 a simple explicit formula exist. A189912(n) = n*M(n-1) + M(n), M(n) the Motzkin numbers. Since the row sums of the Scambler triangle are the Catalan numbers C(n), we have A217539(n) = C(n) - (n-1)*M(n-2) - M(n-1), where the shift in the index takes the different offsets into account.

Pi, Euler and Bernoulli

Nov 02, 2012

Some say "Converting Pi to binary: DON'T DO IT". (Google will find the reason for you.) But recently I did and found out:

Dividing the Bernoulli number B(n) by the Euler number E(n) (n even of course) is the same as binary shifting Pi (i.e. dividing Pi by some power of 2). Well, not exactly, but very nearly. Here an example (Maple):p>

convert(evalf(-bernoulli(50)/euler(50),100),binary,50);
convert(evalf(Pi/2^101,100),binary,50);

1.1001001000011111101101010100010001000010110100011 x 10^(-100)    
1.1001001000011111101101010100010001000010110100011 x 10^(-100)    

<p>I have added a comment to the database. Riddle: To what sequence?

How to compute the 1000-th derivative of Lambert W(exp(x))?

Nov 13, 2012

Admittedly this is a somewhat esoteric question. But it entails some nice relations to the second-order Eulerian numbers and polynomials.

AsAssuming the definitions of the Eulerian numbers (as defined by Knuth) of first and second known order (they are given for the first order A(n, k) in A173018 and for the second order A2(n, k) in A201637) we can introduce the alternating second-order Eulerian polynomials as

A2(n,w) = Sum_{k=0..n} (-1)^k A2(n,k) w^k.

With Maple:

A2poly := proc(n,w) 
add((-1)^k*eulerian2(n,k)*w^k, k=0..n) end:

The first few polynomials are:

                        [0] 1
                        [1] 1
                     [2] 1 - 2*w
                 [3] 1 - 8*w + 6*w^2 
            [4] 1 - 22*w + 58*w^2  - 24*w^3 
       [5] 1 - 52*w + 328*w^2  - 444*w^3  + 120*w^4 
[6] 1 - 114*w + 1452*w^2 - 4400*w^3  + 3708*w^4  - 720*w^5 

Evaluated at w = 1 we find:

1, 1, -1, -1, 13, -47, -73, 2447, -16811, -15551, 1726511, 

This is a (nicer) version of A217538.

Now let's look at the n-th derivate of LambertW(exp(x)) at x=1.

L := n -> (D@@n)(LambertW@exp)(1);

FoFor n = 0...7 we find:

1, 1/2, 1/8, -1/32, -1/128, 13/512, -47/2048, -73/8192, 

To get rid of the denominator we multiply with 2^(2*n-1).

L2 := n -> 2^(2*n-1)*(D@@n)(LambertW@exp)(1);

1/2, 1, 1, -1, -1, 13, -47, -73, 2447, -16811, ...

Thus we suspect that the (n+1)-th derivative of LambertW(exp(w)) at x=1 for n>0 is the value of the n-th alternating Eulerian polynomial of second-order at x=1, up to an factor of 2^(2*n-1).

And this is indeed true, see for example the paper of  Corless, Gonnet, Hare, Jeffrey, and Knuth linked to in A001662.

NoNow the question was how to compute the 1000-th derivative of LambertW(exp(x)) at x=1. It was beyond my patience to compute this value via (D@@n)(LambertW@exp)(1) for any n>100 with Maple.

On the other hand using the above relation

n := 1000; A2poly(n-1,1)/2^(2*n-1); evalf(%);

Maple almost immediately answered

-979[..]313 / 574[..]688 = -1.70..*10^1992

Here [..] indicates many suppressed digits.

What does 'frac' mean?

Nov 17, 2012

Recently I stumbled upon this triangle:

 1
-1   2 1  -2   8
-1   2  -8   16
 1  -2   8  -16   128
-1   2  -8   16  -128   256
 1  -2   8  -16   128  -256  1024

def T(n,k):
    A005187 = lambda n: A005187(n//2) + n if n > 0 else 0
    return (-1)^(n-k)*2^A005187(k)

It has the keyword 'frac'. But why? So I looked up some 'offical' explanations of this keyword in the OEIS demonstration page.

Rational numbers (or fractions) are given as a pair of sequences, listing the numerators and denominators separately.

This is indicated by the keyword "frac", and the two sequences are cross-referenced to each other.

I think that the author means that sequences of rational numbers are given as a pair of sequences, not rational numbers.

But this definition did not help me much. So I looked up some 'official' examples in the index to fractions.

1, 1/2, 1/3, 1/4, 1/5, 1/6, 1/7, 1/8, 1/9, 1/10, 1/11, 
1/12, 1/13, 1/14, 1/15, 1/16 ... = 1/A000027

0, 1/2, 2/3, 3/4, 4/5, 5/6, 6/7, 7/8, 8/9, 9/10, 10/11,
1111/12, 12/13, 13/14, 14/15, 15/16, ... = A001477/A000027

Now I was really stunned. I never would have expected these pairs as 'frac' sequences. And of course neither A000027 (the natural numbers) nor A001477 (the nonnegative integers) have the keyword 'frac'.

I gave up, as I often do, to try to understand the inconsistent and to me incomprehensible definitions and conventions of the OEIS.

So I had to make up my mind when I am inclined to give the keyword 'frac' to a pair of integer sequences.

Given two sequences (a0,a1,a2,...) and (b0,b1,b2,...) we can consider the two mapsp>

( I) ((a0,a1,a2,..),(b0,b1,b2,..)) -> (a0/b0,a1/b1,a2/b2,..)
(II) (a0/b0,a1/b1,a2/b2,..) -> ((a0,a1,a2,..),(b0,b1,b2,..)).

<p>The first case, NxN -> Q, makes sense whenever bn <> 0 for all n. But even in this case I would not consider the sequence of rational numbers (a0/b0, a1/b1, a2/b2,...) a sequence which has to be distinguished by the keyword 'frac', and

the above pair A001477/A000027 is an example for this case.

The second case, Q -> NxN, is much more interesting. From the demonstration page:

The Bernoulli numbers are a good example. This is probably
the most important sequence of fractions in number theory.
It begins:

1, -1/2, 1/6, 0, -1/30, 0, 1/42, 0, -1/30, 0, 5/66, ...
 
The numerators and denominators are listed separately as 
sequences A027641 and A027642 as shown here.

Yes, I agree that this sequence deserves the attribute 'frac'. But why? Now I try to give my definition:

A A pair of integer sequences (a, b) has the attribute 'frac' if and only if at least one of the two sequences cannot be defined independently from the other.

In this definition the difficulty is shifted to the phrase "defined independently from the other". Admittedly this phrase is vague. So let me make it more concrete by a rule of thumb: The sequences a and b are dependent on each other if a or b cannot be derived from a/b without using the operators 'numerator' or 'denominator'.

As a practical rule of thumb I recommend avoiding the use of the operators 'numerator' or 'denominator' whenever this is possible. This will always lead to a better description of the sequence, 'better' here means a more intrinsic description.

To come back to the introductory example. It is called: "Table of the numerators of the higher order differences of the binomial transform of the dhava-Gregory-Leibniz series for Pi/4." A nice observation, no doubt worth a comment, but by no means a definition of A173755.

---

Let's look at a second, more interesting example. It is defined as

• a(n) = denominator of Π²ⁿ/(Γ(2n)*(1-2^-2n)*ζ(2n)).

This sequence has some nice characterizations. For example, as noticed by Sloane,

• (-1)ⁿ a(n) = numerator of Euler(2n+1, 1).

This looks like a genuine 'frac' sequence, both definitions refer to a projection operator, to the denominator and to the numerator operator, respectively. But is it really a 'frac' sequence?

Trying to compute this sequence by elementary means I came up with:

def A_list(n):
    T = [0]*n; T[0] = 1; S = [0]*n; k2 = 0
    for k in (1..n-1): T[k] = k*T[k-1]
    for k in (1..n):
        if is_odd(k): S[k-1] = 4*k2; k2 += 1
        else: S[k-1] = S[k2-1]+2*k2-1
        for j in (k..n-1): T[j] = (j-k)*T[j-1]+(j-k+2)*T[j]
      return [T[j]>>S[j] for j in (0..n-1)]

The clue was given by Ralf Stephan who observed that the sequence is the odd part of the tangent numbers (which is by far the best definition of this sequence).

However, the tangent numbers A000182 do not have the attribute 'frac'; does taking the odd part of these numbers turn the tangent numbers into a 'frac' sequence? Also notice that no 'numerator' or 'denominator' operator is used in the above script. Thus is this sequence rightly categorized by 'frac'?

Well, I think it is. The reason for me is the innocent looking '>>' in the script, which means division by some power of 2. So it is time to give a second, sharper definition of what I think the attribute 'frac' stands for.

A sequence of integers has the attribute 'ac' iff it cannot be computed by the operations '+', '-' and '*' alone. In other words: One cannot compute this sequence without using divisions.

A great advantage of this definition compared to the earlier one is that no reference to the existence of a second sequence is made.

A sequence of fractions f(n) has the attribute 'frac' iff at least one of the sequences numerator(f(n)) and denominator(f(n)) has the attribute 'frac'.

For example A001477/A000027 from the index of fractions is not a 'frac' sequence by this definition.

A consistent choice of some classical numbers.

I updated the links of Eulerian and Stirling numbers, a page where I discuss which variants of these classical numbers I prefer and why. I also like the intuitive notations invented by D. Knuth.

Power of primes versus 'prime powers'.

We all know what the primes are: the infinite sequence of integers 2, 3, 5, 7,...

ThThe integer power of a number n is also a common concept: Given a nonnegative integer m then n^m denotes n*n*...*n, where ... indicates that n is taken m times.

Thus we also know what the integer power of primes are, the set

IP = {p^m : p a prime and m a nonnegative integer}.

If we enumerate this set in the usual order we have

IP = {1, 2, 3, 4, 5, 9,...}.

The '1' is an element of IP because p^0 = 1 for all primes p. All this is elementary and needs no further comment.

Except one thing: There is also the notion of a 'prime power'. Now this invites confusion because the names 'power of primes' and 'prime powers' look so similar. But 'prime powers' is a different concept.

What 'prime powers' are is defined by the fundamental theorem of arithmetic. They are the building blocks of all the integers (except for 0).p> <p>Every integer n different from 0 can be written in the form

n = u * p_1^e_1 * p_2^e_2 *...* p_m^e_m

u is a unit (i. e. is either -1 or 1), the p_k are distinct prime numbers and e_k positive numbers. Up to a rearrangement of the factors this representation is unique. And it is this uniqueness property which gives this theorem the power as a fundamental tool of arithmetic.

The e_k are positive numbers because otherwise for example 10 = 1*2^1*3^0*5^1 and 10 = 1*2^1*5^1*7^0 would be two different prime factorizations of 10.

ThThe distinct "fundamental building blocks" p_k^e_k of all integers showing up in the theorem are the 'prime powers'. Although they obviously differ only by the exclusion of '1' from the integer powers of primes this exclusion is conceptually important. This is similar to the now generally accepted exclusion of the '1' from the prime numbers and the exclusion is for the same reason: to ensure the uniqueness property of the prime factorization.

Finally note that all major computer algebra systems have a query function 'PrimePower(n)'. And they understand this question in the sense of the fundamental theorem, not in the trivial sense. This means Pari, Magma, GAP and Mathematica return for 'isPrimePower(1)' the value 'False'.

And they certainly do not do so because the programmers of these software don't understand that n^0 = 1 for any nonzero number n as has been foolishly suspected recently.

To understand the term 'prime power' in the fundamental and not in the trivial sense is not a new invention. It is has been used in mathematics probably since L. Euler investigated the number phi(n) of positive integers not exceeding n which are relatively prime to n.

Since then this formula is part of mathematical folklore:

phi(p^n) = p^n (1 - 1/p)

This formula holds for p prime since p, 2p, ..., p(n-1)p are the only ones of the p^n positive integers ≤ p^n not prime to p^n.

Here p^n is assumed exactly in the sense of a 'prime power' not in the sense of 'power of a prime'. In fact this formula is false if interpreted in the latter sense. For instance phi(3^0) = 1 whereas 3^0*(1-1/3) = 2/3.

This is only the beginning of a tower of relations which become invalid if the exponent n=0 is not excluded. For example John Cremona gave the following case in a recent thread at Sage-devel:

"A positive integer n is a prime power if and only if nZ is a primary ideal of Z, and the definition of primary ideal (as with prime ideal) explicitly requires the ideal to be proper."

So it is not surprising that Cremona exclaimed: "No! No! No! please don't let Sage match OEIS here!"

And indeed the definition of OEIS A000961 follows the naive definition and not the mathematical one, (or discloses a lack of understanding of the latter). Even the web pages which explain the term 'prime power' and refer to A000961 like Wikipedia or MathWorld actually display the sequence [2,3,4,5,7,8, ...] and not the data of A000961.

Meanwhile, and mainly due to the efforts of Ralf Stephan, the name and the comments of A000961 are vastly improved. It now correctly says: "Powers of primes. Alternatively, 1 and the prime powers (p^k, p prime, k ≥ 1)." Nevertheless the fact remains that the OEIS to the day does not have the sequence of the prime powers 2, 3, 4, ..., which is a quite remarkable fact.