I made a function of 116 lines in assembler:

isprime(n: cardinal): boolean;

which returns true if a cardinal is prime, false otherwise. The routine is fast and uses only an array of 18 cardinals and six cardinals as local variables. Obviously it may be used by an external unit, anyhow let me know if it might not be bad the idea to implement the function on the compiler.

Franco Milani

Hello, List!

I made a function of 116 lines in assembler: isprime(n: cardinal): boolean;

Sounds good.

The RTL?

Does the function return false for 0 and 1?

How easy isit to upgrade for > 64 bit lengths?

MfG,

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:

----------

Hi Anton :)

Does the function return false for 0 and 1?

Yes

How easy isit to upgrade for > 64 bit lengths?

Indeed, isprime() is part of a routine working with the 'nat' big numbers of the unit bigint10. Its purpose is to perform a primality check by a strong test to the base 2 and a Lucas test. But up to 32 bits it is enough a Miller-Rabin test to bases 2, 7 and 61, which covers up to 4,759,123,141. To do this, isprime() must evaluate some powers and squares modulo a cardinal involving 64 bit operands. I write here the function, so it may be used and tested. If there are problems let me know.

Franco Milani
-------------------------------------------------------------------------

CONST
arlc_ispv: array [1..18] of cardinal =
(61,59,53,47,43,41,37,31,29,23,19,17,13,11,7,5,3,2);
VAR
vloc_ispa,vloc_ispb,vloc_ispc,vloc_ispd,vloc_ispe,vloc_ispf: cardinal;

FUNCTION isprime (n: cardinal): boolean; assembler; asm
movl 8(%ebp),%ebx
leal arlc_ispv,%edi
movl $17,%ecx
cmpl $61,%ebx
ja .L5
.L6: cmpl (%edi,%ecx,4),%ebx
je .L7
decl %ecx
jns .L6
jmp .L8
.L5: movl %ebx,%eax
xorl %edx,%edx
movl (%edi,%ecx,4),%esi
divl %esi
cmpl $0,%edx
je .L8
decl %ecx
jns .L5
decl %ebx
movl %ebx,vloc_ispa
xorl %ecx,%ecx
testl $1,%ebx
jnz .L9
.L10: shrl $1,%ebx
incl %ecx
testl $1,%ebx
jz .L10
.L9: movl %ebx,vloc_ispb
cmpl $0,%ecx
je .L11
decl %ecx
.L11: movl %ecx,vloc_ispc
movl $2,vloc_ispd
movl $2,%eax
jmp .L12
.L13: movl $7,%eax
jmp .L12
.L14: movl $61,%eax
.L12: movl vloc_ispb,%ebx
movl 8(%ebp),%ecx
cmpl $1,%ebx
jne .L15
xorl %edx,%edx
divl %ecx
movl %edx,%eax
jmp .L16
.L15: movl %eax,%edi
bsrl %ebx,%edx
movl %edx,vloc_ispe
btl $0,%ebx
jnc .L17
movl %eax,%esi
jmp .L18
.L17: movl $1,%esi
.L18: movl $1,vloc_ispf
.L19: movl %edi,%eax
mull %edi
cmpl %ecx,%edx
jb .L20
movl %eax,%ebx
movl %edx,%eax
xorl %edx,%edx
divl %ecx
movl %ebx,%eax
.L20: divl %ecx
movl %edx,%edi
movl vloc_ispb,%eax
movl vloc_ispf,%ebx
btl %ebx,%eax
jnc .L21
movl %esi,%eax
mull %edi
cmpl %ecx,%edx
jb .L22
movl %eax,%ebx
movl %edx,%eax
xorl %edx,%edx
divl %ecx
movl %ebx,%eax
.L22: divl %ecx
movl %edx,%esi
.L21: incl vloc_ispf
decl vloc_ispe
jnz .L19
movl %esi,%eax
.L16: cmpl $1,%eax
je .L23
cmpl vloc_ispa,%eax
je .L23
cmpl $0,vloc_ispc
je .L8
movl vloc_ispc,%esi
.L24: mull %eax
cmpl %ecx,%edx
jb .L25
movl %eax,%ebx
movl %edx,%eax
xorl %edx,%edx
divl %ecx
movl %ebx,%eax
.L25: divl %ecx
movl %edx,%eax
cmpl vloc_ispa,%eax
je .L23
cmpl $1,%eax
je .L8
.L26: decl %esi
jnz .L24
jmp .L8
.L23: decl vloc_ispd
jz .L14
jns .L13
.L7: movl $1,%eax
jmp .L27
.L8: xorl %eax,%eax
.L27: end;

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:

----------

whereas you analytically check one single number for primality.

Just out of curiosity, is Franco's algorithm based upon the recent discovery by mathematicians in India?

Mark

Hi Anton

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.

I'm happy to read that the function works.

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?

Also 1 may be an exception. 0 and 1 have no proper factors so maybe it is more simple to return false for both.

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?

I don't know this, it is interesting and I'll look about it.

I agree, the sieve of Erathostenes is very fast so for some applications may be suitable. But it requires a lot of memory and large arrays.

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.

Good suggestion. 61 is the minimum, because the function must check the primes up to the bigger base. We can increase the array if there are no problems of memory or implementation.

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

Sorry, some tabs of the code end up in the message :). I'll pay more attention.
Thanks for your test.

Franco

Hi Mark

Just out of curiosity, is Franco's algorithm based upon the recent discovery by mathematicians in India?

No, the method of Agrawal, Kayal and Saxena is still in a theoretical state and we have to wait for an efficient algorithm to apply it. I used an old Miller-Rabin method. By testing the strong probable primality to some bases, in particular ranges we can have the certainty that the number is prime.
See for example http://www.utm.edu/research/primes/prove/prove2_3.html

Franco

I've revised isprime(), now it is a bit faster. The new version embedded in a small unit can be found at:
http://www.polarhome.com/~franco/isprime.pas

Franco