### Bookstore

 Computer Math and Games in Pascal (preview) Lazarus Handbook (preview only)

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

#### julkas

• Sr. Member
• Posts: 306
• KISS principle / Lazarus 2.0.0 / FPC 3.0.4
##### 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
procedure mulu64(a, b: QWORD; out clo, chi: QWORD); assembler;
asm
mov rax, a
mov rdx, b
mul rdx
mov [clo], rax
mov [chi], rdx
end;
(* Pointer game *) Inc(ptr, 1); (* vs *) ptr := ptr + 1;

#### avk

• Jr. Member
• Posts: 95
##### 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.
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
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

• Sr. Member
• Posts: 306
• KISS principle / Lazarus 2.0.0 / FPC 3.0.4
##### 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()
procedure mulu64(a, b: QWORD; out clo, chi: QWORD); assembler;
asm
mov rax, a
mov rdx, b
mul rdx
mov [clo], rax
mov [chi], rdx
end;
(* Pointer game *) Inc(ptr, 1); (* vs *) ptr := ptr + 1;

#### avk

• Jr. Member
• Posts: 95
##### 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

• Jr. Member
• Posts: 95
##### Re: Pascal programs in the Rosetta Code collection
« Reply #34 on: August 15, 2019, 10:27:23 am »
Code: Pascal  [Select]
1. unit sieves;
2. {\$mode objfpc}{\$H+}
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

• Full Member
• Posts: 158
• 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+}
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');
121. end.
122.
Lazarus trunk r. 59978/03.01.2019 (+/- patches regarding enabled, TScrollBar, TCursorImage). FPC 3.0.4 32 bits. (+heaptrc with leaked ClassName+Revisited TList) , Windows 10 Pro x64 (v. 1803)

• Hero Member
• Posts: 8484
##### 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.... (ok, shoot... )

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    >
This goes for many other languages entries too, so don't be offended.
« Last Edit: August 15, 2019, 09:51:08 pm by Thaddy »
Read the manuals and if you are a professional get a proper education in computer science. Makes the forum a lot cleaner.

#### avk

• Jr. Member
• Posts: 95
##### 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

• Full Member
• Posts: 170
• Whasup!?
##### 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

• Hero Member
• Posts: 8484
##### Re: Pascal programs in the Rosetta Code collection
« Reply #39 on: August 16, 2019, 07:20:24 pm »

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.
Read the manuals and if you are a professional get a proper education in computer science. Makes the forum a lot cleaner.

#### avk

• Jr. Member
• Posts: 95
##### Re: Pascal programs in the Rosetta Code collection
« Reply #40 on: August 17, 2019, 04:33:55 pm »
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?

• Hero Member
• Posts: 8484
##### 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.
Read the manuals and if you are a professional get a proper education in computer science. Makes the forum a lot cleaner.

#### avk

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