Calculating Value of Zeta function using Julia - Part1

Author: Kazi Abu Rousan

Where are the zeros of zeta of s?

G.F.B. Riemann has made a good guess;

They're all on the critical line, saith he,

And their density's one over 2 p log t.

Source https://www.physicsforums.com/threads/a-poem-on-the-zeta-function.16280/

If you are a person who loves to read maths related stuff then sure you have came across the words Riemann Zeta function and Riemann Hypothesis at least once. But today we are not going into the Riemann Hypothesis, rather we are going into the Zeta Function. To be more specific, we will see how to calculate the value of zeta function using a simple Julia Program only for $Re(input)>1$.

What is Zeta Function?

What is the Zeta Function?

The Riemann zeta function $\zeta (s)$ is a function of a complex variables $z = \sigma + i t $. When $Re(z) = \sigma >1$, we can define this as a converging summation given by,

$$ \zeta(z) = \sum_{n = 1}^{\infty} \frac{1}{n^z} = \frac{1}{1^z} + \frac{1}{2^z} + \frac{1}{3^z} +\cdots $$

This definition is so simple right?, Here $z$ is a complex number. If it is taken as real, then we will get many famous series like,

Harmonic Series(Diverge as $Re(z)=1$): $H = \zeta(1) = \frac{1}{1} + \frac{1}{2} + \frac{1}{3} + \frac{1}{4} + \cdots$

Basel problem Series(Converges): $ \zeta(2) = \frac{1}{1^2} + \frac{1}{2^2} + \frac{1}{3^2} + \frac{1}{4^2} + \cdots = \frac{\pi^2}{6}$

and many more.

There are many methods to calculate the values of this function for each value of $n$. You can surely do this by hand although I can't suggest that. The easiest way will be to utilize program.

Code to find the value of zeta for Re(z)>1

Before writing the code, you should remember that the sum is actually infinite. Hence, our function will give us more and more fine result as we include more and more terms.

function ζ(z,limit=1000)
    result = 0
    for i in 1:limit
        result += (1/i)^z
    end
    return result
end

This is the function. How simple!!.. and If we use julia as you can see, we can actually use $\zeta $ symbol. Pretty cool right?

ζ(2,10_0000_000)
#Output: 1.644934057834575
#where pi^2/6 = 1.6449340668482264
#pretty close

It's not all. We can actually apply this function to find the value of zeta function for complex inputs too.

ζ(2+5im)
#Output: 0.8510045028264933 + 0.09880538410302253im

When we use complex inputs, there is a beautiful hidden beauty. Let's see that using a plotting library called Plots and we also have to rewrite our function a little bit.

function ζ(z,limit=1000;point_ar = false)
    points = ComplexF64[0]
    result = 0
    for i in 1:limit
        result += (1/i)^z
        if point_ar
            push!(points,result)
        end
    end
    return result, points
end

Now, we can use this to plot.

z = 2+5im
result, points = ζ(z;point_ar = true)
plot((points),color=:blue,width=3, title="Zeta function spiral",framestyle= :origin,label="$z")
scatter!((points),color=:red, label="Points")S

The output is:

The complex power rotates each point.
Zoomed Image of the spiral

Who would have thought that there will be something like this hidden.

Now, If we apply this function to all possible points of the NumberPlane, then we will get a mesmerizing pattern.

Beautiful!!! But why is it feels like incomplete?

Looking at the image it feels like It's begging to be extended to the other portion. This extension is done using something we called Analytic Continuation. We will not go into much detail here.

If you want to know about the remaining story visit by lecture here:

This is all for today. Why not try to extend the idea to the other side?

Hope you learnt something new.

A Probability Birthday problem along with Julia Programming

Probability theory is nothing but common sense reduced to calculation.

Pierre-Simon Laplace

Today we will be discussing a problem from the second chapter of A First Course in Probability(Eighth Edition) by Sheldon Ross.

Let's see what the problem says:

Describing the Problem

The problem(prob-48) says:

Given 20 people, what is the probability that among the 12 months in the year, there are 4 months containing exactly 2 birthdays and 4 containing exactly 3 birthdays?

So, what do we mean by this? Let's suppose we have $20$ people, named $p_1,p_2,\cdots , p_{20}$ and we will be representing the months by numbers like, January$\to 1$, February $\to 2$ , March $\to 3$ ,$\cdots $,December $\to 12$.

Then, One possible combination of birthday can be,

People Name Birthday Month
$p_1$$01$
$p_2$$05$
$p_3$$09$
$p_4$$04$
$p_5$$07$
$p_6$$12$
$p_7$$03$
$p_8$$09$
$p_9$$06$
$p_{10}$$11$
$p_{11}$$12$
$p_{12}$$01$
$p_{13}$ $07$
$p_{14}$ $02$
$p_{15}$ $09$
$p_{16}$ $12$
$p_{17}$ $02$
$p_{18}$ $11$
$p_{19}$ $06$
$p_{20}$ $08$

