Blum integer: Difference between revisions
m
→{{header|Phix}}: use pygments
(Haskell implementation) |
m (→{{header|Phix}}: use pygments) |
||
| (10 intermediate revisions by 5 users not shown) | |||
Line 139:
24.985% end with 9
</pre>
=={{header|BASIC}}==
==={{header|BASIC256}}===
{{trans|FreeBASIC}}
<syntaxhighlight lang="vb">global Prime1
n = 3
c = 0
print "The first 50 Blum integers:"
while True
if isSemiprime(n) then
if Prime1 % 4 = 3 then
Prime2 = n / Prime1
if (Prime2 <> Prime1) and (Prime2 % 4 = 3) then
c += 1
if c <= 50 then
print rjust(string(n), 4);
if c % 10 = 0 then print
end if
if c >= 26828 then
print : print "The 26828th Blum integer is: "; n
exit while
end if
end if
end if
end if
n += 2
end while
end
function isSemiprime(n)
d = 3
c = 0
while d*d <= n
while n % d = 0
if c = 2 then return false
n /= d
c += 1
end while
d += 2
end while
Prime1 = n
return c = 1
end function</syntaxhighlight>
==={{header|FreeBASIC}}===
<syntaxhighlight lang="vb">Dim Shared As Uinteger Prime1
Dim As Uinteger n = 3, c = 0, Prime2
Function isSemiprime(n As Uinteger) As Boolean
Dim As Uinteger d = 3, c = 0
While d*d <= n
While n Mod d = 0
If c = 2 Then Return False
n /= d
c += 1
Wend
d += 2
Wend
Prime1 = n
Return c = 1
End Function
Print "The first 50 Blum integers:"
Do
If isSemiprime(n) Then
If Prime1 Mod 4 = 3 Then
Prime2 = n / Prime1
If (Prime2 <> Prime1) And (Prime2 Mod 4 = 3) Then
c += 1
If c <= 50 Then
Print Using "####"; n;
If c Mod 10 = 0 Then Print
End If
If c >= 26828 Then
Print !"\nThe 26828th Blum integer is: " ; n
Exit Do
End If
End If
End If
End If
n += 2
Loop
Sleep</syntaxhighlight>
{{out}}
<pre>The 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 26828th Blum integer is: 524273</pre>
==={{header|Gambas}}===
{{trans|FreeBASIC}}
<syntaxhighlight lang="vbnet">Public Prime1 As Integer
Public Sub Main()
Dim n As Integer = 3, c As Integer = 0, Prime2 As Integer
Print "The first 50 Blum integers:"
Do
If isSemiprime(n) Then
If Prime1 Mod 4 = 3 Then
Prime2 = n / Prime1
If (Prime2 <> Prime1) And (Prime2 Mod 4 = 3) Then
c += 1
If c <= 50 Then
Print Format$(n, "####");
If c Mod 10 = 0 Then Print
End If
If c >= 26828 Then
Print "\nThe 26828th Blum integer is: "; n
Break
End If
End If
End If
End If
n += 2
Loop
End
Function isSemiprime(n As Integer) As Boolean
Dim d As Integer = 3, c As Integer = 0
While d * d <= n
While n Mod d = 0
If c = 2 Then Return False
n /= d
c += 1
Wend
d += 2
Wend
Prime1 = n
Return c = 1
End Function </syntaxhighlight>
=={{header|C}}==
Line 218 ⟶ 359:
Same as Wren example.
</pre>
=={{header|C#}}==
{{trans|Java}}
<syntaxhighlight lang="C#">
using System;
using System.Collections.Generic;
public class BlumInteger
{
public static void Main(string[] args)
{
int[] blums = new int[50];
int blumCount = 0;
Dictionary<int, int> lastDigitCounts = new Dictionary<int, int>();
int number = 1;
while (blumCount < 400000)
{
int prime = LeastPrimeFactor(number);
if (prime % 4 == 3)
{
int quotient = number / prime;
if (quotient != prime && IsPrimeType3(quotient))
{
if (blumCount < 50)
{
blums[blumCount] = number;
}
if (!lastDigitCounts.ContainsKey(number % 10))
{
lastDigitCounts[number % 10] = 0;
}
lastDigitCounts[number % 10]++;
blumCount++;
if (blumCount == 50)
{
Console.WriteLine("The first 50 Blum integers:");
for (int i = 0; i < 50; i++)
{
Console.Write($"{blums[i],3}");
Console.Write((i % 10 == 9) ? Environment.NewLine : " ");
}
Console.WriteLine();
}
else if (blumCount == 26828 || blumCount % 100000 == 0)
{
Console.WriteLine($"The {blumCount}th Blum integer is: {number}");
if (blumCount == 400000)
{
Console.WriteLine();
Console.WriteLine("Percent distribution of the first 400000 Blum integers:");
foreach (var key in lastDigitCounts.Keys)
{
Console.WriteLine($" {((double)lastDigitCounts[key] / 4000):0.000}% end in {key}");
}
}
}
}
}
number += (number % 5 == 3) ? 4 : 2;
}
}
private static bool IsPrimeType3(int number)
{
if (number < 2) return false;
if (number % 2 == 0) return number == 2;
if (number % 3 == 0) return number == 3;
for (int divisor = 5; divisor * divisor <= number; divisor += 2)
{
if (number % divisor == 0) return false;
}
return number % 4 == 3;
}
private static int LeastPrimeFactor(int number)
{
if (number == 1) return 1;
if (number % 2 == 0) return 2;
if (number % 3 == 0) return 3;
if (number % 5 == 0) return 5;
for (int divisor = 7; divisor * divisor <= number; divisor += 2)
{
if (number % divisor == 0) return divisor;
}
return number;
}
}
</syntaxhighlight>
{{out}}
<pre>
The 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 26828th Blum integer is: 524273
The 100000th Blum integer is: 2075217
The 200000th Blum integer is: 4275533
The 300000th Blum integer is: 6521629
The 400000th Blum integer is: 8802377
Percent distribution of the first 400000 Blum integers:
25.001% end in 1
25.017% end in 3
24.997% end in 7
24.985% end in 9
</pre>
=={{header|C++}}==
Line 315 ⟶ 573:
=={{header|EasyLang}}==
{{trans|FreeBASIC}}
<syntaxhighlight lang=easylang>fastfunc semiprim n .
d = 3
while d * d <= n
Line 353 ⟶ 610:
.
print ""
print "The 26828th Blum integer is: " & n</syntaxhighlight>
=={{header|Go}}==
Line 975 ⟶ 1,181:
end
</syntaxhighlight>{{out}} Same as Wren, Go, etc
=={{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,278 ⟶ 1,585:
{{trans|Raku}}
{{libheader|ntheory}}
<syntaxhighlight lang="perl" line>use v5.36;
use ntheory qw(is_square is_semiprime factor vecall);
sub comma { reverse
sub table ($c, @V) { my $t = $c * (my $w = 5); (
sub is_blum ($n) {
($n % 4) == 1 && is_semiprime($n) && !is_square($n) && vecall { ($_ % 4) == 3 } factor($n);
}
my @nth = (26828, 1e5, 2e5, 3e5, 4e5);
my (@blum, %C);
for (my $i = 1 ; ; ++$i) {
push @blum, $i if is_blum $i;
last if $nth[-1] ==
}
$C{
say "The first fifty Blum integers:\n" . table 10, @blum[0 .. 49];
printf "The %7sth Blum integer: %9s\n", comma($_), comma $blum[$_ - 1] for @nth;
say '';
printf "$_: %6d (%.3f%%)\n", $C{$_}, 100 * $C{$_} / scalar @blum for <1 3 7 9>;</syntaxhighlight>
{{out}}
<pre>
Line 1,348 ⟶ 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>
Line 1,698 ⟶ 1,983:
26828th Blum number: 524273
</pre>
=={{header|Scala}}==
{{trans|Java}}
<syntaxhighlight lang="Scala">
import scala.collection.mutable
object BlumInteger extends App {
var blums = new Array[Int](50)
var blumCount = 0
val lastDigitCounts = mutable.Map[Int, Int]()
var number = 1
while (blumCount < 400_000) {
val prime = leastPrimeFactor(number)
if (prime % 4 == 3) {
val quotient = number / prime
if (quotient != prime && isPrimeType3(quotient)) {
if (blumCount < 50) {
blums(blumCount) = number
}
lastDigitCounts(number % 10) = lastDigitCounts.getOrElse(number % 10, 0) + 1
blumCount += 1
if (blumCount == 50) {
println("The first 50 Blum integers:")
blums.grouped(10).foreach(group => println(group.map(i => f"$i%3d").mkString(" ")))
println("")
} else if (blumCount == 26828 || blumCount % 100_000 == 0) {
println(f"The ${blumCount}th Blum integer is: $number%7d")
if (blumCount == 400_000) {
println("\nPercent distribution of the first 400000 Blum integers:")
lastDigitCounts.foreach { case (key, count) =>
println(f" ${count.toDouble / 4000}%6.3f%% end in $key")
}
}
}
}
}
number += (if (number % 5 == 3) 4 else 2)
}
def isPrimeType3(aNumber: Int): Boolean = {
if (aNumber < 2) return false
if (aNumber % 2 == 0) return aNumber == 2
if (aNumber % 3 == 0) return aNumber == 3
var divisor = 5
while (divisor * divisor <= aNumber) {
if (aNumber % divisor == 0) return false
divisor += 2
}
aNumber % 4 == 3
}
def leastPrimeFactor(aNumber: Int): Int = {
if (aNumber == 1) return 1
if (aNumber % 2 == 0) return 2
if (aNumber % 3 == 0) return 3
var divisor = 5
while (divisor * divisor <= aNumber) {
if (aNumber % divisor == 0) return divisor
divisor += 2
}
aNumber
}
}
</syntaxhighlight>
{{out}}
<pre>
The 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 26828th Blum integer is: 524273
The 100000th Blum integer is: 2075217
The 200000th Blum integer is: 4275533
The 300000th Blum integer is: 6521629
The 400000th Blum integer is: 8802377
Percent distribution of the first 400000 Blum integers:
25.001% end in 1
25.017% end in 3
24.997% end in 7
24.985% end in 9
</pre>
=={{header|Sidef}}==
Takes about 30 seconds:
<syntaxhighlight lang="ruby">func blum_integers(upto) {
var L = []
var P = idiv(upto, 3).primes.grep{ .is_congruent(3, 4) }
for i in (1..P.end) {
var p = P[i]
for j in (^i) {
var t = p*P[j]
break if (t > upto)
L << t
}
}
L.sort
}
func blum_first(n) {
var upto = int(4.5*n*log(n) / log(log(n)))
loop {
var B = blum_integers(upto)
if (B.len >= n) {
return B.first(n)
}
upto *= 2
}
}
with (50) {|n|
say "The first #{n} Blum integers:"
blum_first(n).slices(10).each { .map{ "%4s" % _ }.join.say }
}
say ''
for n in (26828, 1e5, 2e5, 3e5, 4e5) {
var B = blum_first(n)
say "#{n.commify}th Blum integer: #{B.last}"
if (n == 4e5) {
say ''
for k in (1,3,7,9) {
var T = B.grep { .is_congruent(k, 10) }
say "#{k}: #{'%6s' % T.len} (#{T.len / B.len * 100}%)"
}
}
}</syntaxhighlight>
{{out}}
<pre>
The 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
26,828th Blum integer: 524273
100,000th Blum integer: 2075217
200,000th Blum integer: 4275533
300,000th Blum integer: 6521629
400,000th Blum integer: 8802377
1: 100005 (25.00125%)
3: 100067 (25.01675%)
7: 99989 (24.99725%)
9: 99939 (24.98475%)
</pre>
Line 1,703 ⟶ 2,148:
{{libheader|Wren-math}}
{{libheader|Wren-fmt}}
<syntaxhighlight lang="
import "./fmt" for Fmt
| |||