Coding Challenge 161: Estimating π from Random Numbers with Euclid's Algorithm

that's right i'm here standing on platform 3.14 the pie train is rolling into the station full of delicious pies and i'm here to celebrate pie with me [Music] that you can eat my favorite treat circumference divided by diameter this year i am going to look at a technique for estimating the number of pi with random numbers this is a direct result of me re-watching one of my favorite internet videos from stand-up maths called generating pi from 1 000 random numbers now in that video matt parker the host of stand-up maths physically rolls dice picking a thousand numbers to estimate the value of pi i am going to write some code it'll probably take me actually longer than physically rolling dice but that's what i'm going to do however i'm going to give myself a little restriction i will not be using the random function in p5.js no rather i'm going to pick random numbers from my favorite book 1 million random digits and with 100 000 normal deviants and see how close can i get to the number pi from all of the digits in that book so how's this all work let's say i pick two random integers call them a and b while in an ideal world i would pick them between one and infinity i need to establish some kind of range so that i can actually work with this stuff in code i will choose the upper bound of a hundred thousand because in the random numbers book the digits are grouped in groupings of five so between zero and ninety nine thousand nine hundred ninety nine i don't want to pick the number zero you'll see why in a second so i my range will be between one and one hundred thousand the probability that these two numbers are co prime is six divided by pi squared huh why is that that's crazy well this is one of the wonderful things about the number pi it shows up all over the place in mathematics and it's just delightful if you want to understand a bit more about the derivation and proof of this number this is where i would suggest you check out matt parker's video he goes through this in detail for me in this video what i want to focus on is this concept of coprime and how do i write some javascript code or really in any programming language to determine if two integers are co-prime what does it mean for two numbers to be co-prime in the first place let's say i pick two numbers 15 and 16. well for any number i could list all of the factors of that number what other integers that are small in that number divide into that number with a remainder of 0. with 15 i've got 5 3 and 1. with 16 i've got 8 four two and one looking at these two numbers the only factor that they share is one if two numbers don't share any factors other than one they are co-prime another way to say co-prime is greatest common divisor so the greatest common divisor of 15 and 16 is one that's the highest number that divides into both of them let's change the 16 now to 21. the factors of 21 are 7 3 and 1. oh look at this they share a greatest common divisor of three i should not be erasing this with my fingers this means they are not co-prime and in fact matt parker calls this cofactor so if i do have a sequence of random numbers which i do if i could just look at pairs of those random numbers determine whether they're co prime or cofactor look at the ratio of co prime to cofactor solve for pi then i could estimate the value of pi which begs the question how do i determine if two integers are co-prime or not i've got some wonderful news for you there's a well-known algorithm known as euclid's algorithm to determine the greatest common divisor of two integers and you can implement that algorithm with recursion i can't wait to look at the youtube retention stats in this part of the video unfortunately something terrible has happened i have to erase this so i can have some room to explain the algorithm i'll provide links in the description where you can read up more on both the derivation of this probability as well as euclid's algorithm or the euclidean algorithm but i'm going to try to describe it to you in my own words here in this video the first thing i need to do is determine which number is larger and for the sake of this moment i'm just going to assume that a is greater than b in my code i can figure something out to determine which is the bigger number i then need to take a divided by b and determine what the remainder is i'll call that r if by the way this remainder is zero then those numbers are not co-prime if i take 15 divided by 3 i'm going to get 5 remainder 0. 15 and 3 are obviously not co-prime because they both share the factor 3. oh i wasn't including the number itself in my earlier list of factors which clearly i should have if the remainder is one they're definitely co-prime take the numbers 15 and 16 again 16 divided by 15 is one remainder one those two numbers are co prime what if my numbers however are say 22 and five if i take 22 divided by five i get four remainder 2. as i'm moving through this you're going to notice that i'm not bothering to write the 4 here it's because this division actually doesn't matter it's only the remainder that matters and the term the operator that you use in javascript for remainder of division is the percent sign is that modulo or modulus oh the comments [Music] i got it the percent sign is the modulo operator the two the result is the modulus of 22 modulo 5. will i remember this the next time i come to it in a video most definitely not i should also mention that i'm kind of mixing things up here i want to determine whether two numbers are co-prime by calculating the greatest common divisor and if the greatest common divisor is one then they are co-prime and so to be a bit more clear if a modulo b is remainder zero then that greatest common divisor is b itself so this should say that the gcd the greatest common divisor is b if it's not zero i need to keep going and the way that i keep going is by then taking b modulus that remainder and i could call this r zero and then call this r1 and so oh sorry call this r0 b modulus r0 is r1 and i'm going to do this over and over continuously again so now r 0 modulus r1 equals r2 and i keep going and going and going if ever this becomes zero then whatever the value is here is the greatest common divisor and by the way if that remainder is ever one any number modulus one is going to be i said it wrong any number modulo one [Music] the modulus of any number and one will be zero so you that in that case the greatest common divisor will be that number one i wonder if looking at um two example numbers will be helpful in explaining this because boy am i doing a terrible job let's try 93 and 10. 93 modulo 10 is 9 remainder 3. that's not zero so keep going ten modulo three is three remainder one so that's not zero keep going three modulo one is zero hey that's zero this is now the greatest common divisor one guess what 93 and 10 co prime let's change this to 95 95 divided by 10 is 9 remainder 5. so 95 modulo 10 is 5 10 modulo 5 is oh wait a second zero so guess what the greatest common divisor is five these are cofactor 95 and 10 not co-prime i feel ready do you feel ready it's okay if you don't feel ready but i feel ready to go and write this algorithm in the code now here i am in the p5 web editor i'm going to write a new function i'm going to call it gcd for greatest common divisor i'm going to pass it two numbers a and b and i'm going to implement euclid's algorithm in that function the first thing that i want to do is figure out which is the bigger number by the way if the numbers are the same they're not co-prime because a modulo a would be zero and then the greatest common divisor would be itself okay got that so first i just need to determine is b greater than a if b is greater than a i just want to like swap those numbers i look forward to all your suggestions about how to handle this more elegantly but this now is keeping both numbers but always guaranteeing that a is greater than b the first thing i want to do is calculate the modulus of a modulo b i'm going to store that in a variable called r if r is 0 return what b that's the greatest common divisor otherwise what do i do well now i want to check b modulo the remainder but ultimately why am i what is that doing that's checking the greatest common divisor of b comma r this is where the recursion comes in return g c d of b comma r get rid of draw let's confirm that it's working correctly what were the numbers that i tried 93 and 10 and 95 and 10 so i should see in the console a greatest common divisor of one and of five it worked it's like very rare that i'll write something it actually worked the first time all right a little birdie just said something to me in my ear which is like i forgot that this sort of like fancy javascript destructuring thingymabobber allows you to swap two variables like this let's make sure it still works indeed it does we're good to go i guess if i really want to be sure that this works i better also check 10 comma 93 and 10 comma 95. okay i'm really good to go here what's next so i need my random numbers i need those 100 000 pairs of random numbers between 1 and 100 000 and i want to get them from the book that i don't actually have with me but that will somehow be appearing like as if i'm holding it oh look at that post-production effect there luckily for us you can download all the digits i'll put a link in the description where you can get them yourselves but i've already done that and i've put them into a file called milliondigits.txt it looks like this it has that full sequence and in p5 i can use the load strings function to get them a little oddity about load strings not that strange actually it loads everything into an array i just want the first element of that array which has all billion digits in it and then i want to split those up into chunks of five all right i'm looking through the string documentation to try to find a function that might work for me pad and repeat replay slice yeah all right all right maybe there's a better way to do this i'm going to try using slice which allows me to slice out a section of the string so let me make an array called randoms for all my random numbers while i have some digits left the number i want is a slice starting with index zero up until but not including index five zero one two three four that's five things and then digits is just the rest of it strings are immutable so this slice doesn't actually change the string digits so i'm pulling the slice out and then i'm slicing digits the beginning part off of the string digits add that number to the array and if i've done everything correctly i should have a array of 200 000 numbers perfect 200 000. now actually those are not numbers i've created little string slices but i want to work with them as numbers if i'm going to be calculating the greatest common divisor so let's add a parseint of num and just to be sure i'm going to take a look at that full array itself you can see there is my random number sequence now i think it would be fun to see fun is relative here the estimation of the value pi converging as i go through that full list of numbers so i'm going to make use of the p5 draw loop and i need two values i need the coprime total our co-prime count and the total that's going to allow me to compute that probability of any two numbers being co-prime and then solve for an estimated value of pi just to sort of know that it's working i'm going to add no loop to do draw just once this time i'm going to have a global variable called index which is keeping track of where i am in the random sequence so a is random's index b is the next number random's index plus one and if a and b have a greatest common divisor of one increase the coprime count and always increase the total count and let me just console log everything so i can see if this works i love p5 so much it seems that you may have accidentally written randoms instead of random just letting you know maybe you wanted to use random no i didn't i wanted to use randoms but i forgot to make randoms a global variable and i probably shouldn't have called it randoms because it's very confusing because randoms is very similar to random so i really really appreciate you p5 error message i'm going to call it random sequence make it a global variable and try this one more time here are the two numbers i got 10 097 32 533 their greatest common divisor is one therefore they are co-prime dare i take out this no loop i'm getting the same number over and over again i forgot that i need to increase the index so instead of increasing the index by 1 i'm taking them as pairs so i'm going to go up by 2 and if index equals the end of that array no loop console.log complete oh look i saw one that had 11. that was exciting so i it appears to be working you know i have to hope that all of my code and math is right but i'm seeing all the numbers coming out in pairs i'm seeing their greatest common divisor i think i'm good i'll wait just to see until it gets to the end so that i can see that it says complete or if i get an error message this seems to be taking an awfully long time i'm doing this a hundred thousand times a hundred thousand at 60 frames per second i don't know a lot of minutes oh this is going to take a while i'm not going to wait till the end i'm going to hit stop now i could do more than one per frame we'll get to that the next step though is that i need to figure out how to solve for pi so this is now a variable in my system the probability of the greatest common divisor of two integers a and b equaling one is six divided by pi squared let's say that this probability is x and in my code x is coprime count divided by total so coprime count divided by total equals 6 over pi squared and this is not pi squared actually it's my variable that i want to solve for it's my pi estimate i'll just call that p-i-e for pi estimate i might just have pi for lunch today and i could rewrite this thing to say pi squared equals 6 times the total divided by coprime count or pi the pi estimate is the square root of all of this i think i did that right let's add that into the code let pi for pi estimate square root of 6 times total divided by co prime count let's be really really safe here and add some extra parentheses although i don't think they're technically necessary let's get rid of this console log and let's just watch console log pi oh let's get in there come on 3.1 let's show me a 3.1 yeah 3.1 can i get a 4 can i get a four let's let's actually draw this into the canvas oop erasing the background might be nice and here we go let's also make a little progress bar so we can see how things are going i haven't tested this yet but i'm just drawing a rectangle that's the full width of the window and then drawing one that's the percentage along the way with a fill wow this is really really really slow let's do more than one pair per frame so that index increasing by two has got to be up here then i can do let's try five per frame and i should do this check here so if i get to the end no loop complete break out of the loop and let's do 500 per frame ooh maximum call stack size exceeded let's take a look at what a and b are it really shouldn't like the numbers can't be so big that it can't possibly do the recursion to find the greatest common divisor i'm going to console log them oh look there's a zero oh my goodness of course i said this earlier in the video let's think back to that time i don't want to pick the number zero you'll see why in a second so i my range will be between 1 and 100 000. the range has to be between one and a hundred thousand i can't use a zero in this algorithm you can't have a greatest common divisor with zero so i can fix that by just saying plus one and plus one and now let's run it and there we go can i get a four yes ah look at this there's my estimated value of pi with all a hundred thousand pairs of five digit numbers in the one million random digits book i would encourage you to try it different ways certainly i could start from a different place in the book i could look at every page of the book and see which page of the book has the closest approximation of pi but i have another thought experiment i would like to try what if instead of using those million random digits what if i used the digits of pi itself here's a text file with the first 1 million digits of pi let's get rid of the 3 and the decimal place save that and upload it to my sketch now all i need to do is replace the file name from million digits to pi 1 million and let me comment out this console log i don't want to be doing that anymore and there we go i'm gonna stop here i'm very pleased with the fact that i have taken the digits of pi divided them into groups of five found out whether those groups of fives are co-prime or not and then use that ratio of co-prime to not co-prime to estimate the digits of pi itself thank you thank you acting thanksgiving i hope you enjoyed this pi day video i feel like there's a lot of things that you could do to improve on what i've done certainly you could try to make the code more efficient use more random numbers pick your random numbers from some type of like a quantum thermodynamic machine that's true randomness is that a thing i don't know also i think there's some visualization possibilities is there a way to visualize the process of euclid's algorithm looking on the wikipedia page you can see some nice diagrams there that might be a starting point i would love for you to make something please share it with me there's a link in the video description to the coding train website where you can add your twist onto this uh coding challenge and i will see you i'll hopefully see you on this internet thing before next pie day but i'll be back every year on the dot march 14th it's pi day i'm gonna go and have blueberry is my favorite kind of pie just putting that out there and i'll see you later i can't stop these videos i got to get on the train uh it was track 3.14228 apparently circumference divided by diameter day

Channel: The Coding Train

Views: 32,203

Keywords:

Id: EvS_a921dBo

Length: 24min 19sec (1459 seconds)

Published: Sun Mar 14 2021