Re: [OctDev] list of primes using D. Bernsteins fast primegen library

[David Bateman]

Torsten,

I'm bring this discussion back to the octave maintainers list. So I've
left most of the thread in place so others can follow our discussion.

Dr.-Ing. Torsten Finke wrote:

>Hello David, 
>
>  
>
>>>excuse me, but primes() seems not to be an octave core function. There is 
>>>only a
>>>list_primes(), that is a quite slow one and it differs from primes(), since 
>>>it
>>>lists N primes not primes up to N (as primes() does). As I understand, 
>>>primes()
>>>is an octave-forge function.
>>>
>>> 
>>>
>>>      
>>>
>>I ported a large number of octave-forge functions to octave core
>>recently primes.m included. Please find the code in octave 2.9.6, or the
>>attached. Your code would be replacing this function.
>>    
>>
>
>sorry, I have only looked at 2.9.5. 
>
>  
>
>>>To be compatible with matlab it would be useful to have primes() as core
>>>function. Matlab too provides primes() as an m-file function - quite slow and
>>>memory consuming. For higher bounds it fails due to memory overflow. Here
>>>primegen() works much further. 
>>> 
>>>
>>>      
>>>
>>Please check your code against this version of primes.m
>>    
>>
>
>I've just done some tests on my machine (P4, 1GB Ram, 1GB swap, Linux 2.6.16,
>Octave 2.1.72). Here are the results for cputime()/sec: 
>
>    k     |  primes(10^k)   | primegen(10^k) |   ratio
>----------+-----------------+----------------+------------
>   2     |     0.0160010   |   0.028002     |    0.57142
>   3     |     0.0080000   |   0.024002     |    0.33331
>   4     |     0.0040010   |   0.024001     |    0.16670
>   5     |     0.0360020   |   0.024002     |    1.49996
>   6     |     0.2560150   |   0.036002     |    7.11113
>   7     |     2.6241640   |   0.180011     |   14.57780
>   8     |    26.6896680   |   1.616101     |   16.51485
>   9     |       (M)       |  16.073004     |     ---
>  10      |       (M)      |  43.598725     |     ---
>
>(M) memory exhausted. 
>For k=10 my machine was extremely using swap. 
>
>test routine: 
>
>     d = []; 
>     for k = 2:8, 
>         t = cputime(); 
>         primes(10^k); % resp. primegen(10^k)
>         d = [d;  cputime() - t]; 
>     end
>
>
>also I have checked correctness of the result: 
>
>  all(primes(1e8) == primegen(1e8))
>
>  ans = 1
>  
>
The reason that the existing primes code is faster upto 10^4 is that
there is the special case code

    a = [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, 101, 103, 107, ...
         109, 113, 127, 131, 137, 139, 149, 151, 157, 163, 167, ...
         173, 179, 181, 191, 193, 197, 199, 211, 223, 227, 229, ...
         233, 239, 241, 251, 257, 263, 269, 271, 277, 281, 283, ...
         293, 307, 311, 313, 317, 331, 337, 347, 349];
    x = a(a <= p);

for p <= 352. The same might easily be done for your code to obtain the
same speed for smaller p as well.

>
>  
>
>>>May I suggest to add primegen() as an additional option, not as a replacement
>>>for primes()? 
>>> 
>>>
>>>      
>>>
>>Its the dependency on a compile code that makes including your code
>>questionable. The policy has been that we should try and keep large
>>parts of octave in its own scripting language even at the cost of slight
>>slow-downs to reduce the maintaince cost of the code. So, to include
>>such compiled code there needs to be a demonstration of large speed-ups
>>in a function that will be used by a large group of people... This is
>>not something I'll answer myself, but will almost always be raised when
>>someone proposes replacing a script file by compiled code..
>>    
>>
>
>May I emphasize, that it is not actually my code but it is Dan Bernsteins's
>(the person who wrote qmail, but I'm sure you know him). I fully understand
>your concerns about C code - especially since it is quite complex and only few
>tests have been made.  So for the time being it would be absolutely ok to me,
>seeing that powerful function on octave-forge rather than among the octave
>core functions, but I think the code is worth to be introduced to octave.
>  
>
Then there is another policy that has been discussed is that we should
try hard not to have an overlap between the code in octave and
octave-forge to avoid issues of knowing which code is being used and
what bug fixes have been applied where. In the octave-forge CVS the
overlap with octave 2.9.6 has been completely eliminated. So I'd
hesitate to include your code in octave-forge if it isn't accepted in
octave core.

With a change to use a lookup table for small values of p, I think
you've demonstrated the case for the speed-up of your code. The next
question is how many people will really use this function and whether
the additional support costs of C code is offset by the number of
potential users who are impacted by this improvement. The final say is
of course up to John.

Regards
David


-- 
David Bateman                                address@hidden
Motorola Labs - Paris                        +33 1 69 35 48 04 (Ph) 
Parc Les Algorithmes, Commune de St Aubin    +33 6 72 01 06 33 (Mob) 
91193 Gif-Sur-Yvette FRANCE                  +33 1 69 35 77 01 (Fax) 

The information contained in this communication has been classified as: 

[x] General Business Information 
[ ] Motorola Internal Use Only 
[ ] Motorola Confidential Proprietary