Recent

Author Topic: Pascal programs in the Rosetta Code collection  (Read 5822 times)

julkas

  • Guest
Re: Pascal programs in the Rosetta Code collection
« Reply #30 on: August 12, 2019, 01:31:17 pm »
Do you think this is good code? <bemused> https://www.nayuki.io/res/the-versatile-sieve-of-eratosthenes/eratosthenes-sieves.c
If you want simple classic sieve of Eratosthenes - yes.
If you want advanced - https://github.com/kimwalisch/primesieve

avk

  • Hero Member
  • *****
  • Posts: 752
Re: Pascal programs in the Rosetta Code collection
« Reply #31 on: August 12, 2019, 02:17:38 pm »
Hmm, I've eliminated the extra pass through the array (although I doubt it will be faster). I will try to consider advanced variants later.
Code: Pascal  [Select][+][-]
  1. program prime_sieve;
  2. {$mode objfpc}{$H+}
  3. uses
  4.   SysUtils;
  5. const
  6.   MAX_PRIME_LIMIT = {$ifdef cpu32}1024*1024*1024{$else}High(Integer){$endif};
  7. type
  8.   TIntArray = array of Integer;
  9.  
  10. function EstimatePrimeCount(aLimit: Integer): Integer; inline;
  11. begin
  12.   if aLimit <= 200 then
  13.     Result := Trunc((1.6 * aLimit)/Ln(aLimit)) + 1
  14.   else
  15.     Result := Trunc(aLimit/(Ln(aLimit) - 2)) + 1;
  16. end;
  17.  
  18. function SelectPrimes(aLimit: Integer): TIntArray;
  19. var
  20.   IsPrime: array of Boolean;
  21.   I, SqrtBound: Integer;
  22.   Counter: Integer = 0;
  23.   J: Int64;
  24. begin
  25.   SetLength(IsPrime, Succ(aLimit));
  26.   FillChar(Pointer(IsPrime)^, Succ(aLimit), Byte(True));
  27.   SqrtBound := Trunc(Sqrt(aLimit));
  28.   SetLength(Result, EstimatePrimeCount(aLimit));
  29.   for I := 2 to aLimit do
  30.     if IsPrime[I] then
  31.       begin
  32.         Result[Counter] := I;
  33.         Inc(Counter);
  34.         if I <= SqrtBound then
  35.           begin
  36.             J := Int64(I) * Int64(I);
  37.             while J <= aLimit do
  38.               begin
  39.                 IsPrime[J] := False;
  40.                 Inc(J, I);
  41.               end;
  42.           end;
  43.       end;
  44.   SetLength(Result, Counter);
  45. end;
  46.  
  47. var
  48.   Primes: TIntArray;
  49.   Limit: Integer = -1;
  50.   I: Integer;
  51.  
  52. function ReadLimit: Boolean;
  53. var
  54.   Lim: DWord;
  55. begin
  56.   if (ParamCount <> 1) or not DWord.TryParse(ParamStr(1), Lim) then
  57.     exit(False);
  58.   if (Lim < 2) or (Lim > MAX_PRIME_LIMIT) then
  59.     exit(False);
  60.   Limit := Integer(Lim);
  61.   Result := True;
  62. end;
  63.  
  64. begin
  65.   if not ReadLimit then
  66.     begin
  67.       WriteLn('Usage: prime_sieve Limit');
  68.       WriteLn('  where Limit in the range [2, ', MAX_PRIME_LIMIT, ']');
  69.       Halt;
  70.     end;
  71.   Primes := SelectPrimes(Limit);
  72.   for I := 0 to High(Primes) - 1 do
  73.     Write(Primes[I], ', ');
  74.   WriteLn(Primes[High(Primes)]);
  75. end.
  76.  

julkas

  • Guest
Re: Pascal programs in the Rosetta Code collection
« Reply #32 on: August 12, 2019, 02:36:19 pm »
@avk In my Python code I always use bytearray. May be for Pascal bool and byte are same (in perf. terms).
Also I don't precompute prime count - simple primes.append()