In this case, January ($01$) has $2$ birthdays, February($02$) has $2$ birthdays. Similarly we can find the number of birthday in each month. But in our question it is said that one possible arrangement exist where any four months (suppose January, February, March and April) each holds $2$ birthday (in total in this group there are $4\times 2 = 8$ birthday) and another four months (suppose June, July, September and December) each holds $3$ birthday (in total this group has $4\times 3 = 12$ birthdays). We have to find the probability of existence of such groups.

I hope the problem is clear. Let's see how to solve this.

Solution - 1

To solve this problem one must notice the most important point. This problem treats each people in same footing, means they are indistinguishable. Normally, we humans are distinguishable as you can see people's face or maybe use their name to make difference.

So, The number of all possible combinations (outcome) = $12^{20}$.

Now, for our particular case, we must group $4$ months from $12$ months where each contains $2$ birthdays. This can be done in ${12 \choose 4}$ ways. After this again we have to choose $4$ months from remaining $8$, where each contains $3$ birthdays. This can be done in ${8 \choose 4}$ ways.

So, we can choose months (or the month combination is ) in ${12 \choose 4}{8 \choose 4}$ ways. Using multinomial theorem, we can write this number equals to ${12 \choose 4 4 4}$ (we will not use this).

Now, we have $20$ people and we have to to assign months with them.

Remember we have to take $3$ people and give them one single month. This can be done in ${20 \choose 3}$ ways. Again we have to take $3$ people and group them. This time the possible ways are ${17 \choose 3}$. Again doing the same gives us ${14 \choose 3}$ ways (this is done 4 times as we have $4$ such months which contains $3$ birthdays). So, for this group the number of possible combinations is ${20 \choose 3} {17 \choose 3} {14 \choose 3} {11 \choose 3} $.

Now, we still have $4$ months where each of those contains 2 birthdays. This gives a total of $ {8 \choose 2} {6 \choose 2} {4 \choose 2} {2 \choose 2} $ ways.

So, the total numbers of ways of assigning $20$ people to those months is $$ {20 \choose 3} {17 \choose 3} {14 \choose 3 } {11 \choose 3} {8 \choose 2} {6 \choose 2} {4 \choose 2} {2 \choose 2} = \frac{20!}{3!\cdot 17!} \frac{17!}{3!\cdot 14!} \frac{14!}{3!\cdot 11!} \frac{11!}{3!\cdot 8!} \frac{8!}{2!\cdot 6!} \frac{6!}{2!\cdot 4!} \frac{4!}{2!\cdot 2!} \frac{2!}{2!\cdot 0!} = \frac{20!}{12^4} $$

So, The probability = ${{12\choose 4}{8\choose 4}\frac{20!}{12^4}\over 12^{20}} = 0.00106042009902$

Easy right? This can also be solved using a analogy of dice with $12$ sides.

Solution-2

In this method we consider a dice with $12$ sides. Each person's birthday is within 12 months and each of them are independent. So, If that dice is rolled 20 times, then what is the probability to get same number $3$ times for $4$ different number and same number $2$ times for other $4$ different number, i.e., what is the probability to get a state like $11223344666777888999$?

This can be solved in $4$ simple steps.

  1. Among $12$ faces, $8$ of them can be taken in $a={12 \choose 8} = 495$ ways.
  2. Among those $8$ faces, $4$ are groups of $2$'s and $4$ are groups of $3$'s. The number of ways = $b = {8 \choose 4}{4 \choose 4} = 70$.
  3. If we roll it $20$ times or take $20$ dices, then there are $4$ sets of $3$ and $4$ sets of $2$. In this case, the possibilities $c = \frac{20!}{(3!)^4 (2!)^4}$.
  4. The sample space's element number $d = 12^{20}$

Hence, probability $=\frac{a\times b \times c}{d} = 0.0010604200990$

Whenever I solve this type of problems, It always feel like may be !! maybe!! I have made some mistake. To be sure that we are correct, we should perform experiment with a dice or may be take many groups of $20$ people and calculate probability. But It is practically impossible!!

Not quite... We can always use programming language to do the experiments and verify.

Verifying using Julia

First we have to know few tricks.

arr = [2,6,2,10]
#Creates an array

Now, how can we know if there are two 2's inside of this array?, an easy solution will be

arr.==2
#Check if any element equal to 2
#This gives an output as:
#BitVector: [true,false,true,false]

And then we use sum

sum(arr.==2)
#return number of times 2 is there
#output: 2

Using this we write a simple code:

begin
    function prob(n::Int)
	favour = 0
	for ex_num=1:n
		months = rand(1:12,20)
		nums = [sum(months.==i) for i=1:12]
		if sum(nums.==3)==4 && sum(nums.==2)==4
			favour += 1
		end
	end
	probability = favour/n
    end
end

