Blum integer: Difference between revisions
m
→{{header|Phix}}: use pygments
Pwmiller74 (talk | contribs) mNo edit summary |
m (→{{header|Phix}}: use pygments) |
||
| (3 intermediate revisions by one other user not shown) | |||
Line 1,183:
=={{header|Mathematica}} / {{header|Wolfram Language}}==
<syntaxhighlight lang="mathematica">
ClearAll[BlumIntegerQ, BlumIntegersInRange, PrimePi2, BlumCount, binarySearch, BlumInts, timing, upperLimitEstimate, lastDigit, lastDigitDistributionPercents];
BlumIntegerQ[n_Integer] := With[{factors = FactorInteger[n]},
n ~ Mod ~ 4 == 1 &&
Length[factors] == 2 &&
factors[[1, 1]] ~ Mod ~ 4 == 3 &&
Last@Total@factors == 2
];
SetAttributes[BlumIntegerQ, Listable];
BlumIntegersInRange[n_Integer] := BlumIntegersInRange[1, n];
BlumIntegersInRange[start_Integer, end_Integer] :=
Select[Range[start + (4 - start) ~ Mod ~ 4, end, 4] + 1, BlumIntegerQ];
(* Counts semiprimes. See https://people.maths.ox.ac.uk/erban/papers/paperDCRE.pdf *)
PrimePi2[x_] := (PrimePi[Sqrt[x]] - PrimePi[Sqrt[x]]^2)/2 + Sum[PrimePi[x/Prime[p]], {p, 1, PrimePi[Sqrt[x]]}];
SetAttributes[PrimePi2, Listable];
(* Blum integers are semiprimes that are 1 mod 4, with two distinct factors where both factors are 3 mod 4. The following function gives an approximation of the number of Blum integers <= x.
According to my tests, this function tends to overestimate by approximately 5% in the range we're interested in.
*)
BlumCount[x_] := Ceiling[(PrimePi2[x] - PrimePi[Sqrt[x]]) / 4];
SetAttributes[BlumCount, Listable];
binarySearch[f_, targetValue_] :=
Module[{lo = 1, mid, hi = 1, currentValue},
While[f[hi] < targetValue,
hi *= 2;];
While[lo <= hi,
mid = Ceiling[(lo + hi) / 2];
currentValue = f[mid];
If[currentValue < targetValue,
lo = mid + 1;];
If[currentValue > targetValue,
hi = mid - 1;];
If[currentValue == targetValue,
While[f[mid] == targetValue,
mid++;
];
Return[mid - 1];
];
];
];
lastDigit[n_Integer] := n ~ Mod ~ 10;
SetAttributes[lastDigit, Listable];
upperLimitEstimate = Ceiling[binarySearch[BlumCount, 400000]* 1.1];
timing = First@AbsoluteTiming[BlumInts = BlumIntegersInRange[upperLimitEstimate];];
lastDigitDistributionPercents = N[Counts@lastDigit@BlumInts[[;; 400000]] / 4000, 5];
Print["Calculated the first ", Length[BlumInts],
" Blum integers in ", timing, " seconds."];
Print[];
Print["First 50 Blum integers:"];
Print[TableForm[Partition[BlumInts[[;; 50]], 10],
TableAlignments -> Right]];
Print[];
Print[Grid[
Table[{"The ", n , "th Blum integer is: ",
BlumInts[[n]]}, {n, {26828, 100000, 200000, 300000, 400000}}] ,
Alignment -> Right]]
Print[];
Print["% distribution of the first 400,000 Blum integers:"];
Print[Grid[
Table[{lastDigitDistributionPercents[n], "% end in ",
n}, {n, {1, 3, 7, 9}} ], Alignment -> Right]];
</syntaxhighlight>
{{out}}
<pre>
Calculated the first 416420 Blum integers in 15.1913 seconds.
First 50 Blum integers:
21 33 57 69 77 93 129 133 141 161
177 201 209 213 217 237 249 253 301 309
321 329 341 381 393 413 417 437 453 469
473 489 497 501 517 537 553 573 581 589
597 633 649 669 681 713 717 721 737 749
The 26828 th Blum integer is: 524273
The 100000 th Blum integer is: 2075217
The 200000 th Blum integer is: 4275533
The 300000 th Blum integer is: 6521629
The 400000 th Blum integer is: 8802377
% distribution of the first 400,000 Blum integers:
25.001% end in 1
25.017% end in 3
24.997% end in 7
24.985% end in 9
</pre>
=={{header|Maxima}}==
Line 1,533 ⟶ 1,632:
{{trans|Pascal}}
You can run this online [http://phix.x10.mx/p2js/Blum.htm here].
<!--
<syntaxhighlight lang="phix">
with javascript_semantics
constant LIMIT = 1e7, N = floor((floor(LIMIT/3)-1)/4)+1
function Sieve4n_3_Primes()
sequence sieve = repeat(0,N), p4n3 = {}
for idx=1 to N do
if sieve[idx]=0 then
integer n = idx*4-1
p4n3 &= n
if idx+n>N then
// collect the rest
for j=idx+1 to N do
if sieve[j]=0 then
end
exit
end if
for j=idx+n to N by n do
sieve[j] = 1
end for
end if
end for
return p4n3
end function
sequence p4n3 = Sieve4n_3_Primes(),
BlumField = repeat(false,LIMIT),
Blum50 = {}, counts = repeat(0,10)
for idx,n in p4n3 do
for bj in p4n3 from idx+1 do
atom k = n*bj
if k>LIMIT then exit end if
BlumField[k] = true
end for
end for
integer count = 0
for n,k in BlumField do
if k then
if count<50 then Blum50 &= n end if
counts[remainder(n,10)] += 1
count += 1
if count=50 then
printf(1,"First 50 Blum integers:\n%s\n",{join_by(Blum50,1,10," ",fmt:="%3d")})
elsif count=26828 or remainder(count,1e5)=0 then
printf(1,"The %,7d%s Blum integer is: %,9d\n", {count,ord(count),n})
if count=4e5 then exit end if
end if
end if
end for
printf(1,"\nPercentage distribution of the first 400,000 Blum integers:\n")
for i,n in counts do
if n then
printf(1," %6.3f%% end in %d\n", {n/4000, i})
end if
end for
</syntaxhighlight>
{{out}}
<pre>
| |||