Wednesday
Jun202012

Testing for Prime Numbers Using Powershell

 So what is a prime number anyway. Lets see what Wikipedia has to say:

A prime number (or a prime) is a natural number greater than 1 that has no positive divisors other than 1 and itself. A natural number greater than 1 that is not a prime number is called a composite number. For example, 5 is prime, as only 1 and 5 divide it, whereas 6 is composite, since it has the divisors 2 and 3 in addition to 1 and 6. The fundamental theorem of arithmetic establishes the central role of primes in number theory: any integer greater than 1 can be expressed as a product of primes that is unique up to ordering. This theorem requires excluding 1 as a prime.
Since there are infintely number of primes, it's always fun to see how fast and how many we can find. The largest know prime number has nearly 13 million decimal digits!

 

Now thats out of the way lets see how we can test to see whether a number is prime or not using Powershell. Today we will look at a few different methods, including my favorite, Sieve of Eratosthenes. The first method we will look at is using a brute force method which is by far the slowest. Now lets look at some code:
function Test-isPrime            
{            
    <#
    .Synopsis
       Tests if a number is a prime using the brute force method.
    .EXAMPLE
       Test-isPrime (1..1000)
    .EXAMPLE
       1..1000 | Test-isPrime
    #>            
    [CmdletBinding()]          
            
    Param            
    (            
        [Parameter(Mandatory=$true,            
                   ValueFromPipeLine=$true,            
                   ValueFromPipelineByPropertyName=$true,            
                   Position=0)]            
        [int[]]$number            
    )            
            
    Begin            
    {            
        $properties = @{'Number'=$null;'Prime'=$null}            
    }            
            
    Process            
    {            
        foreach ($n in $number)            
        {            
            $i = 2            
            $properties.Number = $n            
            $properties.Prime = $true            
               
            if ($n -lt 2)            
            {            
                 $properties.Prime = $false            
            }            
               
            if ($n -eq 2)            
            {            
                $properties.Prime = $true            
            }            
            
            while ($i -lt $n)            
            {            
                if ($n % $i -eq 0)            
                {            
                    $properties.Prime = $false            
                }            
            
                $i++            
            
            }
            $object = New-Object PSObject –Prop $properties            
            Write-Output $object            
        }            
    }            
    End            
    {            
    }            
}
Now lets measure how long it takes to test 1 through 10,000.
PS C:\Windows\system32> Measure-Command {1..10000 | Test-isPrime}            
            
Days              	: 0            
Hours             	: 0            
Minutes           	: 3            
Seconds           	: 20            
Milliseconds      	: 93            
Ticks             	: 2000930921            
TotalDays         	: 0.00231589226967593            
TotalHours        	: 0.0555814144722222            
TotalMinutes      	: 3.33488486833333            
TotalSeconds      	: 200.0930921            
TotalMilliseconds 	: 200093.0921
This certainly isn't the most efficient way to check whether a number is a prime. So lets spice it up a bit and test if a number is a prime by checking only the numbers between 1 and n/2-1. If we find such a divisor that will be enough to say that n isn't prime:
function Test-isPrime2            
{            
    <#
    .Synopsis
       Tests if a number is a prime by checking only the numbers between 1 and n/2-1. If 
       we find such a divisor that will be enough to say that n isn’t prime.
    .EXAMPLE
       Test-isPrime2 (1..1000)
    .EXAMPLE
       1..1000 | Test-isPrime2
    #>            
            
    [CmdletBinding()]           
            
    Param            
    (            
        [Parameter(Mandatory=$true,            
                   ValueFromPipeLine=$true,            
                   ValueFromPipelineByPropertyName=$true,            
                   Position=0)]            
        [int[]]$number            
    )            
            
    Begin            
    {            
        $properties = @{'Number'=$null;'Prime'=$null}            
    }            
            
    Process            
    {            
        foreach ($n in $number)            
        {            
            $i = 2            
            $properties.Number = $n            
            $properties.Prime = $true            
            
            if ($n -lt 2)            
            {            
                $properties.Prime = $false            
            }            
                
            if ($n -eq 2)            
            {            
                $properties.Prime = $true            
            }            
            
            while ($i -le $n/2)            
            {            
                if ($n % $i -eq 0)            
                {            
                    $properties.Prime = $false            
                }            
            
                $i++            
            
            }
            $object = New-Object PSObject –Prop $properties            
            Write-Output $object            
        }            
    }            
    End            
    {            
    }            
}
Now lets measure how long it takes to test 1 through 10,000.
PS C:\Windows\system32> Measure-Command {1..10000 | Test-isPrime2}            
            
