Amazon.com Widgets Fun Code of the Week #6: Primality Checking

Fun Code of the Week #6: Primality Checking

By Nick at October 27, 2012 18:02
Filed Under: Delphi, Fun Code

 

I decided to give KhanAcademy a spin, and watched this set of videos:

Level 1: Primality Test

And that let to me writing this code:

function IsPrime(const x: integer): Boolean;
var
  i: integer;
begin
  i := 2;
  repeat
    if X mod i = 0 then
    begin
      Result := False;
      Exit;
    end;
    Inc(i);
  until i > Sqrt(x);
  Result := True;
end;

I know that this can be optimized (which I'll do if I end up watching the next video :-) ), but I don't often write really geeky code like this, so I thought I'd post it so you all can write lots of comments like you always do when I post code. ;-)

UPDATE #1: There was a bug if you passed in a negative number. 

function IsPrime(const x: integer): Boolean;
var
  i: integer;
begin
  if (x <= 2) then
  begin
    Result := False;
    Exit;
  end;
  i := 3;
  repeat
    if X mod i = 0 then
    begin
      Result := False;
      Exit;
    end;
    Inc(i);
  until i > Sqrt(x);
  Result := True;
end;

UPDATE #2:  Took the Sqrt() call out of the loop.

function IsPrime(const x: integer): Boolean;
var
  i: integer;
  Wall: double;
begin
  if (x <= 2) then
  begin
    Result := False;
    Exit;
  end;
  i := 3;
  Wall := Sqrt(x);
  repeat
    if X mod i = 0 then
    begin
      Result := False;
      Exit;
    end;
    Inc(i);
  until i > Wall;
  Result := True;
end;

UPDATE #3:   The learning and optimizing continues!  I updated the code based on more suggestions, but the special case of 2 is irritating me.  I confess I thought a while (I’m not cheating by looking up implementations on Google…) and so I came up with something of a hack.  Also, just for Julian, all the x’s are lowercase now.  Smile with tongue out  Also, I learned/re-learned that 1 is not a prime either, by definition.  Interesting.  Any more optimizations out there?  Is there a “correct” way to handle 2?  The end goal here is to have a really nice, really efficient IsPrime function. 

function IsPrime(const x: integer): Boolean;
var
  i: integer;
  Wall: integer;
begin
  if x = 2 then
  begin
    Result := True;
    Exit;
  end;
  i := 3;
  Result := not ((x < i) or (x mod 2 = 0));
  if not Result then Exit;
  Wall := Trunc(Sqrt(x));
  repeat
    Result := not (X mod i = 0);
    if not Result then Exit;
    Inc(i, 2);
  until i > Wall;
  Result := True;
end;

Comments (13) -

10/27/2012 7:15:23 PM #

Warren A Postma

I would store sqrt(x) instead of repeatedly calculating it in real world code. Smile

For example, check once for X mod 2 outside the loop using ((x and 1)=0), then loop starting with I=3, adding 2 to I each time, only use odd numbers in the loop, you would roughly double the speed of your primality testing.

W

Warren A Postma Canada |

10/27/2012 9:14:50 PM #

Nick Hodges

Warren --

Good point.  Wink  I'll update the code.

Nick Hodges United States |

10/27/2012 10:01:31 PM #

Julian M Bucknall

I notice that you only took one of Warren's recommendations to heart (but in doing so, introduced another bug: 2 is prime). I would also calculate the root of x as an integer to slightly improve the speed of the comparison at the end of the repeat loop (you're comparing an integer value to a double, so one value has to be converted to the other type).

Also, one reason I now don't use Pascal (er, Delphi) much: you have x, x, X, x. Meh. ;)

Cheers, Julian

Julian M Bucknall United States |

10/28/2012 5:18:47 AM #

&#201;ric

Check out the Euler Project Smile

Éric France |

10/28/2012 6:57:01 AM #

Just Curious

I was going to describe in words, but copy and paste works fine ;)

Thanks for previous contributors

function IsPrime(const x: integer): Boolean;
var
  i: integer;
  Wall: integer;
begin
  i := 3;
  Result := not ((x < i) OR (X mod 2 = 0));
  if not Result then Exit;
  Wall := Trunc(Sqrt(x));
  repeat
    Result := not (X mod i = 0);
    if not Result then Exit;
    Inc(i,2);
  until i > Wall;
  Result := True;
end;

Just Curious United Kingdom |

10/28/2012 3:00:13 PM #