avk

  • Hero Member
  • *****
  • Posts: 752
Re: Pascal programs in the Rosetta Code collection
« Reply #33 on: August 12, 2019, 02:50:24 pm »
@avk In my Python code I always use bytearray. May be for Pascal bool and byte are same (in perf. terms)...
I beleave, it is.

...Also I don't precompute prime count - simple primes.append()
It would be possible to use TVector<Integer>, but I suppose it would be slower.

avk

  • Hero Member
  • *****
  • Posts: 752
Re: Pascal programs in the Rosetta Code collection
« Reply #34 on: August 15, 2019, 10:27:23 am »
Some advanced sieve versions:
Code: Pascal  [Select][+][-]
  1. unit sieves;
  2. {$mode objfpc}{$H+}
  3. {$modeswitch advancedrecords}
  4. {$inline on}
  5. interface
  6. uses
  7.   SysUtils;
  8.  
  9. const
  10.   MAX_LIMIT = High(DWord) - 131072; //to avoid overflow
  11.  
  12. type
  13.   TDWordArray = array of DWord;
  14.  
  15. { uses bit-packed array }
  16.   function BpSieve(aLimit: DWord): TDWordArray;
  17. { uses bit-packed array and skips even numbers}
  18.   function OddSieve(aLimit: DWord): TDWordArray;
  19. { segmented sieve version }
  20.   function SegmentSieve(aLimit: DWord): TDWordArray;
  21.  
  22. implementation
  23.  
  24. type
  25.   TBitArray   = record
  26.   private
  27.     FBits: array of Byte;
  28.     function  GetBit(aIndex: DWord): Boolean; inline;
  29.     procedure SetBit(aIndex: DWord; aValue: Boolean); inline;
  30.   public
  31.     constructor Create(aSize: DWord);
  32.     property Bits[aIndex: DWord]: Boolean read GetBit write SetBit; default;
  33.   end;
  34.  
  35. function TBitArray.GetBit(aIndex: DWord): Boolean;
  36. begin
  37.   Result := Boolean(FBits[aIndex shr 3] and (1 shl (aIndex and 7)));
  38. end;
  39.  
  40. procedure TBitArray.SetBit(aIndex: DWord; aValue: Boolean);
  41. begin
  42.   if aValue then
  43.     FBits[aIndex shr 3] := FBits[aIndex shr 3] or (1 shl (aIndex and 7))
  44.   else
  45.     FBits[aIndex shr 3] := FBits[aIndex shr 3] and not(1 shl (aIndex and 7));
  46. end;
  47.  
  48. constructor TBitArray.Create(aSize: DWord);
  49. begin
  50.   SetLength(FBits, aSize shr 3 + Ord(aSize and 7 <> 0));
  51.   FillChar(Pointer(FBits)^, Length(FBits), $ff);
  52. end;
  53.  
  54. function EstimatePrimeCount(aLimit: DWord): DWord;
  55. begin
  56.   if aLimit < 2 then
  57.     exit(0);
  58.   if aLimit <= 200 then
  59.     Result := Trunc((1.6 * aLimit)/Ln(aLimit)) + 1
  60.   else
  61.     Result := Trunc(aLimit/(Ln(aLimit) - 2)) + 1;
  62. end;
  63.  
  64. function BpSieve(aLimit: DWord): TDWordArray;
  65. var
  66.   IsPrime: TBitArray;
  67.   I, J, SqrtBound: DWord;
  68.   Count: Integer = 0;
  69. begin
  70.   if aLimit > MAX_LIMIT then
  71.     raise Exception.Create('Prime limit exceeded');
  72.   if aLimit < 2 then
  73.     exit(nil);
  74.   IsPrime := TBitArray.Create(Succ(aLimit));
  75.   SqrtBound := Trunc(Sqrt(aLimit));
  76.   SetLength(Result, EstimatePrimeCount(aLimit));
  77.   for I := 2 to aLimit do
  78.     if IsPrime[I] then
  79.       begin
  80.         Result[Count] := I;
  81.         Inc(Count);
  82.         if I <= SqrtBound then
  83.           begin
  84.             J := I * I;
  85.             repeat
  86.               IsPrime[J] := False;
  87.               Inc(J, I);
  88.             until J > aLimit;
  89.           end;
  90.       end;
  91.   SetLength(Result, Count);
  92. end;
  93.  
  94. function OddSieve(aLimit: DWord): TDWordArray;
  95. var
  96.   IsPrime: TBitArray;
  97.   I: DWord = 3;
  98.   J, SqrtBound: DWord;
  99.   Count: Integer = 1;
  100. begin
  101.   if aLimit > MAX_LIMIT then
  102.     raise Exception.Create('Prime limit exceeded');
  103.   if aLimit < 2 then
  104.     exit(nil);
  105.   IsPrime := TBitArray.Create(aLimit div 2);
  106.   SqrtBound := Trunc(Sqrt(aLimit));
  107.   SetLength(Result, EstimatePrimeCount(aLimit));
  108.   Result[0] := 2;
  109.   while I <= aLimit do
  110.     begin
  111.       if IsPrime[(I - 3) shr 1] then
  112.         begin
  113.           Result[Count] := I;
  114.           Inc(Count);
  115.           if I <= SqrtBound then
  116.             begin
  117.               J := I * I;
  118.               repeat
  119.                 IsPrime[(J - 3) shr 1] := False;
  120.                 Inc(J, I shl 1);
  121.               until J > aLimit;
  122.             end;
  123.         end;
  124.       Inc(I, 2);
  125.     end;
  126.   SetLength(Result, Count);
  127. end;
  128.  
  129. function SegmentSieve(aLimit: DWord): TDWordArray;
  130. const
  131.   SEG_SIZE = $8000;
  132. var
  133.   Segment: array[0..Pred(SEG_SIZE)] of Boolean;
  134.   FirstPrimes, Primes: TDWordArray;
  135.   I, J, Prime: DWord;
  136.   Count: Integer = 0;
  137. begin
  138.   if aLimit > MAX_LIMIT then
  139.     raise Exception.Create('Prime limit exceeded');
  140.   if aLimit <= SEG_SIZE * 2 then
  141.     exit(OddSieve(aLimit));
  142.   I := Trunc(Sqrt(aLimit)) + 1;
  143.   FirstPrimes := OddSieve(I);
  144.   SetLength(Primes, EstimatePrimeCount(aLimit) - Length(FirstPrimes));
  145.   Dec(I);
  146.   while I < aLimit do
  147.     begin
  148.       FillChar(Segment, SEG_SIZE, Byte(True));
  149.       for Prime in FirstPrimes do
  150.         begin
  151.           J := I mod Prime;
  152.           if J <> 0 then
  153.             J := Prime - J;
  154.           while J < SEG_SIZE do
  155.             begin
  156.               Segment[J] := False;
  157.               Inc(J, Prime);
  158.             end;
  159.         end;
  160.       for J := 0 to Pred(SEG_SIZE) do
  161.         if (J + I <= aLimit) and Segment[J] then
  162.           begin
  163.             Primes[Count] := J + I;
  164.             Inc(Count);
  165.           end;
  166.       Inc(I, SEG_SIZE);
  167.     end;
  168.   SetLength(FirstPrimes, Length(FirstPrimes) + Count);
  169.   Move(Primes[0], FirstPrimes[Length(FirstPrimes) - Count], Count * SizeOf(DWord));
  170.   Result := FirstPrimes;
  171. end;
  172.  
  173. end.
  174.  