This function beautifully give you the value. We need more sample, so I have used $n=1,00,00,000$.

This gives the result:3

prob_10_000_000 = prob(10_000_000)
#output: 0.0010616

Beautiful isn't it?, just using few lines of code we have verified our answer.

I encourage you to play with the code.

This is all for today. I hope you have learnt something new.

A code to find Primes - Sieve of Eratoshenes

To some extent the beauty of number theory seems to be related to the contradiction between the simplicity of the integers and the complicated structure of the primes, their building blocks. This has always attracted people.

A. Knauf from "Number theory, dynamical systems and statistical mechanics" 

This quote is indeed true. If you just think about the simplicity of integers and complexity of primes, it's really horrifying. It's almost like how easy and intuitive it is to analyse big "classical objects" but extremely un-intuitive to discuss about "quatum objects".

Today's blog will be a bit special. Today we will discuss about a simple algorithm which is one of the most easiest algorithm which find primes in a certain range. This is called The Sieve of Eratosthenes.

Main Idea of Eratosthenes

The main idea behind this algorithm is very simple. Let's try to understand this using an example.

Suppose, you want to find all primes within $20$.

This is basically our algorithm. Below you can find the algorithm for $30$.

Sieve of Eratosthenes
Code to perform this

The code to perform this is quite easy, specially in python, because of the filter function. Let's see the code;

n = 20
all_list = [i for i in range(2,n+1)]
factor_2_remove = filter(lambda m: m %2 !=0, all_list)
factor_2_remove = list(factor_2_remove)
factor_2_remove.append(2)
print(factor_2_remove)
#output:[3, 5, 7, 9, 11, 13, 15, 17, 19, 2]

As you can clearly see, first we have created a list containing a list from $2$ to $n$, then we use filter. This function filter out the numbers which doesn't satisfy the given. Here the condition is $m \not\equiv 0$ mod $2$. This is written as $m%2 !=0$ in python. After using this we again have to add $2$ as it also doesn't satisfy the condition so, it also is removed.

So, to find all primes within any range, we simply do this for all number, not just $2$. You can use a simple for loop for this.

n = 100
all_list = [i for i in range(2,n+1)]
for i in len(all_list):
    first_ele = all_list[0]
    all_list = filter(lambda m: m %first_ele !=0, all_list)
    all_list = list(all_list)
    all_list.append(first_ele)
print(all_list)
#Output: [2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97]

See!, how easy it is.. But we can do even better. Notice, here we have to run the loop $\pi(n)$ times. But we don't have to do that.

Note: $\pi(n)$ is a function, which represent number of primes in between $2$ to $n$.

Remember: If an integer $a>1$ is not divisible by any prime $p\leq \sqrt{a}$, then $a$ is of necessity a prime.

Using this we can increase the efficiency of the program. Let's try to understand this using an example.

Suppose, we want to find primes up to $100$. Then, we create a list up to $100$, i.e., $2,3,4,5,6,7,\cdots ,98,99,100$. Now, we take $2$ and remove all multiples of $2$ from the list, except $2$ itself. We again do the same for $3$, $5$,$\cdots $. But notice $\sqrt{100}=10$, so we just do this up to $10$, i.e., up to $7$. And this will generate all primes within $100$.

If we do this, now the we only have to run the loop for $\pi(n)$ times. Let's see how much the speed has increased. I will now make two functions, one which runs the loop $\pi(n)$ times and one which runs $\pi(\sqrt{n})$ times.

import math as mp
def with_just_sqrtn(n):#uses sqrt n
    all_list = [i for i in range(2,n+1)]
    root_n = mp.ceil(mp.sqrt(n))
    for i in range(n):
        first_ele = all_list[0]
        all_list = filter(lambda m: m %first_ele !=0, all_list)
        all_list = list(all_list)
        all_list.append(first_ele)
        if first_ele == root_n:
            break
    return all_list
print(with_just_sqrtn(100))
#output: [2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97]

See, the output is same as before.

The previous one can also be written as a function.

def with_all(n):
    all_list = [i for i in range(2,n+1)]
    for i in range(n):
        first_ele = all_list[0]
        all_list = filter(lambda m: m %first_ele !=0, all_list)
        all_list = list(all_list)
        all_list.append(first_ele)
    return all_list

Now, let's find which one takes much time. We will use $n=1,00,000$

n = 1_00_000
import time
start_t = time.time()
with_all(n)
print("Time = %s seconds"%(time.time()-start_t))
#Output: Time = 131.59189105033875 seconds

Just close to $2$ min $11$ sec. Now, let's run the one with small loop number.

start_t = time.time()
with_just_sqrtn(n)
print("Time = %s seconds"%(time.time()-start_t))
#Output: Time = 0.1107027530670166 seconds

What the !!, just 0.11 sec?

So, i hope now you can see, how much faster it is compared to the previous one.

This is all for today, I hope you have learnt something new.