Rosetta Code

Prime reciprocal sum: Difference between revisions

Page Discussion

Prime reciprocal sum (view source)

Revision as of 07:20, 13 May 2023

2,966 bytes added ,  1 year ago
m
→‎{{header|Free Pascal}}: using presieving with small primes ( 2..65519)
m (→‎{{header|Python}}: paste error)
m (→‎{{header|Free Pascal}}: using presieving with small primes ( 2..65519))
Line 196:
=={{header|Pascal}}==
==={{header|Free Pascal}}===
Most time consuming is finding the next prime.<br>
Now pre-sieving with the primes til 65535, to reduce tests.<br>
<syntaxhighlight lang=pascal>
That is the same as checking all numbers in sieve as divisible by the small primes.Since the prime are in big distance in that region, that's an improvement.
program primeRezSum;
<syntaxhighlight lang=pascal>program primeRezSum;
{$IFDEF FPC} {$MODE DELPHI}{$Optimization ON,ALL} {$ENDIF}
{$IFDEF WINDOWS}{$APPTYPE CONSOLE}{$ENDIF}
Line 204 ⟶ 205:
sysutils,
gmp;
const
PrimeCount = 6542;
PRIMEMAXVAL = 65535;
SIEVESIZE = 65536;
type
Tsieveprime = record
prime,
offset : Uint16;
end;
tSievePrimes = array[0..PrimeCount-1] of Tsieveprime;
tSieve = array[0..SIEVESIZE-1] of byte;
var
s : AnsiString;
MyPrimes : tSievePrimes;
sieve : tSieve;
procedure OutStr(const s: AnsiString);
var
Line 220 ⟶ 236:
writeln(pChar(@myString[1]),'...',pChar(@myString[l-19]),' (',l:6,' digits)');
end;
end;
 