BrunoK

  • Sr. Member
  • ****
  • Posts: 452
  • Retired programmer
Re: Pascal programs in the Rosetta Code collection
« Reply #35 on: August 15, 2019, 04:28:10 pm »
And another entry that goes up to High(DWORD) values.
Code: Pascal  [Select][+][-]
  1. program pgmSieve;
  2.  
  3. {$mode objfpc}{$H+}
  4. {$modeswitch advancedrecords}
  5. {$inline on}
  6.  
  7. // Checking primes : https://primes.utm.edu/nthprime/index.php#nth
  8.  
  9. // Compiled -O1 -OoREGVAR so one can debug
  10.  
  11. uses
  12.   Classes,
  13.   SysUtils;
  14.  
  15. type
  16.  
  17.   { RSieveBits }
  18.  
  19.   RSieveBits = record
  20.   private
  21.     FBits : array of Byte;
  22.     FMaxPrime : DWord;
  23.  
  24.     FPrimeCount : integer; // Number of primes found
  25.     FMaxCount : integer;   // Maximum primes that can be stored
  26.     procedure FlagNonPrimes(aPrime : DWord);
  27.   public
  28.     procedure Init(aSize: DWord);
  29.     procedure BpSieve(var aPrimes : array of DWord; aSize : DWord);
  30.   end;
  31.  
  32. { Flag multiples as non primes }
  33. procedure RSieveBits.FlagNonPrimes(aPrime: DWord);
  34. var
  35.   lByteIx: integer;
  36.   lpByte: PByte;
  37.   lPrime : DWord;
  38.   lLastPrime : DWord;
  39. begin
  40.   lPrime := aPrime + aPrime;
  41.   while lPrime <= FMaxPrime do begin
  42.     if lPrime<=aPrime then // Wrap around -> exit
  43.       exit;
  44.     if (lPrime and 1) <> 0 then begin // Only odd values
  45.       lByteIx := lPrime shr 4;
  46.       lpByte:=@FBits[lByteIx];
  47.       if lpByte^ <> 0 then begin
  48.         lpByte^ := lpByte^ and not (1 shl ((lPrime and 15) shr 1));
  49.       end;
  50.     end;
  51.     inc(lPrime, aPrime);
  52.   end;
  53. end;
  54.  
  55. procedure RSieveBits.Init(aSize: DWord);
  56. var
  57.   lQWord : QWord;
  58. begin
  59.   lQWord := (QWord(aSize) shr 4) + 1;
  60.   SetLength(FBits, lQWord);  // Only odd bits
  61.   FillChar(FBits[0], Length(FBits), $ff);
  62.   FMaxPrime := aSize;
  63. end;
  64.  
  65. procedure RSieveBits.BpSieve(var aPrimes: array of DWord; aSize: DWord);
  66. var
  67.   lPrime: DWORD;
  68.   lByteIx: integer;
  69.   lByte: byte;
  70. begin
  71.   if aSize<2 then
  72.     exit;
  73.   Init(aSize);
  74.   FMaxCount := Length(aPrimes);
  75.   FBits[0] := 254; // One is NOT prime
  76.   aPrimes[FPrimeCount] := 2;
  77.   inc(FPrimeCount);
  78.   lPrime := 3;
  79.   while lPrime <= FMaxPrime do begin
  80.     lByteIx := lPrime shr 4;
  81.     lByte := FBits[lPrime shr 4];
  82.     if lByte<>0 then begin
  83.       if (lByte and (1 shl ((lPrime and 15) shr 1))) <> 0 then begin
  84.         if FPrimeCount>=FMaxCount then begin
  85.           WriteLn('Stopped after ', FPrimeCount, ' primes. Next prime to store = ', lPrime,
  86.                   ' due to result set array full');
  87.           Break;
  88.         end;
  89.         aPrimes[FPrimeCount] := lPrime;
  90.         inc(FPrimeCount);
  91.         FlagNonPrimes(lPrime);
  92.       end;
  93.       inc(lPrime, 2);   // Only odd numbers
  94.     end
  95.     else
  96.       lPrime := ((lByteIx + 1) shl 4) + 1;
  97.     if lPrime<=3 then // Wrap around -> done
  98.       break;
  99.   end;
  100. end;
  101.  
  102. const             // 1'000'000'000
  103.   cMaxPrime: DWord = 1000000000; // ~60 secs for high(DWORD);
  104. var
  105.   vPrimes: Array[0..Maxint div 8] of DWord;  // Array of primes found
  106.   i : integer;
  107.   lSieveBits : RSieveBits;
  108.   lTickStart : int64;
  109.   lFrom : integer;
  110. begin
  111.   lTickStart := GetTickCount64;
  112.   lSieveBits.BpSieve(vPrimes, cMaxPrime);
  113.   WriteLn('Max Prime=',cMaxPrime, ' Nb of primes found=',lSieveBits.FPrimeCount);
  114.   lFrom := lSieveBits.FPrimeCount-20; // Write last 20 primes in ranges
  115.   if lFrom<0 then
  116.     lFrom := 0;
  117.   for i := lFrom to lSieveBits.FPrimeCount-1 do
  118.     WriteLn(DWord(vPrimes[i]));
  119.   WriteLn('Done in ', GetTickCount64 - lTickStart, ' milliseconds');
  120.   ReadLn;
  121. end.
  122.  