Bunny

Just to be complete ... Sieve of Eratosthenes and/or store the prime numbers found into a file.

Bunny Austria |

10/28/2012 3:55:19 PM #

David M

Also you can change
  if x = 2 then
  begin
    Result := True;
    Exit;
  end;
to
  if x = 2 then
    Exit(True);

David M New Zealand |

10/29/2012 4:18:18 AM #

Uwe Raabe

The new exit syntax can be quite a line saver here. It also spares setting and checking for the result variable.

Another optimization is to replace the x mod 2 = 0 with a not odd(x). The compiled code looks much shorter (XE2 Win32):

Project191.dpr.24: if x mod 2 = 0 then Exit(false);
0040512C 8B45FC           mov eax,[ebp-$04]
0040512F 2501000080       and eax,$80000001
00405134 7905             jns $0040513b
00405136 48               dec eax
00405137 83C8FE           or eax,-$02
0040513A 40               inc eax
0040513B 85C0             test eax,eax
0040513D 7506             jnz $00405145
0040513F C645FB00         mov byte ptr [ebp-$05],$00
00405143 EB0A             jmp $0040514f
Project191.dpr.25: if not Odd(x) then Exit(false);
00405145 F645FC01         test byte ptr [ebp-$04],$01
00405149 7504             jnz $0040514f
0040514B C645FB00         mov byte ptr [ebp-$05],$00


Putting all together and shuffling a bit, this is my small contribution (keeping the algorithm):

function IsPrime(const x: integer): Boolean;
var
  i: integer;
  Wall: integer;
begin
  if Odd(x) then begin
    if x < 3 then Exit(false);
    i := 3;
    Wall := Trunc(Sqrt(x));
    repeat
      if (x mod i = 0) then Exit(false);
      Inc(i, 2);
    until i > Wall;
    Result := True;
  end
  else begin
    Result := (x = 2);
  end;
end;

Uwe Raabe Germany |

10/29/2012 10:46:40 AM #

Nick Hodges

Uwe -- nice!

Nick Hodges United States |

11/3/2012 7:52:40 PM #

Thomas Vedel

In fact the

repeat
until i > Wall;

should be replaced by a

while i <= Wall do
begin
end;

which also prevents that 3 is mistakenly rejected as a prime.

To furthermore minimize the calculations, it is actually only necessary to check division by prime numbers instead of all odd numbers. However I guess that the housekeeping needed to only divide by primes may actually take more computing power than dividing by all odd numbers.

Thomas Vedel Denmark |

10/29/2012 12:13:42 PM #

Eric

> The end goal here is to have a really nice, really efficient IsPrime function.

Check the one in DWScript's dwsMathFunctions, it uses a slightly improved algorithm, and when benchmarking on

for i:=1 to 1000000 do
   IsPrime(i);

using your Update #3 (which btw fails for 3) I see 0.49 sec in Delphi using Integer
The same loop runs in 0.39 sec in DWScript (using Int64 in a 32bit build).

Eric France |

10/30/2012 5:56:23 PM #

Nick Hodges

Eric --

  Will do.  I should have thought to look there.  You don't let inefficient code into your libraries.  Wink

Nick Hodges United States |

11/10/2012 3:57:26 PM #

Ulf Ravn

The DWScript code skips every third odd number, avoiding all odd numbers divisible by 3. This males the routine about a third faster. Applying this idea to Uwe's code gets:
function TForm1.IsPrime(const x: int64): Boolean;
var
  i: Int64;
  Wall: Int64;
begin
  if Odd(x) then begin
    if x < 3 then Exit(false);
    if ((x mod 3) = 0) then Exit(false);
    i := 5;
    Wall := Trunc(Sqrt(x));
    repeat
      if (x mod i = 0) then Exit(false);
      Inc(i, 2);
      if (x mod i = 0) then Exit(false);
      Inc(i, 4);
    until i > Wall;
    Result := True;
  end
  else begin
    Result := (x = 2);
  end;
end;

Ulf Ravn Denmark |

Comments are closed

My Book

A Pithy Quote for You

"Luck is what happens when preparation meets opportunity."    –  Seneca

Amazon Gift Cards

General Disclaimer

The views I express here are entirely my own and not necessarily those of any other rational person or organization.  However, I strongly recommend that you agree with pretty much everything I say because, well, I'm right.  Most of the time. Except when I'm not, in which case, you shouldn't agree with me.

Month List