function InitPrimes:Uint32;
var
f : extended;
idx,p,pr_cnt : Uint32;
Begin
fillchar(sieve,Sizeof(sieve),#0);
pr_cnt := 0;
p := 2;
f := 1.0;
repeat
while Sieve[p]<> 0 do
inc(p);
MyPrimes[pr_cnt].prime := p;
f := f*(p-1)/p;
inc(pr_cnt);
idx := p*p;
if idx > PRIMEMAXVAL then
Break;
repeat
sieve[idx] := 1;
inc(idx,p);
until idx > high(sieve);
inc(p);
until sqr(p)>PRIMEMAXVAL;
 
while (pr_cnt<= High(MyPrimes)) AND (p<PRIMEMAXVAL) do
begin
while Sieve[p]<> 0 do
inc(p);
MyPrimes[pr_cnt].prime := p;
f := f*(p-1)/p;
inc(p);
inc(pr_cnt);
end;
Writeln ('reducing factor ',f:10:8);
result := pr_cnt-1;
end;
 
procedure DoSieveOffsetInit(var tmp:mpz_t);
var
dummy :mpz_t;
i,j,p : Uint32;
Begin
mpz_init(dummy);
for i := 0 to High(MyPrimes) do
with MyPrimes[i] do
Begin
if prime = 0 then Begin writeln(i);halt;end;
offset := prime-mpz_tdiv_q_ui(dummy,tmp,prime);
end;
mpz_set(dummy,tmp);
repeat
//one sieve
fillchar(sieve,Sizeof(sieve),#0);
//sieve
For i := 0 to High(MyPrimes) do
begin
with MyPrimes[i] do
begin
p := prime;
j := offset;
end;
repeat
sieve[j] := 1;
j += p;
until j >= SIEVESIZE;
MyPrimes[i].offset := j-SIEVESIZE;
end;
 
j := 0;
For i := 0 to High(sieve) do
begin
// function mpz_probab_prime_p(var n: mpz_t; reps: longint): longint;1 = prime
if (sieve[i]= 0) then
begin
mpz_add_ui(dummy,dummy,i-j);
j := i;
if (mpz_probab_prime_p(dummy,1) >0) then
begin
mpz_set(tmp,dummy);
mpz_clear(dummy);
EXIT;
end;
end;
end;
until false;
end;
 
var
nominator,denominator,tmp,tmpDemP,p : mpz_t;
T1,T0:Int64;
s : AnsiString;
cnt : NativeUint;
begin
InitPrimes;
setlength(s,10000);
setlength(s,100000);
cnt := 1;
mpz_init(nominator);
Line 233 ⟶ 338:
mpz_init(tmpDemP);
mpz_init_set_ui(denominator,1);
mpz_init_set_ui(p,21);
 
repeat
mpz_sub_uimpz_set(ptmpDemP,p,1);
T0 := GetTickCount64;
mpz_nextprime(p,p);
write(if cnt:3,' > 9 ');then
DoSieveOffsetInit(p)
else
mpz_nextprime(p,p);
T1 := GetTickCount64;
write(cnt:3,' ',T1-T0,' ms ,delta ');
mpz_sub(tmpDemP,p,tmpDemP);
mpz_get_str(pChar(@s[1]),10, tmpDemP); OutStr(s);
mpz_get_str(pChar(@s[1]),10, p); OutStr(s);
if cnt >=15 then
Break;
repeat
mpz_mul(tmp,nominator,p);
Line 244 ⟶ 360:
if mpz_cmp(tmp,tmpDemP)< 0 then
BREAK;
mpz_nextprime mpz_add_ui(p,p,1);
until false;
mpz_get_str(pChar(@s[1]),10, p);
OutStr(s);
mpz_set(nominator,tmp);
mpz_mul(denominator,denominator,p);
Line 253 ⟶ 367:
//next smallest possible number denominator/delta
mpz_sub(tmp,denominator,nominator);
 
mpz_fdiv_q(p,denominator,tmp);
 
inc(cnt);
until cnt> 1417;
end.</syntaxhighlight>
{{out|@TIO.RUN}}
<pre>
 
1 2
reducing factor 0.05041709
2 3
31 70 ms ,delta 1
2
4 43
52 18110 ms ,delta 1
3
6 654149
73 270823151090 ms ,delta 1
7
8 153694141992520880899
4 0 ms ,delta 1
9 337110658273917297268061074384231117039
43
10 8424197597064114319...13803869133407474043 ( 76 digits)
5 0 ms ,delta 5
11 2030075381384823476...91313959045797597991 ( 150 digits)
1811
12 2032370538147127284...21649394434192763213 ( 297 digits)
6 0 ms ,delta 16
13 1274824659267207819...20708715953110886963 ( 592 digits)
654149
14 4674902516513883824...65355869250350888941 ( 1180 digits)
7 0 ms ,delta 5
Real time: 1.147 s User time: 1.093 s Sys. time: 0.046 s CPU share: 99.28 %</pre>
27082315109
8 0 ms ,delta 71
153694141992520880899
9 1 ms ,delta 14
337110658273917297268061074384231117039
10 1 ms ,delta 350
8424197597064114319...13803869133407474043 ( 76 digits)
11 1 ms ,delta 203
2030075381384823476...91313959045797597991 ( 150 digits)
12 3 ms ,delta 33
2032370538147127284...21649394434192763213 ( 297 digits)
13 69 ms ,delta 348
1274824659267207819...20708715953110886963 ( 592 digits)
14 234 ms ,delta 192
4674902516513883824...65355869250350888941 ( 1180 digits)
15 20391 ms ,delta 3510
1139012563947167462...31060548964273180103 ( 2358 digits)
 
Real time: 21.024 s User time: 20.768 s Sys. time: 0.067 s CPU share: 99.10 %
</pre>
 
=={{header|Python}}==
Horst
132

edits

Retrieved from "https://rosettacode.org/wiki/Special:MobileDiff/342241"