Days              	: 0            
Hours             	: 0            
Minutes           	: 2            
Seconds           	: 30            
Milliseconds     	: 357            
Ticks             	: 1503579593            
TotalDays         	: 0.00174025415856481            
TotalHours        	: 0.0417660998055556            
TotalMinutes      	: 2.50596598833333            
TotalSeconds      	: 150.3579593            
TotalMilliseconds 	: 150357.9593
Slightly better, but still not very effective for larger numbers. Lets try to improve our code by checking our number against [2, sqrt(n)]. This is correct because if n isn't prime it can be represented as p*q = n. Of course if p > sqrt(n), which we assume can't be true, that will mean that q < sqrt(n).
function Test-isPrime3            
{            
    <#
    .Synopsis
       Tests if a number is a prime by checking against [2, sqrt(n)]. This is correct, 
       because if n isn’t prime it can be represented as p*q = n. Of course if 
       p > sqrt(n), which we assume can’t be true, that will mean that q < sqrt(n).
    .EXAMPLE
       Test-isPrime3 (1..1000)
    .EXAMPLE
       1..1000 | Test-isPrime3
    #>            
            
    [CmdletBinding()]           
            
    Param            
    (            
        [Parameter(Mandatory=$true,            
                   ValueFromPipeLine=$true,            
                   ValueFromPipelineByPropertyName=$true,            
                   Position=0)]            
        [int[]]$number            
    )            
            
    Begin            
    {            
        $properties = @{'Number'=$null;'Prime'=$null}            
    }            
            
    Process            
    {            
        foreach ($n in $number)            
        {            
            $i = 2            
            $sqrtN = [math]::sqrt($n)            
            $properties.Number = $n            
            $properties.Prime = $true            
            
            if ($n -lt 2)            
            {            
                $properties.Prime = $false            
            }            
               
            if ($n -eq 2)            
            {            
                $properties.Prime = $true            
            }            
            
            while ($i -le $sqrtN)            
            {            
                if ($n % $i -eq 0)            
                {            
                    $properties.Prime = $false            
                }            
            
                $i++            
            
            }
            $object = New-Object PSObject –Prop $properties            
            Write-Output $object            
        }            
    }            
    End            
    {            
    }            
}
Now lets measure how long it takes to test 1 through 10,000.
PS C:\Windows\system32> Measure-Command {1..10000 | Test-isPrime3}            
            