Thaddy

  • Hero Member
  • *****
  • Posts: 14197
  • Probably until I exterminate Putin.
Re: Pascal programs in the Rosetta Code collection
« Reply #36 on: August 15, 2019, 07:51:13 pm »
Please note Rosetta Code is supposed to be algorithmically fast, but clean code. Not optimized to death.
I would suggest avk's odd solution as a comprehensible and fast solution compared to the current (elegant and short) one.
It is more about the expressiveness of a language than its speed. It is not a compiler shoot-out. Remember that.
Keep code simple and elegant on Rosetta code and use the expressiveness of the language.

You impress more with few lines. That's what it is about. (everything else can theoretically be handled by the compiler)

Alternatively write code for the compiler shout-out.... :D (ok, shoot... O:-) )

Btw: you are only allowed standard libraries from FreePascal, note that too. I always try to limit it to system and classes.

And there will be always noise on the line by speed freaks. Then they do not understand the purpose of Rosetta code and have zilch.none, whatsoever knowledge....about computer science. <grumpy   >:D >
This goes for many other languages entries too, so don't be offended.
« Last Edit: August 15, 2019, 09:51:08 pm by Thaddy »
Specialize a type, not a var.

avk

  • Hero Member
  • *****
  • Posts: 752
Re: Pascal programs in the Rosetta Code collection
« Reply #37 on: August 16, 2019, 08:05:47 am »
As for performance, on my machine for the limit of 1,000,000,000, the results look like this:
Code: Text  [Select][+][-]
  1.   basic sieve:      16.7s
  2.   BrunoK's BpSieve: 13s
  3.   BpSieve:          12.1s
  4.   OddSieve:         6.1s
  5.   wheel sieve:      4.2s
  6.   SegmentSieve:     3.3s  
  7.  
