tag:blogger.com,1999:blog-6272748843257835470Sat, 29 Feb 2020 07:40:25 +0000pythonfsharpgoogleaimathprogrammingjavascriptlispA Mingled Cuphttp://mingledcup.dkjones.org/search/label/mathnoreply@blogger.com (Dave)Blogger5125tag:blogger.com,1999:blog-6272748843257835470.post-8209006214425050273Sun, 22 Aug 2010 01:14:00 +00002010-08-21T22:21:21.406-04:00fsharpmathGCD in F#<p>
Now we will write a function to calculate the greatest common divisor
of two numbers. The <a href="http://en.wikipedia.org/wiki/Euclidean_algorithm">Euclidean Algorithm</a> is an efficient, and ancient, method of finding the GCD.
<p>
Euclid's algorithm takes advantage of some properties of the gcd function and the modulus function. Here is a loose explanation of why the algorithm works
<pre>
gcd(a,b) = gcd(b, a - q*b) for any integer q
a % b = a - q*b from some integer q,
so gcd(a,b) = gcd(b, a % b)
</pre>
<p>
Then, we just stop the algorithm when b is zero. And now for the F# code.
<pre>
let rec gcd a b =
if b = 0 then a
else
gcd b (a % b)
</pre>
<p>
Aside from the syntax, and checking for the end condition, this looks like the math.
<p>
And now, a few tests.
<pre class="code">
let f a b =
printfn "gcd(%d,%d) = %d" a b (gcd a b)
f 10 9
f (13*61*2*2*2*3) (13*2*3*5*7)
f -10 12345
</pre>
<p>
With the output,
<pre>
gcd(10,9) = 1
gcd(19032,2730) = 78
gcd(-10,12345) = 5
</pre>
<a href="http://ideone.com/hHN7j">ideone.com/hHN7j</a>http://mingledcup.dkjones.org/2010/08/gcd-in-f.htmlnoreply@blogger.com (Dave)0tag:blogger.com,1999:blog-6272748843257835470.post-210046780364552898Mon, 16 Aug 2010 19:46:00 +00002010-08-16T15:52:57.861-04:00mathpythonUsing the sieve to factor numbersPreviously, we saw an overly simple implementation of the Sieve of Eratosthenes. Without any optimizations, that code is pretty much only good for puzzles and learning. Even for puzzles, it may be too slow.
Some of the puzzle problems out there require code to find the prime factorization of numbers. For small numbers, trial division might be ok. Often, the problems require the factorization for every number less than one million. Then... it gets slow.
Using the same idea as the Sieve, we can prepare an array that contains a prime factor for every number, or one if the number is prime. Then, to get the full factorization of the number, we can walk the array.
Here is the code to create the array:
<pre class="code">
def divisors(N):
"""return the biggest prime divisor for numbers less than N"""
result = [1 for x in xrange(N)]
for d in xrange(2,N):
if result[d] > 1: continue
for j in xrange(d, N, d):
result[j] = d
return result
</pre>
Then, to find the factorization of a number, we can look it up in the array to find a factor, then look up (n/p) in the array to find the next factor, and keep going until we are down to a prime number.
<pre class="code">
divs = divisors(2000)
def factor(n, divisors=divs):
result = []
while n > 1:
p = divisors[n]
if p == 1:
result.append(n)
return result
while n % p == 0:
result.append(p)
n /= p
return result
print 1000, factor(1000)
print 24, factor(24)
print 17*31, factor(17*31)
</pre>
Often, the factorization is used for things like computer Euler's totient function, and then we don't care to return the factors, just the result of using them.
I used ideone.com to work on this code:
<a href="http://ideone.com/RLRZt"> Ideone.com | RLRZt </a>http://mingledcup.dkjones.org/2010/08/using-sieve-to-factor-numbers.htmlnoreply@blogger.com (Dave)0tag:blogger.com,1999:blog-6272748843257835470.post-7744328213141706733Sun, 15 Aug 2010 00:20:00 +00002010-08-14T20:33:53.978-04:00mathpythonSieve of EratosthenesHere is a simple implementation of the <a href="http://en.wikipedia.org/wiki/Sieve_of_Eratosthenes">Sieve of Eratosthenes</a> in Python.
<pre class="code">
from math import sqrt
def sieve(N):
"""Sieve of Eratosthenes. Returns primes less than N
http://en.wikipedia.org/wiki/Sieve_of_Eratosthenes
"""
# Set up the array. 2 is the first prime, not 1.
ar = range(N)
ar[1] = 0
limit = 1 + int(sqrt( N))
for d in xrange(2,limit):
if ar[d]:
# d is prime! mark off the multiples of d.
# note that we can start with d*d instead of 2*d.
for j in xrange(d*d,N,d):
ar[j] = 0
return [x for x in ar if x]
</pre>
There are several optimizations that people tend to do when they write this algorithm. But, this gets us most of the way. I keep rewriting this function for <a href="http://projecteuler.net/">Project Euler</a> problems, but I should really keep a library of the functions I use, or better yet, use <a href="http://sagemath.com/">Sage</a>.http://mingledcup.dkjones.org/2010/08/sieve-of-eratosthenes.htmlnoreply@blogger.com (Dave)0tag:blogger.com,1999:blog-6272748843257835470.post-6786774449957375128Sat, 14 Aug 2010 01:23:00 +00002010-08-13T21:24:32.944-04:00fsharpmathSimple prime test in F#Here is the same algorithm as the last <a href="2010/08/simple-prime-test-in-python.html"> post </a>, but in F#. There is no way to break out of a for loop in F#, besides throwing an exception, so this code uses recursion instead.
<pre class="code">
let is_prime n =
if n < 2 then false
else
let limit = float n |> sqrt |> fun x -> x + 1. |> int
let rec loop d =
if d = limit then true
elif n % d = 0 then false
else loop (d + 1)
loop 2
> [1..100] |> List.filter is_prime;;
val it : int list =
[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]
</pre>http://mingledcup.dkjones.org/2010/08/simple-prime-test-in-f.htmlnoreply@blogger.com (Dave)0tag:blogger.com,1999:blog-6272748843257835470.post-2287714993431006418Fri, 13 Aug 2010 02:51:00 +00002010-08-16T06:48:57.558-04:00mathpythonSimple prime test in PythonA simple way to test if a number is prime is to check if the number is divisible by any other number besides one and itself. This can take a lot of time, but for small numbers, it might be ok.
<pre class="code">
def is_prime(n):
# assumes n > 0
if n < 2: return True
for d in range(2,n):
if n % d == 0: return False
return True
</pre>
This simple function can be sped up by not checking as many numbers. We don't need to test all divisors up to n-1. The largest possible divisor is n/2, so we could stop there. But, even better, we can show that if there are any divisors, there must be divisor less than the sqrt(n). If p*q = n, then either p or q must be less than or equal to sqrt(n). Otherwise, p*q would be bigger than n. So, we can just loop up to sqrt(n) because we only need to find a divisor of n, not every divisor.
Here is the updated code:
<pre class="code">
import math
def is_prime2(n):
if n < 2: return False
for d in xrange(2, int(math.sqrt(n)) + 1):
if n % d == 0: return False
return True
>>> [x for x in xrange(40) if is_prime2(x)]
[2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37]
</pre>
Finally, we could use the same code to find a divisor of a number. In fact, it finds the smallest divisor. Some of the more advanced primality tests do not determine any of the factors of the number.
<pre>
def is_prime3(n):
"""returns the smallest divisor of n.
returns 1 if the number is prime, or less than 2"""
if n < 2: return 1
for d in xrange(2, int(math.sqrt(n)) + 1):
if n % d == 0: return d
return 1
</pre>http://mingledcup.dkjones.org/2010/08/simple-prime-test-in-python.htmlnoreply@blogger.com (Dave)2