Blum integer: Difference between revisions
Mathematica / Wolfram language implementation of Blum integer task.
Pwmiller74 (talk | contribs) mNo edit summary |
Pwmiller74 (talk | contribs) (Mathematica / Wolfram language implementation of Blum integer task.) |
||
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}}==
| |||