« Last Edit: August 16, 2019, 08:18:35 am by avk »

Kays

  • Hero Member
  • *****
  • Posts: 569
  • Whasup!?
    • KaiBurghardt.de
Re: Pascal programs in the Rosetta Code collection
« Reply #38 on: August 16, 2019, 12:39:19 pm »
[…] There is a Object Pascal category. […]
There is even a FreePascal category. And FreePascal/Lazarus.

I second Thaddy’s latest reasoning. I’d use RC to “lure” programmers into using FPC. If you wanna show (off your language skills, and) how to tweak code, write an article in our own wiki.
Yours Sincerely
Kai Burghardt

Thaddy

  • Hero Member
  • *****
  • Posts: 14197
  • Probably until I exterminate Putin.
Re: Pascal programs in the Rosetta Code collection
« Reply #39 on: August 16, 2019, 07:20:24 pm »
second Thaddy’s latest reasoning.

Thanks. The Rosetta stone history and its meaning gets often lost on Rosetta code.
It is really about the expressiveness, not about speed. Many "slow" languages have entries too, even scripting.
Specialize a type, not a var.

avk

  • Hero Member
  • *****
  • Posts: 752
Re: Pascal programs in the Rosetta Code collection
« Reply #40 on: August 17, 2019, 04:33:55 pm »
...I second Thaddy’s latest reasoning...
Any critique appreciated. :)

