Ormiston triples: Difference between revisions
Content added Content deleted
SqrtNegInf (talk | contribs) (Added Perl) |
No edit summary |
||
Line 216: | Line 216: | ||
Number of Ormiston triples < 10000000000: 4925 |
Number of Ormiston triples < 10000000000: 4925 |
||
</pre> |
</pre> |
||
=={{header|Delphi}}== |
|||
{{works with|Delphi|6.0}} |
|||
{{libheader|SysUtils,StdCtrls}} |
|||
I went a little bit overboard on this problem. Delphi-6 only has a 32-bit memory space, so it would be hard to find triples to 10-billion when the memory only goes up to 2 billion. To solve the problem, I created a bit-boolean array object that stores 8 booleans per byte. I made this into a an object that simulates the functioning of an array, so access to array items is identical accessing a regular boolean array: |
|||
PrimeArray[I]:=True; |
|||
The whole problem is built as a series of hierchical objects: |
|||
1. Bit-Boolean Array Object. This object is responsilbe for accessing bits in each byte. It uses Delphi's array property features, which allows objects to be treated as arrays accessed with the same syntax as normal arrays: |
|||
PrimeArray[I]:=PrimeArray[I+1]; |
|||
2. Sieve Prime Object. This object builds and maintains a list of booliean flags indicating whether a number is prime or not. Since all even numbers are guaranteed not prime, the array only stores odd values. This cuts the space requirements by 50%. However, to the outside world the, array of primes appears to have all the even and odd values. The object uses Delphi's array property feature that allows array syntax to be used, but calls a subroutine to actually actuall read or write the data. The subroutine returns false for every even number and handles special cases such 0, 1 and 2. |
|||
3. Ormiston Triple Iterator Object. This object creats a prime table and then allows the user to iterate through all the Ormiston Triples. This makes it easy to perform different task such reading the first 25 triples or rapidly scanning through 1-billion or 10-billion triples to count the one that satisfy the Ormiston criteria. |
|||
The program takes 17 seconds to scan 1-billion primes and 10 minutes to test 10 billion. |
|||
<syntaxhighlight lang="Delphi"> |
|||
{Bit boolean - because it stores 8 bools per bytye, it will} |
|||
{handle up to 16 gigabyte in a 32 bit programming environment} |
|||
type TBitBoolArray = class(TObject) |
|||
private |
|||
FSize: int64; |
|||
ByteArray: array of Byte; |
|||
function GetValue(Index: int64): boolean; |
|||
procedure WriteValue(Index: int64; const Value: boolean); |
|||
function GetSize: int64; |
|||
procedure SetSize(const Value: int64); |
|||
protected |
|||
public |
|||
property Value[Index: int64]: boolean read GetValue write WriteValue; default; |
|||
constructor Create; |
|||
property Count: int64 read GetSize write SetSize; |
|||
procedure Clear(Value: boolean); |
|||
end; |
|||
{ TBitBoolArray } |
|||
const BitArray: array [0..7] of byte = ($01, $02, $04, $08, $10, $20, $40, $80); |
|||
function TBitBoolArray.GetValue(Index: int64): boolean; |
|||
begin |
|||
{Note: (Index and 7) is faster than (Index mod 8)} |
|||
Result:=(ByteArray[Index shr 3] and BitArray[Index and 7])<>0; |
|||
end; |
|||
procedure TBitBoolArray.WriteValue(Index: int64; const Value: boolean); |
|||
var Inx: int64; |
|||
begin |
|||
Inx:=Index shr 3; |
|||
{Note: (Index and 7) is faster than (Index mod 8)} |
|||
if Value then ByteArray[Inx]:=ByteArray[Inx] or BitArray[Index and 7] |
|||
else ByteArray[Inx]:=ByteArray[Inx] and not BitArray[Index and 7] |
|||
end; |
|||
constructor TBitBoolArray.Create; |
|||
begin |
|||
SetLength(ByteArray,0); |
|||
end; |
|||
function TBitBoolArray.GetSize: int64; |
|||
begin |
|||
Result:=FSize; |
|||
end; |
|||
procedure TBitBoolArray.SetSize(const Value: int64); |
|||
var Len: int64; |
|||
begin |
|||
FSize:=Value; |
|||
{Storing 8 items per byte} |
|||
Len:=Value div 8; |
|||
{We need one more to fill partial bits} |
|||
if (Value mod 8)<>0 then Inc(Len); |
|||
SetLength(ByteArray,Len); |
|||
end; |
|||
procedure TBitBoolArray.Clear(Value: boolean); |
|||
var Fill: byte; |
|||
begin |
|||
if Value then Fill:=$FF else Fill:=0; |
|||
FillChar(ByteArray[0],Length(ByteArray),Fill); |
|||
end; |
|||
+++++++++++++++++++++++++++++++++++++++++++++++++++++++++ |
|||
{Sieve object the generates and holds prime values} |
|||
{Enable this flag if you need primes past 2 billion. |
|||
The flag signals the code to use bit-booleans arrays |
|||
which can contain up to 8 x 4 gigabytes = 32 gig booleans.} |
|||
{$define BITBOOL} |
|||
type TPrimeSieve = class(TObject) |
|||
private |
|||
{$ifdef BITBOOL} |
|||
PrimeArray: TBitBoolArray; |
|||
{$else} |
|||
PrimeArray: array of boolean; |
|||
{$endif} |
|||
FArraySize: int64; |
|||
function GetPrime(Index: int64): boolean; |
|||
procedure Clear; |
|||
protected |
|||
procedure DoSieve; |
|||
property ArraySize: int64 read FArraySize; |
|||
public |
|||
BitBoolean: boolean; |
|||
constructor Create; |
|||
destructor Destroy; override; |
|||
procedure Intialize(Size: int64); |
|||
property Prime[Index: int64]: boolean read GetPrime; default; |
|||
function NextPrime(Start: int64): int64; |
|||
end; |
|||
procedure TPrimeSieve.Clear; |
|||
begin |
|||
{$ifdef BITBOOL} |
|||
PrimeArray.Clear(True); |
|||
{$else} |
|||
FillChar(PrimeArray[0],Length(PrimeArray),True); |
|||
{$endif} |
|||
end; |
|||
constructor TPrimeSieve.Create; |
|||
begin |
|||
{$ifdef BITBOOL} |
|||
PrimeArray:=TBitBoolArray.Create; |
|||
BitBoolean:=True; |
|||
{$else} |
|||
BitBoolean:=False; |
|||
{$endif} |
|||
end; |
|||
destructor TPrimeSieve.Destroy; |
|||
begin |
|||
{$ifdef BITBOOL} |
|||
PrimeArray.Free; |
|||
{$endif} |
|||
inherited; |
|||
end; |
|||
procedure TPrimeSieve.DoSieve; |
|||
{Load flags with true/false to flag that number is prime} |
|||
{Note: does not store even values, because except for 2, all primes are even} |
|||
{Starts storing flags at Index=3, so reading/writing routines compensate} |
|||
{Uses for-loops for boolean arrays and while-loops for bitbooleans arrays} |
|||
{$ifdef BITBOOL} |
|||
var Offset, I, K: int64; |
|||
{$else} |
|||
var Offset, I, K: cardinal; |
|||
{$endif} |
|||
begin |
|||
Clear; |
|||
{$ifdef BITBOOL} |
|||
I:=0; |
|||
while I<ArraySize do |
|||
{$else} |
|||
for I:=0 to ArraySize-1 do |
|||
{$endif} |
|||
begin |
|||
if PrimeArray[I] then |
|||
begin |
|||
Offset:= I + I + 3; |
|||
K:= I + Offset; |
|||
while K <=(ArraySize-1) do |
|||
begin |
|||
PrimeArray[K]:= False; |
|||
K:= K + Offset; |
|||
end; |
|||
end; |
|||
{$ifdef BITBOOL} Inc(I); {$endif} |
|||
end; |
|||
end; |
|||
function TPrimeSieve.GetPrime(Index: int64): boolean; |
|||
{Get a prime flag from array - compensates} |
|||
{ for 0,1,2 and even numbers not being stored} |
|||
begin |
|||
if Index in [0,1,2] then Result:=True |
|||
else if (Index and 1)=0 then Result:=false |
|||
else Result:=PrimeArray[(Index div 2)-1]; |
|||
end; |
|||
function TPrimeSieve.NextPrime(Start: int64): int64; |
|||
{Get next prime after Start} |
|||
begin |
|||
Result:=Start+1; |
|||
while Result<=((ArraySize-1) * 2) do |
|||
begin |
|||
if Self.Prime[Result] then break; |
|||
Inc(Result); |
|||
end; |
|||
end; |
|||
procedure TPrimeSieve.Intialize(Size: int64); |
|||
{Set array size and do Sieve to load flag array with} |
|||
begin |
|||
FArraySize:=Size div 2; |
|||
{$ifdef BITBOOL} |
|||
PrimeArray.Count:=FArraySize; |
|||
{$else} |
|||
SetLength(PrimeArray,FArraySize); |
|||
{$endif} |
|||
DoSieve; |
|||
end; |
|||
{-------------------------------------------------------------------------------} |
|||
type TTripleInfo = record |
|||
Prime1,Prime2,Prime3: int64; |
|||
Count: int64; |
|||
end; |
|||
{Iterator for Ormiston Triple} |
|||
type TOrm3Iterator = class(TObject) |
|||
FInfo: TTripleInfo; |
|||
PS: TPrimeSieve; |
|||
private |
|||
function IsOrmistonTriple(P1, P2, P3: int64): boolean; |
|||
function EncodeNumber(N: int64): int64; |
|||
protected |
|||
public |
|||
procedure Reset; |
|||
procedure SetSize(Size: int64); |
|||
function GetNext(Limit: int64; var Info: TTripleInfo): boolean; |
|||
constructor Create; |
|||
destructor Destroy; override; |
|||
end; |
|||
procedure TOrm3Iterator.Reset; |
|||
{Restart iterator} |
|||
begin |
|||
FInfo.Count:=0; |
|||
FInfo.Prime1:=1; FInfo.Prime2:=3; FInfo.Prime3:=5; |
|||
end; |
|||
procedure TOrm3Iterator.SetSize(Size: int64); |
|||
begin |
|||
PS.Intialize(Size); |
|||
Reset; |
|||
end; |
|||
constructor TOrm3Iterator.Create; |
|||
begin |
|||
PS:=TPrimeSieve.Create; |
|||
{Start with trivial prime set} |
|||
SetSize(100); |
|||
end; |
|||
destructor TOrm3Iterator.Destroy; |
|||
begin |
|||
PS.Free; |
|||
inherited; |
|||
end; |
|||
function TOrm3Iterator.EncodeNumber(N: int64): int64; |
|||
{Encode N by counting digits 0..9 into nibbles} |
|||
{Get the product of the integers in a number} |
|||
var T: integer; |
|||
const NibMap: array [0..9] of int64 = ( |
|||
{0} $1, {1} $10, {2} $100, {3} $1000, {4} $10000, {5} $100000, |
|||
{6} $1000000, {7} $10000000, {8} $100000000, {9} $1000000000); |
|||
begin |
|||
Result:=0; |
|||
repeat |
|||
begin |
|||
T:=N mod 10; |
|||
N:=N div 10; |
|||
Result:=Result + NibMap[T]; |
|||
end |
|||
until N<1; |
|||
end; |
|||
function TOrm3Iterator.IsOrmistonTriple(P1,P2,P3: int64): boolean; |
|||
var Pd1,Pd2,Pd3: int64; |
|||
begin |
|||
Result:=False; |
|||
{Optimization - difference in primes should be multiple of 18} |
|||
if (((P2 - P1) mod 18)<>0) or (((p3 - p2) mod 18)<>0) then exit; |
|||
Pd1:=EncodeNumber(P1); |
|||
Pd2:=EncodeNumber(P2); |
|||
if Pd1<>Pd2 then exit; |
|||
Pd3:=EncodeNumber(P3); |
|||
Result:=Pd2=Pd3; |
|||
end; |
|||
function TOrm3Iterator.GetNext(Limit: int64; var Info: TTripleInfo): boolean; |
|||
{Iterate to next Ormiston Pair - automatically stop at prime>Limit} |
|||
{Returns false if it hits limit - true if it found Next Ormiston Pair} |
|||
begin |
|||
Result:=False; |
|||
while true do |
|||
begin |
|||
{Get next set of primes} |
|||
FInfo.Prime1:=FInfo.Prime2; |
|||
FInfo.Prime2:=FInfo.Prime3; |
|||
FInfo.Prime3:=PS.NextPrime(FInfo.Prime3); |
|||
{Abort if 3rd prime is ove limit} |
|||
if FInfo.Prime3>=Limit then break; |
|||
{Test if it is an Ormiston triple} |
|||
if IsOrmistonTriple(FInfo.Prime1,FInfo.Prime2,FInfo.Prime3) then |
|||
begin |
|||
{Return info on triple} |
|||
Inc(FInfo.Count); |
|||
Info:=FInfo; |
|||
Result:=True; |
|||
break; |
|||
end; |
|||
end; |
|||
end; |
|||
procedure ShowOrmistonTriple(Memo: TMemo); |
|||
var I: integer; |
|||
var S: string; |
|||
var OI: TOrm3Iterator; |
|||
var Info: TTripleInfo; |
|||
const Limit = 10000000000; |
|||
const DisMod = Limit div 10000000; |
|||
const Bill = Limit div 1000000000; |
|||
var NS: string; |
|||
procedure DisplayTitle; |
|||
begin |
|||
NS:=IntToStr(Bill)+'-Billion'; |
|||
Memo.Lines.Add('====== Find Ormiston Triples ======'); |
|||
Memo.Lines.Add('First 25 plus number to '+NS); |
|||
S:='Bit-Boolean Array: '; |
|||
if OI.PS.BitBoolean then S:=S+'Yes' else S:=S+'No'; |
|||
Memo.Lines.Add(S); |
|||
end; |
|||
begin |
|||
{Create iterator} |
|||
OI:=TOrm3Iterator.Create; |
|||
try |
|||
DisplayTitle; |
|||
Memo.Lines.Add('Sieving Primes'); |
|||
OI.SetSize(Limit); |
|||
Memo.Lines.Add('Finding first 25 Triples'); |
|||
{Iterate throug 1st 25 tuples} |
|||
for I:=1 to 25 do |
|||
begin |
|||
OI.GetNext(High(Int64),Info); |
|||
Memo.Lines.Add(Format('%3d - (%6D %6D %6D) ',[I,Info.Prime1,Info.Prime2,Info.Prime3])); |
|||
end; |
|||
Memo.Lines.Add('Count='+IntToStr(Info.Count)); |
|||
Memo.Lines.Add('Counting triples to '+NS); |
|||
{Iterate to limit number of triples} |
|||
while OI.GetNext(Limit,Info) do |
|||
if (Info.Count mod DisMod)=0 then Memo.Lines.Add('Count='+IntToStr(Info.Count)); |
|||
Memo.Lines.Add(NS+'='+IntToStr(Info.Count)); |
|||
finally OI.Free; end; |
|||
end; |
|||
</syntaxhighlight> |
|||
{{out}} |
|||
<pre> |
|||
====== Find Ormiston Triples ====== |
|||
First 25 plus number to 1-Billion |
|||
Bit-Boolean Array: No |
|||
Sieving Primes |
|||
Finding first 25 Triples |
|||
1 - (11117123 11117213 11117321) |
|||
2 - (12980783 12980837 12980873) |
|||
3 - (14964017 14964071 14964107) |
|||
4 - (32638213 32638231 32638321) |
|||
5 - (32964341 32964413 32964431) |
|||
6 - (33539783 33539837 33539873) |
|||
7 - (35868013 35868031 35868103) |
|||
8 - (44058013 44058031 44058103) |
|||
9 - (46103237 46103273 46103327) |
|||
10 - (48015013 48015031 48015103) |
|||
11 - (50324237 50324273 50324327) |
|||
12 - (52402783 52402837 52402873) |
|||
13 - (58005239 58005293 58005329) |
|||
14 - (60601237 60601273 60601327) |
|||
15 - (61395239 61395293 61395329) |
|||
16 - (74699789 74699879 74699897) |
|||
17 - (76012879 76012897 76012987) |
|||
18 - (78163123 78163213 78163231) |
|||
19 - (80905879 80905897 80905987) |
|||
20 - (81966341 81966413 81966431) |
|||
21 - (82324237 82324273 82324327) |
|||
22 - (82523017 82523071 82523107) |
|||
23 - (83279783 83279837 83279873) |
|||
24 - (86050781 86050817 86050871) |
|||
25 - (92514341 92514413 92514431) |
|||
Count=25 |
|||
Counting triples to 1-Billion |
|||
Count=100 |
|||
Count=200 |
|||
Count=300 |
|||
1-Billion=368 |
|||
Elapsed Time: 17.303 Sec. |
|||
====== Find Ormiston Triples ====== |
|||
First 25 plus number to 10-Billion |
|||
Bit-Boolean Array: Yes |
|||
Sieving Primes |
|||
Finding first 25 Triples |
|||
1 - (11117123 11117213 11117321) |
|||
2 - (12980783 12980837 12980873) |
|||
3 - (14964017 14964071 14964107) |
|||
4 - (32638213 32638231 32638321) |
|||
5 - (32964341 32964413 32964431) |
|||
6 - (33539783 33539837 33539873) |
|||
7 - (35868013 35868031 35868103) |
|||
8 - (44058013 44058031 44058103) |
|||
9 - (46103237 46103273 46103327) |
|||
10 - (48015013 48015031 48015103) |
|||
11 - (50324237 50324273 50324327) |
|||
12 - (52402783 52402837 52402873) |
|||
13 - (58005239 58005293 58005329) |
|||
14 - (60601237 60601273 60601327) |
|||
15 - (61395239 61395293 61395329) |
|||
16 - (74699789 74699879 74699897) |
|||
17 - (76012879 76012897 76012987) |
|||
18 - (78163123 78163213 78163231) |
|||
19 - (80905879 80905897 80905987) |
|||
20 - (81966341 81966413 81966431) |
|||
21 - (82324237 82324273 82324327) |
|||
22 - (82523017 82523071 82523107) |
|||
23 - (83279783 83279837 83279873) |
|||
24 - (86050781 86050817 86050871) |
|||
25 - (92514341 92514413 92514431) |
|||
Count=25 |
|||
Counting triples to 10-Billion |
|||
Count=1000 |
|||
Count=2000 |
|||
Count=3000 |
|||
Count=4000 |
|||
10-Billion=4925 |
|||
Elapsed Time: 09:52.882 min |
|||
</pre> |
|||
=={{header|F_Sharp|F#}}== |
=={{header|F_Sharp|F#}}== |