Trying Julia for Analytic combinatorics

General Usage
, , ,
1

Hello, fellow colleagues in the Julia community.

This question may not be a good fit for the Julia Discourse, but any help is appreciated.

I was exploring Analytic Combinatorics (by Flajolet & Sedgewick). The part about binary words (words over a binary alphabet). contains the following generating function. in eq. 54 (pg. 52)

image

The following part in page 53 (equation without number) describes the probabilities associated with the maximum length of runs of consecutive letters:

image
image739×326 50.2 KB

I used SymPy.jl to implement it:

using SymPy
#using FastRationals

@vars z r n

# Binary words that never have more than r consecutive identical letters is found to be (set α = β = r)
# Flajolet pag. 52
w_rr(r,z) = (1-z^(r+1))/(1-2z+z^(r+1)) # OGF
w_rr(r,z) = sum(z^x for x in 0:r)/(1 - sum(z^x for x in 1:r)) # Alternate form

function n_words(r)
    coefs = collect(series(w_rr(r,z),z,0,r+1),z)
    coefs.coeff(z,r)
end

function p_k(k,n)
    #a =  FastRational{Int128}(1/(2^n))
    #a = Rational{BigInt}(1/(2^n))
    a = 1/(2^n)
    result_sym = a*(n_words(k) - n_words(k-1))
    return(result_sym)
end

Testing for k=6 and verifying the Taylor series for r=15:

julia> p_k(6,n)
    -n
32⋅2
julia>series(w_rr(15,z),z,0,15+1)
             2      3       4       5       6        7        8        9         10         11         12         13          14        
1 + 2⋅z + 4⋅z  + 8⋅z  + 16⋅z  + 32⋅z  + 64⋅z  + 128⋅z  + 256⋅z  + 512⋅z  + 1024⋅z   + 2048⋅z   + 4096⋅z   + 8192⋅z   + 16384⋅z   + 32768

  15    ⎛ 16⎞
⋅z   + O⎝z  ⎠

My problems:

1 - Running p_k(6,200) directly returns ∞ or another undefined value, even in with FastRational{Int128}(1/(2^n)) or Rational{BigInt}(1/(2^n)) that are commented out. I can run p_k(6,n) passing symbolic n and then evaluate the output directly in the REPL (julia>32*2^(-200)) though.

2 - The result seems to be in accordance with the formula (1/2^n)*([z^6]W<6,6> - [z^5]W<5,5>), but it is quite different from ‘template’ for k=6. The coefs. in the Taylor expansion seem to be just [z^n]W<k,k> = 2^n. Something must be wrong, since monotonically increasing values for k would not give the P(K) distribution presented in the table for integers k in [1,12].

3 - I intend to calculate probabilities and test hypotheses about structures in a recurrence matrix (trapping time, vertical, diagonal), which are also binary words. Any thoughts on that?

2

Maybe Rational{BigInt}(1/(2^n)) doesn’t give the expected result, as the 2 in 2^n is parsed as an Int64 and when raised to 200 gives 0 (i.e. zero(Int)) making the whole expression 1//0. The intended value can be obtained with 1/BigInt(2)^200:

julia> Rational{BigInt}(1/(2^200))  # bad result
1//0

julia> 1//(BigInt(2)^200) # good result
1//1606938044258990275541962092341162602522202993782792835301376
1 Like
3

Thank you!
It worked like a charm.

Do you have any idea about my 2nd and 3rd questions?

4

I fixed it thanks to RicBit.

W_n<k,k> notation stands for [z^n] Taylor(W<k,k>). That is, we need the 200th coefficient of the Taylor expansion of W<k,k>.

This version works:

function n_words(r;n_tot=200)
    coefs = collect(series(w_rr(r,z),z,0,n_tot+1),z)
    coefs.coeff(z,n_tot)
end

function p_k(k,n)
    #a =  FastRational{Int128}(1/(2^n))
    #a = Rational{BigInt}(1/(2^n))
    a = 1//(BigInt(2)^n)
    #a = 1/(2^n)
    result_sym = a*(n_words(k) - n_words(k-1))
    return(Float64(result_sym))
end

julia>for i in 3:12
    print(p_k(i,200))
    print("\n")
end

# Correct results
6.549422864002916e-8
0.0007075224314585063
0.033979610958548824
0.1660282049825555
0.2574174128392719
0.22352195216110582
0.14594875198080856
0.08292457405578264
0.0440234813776148
0.022607888203734293
5

This is nice, but I don’t think using SymPy is really using Julia :wink:

Maybe GitHub - JuliaDiff/TaylorSeries.jl: Taylor polynomial expansions in one and several independent variables. can help. (Disclaimer: I am one of the co-authors.)

2 Likes
6

I gave it a shot, but it seems that there’s something wrong with the types again.
The results are blowing up.

using TaylorSeries

w_rr(r,z) = (1-z^(r+1))/(1-2z+z^(r+1))
n_words(r;n_tot=200) = getcoeff(taylor_expand(z -> w_rr(z,r),order=n_tot+1),n_tot)

function p_k(k,n)
    #a =  FastRational{Int128}(1/(2^n))
    #a = Rational{BigInt}(1/(2^n))
    a = 1//(BigInt(2)^n)
    #a = 1/(2^n)
    result_sym = a*(n_words(k;n_tot=n) - n_words(k-1;n_tot=n))
    return(Float64(result_sym))
end

julia>for i in 3:12
    print(p_k(i,200))
    print("\n")
end

3.1854727641374944e6
5.919210583485605e18
2.9142600879986195e27
1.2941726894119323e34
2.808302147070106e39
7.566803340725184e43
4.4580962403509384e47
8.293787625532763e50
6.26377090076716e53
2.29651524957329e56