{"id":117,"date":"2023-04-18T02:00:00","date_gmt":"2023-04-18T02:00:00","guid":{"rendered":"https:\/\/renegrothmann.de\/?p=117"},"modified":"2023-05-09T20:08:19","modified_gmt":"2023-05-09T20:08:19","slug":"computing-primes-with-python-and-emt","status":"publish","type":"post","link":"https:\/\/renegrothmann.de\/?p=117","title":{"rendered":"Computing Primes with Python and EMT"},"content":{"rendered":"\n
Another old posting from my previous blog.<\/p>\n\n\n\n
\n\n\n\n
Continuing to get acquainted with Python, I looked at the book \u201eIntroduction to Data Science\u201c by Joel Grus. I have it in German as an electronic book, which is unfortunately no fun to read, although the content is solid \u2013 with restrictions. The most important problem is the formatting of the Python code, which does not even work in landscape mode. The printed book uses a smaller font and lines are completely in one line. You can and should download the examples. There are also some misprints. I found one in the introduction to the Bayes formula, where the fraction is reversed. But most importantly, I have my doubts if the book is really a help for the starting data scientist. The presentation is very short. You will need a good background in math and programming, or have someone to ask if you do not understand the theory or the code. The code gets pretty advance pretty fast.<\/p>\n\n\n\n
However, this is not a review of the book. But I take an example of the book to show how Python can be used from within Euler Math Toolbox. The following code reads the bible and prints the 10 most common words.<\/p>\n\n\n\n
infile = open(\"bible.txt\",\"r\")\n\nfrom collections import Counter\n\nnum_words = 10\n\ncounter = Counter(word.lower()\n for line in infile\n for word in line.strip().split()\n if word)\n\nfor word, count in counter.most_common(num_words):\n print(str(count)+\"\\t\"+word)\n\ninfile.close()<\/code><\/pre>\n\n\n\n
I downloaded the file \u201ebible.txt\u201c from the net and put it into the folder, where the IPython notebook resides. The output is the following.<\/p>\n\n\n\n
Now, let us do this programming exercise in Euler Math Toolbox via Python.<\/p>\n\n\n\n
>>> from collections import Counter\n>function python testbible (filename) ...\n$\"\"\" 10 most common words in the bible\n$\"\"\"\n$infile = open(filename,\"r\")\n$num_words = 10\n$counter = Counter(word.lower()\n$ for line in infile\n$ for word in line.strip().split()\n$ if word)\n$res = \"\"\n$for word, count in counter.most_common(num_words):\n$ res = res+str(count)+\"\\t\"+word+\"\\n\"\n$infile.close()\n$return res\n$endfunction\n>testbible(\"bible.txt\")\n 64183 the\n 51379 and\n 34746 of\n 13643 to\n 12799 that\n 12559 in\n 10263 he\n 9840 shall\n 8987 unto\n 8835 for\n \n>>> help(testbible)\n Help on function testbible in module __main__:\n \n testbible(filename)\n 10 most common words in the bible\n<\/code><\/pre>\n\n\n\n
The import was done with a direct Python command. The function returns a string with the answer. But it could just as well use \u201eeumat.dump(\u201e\u2026.\u201c)\u201c for a direct output. Note, that you do not need to indent the Python definition. That is done automatically by EMT.<\/p>\n\n\n\n
The \u201e$\u201c-signs in the function definition allow copying and pasting the code into an EMT notebook. They are not visible in the EMT window.<\/p>\n\n\n\n
Here is another example, done in IPython.<\/p>\n\n\n\n
def sieve(N):\n v = [True for _ in range(N+1)]\n for k in range(2,N+1):\n if k*k>N:\n break\n if v[k]:\n m=2*k\n while m<=N:\n v[m]=False\n m=m+k\n return [k for k in range(N+1) if v[k] and k>1]\n\ns=sieve(100_000_000)\n\nfrom math import log\nx=[k for k in range(1000,len(s),1000)]\ny=[(s[k]\/(k*log(k))-1-(log(log(k))-1)\/log(k))*log(k)**2 for k in \n range(1000,len(s),1000)]\n\nimport matplotlib as mpl\nimport matplotlib.pyplot as plt\nfig, ax = plt.subplots(figsize=(8,8))\nax.plot(x,y)\n<\/code><\/pre>\n\n\n\n
This code computes the primes up to 100 million and plots the actual size of the k-th prime, divided by a known asymptotic for the size of the k-th prime.<\/p>\n\n\n\n<\/figure>\n\n\n\n
We only use Python to do the sieve. For the plot, we use Euler Math Toolbox.<\/p>\n\n\n\n
>function python sieve (N) ...\n$v = [True for _ in range(N+1)]\n$for k in range(2,N+1):\n$ if k*k>N:\n$ break\n$ if v[k]:\n$ m=2*k\n$ while m<=N:\n$ v[m]=False\n$ m=m+k\n$return [k for k in range(N+1) if v[k] and k>1]\n$endfunction\n>tic; s = sieve(100 000 000); toc;\n Used 29.407 seconds\n>n = length(s)\n 5761455\n>k = 1000:1000:n;\n>y = (s[k]\/(k*log(k))-1-(log(log(k))-1)\/log(k))*log(k)^2;\n>plot2d(k,y):<\/code><\/pre>\n\n\n\n<\/figure>\n\n\n\n
As you see, the sieve takes about 30 seconds on my computer (Core I7, 4.0GHz). Of course, it is much faster to do that in C.<\/p>\n\n\n\n
The procedures to communicate from EMT to Tiny C are a bit involved. But the result is 15 times faster. The basic code is a direct translation of the Python code.<\/p>\n\n\n\n
Depending on your heap space of Euler Math Toolbox, you can do more primes. If you compress the sieve using bits, you could to a lot more. But for this code and the default size, 100 million is the limit. Increasing the stack size to 2048 MB, I could do all primes up to 1 billion in 15 seconds.<\/p>\n\n\n\n<\/figure>\n\n\n\n
![\"\"](\"https:\/\/renegrothmann.de\/wp-content\/uploads\/2023\/04\/Download.png\")
<\/figure>\n\n\n\n![\"\"](\"https:\/\/renegrothmann.de\/wp-content\/uploads\/2023\/04\/Primes.png\")
<\/figure>\n\n\n\n![\"\"](\"https:\/\/renegrothmann.de\/wp-content\/uploads\/2023\/04\/Primes-1.png\")
<\/figure>\n\n\n\n