[fpc-pascal]isprime()

[Anton Tichawa]

Hi, Franco Milani!

To me it seems that your function works properly. I compared the results for 
all cardinal values, using your algorithm (cprimes_v.create below) against 
the results of my own algorithm (cprimes_m.create below), which I updated to 
32 bits tonight. All I changed in your function (isprime) is that I renamed 
it to 'is_prime' and moved the constants and variables into the function's 
context. The comparison was done file-to-file (about 257 MB).

The report file may shed some light on the timing trend behaviour of the 
algorithms - see report below. Note that due to a bug, my primes_m num_primes 
values shown are too low by 1 for all but the last two runs (bug).

As I use a different algorithm, the function may be considered to be tested 
by diversification (?).

Shouldn't the function behaviour for n = 0 be discussed? Maybe it's worth 
thinking about returning true or even throwing an exception for n = 0?

Compared to yours, my approach to primes is a 'synthetic' one - I build a 
series of primes to work with, whereas you analytically check one single 
number for primality. Did you know that using a series of primes as input to 
'turtle graphics', where prime modulos determine turtle movement, and each 
turtle move changes pixel colors, leads to very impressive images?

A note on efficiency: On a typical PC, with enough memory but performance 
needs, it might be better to use a one-time created table instead of run-time 
calculation. When the the sieve of Erasthotenes (?) is applied with a defined 
order, there is a relation of similar order from a given cardinal number to 
the number (index) within the set of "sieved" numbers. The table then holds a 
primality flag for each "sieved" number. The table is calculated once, and 
allows very fast primality checking: At runtime, the sieve is applied, and 
those numbers "falling through the sieve" are looked up in the table. The 
higher the order, the smaller the table. But that's an idea - it's not coded 
yet. The most simple example is a table that contains at index (n div 2) a 
primality flag for each odd number (sieve order 1), and holds high(cardinal) 
div 2 entries.

In the first lines of assembly, I found a possible optimization, so I dared 
skipping the rest of the assembly text. This optimization would be checking > 
66 instead of > 61, which speeds up the function for the non-prime numbers 62 
.. 66.

here's my procedure testing your algorithm:

constructor cprimes_v.create(aname: string);
var
  i: cardinal;
  s: string;  // -tbo-
  fc: coutputtextfile;  // -tbn- the common report file
  f: coutputtextfile;
  my_num_primes: cardinal;
  duration: cardinal;
begin
  log('* cprimes_v.create');
  duration := getabsticks;
  my_num_primes := 0;
  fc := coutputtextfile.create('/io/txt/primes.txt');
  fc.append;
  f := coutputtextfile.create('/io/txt/primes_v.txt');
  inherited create(aname);
  for i := 0 to max_prime_mask do begin
    if is_prime(i) then begin
      s := intstr(i);
      f.writeln(s);
      inc(my_num_primes);
      // add_string(s);
    end;
  end;
  duration := getabsticks - duration;
  fc.writeln('primes_v:');
  fc.writeln('  max_prime_mask ' + intstr(max_prime_mask, 10));
  fc.writeln('  num_primes     ' + intstr(my_num_primes, 10));
  fc.writeln('  duration (ms)  ' + intstr(duration, 10));
  fc.free;
  f.free;
  log('+ cprimes_v.create');
end;

here's my procedure testing my algorithm:

constructor cprimes_m.create(aname: string);
var fc: coutputtextfile;  // -tbn- common report file
var f: coutputtextfile;
var i: cardinal;
var j: cardinal;
var q: cardinal;
var my_num_roots: cardinal;
var my_num_primes: cardinal;
var my_primes: array [0 .. max_prime_root - 1] of cardinal;
var my_overrun: boolean;
var duration: cardinal;
begin
  log('* cprimes.create');
  inherited create(aname);
  fc := coutputtextfile.create('/io/txt/primes.txt');
  fc.append;
  f := coutputtextfile.create('/io/txt/primes_m.txt');
  duration := getabsticks;
  my_num_primes := 1;
  my_primes [0] := 2;
  my_num_roots := 1;
  f.writeln('2');  // -tbo-
  i := 3;
  my_overrun := false;
  repeat
    for j := 0 to my_num_roots - 1 do begin
      q := my_primes [j];
      // log('checking i ' + intstr(i) + ' q ' + intstr(q));
      if (i div q) * q = i then begin
        // log('  no');
        break;
      end;
      // log('  yes');
      if q * q >= i then begin
        inc(my_num_primes);
        f.writeln(intstr(i));
        if my_num_roots < max_prime_root then begin
          my_primes [my_num_roots] := i;
          inc(my_num_roots);
        end;
        break;
      end;
    end;
    if my_overrun then begin
      break;
    end;
    if i >= max_prime_mask then begin
      break;
    end;
    inc(i);
  until false;
  if my_overrun then begin
    f.writeln('overrun');
  end;
  duration := getabsticks - duration;
  fc.writeln('primes_m:');
  fc.writeln('  max_prime_mask ' + intstr(max_prime_mask, 10));
  fc.writeln('  num_primes     ' + intstr(my_num_primes, 10));
  fc.writeln('  duration (ms)  ' + intstr(duration, 10));
  fc.free;
  f.free;
  log('+ cprimes.create');
end;

here's the report (file '/io/txt/primes.txt'):

primes_m:
  max_prime_mask     999999
  num_primes          78497  // should read 78498
  duration (ms)        1347
primes_v:
  max_prime_mask     999999
  num_primes          78498
  duration (ms)        1526
primes_m:
  max_prime_mask    9999999
  num_primes         664578  // .. '9
  duration (ms)       27956
primes_v:
  max_prime_mask    9999999
  num_primes         664579
  duration (ms)        9619
primes_m:
  max_prime_mask    3999999
  num_primes         283145  // .. '6
  duration (ms)        8757
primes_v:
  max_prime_mask    3999999
  num_primes         283146
  duration (ms)        3764
primes_m:
  max_prime_mask   39999999
  num_primes        2433653  // .. '4
  duration (ms)      172097
primes_v:
  max_prime_mask   39999999
  num_primes        2433654
  duration (ms)       44459
primes_m:
  max_prime_mask   99999999
  num_primes        5761454  // .. '5
  duration (ms)      582844
primes_v:
  max_prime_mask   99999999
  num_primes        5761455
  duration (ms)      114538
primes_m:
  max_prime_mask  499999999
  num_primes       26355867  // ok
  duration (ms)     5146543
primes_v:
  max_prime_mask  499999999
  num_primes       26355867
  duration (ms)      606424
primes_m:
  max_prime_mask       9999
  num_primes           1229  // ok
  duration (ms)         357
primes_v:
  max_prime_mask       9999
  num_primes           1229
  duration (ms)         358

Another note: During inclusion of your sources, I had to convert characters 
of order 160 (0xA0) to spaces - where do those characters come from?

Hope this Helps.

Anton Tichawa.



----------

"Adas Methode war, wie sich zeigen wird, Tagträume in offenbar korrekte 
Berechnungen einzuweben."

Doris Langley Moore: Ada, Countess of Lovelace (London 1977).

----------

Anton Tichawa
Volkertstrasse 19 / 20
A-1020 Wien
mobil: +43 664 52 07 907
email: anton.tichawa at chello.at

----------