Days              	: 0            
Hours             	: 0            
Minutes           	: 0            
Seconds           	: 6            
Milliseconds      	: 32            
Ticks             	: 60323509            
TotalDays         	: 6.98188761574074E-05            
TotalHours        	: 0.00167565302777778            
TotalMinutes      	: 0.100539181666667            
TotalSeconds      	: 6.0323509            
TotalMilliseconds 	: 6032.3509
As you can see this is a dramatic increase in performance, but we can still do better. Time to implement the Sieve of Eratosthenes.
function Test-isPrime4            
{            
    <#
    .Synopsis
       Tests if a number is a prime by using the Eratosthenes sieve.
    .EXAMPLE
       Test-isPrime4 -max 10
    #>            
            
    Param            
    (            
        [Parameter(Mandatory=$true,            
                   ValueFromPipeLine=$true,            
                   ValueFromPipelineByPropertyName=$true,            
                   Position=0)]            
        [int]$maxnumber            
    )            
            
    Begin            
    {            
        $sieve = New-Object 'System.Collections.Generic.SortedDictionary[int, bool]'            
        for ($i=2;$i -lt $maxnumber+1; $i++)            
        {            
            <#Intially set all numbers in the list to true. We will use
	    Eratosthenes sieve to find which integers are not prime.#>            
            $sieve[$i] = $true            
        }            
    }            
             
 Process            
 {             
    foreach ($i in @($sieve.GetEnumerator()))            
    {            
         <#Starting from p, count up in increments of p and mark each of these 
	 numbers greater than p itself in the list. These numbers will be 2p, 3p, 
	 4p, etc.; note that some of them may have already been marked. 
	 Initially, let p equal 2, the first prime number #>            
        $p = 2            
               
        while (($i.Key * $p) -le $maxnumber)            
        {            
            if ($sieve.ContainsKey($i.Key * $p))            
            {            
                $sieve[($i.Key * $p)] = $false            
            }            
                
            $p++            
         }            
    }            
              
    Write-Output $sieve            
 }            
             
 End            
 {            
 }            
}
Now lets measure how long it takes to test 1 through 10,000.
PS C:\Windows\system32> Measure-Command {Test-isPrime4 -maxnumber 10000}            
            
Days              	: 0            
Hours             	: 0            
Minutes           	: 0            
Seconds           	: 2            
Milliseconds      	: 979            
Ticks             	: 29797250            
TotalDays         	: 3.44875578703704E-05            
TotalHours        	: 0.000827701388888889            
TotalMinutes      	: 0.0496620833333333            
TotalSeconds      	: 2.979725            
TotalMilliseconds 	: 2979.725
Implementing the right algortihm, with relatively minor changes to our function, can have substantial benefits. Even if it is an ancient algorithm!

 

Thursday
Oct062011

Back to Basics: Newton-Raphson Method in Powershell

In the last post we looked at the Bisection Method to solve a simple problem, finding the square root of a real number. This week we are going to use the same exact problem, but use a better algorthim.

Newton's method (also known as the Newton–Raphson method), named after Isaac Newton and Joseph Raphson, is a method for finding successively better approximations to the roots (or zeroes) of a real-valued function. The algorithm is first in the class of Householder's methods, succeeded by Halley's method.

The Newton-Raphson method in one variable:

Given a function ƒ(x) and its derivative ƒ '(x), we begin with a first guess x0 for a root of the function. Provided the function is reasonably well-behaved a better approximation x1 is