Well, I’ll try to explain why I personally don’t like this "pretty good" algorithm(first RC entry).
RC for this task says: "It is important that the sieve algorithm be the actual algorithm used to find prime numbers for the task.", but the above algorithm is not very practical. The bottleneck of the classic sieve of Eratosthenes is a rather large space complexity (O(N)),
which in particular leads to the fact that the real time complexity is worse than theoretical due to cache misses.
As for this algorithm, it more aggravates the situation by using a hashset(my guess) for storage, which increases memory consumption by an order of magnitude and is also much slower than plane array.
The snippets that I submitted demonstrate step by step how to reduce memory consumption and how it affects performance.
And yes, I have one more snippet - odd segmented sieve: :)
Code: Pascal  [Select][+][-]
  1. function OddSegmentSieve(aLimit: DWord): TDWordArray;
  2. const
  3.   SEG_SIZE = $8000;
  4. var
  5.   Segment: array[0..Pred(SEG_SIZE)] of Boolean;
  6.   NewPrimes: TDWordArray;
  7.   I, J, K, Prime: DWord;
  8.   Count: Integer = 0;
  9. begin
  10.   if aLimit > MAX_LIMIT then
  11.     raise Exception.Create('Prime limit exceeded');
  12.   if aLimit <= SEG_SIZE * 2 then
  13.     exit(OddSieve(aLimit));
  14.   I := Trunc(Sqrt(aLimit));
  15.   Result := OddSieve(Succ(I));
  16.   I += Ord(not Odd(I));
  17.   SetLength(NewPrimes, EstimatePrimeCount(aLimit) - Length(Result));
  18.   while I <= aLimit do
  19.     begin
  20.       FillChar(Segment, SizeOf(Segment), Byte(True));
  21.       for K := 1 to High(Result) do
  22.         begin
  23.           Prime := Result[K];
  24.           J := I mod Prime;
  25.           if J <> 0 then
  26.             if Boolean(J and 1) then
  27.               J := Prime - J
  28.             else
  29.               J := Prime shl 1 - J;
  30.           while J < SEG_SIZE * 2 do
  31.             begin
  32.               Segment[J shr 1] := False;
  33.               J += Prime shl 1;
  34.             end;
  35.         end;
  36.       K := Min(Pred(SEG_SIZE * 2), aLimit - I);
  37.       J := 0;
  38.       while J <= K do
  39.         begin
  40.           if Segment[J shr 1] then
  41.             begin
  42.               NewPrimes[Count] := J + I;
  43.               Inc(Count);
  44.             end;
  45.           J += 2;
  46.         end;
  47.       I += SEG_SIZE * 2;
  48.     end;
  49.   SetLength(Result, Length(Result) + Count);
  50.   Move(NewPrimes[0], Result[Length(Result) - Count], Count * SizeOf(DWord));
  51. end;
  52.  

The performance comparison table now looks like this:
Code: Text  [Select][+][-]
  1.   basic sieve:      16.7s
  2.   BrunoK's BpSieve: 13s
  3.   BpSieve:          12.1s
  4.   OddSieve:         6.1s
  5.   wheel sieve:      4.2s
  6.   SegmentSieve:     3.3s
  7.   OddSegmentSieve:  2.1s  
  8.  
Also lazy prime iterator seems to be very useful.
Any opinions? 

Thaddy

  • Hero Member
  • *****
  • Posts: 14197
  • Probably until I exterminate Putin.
Re: Pascal programs in the Rosetta Code collection
« Reply #41 on: August 17, 2019, 06:25:45 pm »
@avk
I would solve it like this:
- provide a very clean one
AND
- provide an understandable optimized one with references.
Specialize a type, not a var.

avk

  • Hero Member
  • *****
  • Posts: 752
Re: Pascal programs in the Rosetta Code collection
« Reply #42 on: August 18, 2019, 05:44:14 pm »
@Thaddy, thanks, you are probably right.

 

TinyPortal © 2005-2018