x_{1} = x_0 - \frac{f(x_0)}{f'(x_0)}.\,\!

Geometrically, x1 is the intersection with the x-axis of a line tangent to f at f(x0).

The process is repeated until a sufficiently accurate value is reached:

x_{n+1} = x_n - \frac{f(x_n)}{f'(x_n)}.\,\!

So basically we start with an initial guess, find the tangent where it crosses the x axis as the next guess, and iterate.

So lets define a function to find a square root using this method:

function squareRootNR {            
 Param(            
  [parameter(Position=0,             
  Mandatory=$true,             
  HelpMessage="Enter a non-negative number." )]             
  [ValidateScript({$_ -gt 0})]             
  [float]$x,            
  [parameter(Position=1,             
  Mandatory=$true,             
  HelpMessage="Enter upper bound on the relative error." )]             
  [ValidateScript({$_ -gt 0})]            
  $epsilon            
 )            
             
 $guess = $x / 2.0            
 $diff = [Math]::Pow($guess,2) - $x            
 $ctr = 1            
             
 while ([Math]::Abs($diff) -gt $epsilon -and $ctr -lt 100) {            
  $guess = $guess - $diff / (2.0 * $guess)            
  $diff = [Math]::Pow($guess,2) - $x            
  $ctr += 1            
 }            
             
 if ($ctr -ge 100) {            
  Write-Warning "Iteration count exceeded."            
 }            
             
 Write-Host "NRMethod. Number of Iterations was:", $ctr, " and Estimate is: ", $guess            
 return $guess            
}
PS H:\Development\Powershell> squareRootNr -x 9 -epsilon 0.000001            
NRMethod. Number of Iterations was: 5  and Estimate is:  3.00000000003932            
3.00000000003932     
       

Lets compare this against our last attempt using the Bisection method:

PS H:\Development\Powershell> squareRootBi -x 9 -epsilon 0.000001 BiMethod. Number of Iterations was: 25 and Estimate is: 3.00000008940697 3.00000008940697

As you can see we need much fewer iterations to solve the same problem.

Tuesday
Oct042011

Back to Basics: Successive Approximation in Powershell

In the last few posts we talked about exhaustive enumeration where we try all "possible" values until we find the solution. But what if there is no exact answer? What if we cannot enumerate all of the guesses? What we need is an ability to improve our guesses.

Enter Successive Approximation:

  1. Start with a guess
  2. Iterate
  3. Improve our guess

So how do we do that? We are going to start with something called the Bisection Method. Lets see what Wikipedia has to say:

The bisection method in mathematics is a root-finding method which repeatedly bisects an interval and then selects a subinterval in which a root must lie for further processing. It is a very simple and robust method, but it is also relatively slow. Because of this, it is often used to obtain a rough approximation to a solution which is then used as a starting point for more rapidly converging methods.[1] The method is also called the binary search method[2] and is similar to the computer science Binary Search, where the range of possible solutions is halved each iteration.

The method is applicable when we wish to solve the equation f(x) = 0 for the real variable x, where f is a continuous function defined on an interval [ab] and f(a) and f(b) have opposite signs. In this case aand b are said to bracket a root since, by the intermediate value theorem, the f must have at least one root in the interval (ab).

At each step the method divides the interval in two by computing the midpoint c = (a+b) / 2 of the interval and the value of the function f(c) at that point. Unless c is itself a root (which is very unlikely, but possible) there are now two possibilities: either f(a) and f(c) have opposite signs and bracket a root, or f(c) and f(b) have opposite signs and bracket a root. The method selects the subinterval that is a bracket as a new interval to be used in the next step. In this way the interval that contains a zero of f is reduced in width by 50% at each step. The process is continued until the interval is sufficiently small.

As always lets use a problem we can all understand, finding the square root of a real number. Lets write a Bisection search function and test it.

            
function squareRootBi {            
 Param(            
  [parameter(Position=0,             
  Mandatory=$true,             
  HelpMessage="Enter a non-negative number." )]             
  [ValidateScript({$_ -gt 0})]             
  $x,            
  [parameter(Position=1,             
  Mandatory=$true,             
  HelpMessage="Enter upper bound on the relative error." )]             
  [ValidateScript({$_ -gt 0})]            
  $epsilon            
 )            
             
 $low = 0            
 $high = $x            
 $guess = ($low + $high) / 2.0            
 $ctr = 1            
             
 while ([Math]::Abs([Math]::Pow($guess,2) - $x) -gt $epsilon -and $ctr -le 100) {            
  if ([Math]::Pow($guess,2) -lt $x) {            
   $low = $guess # Set the lower bound bisection value            
  } else {            
   $high = $guess # Set the higher bound bisection value            
  }            
              
  $guess = ($low + $high) / 2.0            
  $ctr += 1            
 }            
             
 if ($ctr -ge 100) {            
  Write-Warning "Iteration count exceeded."            
 }            
             
 Write-Host "BiMethod. Number of Iterations was:", $ctr, " and Estimate is: ", $guess            
 return $guess            
}
PS H:\Development\Powershell> squareRootBi -x 9 -epsilon 0.0001            
BiMethod. Number of Iterations was: 18  and Estimate is:  2.9999885559082            
2.9999885559082

Stay tuned for the next post where we dive into a much better way to solve this, the Newton-Raphson Method.