Blum integer

You are encouraged to solve this task according to the task description, using any language you may know.
- Definition
A positive integer n is a Blum integer if n = p x q is a semi-prime for which p and q are distinct primes congruent to 3 mod 4. In other words, p and q must be of the form 4t + 3 where t is some non-negative integer.
- Example
21 is a Blum integer because it has two prime factors: 3 (= 4 x 0 + 3) and 7 (= 4 x 1 + 3).
- Task
Find and show on this page the first 50 Blum integers.
Also show the 26,828th.
- Stretch
Find and show the 100,000th, 200,000th, 300,000th and 400,000th Blum integers.
For the first 400,000 Blum integers, show the percentage distribution by final decimal digit (to 3 decimal places). Clearly, such integers can only end in 1, 3, 7 or 9.
- Related task
- References
- Wikipedia article Blum integer
- OEIS sequence A016105: Blum integers
ALGOL 68
If running this with ALGOL 68G, you will need to specify a larger heap size with -heap 256M
on the command line.
Builds tables of the unique prime factor counts and the last prime factor to handle the stretch goal, uses prime factorisation to find the first 50 Blum integers.
Note that Blum integers must be 1 modulo 4 and if one prime factor is 3 modulo 4, the other must also be.
BEGIN # find Blum integers, semi-primes where both factors are 3 mod 4 #
# and the factors are distinct #
INT max number = 10 000 000; # maximum possible Blum we will consider #
[ 1 : max number ]INT upfc; # table of unique prime factor counts #
[ 1 : max number ]INT lpf; # table of last prime factors #
FOR i TO UPB upfc DO upfc[ i ] := 0; lpf[ i ] := 1 OD;
FOR i FROM 2 TO UPB upfc OVER 2 DO
IF upfc[ i ] = 0 THEN
FOR j FROM i + i BY i TO UPB upfc DO
upfc[ j ] +:= 1;
lpf[ j ] := i
OD
FI
OD;
# returns TRUE if n is a Blum integer, FALSE otherwise #
PROC is blum = ( INT n )BOOL:
IF n < 21 THEN FALSE # the lowest primes modulo 4 that are 3 are #
# 3 & 7, so 21 is the lowest possible number #
ELIF NOT ODD n THEN FALSE
ELSE
INT v := n;
INT f count := 0;
INT f := 3;
WHILE f <= v AND f count < 3 DO
IF v MOD f = 0 THEN
IF f MOD 4 /= 3 THEN
f count := 9 # the prime factor mod 4 isn't 3 #
ELSE
v OVERAB f;
f count +:= 1
FI;
IF v MOD f = 0 THEN
f count := 9 # f is not a distinct factor of n #
FI
FI;
f +:= 2
OD;
f count = 2
FI # is blum # ;
# show various Blum integers #
print( ( "The first 50 Blum integers:", newline ) );
INT b count := 0;
[ 0 : 9 ]INT last count := []INT( 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 )[ AT 0 ];
FOR i FROM 21 BY 4 WHILE b count < 400 000 DO
# Blum integers must be 1 modulo 4 #
IF b count < 50 THEN
IF is blum( i ) THEN
b count +:= 1;
last count[ i MOD 10 ] +:= 1;
print( ( whole( i, -4 ) ) );
IF b count MOD 10 = 0 THEN print( ( newline ) ) FI
FI
ELIF upfc[ i ] = 2 THEN
# two prime factors - could be a Blum integer #
IF lpf[ i ] MOD 4 = 3 THEN
# the last prime factor mod 4 is three #
IF upfc[ i OVER lpf[ i ] ] = 0 THEN
# and the other prime factor is unique (e.g. not 3^3) #
b count +:= 1;
last count[ i MOD 10 ] +:= 1;
IF b count = 26 828
OR b count = 100 000
OR b count = 200 000
OR b count = 300 000
OR b count = 400 000
THEN
print( ( "The ", whole( b count, -6 )
, "th Blum integer is ", whole( i, -11 )
, newline
)
)
FI
FI
FI
FI
OD;
# show some statistics of the last digits #
print( ( "For Blum integers up to ", whole( b count, 0 ), ":", newline ) );
FOR i FROM LWB last count TO UPB last count DO
IF last count[ i ] /= 0 THEN
print( ( " ", fixed( ( last count[ i ] * 100 ) / b count, -8, 3 )
, "% end with ", whole( i, 0 )
, newline
)
)
FI
OD
END
- Output:
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 For Blum integers up to 400000: 25.001% end with 1 25.017% end with 3 24.997% end with 7 24.985% end with 9
BASIC
BASIC256
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
FreeBASIC
#include "isprime.bas"
Type PrimeHelper
inc(7) As Integer
idx As Integer
End Type
Function initPrimeHelper() As PrimeHelper
Dim helper As PrimeHelper
helper.inc(0) = 4 : helper.inc(1) = 2 : helper.inc(2) = 4
helper.inc(3) = 2 : helper.inc(4) = 4 : helper.inc(5) = 6
helper.inc(6) = 2 : helper.inc(7) = 6
helper.idx = 0
Return helper
End Function
Function firstPrimeFactor(n As Integer) As Integer
If n = 1 Then Return 1
If n Mod 3 = 0 Then Return 3
If n Mod 5 = 0 Then Return 5
Dim helper As PrimeHelper = initPrimeHelper()
Dim k As Integer = 7
While k * k <= n
If n Mod k = 0 Then Return k
k += helper.inc(helper.idx)
helper.idx = (helper.idx + 1) Mod 8
Wend
Return n
End Function
Sub main()
Dim As Integer blum(49), counts(9)
Dim As Integer bc = 0, i = 1, p, q
Dim As Integer j
Dim As Double t0 = Timer
Do
p = firstPrimeFactor(i)
If p Mod 4 = 3 Then
q = i \ p
If q <> p Andalso q Mod 4 = 3 Andalso isPrime(q) Then
If bc < 50 Then blum(bc) = i
counts(i Mod 10) += 1
bc += 1
If bc = 50 Then
Print "First 50 Blum integers:"
For j = 0 To 49
Print Using "####"; blum(j);
If (j + 1) Mod 10 = 0 Then Print
Next
Print
Elseif bc = 26828 Orelse bc Mod 100000 = 0 Then
Print Using "The ###,###th Blum integer is: #,###,###"; bc; i
If bc = 400000 Then
Print !"\n% distribution of the first 400,000 Blum integers:"
For j = 1 To 9 Step 2
If j <> 5 Then Print Using " ##.###% end in #"; (counts(j)/4000); j
Next
Exit Do
End If
End If
End If
End If
i += Iif(i Mod 5 = 3, 4, 2)
Loop
Print Chr(10); Timer - t0; " sec."
End Sub
main()
Sleep
- Output:
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 26,828th Blum integer is: 524,273 The 100,000th Blum integer is: 2,075,217 The 200,000th Blum integer is: 4,275,533 The 300,000th Blum integer is: 6,521,629 The 400,000th Blum integer is: 8,802,377 % 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 6.020220100006554 sec.
Gambas
Use "isprime.bas"
Private inc As Integer[] = [4, 2, 4, 2, 4, 6, 2, 6]
Private Function FirstPrimeFactor(n As Long) As Long
If n = 1 Then Return 1
If n Mod 3 = 0 Then Return 3
If n Mod 5 = 0 Then Return 5
Dim k As Long = 7
Dim i As Integer = 0
While k * k <= n
If n Mod k = 0 Then Return k
k += inc[i]
i = (i + 1) Mod 8
Wend
Return n
End
Public Sub Main()
Dim blum As New Long[50]
Dim counts As New Collection
counts[1] = 0
counts[3] = 0
counts[7] = 0
counts[9] = 0
Dim bc As Long = 0, i As Long = 1
While True
Dim p As Long = FirstPrimeFactor(i)
If p Mod 4 = 3 Then
Dim q As Long = i \ p
If q <> p And q Mod 4 = 3 And IsPrime(q) Then
If bc < 50 Then blum[bc] = i
counts[i Mod 10] += 1
bc += 1
If bc = 50 Then
Print "First 50 Blum integers:"
For j As Integer = 0 To 49
Print Format(blum[j], "####"); " ";
If (j + 1) Mod 10 = 0 Then Print
Next
Print
Else If bc = 26828 Or bc Mod 100000 = 0 Then
Print "The "; Format(bc, " ###,###"); "th Blum integer is: "; Format(i, " #,###,###")
If bc = 400000 Then
Print Chr(10); "% distribution of the first 400,000 Blum integers:"
For Each j As Integer In [1, 3, 7, 9]
Print Format(counts[j] / 4000, " ##.###"); "% end in "; j
Next
Return
Endif
Endif
Endif
Endif
i += If(i Mod 5 = 3, 4, 2)
Wend
End
- Output:
Similar to FreeBASIC entry.
OxygenBasic
#include "isprime.bas"
uses console
dim inc(7) as integer
inc(0) = 4: inc(1) = 2: inc(2) = 4
inc(3) = 2: inc(4) = 4: inc(5) = 6
inc(6) = 2: inc(7) = 6
function firstPrimeFactor(n as long) as long
if n = 1 then return 1
if n mod 3 = 0 then return 3
if n mod 5 = 0 then return 5
long k = 7
int idx = 0
while k * k <= n
if mod(n, k) = 0 then return k
k += inc(idx)
idx = mod((idx + 1), 8)
wend
return n
end function
sub main()
dim as long blum(49), counts(9)
long bc = 0, i = 1, pct, p, q
int j
string s, t
do
p = firstPrimeFactor(i)
if p mod 4 = 3 then
q = i \ p
if q <> p and q mod 4 = 3 and isPrime(q) then
if bc < 50 then blum(bc) = i
counts(i mod 10) = counts(i mod 10) + 1
bc += 1
if bc = 50 then
printl "First 50 Blum integers:"
for j = 0 to 49
s = str(blum(j))
while len(s) < 4: s = " " + s: wend
print s;
if mod((j + 1), 10) = 0 then printl
next
printl
elseif bc = 26828 or bc mod 100000 = 0 then
s = str(bc)
while len(s) < 6: s = " " + s: wend
t = str(i)
while len(t) < 7: t = " " + t: wend
printl "The " + s + "th Blum integer is: " + t
if bc = 400000 then
printl cr "% distribution of the first 400,000 Blum integers:"
for j = 1 to 9 step 2
if j <> 5 then
pct = counts(j)/4000
s = str(pct)
while len(s) < 5: s = " " + s: wend
printl s + "% end in " + str(j)
end if
next
end
end if
end if
end if
end if
if mod(i, 5) = 3 then i += 4 else i += 2
end do
end sub
main()
printl cr "Enter ..."
waitkey
PureBasic
XIncludeFile "isprime.pb"
Structure PrimeHelper
inc.i[8]
index.i
EndStructure
Procedure.i firstPrimeFactor(n.q)
If n = 1 : ProcedureReturn 1 : EndIf
If n % 3 = 0 : ProcedureReturn 3 : EndIf
If n % 5 = 0 : ProcedureReturn 5 : EndIf
Define helper.PrimeHelper
helper\inc[0] = 4 : helper\inc[1] = 2 : helper\inc[2] = 4
helper\inc[3] = 2 : helper\inc[4] = 4 : helper\inc[5] = 6
helper\inc[6] = 2 : helper\inc[7] = 6
Define k.q = 7
While k * k <= n
If n % k = 0 : ProcedureReturn k : EndIf
k + helper\inc[helper\index]
helper\index = (helper\index + 1) % 8
Wend
ProcedureReturn n
EndProcedure
OpenConsole()
Define Dim blum.q(49)
Define.q bc = 0, i = 1
Define Dim counts.q(9)
Repeat
Define p.q = firstPrimeFactor(i)
If p % 4 = 3
Define q.q = i / p
If q <> p And q % 4 = 3 And isPrime(q)
If bc < 50 : blum(bc) = i : EndIf
counts(i % 10) + 1
bc + 1
If bc = 50
PrintN("First 50 Blum integers:")
For j = 0 To 49
Print(" " + RSet(Str(blum(j)), 3))
If (j + 1) % 10 = 0 : PrintN("") : EndIf
Next
PrintN("")
ElseIf bc = 26828 Or bc % 100000 = 0
PrintN("The " + RSet(Str(bc), 6) + "th Blum integer is: " + RSet(Str(i), 7))
If bc = 400000
PrintN(#CRLF$ + "% distribution of the first 400,000 Blum integers:")
For j = 1 To 9 Step 2
If j <> 5
PrintN(RSet(StrF(counts(j)/4000, 3), 5) + "% end in " + Str(j))
EndIf
Next
Break
EndIf
EndIf
EndIf
EndIf
If i % 5 = 3
i + 4
Else
i + 2
EndIf
ForEver
PrintN(#CRLF$ + "Press ENTER to exit"): Input()
- Output:
Similar to FreeBASIC entry.
QB64
Dim Shared inc(7) As Integer
inc(0) = 4: inc(1) = 2: inc(2) = 4
inc(3) = 2: inc(4) = 4: inc(5) = 6
inc(6) = 2: inc(7) = 6
Dim blum(49) As Long
Dim counts(9) As Long
Dim bc As Long, i As Long, p As Long, q As Long, j As Integer
i = 1
Do
p = firstPrimeFactor(i)
If p Mod 4 = 3 Then
q = i \ p
If q <> p And q Mod 4 = 3 And isPrime(q) Then
If bc < 50 Then blum(bc) = i
counts(i Mod 10) = counts(i Mod 10) + 1
bc = bc + 1
If bc = 50 Then
Print "First 50 Blum integers:"
For j = 0 To 49
Print Using "####"; blum(j);
If (j + 1) Mod 10 = 0 Then Print
Next
Print
ElseIf bc = 26828 Or bc Mod 100000 = 0 Then
Print Using "The ###,###th Blum integer is: #,###,###"; bc; i
If bc = 400000 Then
Print Chr$(10); "% distribution of the first 400,000 Blum integers:"
For j = 1 To 9 Step 2
If j <> 5 Then
Print Using " ##.###% end in #"; (counts(j%) / 4000); j%
End If
Next
End
End If
End If
End If
End If
If i Mod 5 = 3 Then i = i + 4 Else i = i + 2
Loop
End
Function firstPrimeFactor& (n As Long)
Dim k As Long, idx As Integer
If n = 1 Then firstPrimeFactor& = 1: Exit Function
If n Mod 3 = 0 Then firstPrimeFactor& = 3: Exit Function
If n Mod 5 = 0 Then firstPrimeFactor& = 5: Exit Function
k = 7: idx = 0
Do While k * k <= n
If n Mod k = 0 Then
firstPrimeFactor& = k
Exit Function
End If
k = k + inc(idx)
idx = (idx + 1) Mod 8
Loop
firstPrimeFactor& = n
End Function
Function isPrime% (n As Long)
Dim i As Long
If n <= 1 Then Exit Function
If n <= 3 Then isPrime% = 1: Exit Function
If n Mod 2 = 0 Or n Mod 3 = 0 Then Exit Function
i = 5
While i * i <= n
If n Mod i = 0 Or n Mod (i + 2) = 0 Then Exit Function
i = i + 6
Wend
isPrime% = 1
End Function
- Output:
Same as FreeBASIC entry.
C
#include <stdio.h>
#include <stdbool.h>
#include <locale.h>
int inc[8] = {4, 2, 4, 2, 4, 6, 2, 6};
bool isPrime(int n) {
if (n < 2) return false;
if (n%2 == 0) return n == 2;
if (n%3 == 0) return n == 3;
int d = 5;
while (d*d <= n) {
if (n%d == 0) return false;
d += 2;
if (n%d == 0) return false;
d += 4;
}
return true;
}
// Assumes n is odd.
int firstPrimeFactor(int n) {
if (n == 1) return 1;
if (!(n%3)) return 3;
if (!(n%5)) return 5;
for (int k = 7, i = 0; k*k <= n; ) {
if (!(n%k)) {
return k;
} else {
k += inc[i];
i = (i + 1) % 8;
}
}
return n;
}
int main() {
int i = 1, j, bc = 0, p, q;
int blum[50], counts[4] = {0}, digits[4] = {1, 3, 5, 7};
setlocale(LC_NUMERIC, "");
while (true) {
p = firstPrimeFactor(i);
if (p % 4 == 3) {
q = i / p;
if (q != p && q % 4 == 3 && isPrime(q)) {
if (bc < 50) blum[bc] = i;
++counts[i % 10 / 3];
++bc;
if (bc == 50) {
printf("First 50 Blum integers:\n");
for (j = 0; j < 50; ++j) {
printf("%3d ", blum[j]);
if (!((j+1) % 10)) printf("\n");
}
printf("\n");
} else if (bc == 26828 || !(bc % 100000)) {
printf("The %'7dth Blum integer is: %'9d\n", bc, i);
if (bc == 400000) {
printf("\n%% distribution of the first 400,000 Blum integers:\n");
for (j = 0; j < 4; ++j) {
printf(" %6.3f%% end in %d\n", counts[j]/4000.0, digits[j]);
}
break;
}
}
}
}
i += (i % 5 == 3) ? 4 : 2;
}
return 0;
}
- Output:
Same as Wren example.
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;
}
}
- Output:
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
C++
#include <algorithm>
#include <cstdint>
#include <iomanip>
#include <iostream>
#include <vector>
bool is_prime_type_3(const uint32_t& number) {
if ( number < 2 ) return false;
if ( number % 2 == 0 ) return false;
if ( number % 3 == 0 ) return number == 3;
for ( uint32_t divisor = 5; divisor * divisor <= number; divisor += 2 ) {
if ( number % divisor == 0 ) { return false; }
}
return number % 4 == 3;
}
uint32_t least_prime_factor(const uint32_t& number) {
if ( number == 1 ) { return 1; }
if ( number % 3 == 0 ) { return 3; }
if ( number % 5 == 0 ) { return 5; }
for ( uint32_t divisor = 7; divisor * divisor <= number; divisor += 2 ) {
if ( number % divisor == 0 ) { return divisor; }
}
return number;
}
int main() {
uint32_t blums[50];
uint32_t blum_count = 0;
uint32_t last_digit_counts[10] = {};
uint32_t number = 1;
while ( blum_count < 400'000 ) {
const uint32_t prime = least_prime_factor(number);
if ( prime % 4 == 3 ) {
const uint32_t quotient = number / prime;
if ( quotient != prime && is_prime_type_3(quotient) ) {
if ( blum_count < 50 ) {
blums[blum_count] = number;
}
last_digit_counts[number % 10] += 1;
blum_count += 1;
if ( blum_count == 50 ) {
std::cout << "The first 50 Blum integers:" << std::endl;
for ( uint32_t i = 0; i < 50; ++i ) {
std::cout << std::setw(3) << blums[i] << ( ( i % 10 == 9 ) ? "\n" : " " );
}
std::cout << std::endl;
} else if ( blum_count == 26'828 || blum_count % 100'000 == 0 ) {
std::cout << "The " << std::setw(6) << blum_count << "th Blum integer is: "
<< std::setw(7) << number << std::endl;
if ( blum_count == 400'000 ) {
std::cout << "\n" << "Percent distribution of the first 400000 Blum integers:"
<< std::endl;
for ( const int32_t& i : { 1, 3, 7, 9 } ) {
std::cout << " " << std::setw(6) << std::setprecision(5)
<< (double) last_digit_counts[i] / 4'000 << "% end in " << i << std::endl;
}
}
}
}
}
number += ( number % 5 == 3 ) ? 4 : 2;
}
}
- Output:
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
Common Lisp
Simple implementation without any performance tuning; execute (main) from REPL.
(defun is-factor (x y)
(zerop (mod x y)))
(defun is-congruent (num a n)
(= (mod num n) a))
(defun is-prime (n)
(cond ((< n 4) (or (= n 2) (= n 3)))
((or (zerop (mod n 2)) (zerop (mod n 3))) nil)
(t (loop for i from 5 to (floor (sqrt n)) by 6
never (or (is-factor n i)
(is-factor n (+ i 2)))))))
(defun first-factor (n)
(loop for i from 2 to (floor (sqrt n))
thereis (if (is-factor n i) i nil)))
(defun is-blum (n)
(let ((factor (first-factor n)))
(and factor
(is-congruent factor 3 4)
(/= factor (/ n factor))
(is-congruent (/ n factor) 3 4)
(is-prime factor)
(is-prime (/ n factor)))))
(defun main ()
(let ((n 0)
(i 1)
(blums '())
(counts '()))
(dotimes (i 10) (push (cons i 0) counts))
(loop while (<= n 400000)
do (setf i (+ i 2))
(when (is-blum i)
(incf n)
(incf (cdr (assoc (mod i 10) counts)))
(when (<= n 50)
(push i blums)
(when (= n 50)
(format t "First 50 Blum integers:~%")
(format t "~{~5d~5d~5d~5d~5d~5d~5d~5d~5d~5d~%~}~%" (reverse blums))))
(if (find n '(26828 100000 200000 300000 400000))
(format t "The ~7:dth Blum integer is: ~9:d~%" n i))))
(format t "~%% distribution of the first 400,000 Blum integers:~%")
(dolist (x '(1 3 7 9))
(format t "~3$% end in ~a~%" (/ (cdr (assoc x counts)) 4000) x))))
- Output:
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 26,828th Blum integer is: 524,273 The 100,000th Blum integer is: 2,075,217 The 200,000th Blum integer is: 4,275,533 The 300,000th Blum integer is: 6,521,629 The 400,000th Blum integer is: 8,802,377 % 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
EasyLang
fastfunc semiprim n .
d = 3
while d * d <= n
while n mod d = 0
if c = 2
return 0
.
n /= d
c += 1
.
d += 2
.
if c = 1
return n
.
.
print "The first 50 Blum integers:"
n = 3
numfmt 0 4
repeat
prim1 = semiprim n
if prim1 <> 0
if prim1 mod 4 = 3
prim2 = n div prim1
if prim2 <> prim1 and prim2 mod 4 = 3
c += 1
if c <= 50
write n
if c mod 10 = 0 ; print "" ; .
.
.
.
.
until c >= 26828
n += 2
.
print ""
print "The 26828th Blum integer is: " & n
Fortran
program BlumInteger
use, intrinsic :: iso_fortran_env, only: int32, int64
implicit none
integer(int32), parameter :: LIMIT = 10*1000*1000
integer(int32), allocatable :: BlumPrimes(:)
integer(int32), allocatable :: BlumPrimes2(:)
logical :: BlumField(0:LIMIT)
integer(int32) :: EndDigit(0:9)
integer(int64) :: k
integer(int32) :: n, idx, j, P4n3Cnt,xx,yy
call system_clock(count=xx)
call Sieve4n_3_Primes(LIMIT, BlumPrimes2)
! allocate(blumprimes(0:size(BlumPrimes2)-1))
! The blumprimes2 array is allocated in the subroutine as a zero based array
! but, the main program doesn't know that and assumes it's 1 based so we allocate
! a zero based array then use move_alloc to correctly resize it a transfer the data
allocate(blumprimes(0:1))
call move_alloc(blumprimes2,blumprimes)
P4n3Cnt = size(BlumPrimes) -1
print *, 'There are ', CommaUint(int(P4n3Cnt, int64)), ' needed primes 4*n+3 to Limit ', CommaUint(int(LIMIT, int64))
P4n3Cnt = P4n3Cnt - 1
print *
! Generate Blum-Integers
BlumField = .false.
do idx = 0, P4n3Cnt
n = BlumPrimes(idx)
do j = idx+1, P4n3Cnt
k = int(n, int64) * int(BlumPrimes(j), int64)
if (k > LIMIT) exit
BlumField(k) = .true.
end do
end do
call system_clock(count=yy)
print *, 'First 50 Blum-Integers '
idx = 0
j = 0
do
do while (idx < LIMIT .and. .not. BlumField(idx))
idx = idx + 1
end do
if (idx == LIMIT) exit
if (mod(j, 10) == 0 .and. j /= 0) print *
write(*, '(I5)', advance='no') idx
j = j + 1
idx = idx + 1
if (j >= 50) exit
end do
print '(//)'
print *, ' relative occurence of digit'
print *, ' n.th |BlumInteger|Digit: 1 3 7 9'
idx = 0
j = 0
n = 0
k = 26828
EndDigit = 0
do
do while (idx < LIMIT .and. .not. BlumField(idx))
idx = idx + 1
end do
if (idx == LIMIT) exit
! Count last decimal digit
EndDigit(mod(idx, 10)) = EndDigit(mod(idx, 10)) + 1
j = j + 1
if (j == k) then
write(*, '(A10,A1,A11,A1)', advance='no') CommaUint(int(j, int64)), '|', CommaUint(int(idx, int64)), '|'
write(*, '(F7.3,A4)', advance='no') real(EndDigit(1))/j*100, '% |'
write(*, '(F7.3,A4)', advance='no') real(EndDigit(3))/j*100, '% |'
write(*, '(F7.3,A4)', advance='no') real(EndDigit(7))/j*100, '% |'
write(*, '(F7.3,A2)') real(EndDigit(9))/j*100, '%'
if (k < 100000) then
k = 100000
else
k = k + 100000
end if
end if
idx = idx + 1
if (j >= 400000) exit
end do
print '(/,a,f8.6,1x,a)', 'Elapsed time = ',(yy-xx)/1000.0,'seconds'
contains
subroutine Sieve4n_3_Primes(Limit, P4n3)
use iso_fortran_env
integer(int32), intent(in) :: Limit
integer(int32), allocatable, intent(out) :: P4n3(:)
integer(kind=1), allocatable :: sieve(:)
integer(int32) :: BlPrCnt, idx, n, j, sieve_size
sieve_size = (Limit / 3 - 3) / 4 + 1
allocate(sieve(0:sieve_size-1))
allocate(P4n3(0:sieve_size-1))
sieve = 0
BlPrCnt = 0
idx = 0
do
if (sieve(idx) == 0) then
n = idx*4 + 3
P4n3(BlPrCnt) = n
BlPrCnt = BlPrCnt + 1
j = idx + n
if (j > ubound(sieve, 1)) exit
do while (j <= ubound(sieve, 1))
sieve(j) = 1
j = j + n
end do
end if
idx = idx + 1
if (idx > ubound(sieve, 1)) exit
end do
! Collect the rest
do idx = idx, ubound(sieve, 1)
if (sieve(idx) == 0) then
P4n3(BlPrCnt) = idx*4 + 3
BlPrCnt = BlPrCnt + 1
end if
end do
P4n3 = P4n3(0:BlPrCnt-1)
end subroutine Sieve4n_3_Primes
function CommaUint(n) result(res)
integer, parameter :: sizer = 30
integer(int64), intent(in) :: n
character(:), allocatable :: res
character(len=sizer) :: temp
integer :: fromIdx, toIdx, i
character :: pRes(sizer)
write(temp, '(I0)') n
fromIdx = len_trim(temp)
toIdx = fromIdx - 1
if (toIdx < 3) then
res = temp(1:fromIdx)
return
end if
allocate(res, mold=repeat(' ',sizer))
toIdx = 4*(toIdx / 3) + mod(toIdx, 3) + 1
pRes = ' '
do i = 1, fromIdx
pRes(toIdx) = temp(fromIdx-i+1:fromIdx-i+1)
toIdx = toIdx - 1
if (mod(i, 3) == 0 .and. i /= fromIdx) then
pRes(toIdx) = ','
toIdx = toIdx - 1
end if
end do
do i = 1,sizer ! Go from character array to string
res(I:I) = pRes(i)
end do
!
res = trim(adjustl(Res))
end function CommaUint
end program BlumInteger
- Output:
There are 119,644 needed primes 4*n+3 to Limit 10,000,000 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 relative occurence of digit n.th |BlumInteger|Digit: 1 3 7 9 26,828| 524,273| 24.933% | 25.071% | 25.030% | 24.966 % 100,000| 2,075,217| 24.973% | 25.026% | 25.005% | 24.996 % 200,000| 4,275,533| 24.990% | 24.986% | 25.032% | 24.992 % 300,000| 6,521,629| 24.982% | 25.014% | 25.033% | 24.970 % 400,000| 8,802,377| 25.001% | 25.017% | 24.997% | 24.985 % Elapsed time = 0.032000 seconds
FutureBasic
begin globals
UInt32 gPrime1
end globals
local fn IsSemiprime( n as UInt32 ) as BOOL
BOOL result = NO
UInt32 d = 3, c = 0
while d * d <= n
while ( n % d = 0 )
if ( c == 2 ) then return NO
n /= d
c += 1
wend
d += 2
wend
gPrime1 = n
result = ( c == 1 )
end fn = result
void local fn BlumIntegers
UInt32 prime2 = 0, n = 3, c = 0
print @"The first 50 Blum integers:"
while ( 1 )
if ( fn IsSemiprime( n ) )
if ( gPrime1 % 4 == 3 )
prime2 = n / gPrime1
if ( prime2 != gPrime1 ) && ( prime2 % 4 == 3 )
c++
if ( c <= 50 )
printf @"%4d\b",n
if ( c % 10 == 0 ) then print
end if
if ( c >= 26828 )
print : print @"The 26828th Blum integer is: "; n
break
end if
end if
end if
end if
n += 2
wend
end fn
fn BlumIntegers
HandleEvents
- Output:
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
Go
package main
import (
"fmt"
"rcu"
)
var inc = []int{4, 2, 4, 2, 4, 6, 2, 6}
// Assumes n is odd.
func firstPrimeFactor(n int) int {
if n == 1 {
return 1
}
if n%3 == 0 {
return 3
}
if n%5 == 0 {
return 5
}
for k, i := 7, 0; k*k <= n; {
if n%k == 0 {
return k
} else {
k += inc[i]
i = (i + 1) % 8
}
}
return n
}
func main() {
blum := make([]int, 50)
bc := 0
counts := make([]int, 4)
i := 1
for {
p := firstPrimeFactor(i)
if p%4 == 3 {
q := i / p
if q != p && q%4 == 3 && rcu.IsPrime(q) {
if bc < 50 {
blum[bc] = i
}
counts[i%10/3]++
bc++
if bc == 50 {
fmt.Println("First 50 Blum integers:")
rcu.PrintTable(blum, 10, 3, false)
fmt.Println()
} else if bc == 26828 || bc%100000 == 0 {
fmt.Printf("The %7sth Blum integer is: %9s\n", rcu.Commatize(bc), rcu.Commatize(i))
if bc == 400000 {
fmt.Println("\n% distribution of the first 400,000 Blum integers:")
digits := []int{1, 3, 7, 9}
for j := 0; j < 4; j++ {
fmt.Printf(" %6.3f%% end in %d\n", float64(counts[j])/4000, digits[j])
}
return
}
}
}
}
if i%5 == 3 {
i += 4
} else {
i += 2
}
}
}
- Output:
Same as Wren example.
Haskell
Works with GHC 9.2.8
and arithmoi-0.13.0.0.
{-# LANGUAGE BangPatterns #-}
{-# LANGUAGE DeriveFunctor #-}
{-# LANGUAGE ScopedTypeVariables #-}
module Rosetta.BlumInteger
( Stream (..)
, Stats (..)
, toList
, blumIntegers
, countLastDigits
) where
import GHC.Natural (Natural)
import Math.NumberTheory.Primes (Prime (..), nextPrime)
-- | A stream is an infinite list.
data Stream a = a :> Stream a
deriving Functor
-- | Converts a stream to the corresponding infinite list.
toList :: Stream a -> [a]
toList (x :> xs) = x : toList xs
unsafeFromList :: [a] -> Stream a
unsafeFromList = foldr (:>) $ error "fromList: finite list"
primes3mod4 :: Stream (Prime Natural)
primes3mod4 = unsafeFromList [nextPrime 3, nextPrime 7 ..]
-- Assume:
-- * All numbers in all the streams are distinct.
-- * Each stream is sorted.
-- * In the stream of streams, the first element of each stream is less than the first element of the next stream.
sortStreams :: forall a. Ord a => Stream (Stream a) -> Stream a
sortStreams ((x :> xs) :> xss) = x :> sortStreams (insert xs xss)
where
insert :: Stream a -> Stream (Stream a) -> Stream (Stream a)
insert ys@(y :> _) zss@(zs@(z :> _) :> zss')
| y < z = ys :> zss
| otherwise = zs :> insert ys zss'
-- | The
blumIntegers :: Stream Natural
blumIntegers = sortStreams $ go $ unPrime <$> primes3mod4
where
go :: Stream Natural -> Stream (Stream Natural)
go (p :> ps) = ((p *) <$> ps) :> go ps
data Stats a = Stats
{ s1 :: !a
, s3 :: !a
, s7 :: !a
, s9 :: !a
} deriving (Show, Eq, Ord, Functor)
lastDigit :: Natural -> Int
lastDigit n = fromIntegral $ n `mod` 10
updateCount :: Stats Int -> Natural -> Stats Int
updateCount !dc n = case lastDigit n of
1 -> dc { s1 = s1 dc + 1 }
3 -> dc { s3 = s3 dc + 1 }
7 -> dc { s7 = s7 dc + 1 }
9 -> dc { s9 = s9 dc + 1 }
_ -> error "updateCount: impossible"
countLastDigits :: forall a. Fractional a => Int -> Stream Natural -> Stats a
countLastDigits n = fmap f . go Stats { s1 = 0, s3 = 0, s7 = 0, s9 = 0 } n
where
go :: Stats Int -> Int -> Stream Natural -> Stats Int
go !dc 0 _ = dc
go !dc m (x :> xs) = go (updateCount dc x) (m - 1) xs
f :: Int -> a
f m = fromIntegral m / fromIntegral n
{-# LANGUAGE NumericUnderscores #-}
{-# LANGUAGE TypeApplications #-}
module Main
( main
) where
import Control.Monad (forM_)
import Text.Printf (printf)
import Numeric.Natural (Natural)
import Rosetta.BlumInteger
main :: IO ()
main = do
let xs = toList blumIntegers
printf "The first 50 Blum integers are:\n\n"
forM_ (take 50 xs) $ \x -> do
printf "%3d\n" x
printf "\n"
nth xs 26_828
forM_ [100_000, 200_000 .. 400_000] $ nth xs
printf "\n"
printf "Distribution by final digit for the first 400000 Blum integers:\n\n"
let Stats r1 r3 r7 r9 = countLastDigits @Double 400_000 blumIntegers
forM_ [(1 :: Int, r1), (3, r3), (7, r7), (9, r9)] $ \(d, r) ->
printf "%d: %6.3f%%\n" d $ r * 100
printf "\n"
where
nth :: [Natural] -> Int -> IO ()
nth xs n = printf "The %6dth Blum integer is %8d.\n" n $ xs !! (n - 1)
- Output:
The first 50 Blum integers are: 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. Distribution by final digit for the first 400000 Blum integers: 1: 25.001% 3: 25.017% 7: 24.997% 9: 24.985%
J
Implementation:
isblum=: {{
ab=. q: y
if. 2= #ab do.
if. </ab do.
*./3=4|ab
else. 0 end.
else. 0 end.
}}"0
blumseq=: {{
r=. (#~ isblum) }.i.b=. 1e4
while. y>#r do.
r=. r, (#~ isblum) b+i.1e4
b=. b+1e4
end.
y{.r
}}
Task examples:
5 10$blumseq 50
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
{: blumseq 26828
524273
Stretch:
B=: blumseq 4e5
{:1e5{.B
2075217
{:2e5{.B
4275533
{:3e5{.B
6521629
{:4e5{.B
8802377
{:B
8802377
(~.,. 4e3 %~ #/.~) 10|B
1 25.0012
3 25.0167
7 24.9973
9 24.9847
Java
import java.util.HashMap;
import java.util.Map;
public final class BlumInteger {
public static void main(String[] aArgs) {
int[] blums = new int[50];
int blumCount = 0;
Map<Integer, Integer> lastDigitCounts = new HashMap<Integer, Integer>();
int number = 1;
while ( blumCount < 400_000 ) {
final int prime = leastPrimeFactor(number);
if ( prime % 4 == 3 ) {
final int quotient = number / prime;
if ( quotient != prime && isPrimeType3(quotient) ) {
if ( blumCount < 50 ) {
blums[blumCount] = number;
}
lastDigitCounts.merge(number % 10, 1, Integer::sum);
blumCount += 1;
if ( blumCount == 50 ) {
System.out.println("The first 50 Blum integers:");
for ( int i = 0; i < 50; i++ ) {
System.out.print(String.format("%3d", blums[i]));
System.out.print(( i % 10 == 9 ) ? System.lineSeparator() : " ");
}
System.out.println();
} else if ( blumCount == 26_828 || blumCount % 100_000 == 0 ) {
System.out.println(String.format("%s%6d%s%7d",
"The ", blumCount, "th Blum integer is: ", number));
if ( blumCount == 400_000 ) {
System.out.println();
System.out.println("Percent distribution of the first 400000 Blum integers:");
for ( int key : lastDigitCounts.keySet() ) {
System.out.println(String.format("%s%6.3f%s%d",
" ", (double) lastDigitCounts.get(key) / 4_000, "% end in ", key));
}
}
}
}
}
number += ( number % 5 == 3 ) ? 4 : 2;
}
}
private static boolean isPrimeType3(int aNumber) {
if ( aNumber < 2 ) { return false; }
if ( aNumber % 2 == 0 ) { return false; }
if ( aNumber % 3 == 0 ) { return aNumber == 3; }
for ( int divisor = 5; divisor * divisor <= aNumber; divisor += 2 ) {
if ( aNumber % divisor == 0 ) { return false; }
}
return aNumber % 4 == 3;
}
private static int leastPrimeFactor(int aNumber) {
if ( aNumber == 1 ) { return 1; }
if ( aNumber % 3 == 0 ) { return 3; }
if ( aNumber % 5 == 0 ) { return 5; }
for ( int divisor = 7; divisor * divisor <= aNumber; divisor += 2 ) {
if ( aNumber % divisor == 0 ) { return divisor; }
}
return aNumber;
}
}
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
jq
Adapted from Wren
### Generic utilities
def lpad($len): tostring | ($len - length) as $l | (" " * $l) + .;
# tabular print
def tprint(columns; wide):
reduce _nwise(columns) as $row ("";
. + ($row|map(lpad(wide)) | join(" ")) + "\n" );
### Primes
def is_prime:
. as $n
| if ($n < 2) then false
elif ($n % 2 == 0) then $n == 2
elif ($n % 3 == 0) then $n == 3
elif ($n % 5 == 0) then $n == 5
elif ($n % 7 == 0) then $n == 7
elif ($n % 11 == 0) then $n == 11
elif ($n % 13 == 0) then $n == 13
elif ($n % 17 == 0) then $n == 17
elif ($n % 19 == 0) then $n == 19
else sqrt as $s
| 23
| until( . > $s or ($n % . == 0); . + 2)
| . > $s
end;
def primes: 2, (range(3;infinite;2) | select(is_prime));
# input: the number to be tested
def isprime($smalls):
if . < 2 then false
else sqrt as $s
| (first( $smalls[] as $p
| if . == $p then 1
elif . % $p == 0 then 0
elif $p > $s then 1
else empty
end) // null) as $result
| if $result then $result == 1
else ($smalls[-1] + 2)
| until( . > $s or ($n % . == 0); . + 2)
| . > $s
end
end;
# Assumes n is odd.
def firstPrimeFactor:
if (. == 1) then 1
elif (. % 3 == 0) then 3
elif (. % 5 == 0) then 5
else . as $n
| [4, 2, 4, 2, 4, 6, 2, 6] as $inc
| { k: 7,
i: 0 }
| ($n | sqrt) as $s
| until (.k > $s or .done;
if $n % .k == 0
then .done = true
else .k += $inc[.i]
| .i = (.i + 1) % 8
end )
| if .done then .k else $n end
end ;
### Blum integers
# Number of small primes to pre-compute
def task($numberOfSmallPrimes):
[limit($numberOfSmallPrimes; primes)] as $smalls
| { blum: [],
bc:0,
counts: { "1": 0, "3": 0, "7": 0, "9": 0 },
i: 1 }
| label $out
| foreach range(0; infinite) as $_ (.;
(.i|firstPrimeFactor) as $p
| .j = null
| if ($p % 4 == 3)
then (.i / $p) as $q
| if $q != $p and ($q % 4 == 3) and ($q | isprime($smalls))
then if (.bc < 50) then .blum[.bc] = .i else . end
| .counts[(.i % 10) | tostring] += 1
| .bc += 1
| .j = .i
else .
end
else .
end
| .i |= if (. % 5 == 3) then . + 4 else . + 2 end;
select(.j)
| if (.bc == 50)
then "First 50 Blum integers:",
(.blum | tprint(10; 3) )
elif .bc == 26828 or .bc % 1e5 == 0
then "The \(.bc) Blum integer is: \(.j)",
if .bc == 400000
then "\n% distribution of the first 400,000 Blum integers:",
((.counts|keys_unsorted[]) as $k
| " \( .counts[$k] / 4000 )% end in \($k)"),
break $out
else empty
end
else empty
end);
task(10000)
- Output:
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 Blum integer is: 524273 The 100000 Blum integer is: 2075217 The 200000 Blum integer is: 4275533 The 300000 Blum integer is: 6521629 The 400000 Blum integer is: 8802377 % distribution of the first 400,000 Blum integers: 25.00125% end in 1 25.01675% end in 3 24.99725% end in 7 24.98475% end in 9
Julia
using Formatting, Primes
function isblum(n)
pe = factor(n).pe
return length(pe) == 2 && all(p -> p[2] == 1 && p[1] % 4 == 3, pe)
end
const blum400k = @view (filter(isblum, 1:9_000_000))[1:400_000]
println("First 50 Blum integers:")
foreach(p -> print(rpad(p[2], 4), p[1] % 10 == 0 ? "\n" : ""), enumerate(blum400k[1:50]))
for idx in [26_828, 100_000, 200_000, 300_000, 400_000]
println("\nThe $(format(idx, commas = true))th Blum number is ",
format(blum400k[idx], commas = true), ".")
end
println("\n% distribution of the first 400,000 Blum integers is:")
for d in [1, 3, 7, 9]
println(lpad(round(count(x -> x % 10 == d, blum400k) / 4000, digits=3), 8), "% end in $d")
end
- Output:
Same as Wren, Go, etc
Mathematica / Wolfram Language
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]];
- Output:
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
Maxima
/* Predicate functions that checks wether an integer is a Blum integer or not */
blum_p(n):=lambda([x],length(ifactors(x))=2 and unique(map(second,ifactors(x)))=[1] and mod(ifactors(x)[1][1],4)=3 and mod(ifactors(x)[2][1],4)=3)(n)$
/* Function that returns a list of the first len Blum integers */
blum_count(len):=block(
[i:1,count:0,result:[]],
while count<len do (if blum_p(i) then (result:endcons(i,result),count:count+1),i:i+1),
result)$
/* Test cases */
/* First 50 Blum integers */
blum_count(50);
/* Blum integer number 26828 */
last(blum_count(26828));
- Output:
[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] 524273
Nim
import std/[strformat, tables]
func isPrime(n: Natural): bool =
## Return "true" is "n" is prime.
if n < 2: return false
if (n and 1) == 0: return n == 2
if n mod 3 == 0: return n == 3
var d = 5
var step = 2
while d * d <= n:
if n mod d == 0:
return false
inc d, step
step = 6 - step
return true
const Inc = [4, 2, 4, 2, 4, 6, 2, 6]
func firstPrimeFactor(n: Positive): int =
## Return the first prime factor.
## Assuming "n" is odd.
if n == 1: return 1
if n mod 3 == 0: return 3
if n mod 5 == 0: return 5
var k = 7
var i = 0
while k * k <= n:
if n mod k == 0:
return k
k += Inc[i]
i = (i + 1) and 7
return n
var blum: array[50, int]
var bc = 0
var counts: CountTable[int]
var n = 1
while true:
var p = n.firstPrimeFactor
if (p and 3) == 3:
let q = n div p
if q != p and (q and 3) == 3 and q.isPrime:
if bc < 50: blum[bc] = n
counts.inc(n mod 10)
inc bc
if bc == 50:
echo "First 50 Blum integers:"
for i, val in blum:
stdout.write &"{val:3}"
stdout.write if i mod 10 == 9: '\n' else: ' '
echo()
elif bc == 26828 or bc mod 100000 == 0:
echo &"The {bc:>6}th Blum integer is: {n:>7}"
if bc == 400000:
echo "\n% distribution of the first 400_000 Blum integers:"
for i in [1, 3, 7, 9]:
echo &" {counts[i]/4000:6.3f} % end in {i}"
break
n = if n mod 5 == 3: n + 4 else: n + 2
- Output:
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 % 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
Pascal
Free Pascal
Generating Blum integer by multiplying primes of form 4*n+3.
Tried to sieve numbers of form 4*n+3.
Could not start at square of prime, since 5,13 are not getting sieved, but 35 =5*7 == 4*8 +3.
Using a simple prime sieve and check for 4*n+3 would be easier and faster.
program BlumInteger;
{$IFDEF FPC} {$MODE DELPHI}{$Optimization ON,ALL} {$ENDIF}
{$IFDEF WINDOWS}{$APPTYPE CONSOLE}{$ENDIF}
{
// for commatize = Numb2USA(IntToStr(n))
uses
sysutils, //IntToStr
strutils; //Numb2USA
}
const
LIMIT = 10*1000*1000;// >750
type
tP4n3 = array of Uint32;
function CommaUint(n : Uint64):AnsiString;
//commatize only positive Integers
var
fromIdx,toIdx :Int32;
pRes : pChar;
Begin
str(n,result);
fromIdx := length(result);
toIdx := fromIdx-1;
if toIdx < 3 then
exit;
toIdx := 4*(toIdx DIV 3)+toIdx MOD 3 +1 ;
setlength(result,toIdx);
pRes := @result[1];
dec(pRes);
repeat
pRes[toIdx] := pRes[FromIdx];
pRes[toIdx-1] := pRes[FromIdx-1];
pRes[toIdx-2] := pRes[FromIdx-2];
pRes[toIdx-3] := ',';
dec(toIdx,4);
dec(FromIdx,3);
until FromIdx<=3;
while fromIdx>=1 do
Begin
pRes[toIdx] := pRes[FromIdx];
dec(toIdx);
dec(fromIdx);
end;
end;
procedure Sieve4n_3_Primes(Limit:Uint32;var P4n3:tP4n3);
var
sieve : array of byte;
BlPrCnt,idx,n,j: Uint32;
begin
//DIV 3 -> smallest factor of Blum Integer
n := (LIMIT DIV 3 -3) DIV 4+ 1;
setlength(sieve,n);
setlength(P4n3,n);
BlPrCnt:= 0;
idx := 0;
repeat
if sieve[idx]= 0 then
begin
n := idx*4+3;
P4n3[BlPrCnt] := n;
inc(BlPrCnt);
j := idx+n;
if j > High(sieve) then
break;
while j <= High(sieve) do
begin
sieve[j] := 1;
inc(j,n);
end;
end;
inc(idx);
until idx>High(sieve);
//collect the rest
for idx := idx+1 to High(sieve) do
if sieve[idx]=0 then
Begin
P4n3[BlPrCnt] := idx*4+3;
inc(BlPrCnt);
end;
setlength(P4n3,BlPrCnt);
setlength(sieve,0);
end;
var
BlumField : array[0..LIMIT] of boolean;
BlumPrimes : tP4n3;
EndDigit : array[0..9] of Uint32;
k : Uint64;
n,idx,j,P4n3Cnt : Uint32;
begin
Sieve4n_3_Primes(Limit,BlumPrimes);
P4n3Cnt := length(BlumPrimes);
writeln('There are ',CommaUint(P4n3Cnt),' needed primes 4*n+3 to Limit ',CommaUint(LIMIT));
dec(P4n3Cnt);
writeln;
//generate Blum-Integers
For idx := 0 to P4n3Cnt do
Begin
n := BlumPrimes[idx];
For j := idx+1 to P4n3Cnt do
Begin
k := n*BlumPrimes[j];
if k>LIMIT then
BREAK;
BlumField[k] := true;
end;
end;
writeln('First 50 Blum-Integers ');
idx :=0;
j := 0 ;
repeat
while (idx<LIMIT) AND Not(BlumField[idx]) do
inc(idx);
if idx = LIMIT then
BREAK;
if j mod 10 = 0 then
writeln;
write(idx:5);
inc(j);
inc(idx);
until j >= 50;
Writeln(#13,#10);
//count and calc and summate decimal digit
writeln(' relative occurence of digit');
writeln(' n.th |Blum-Integer|Digit: 1 3 7 9');
idx :=0;
j := 0 ;
n := 0;
k := 26828;
repeat
while (idx<LIMIT) AND Not(BlumField[idx]) do
inc(idx);
if idx = LIMIT then
BREAK;
//count last decimal digit
inc(EndDigit[idx MOD 10]);
inc(j);
if j = k then
begin
write(CommaUint(j):10,' | ',CommaUint(idx):11,'| ');
write(EndDigit[1]/j*100:7:3,'% |');
write(EndDigit[3]/j*100:7:3,'% |');
write(EndDigit[7]/j*100:7:3,'% |');
writeln(EndDigit[9]/j*100:7:3,'%');
if k < 100000 then
k := 100000
else
k += 100000;
end;
inc(idx);
until j >= 400000;
Writeln;
end.
- @TIO.RUN:
There are 119,644 needed primes 4*n+3 to Limit 10,000,000 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 relative occurence of digit n.th |Blum-Integer| 1 3 7 9 26,828 | 524,273| 24.933% | 25.071% | 25.030% | 24.966% 100,000 | 2,075,217| 24.973% | 25.026% | 25.005% | 24.996% 200,000 | 4,275,533| 24.990% | 24.986% | 25.033% | 24.992% 300,000 | 6,521,629| 24.982% | 25.014% | 25.033% | 24.970% 400,000 | 8,802,377| 25.001% | 25.017% | 24.997% | 24.985% Real time: 0.097 s User time: 0.068 s Sys. time: 0.028 s CPU share: 99.08 % C-Version //-O3 -marchive=native Real time: 1.658 s User time: 1.612 s Sys. time: 0.033 s CPU share: 99.18 %
Perl
use v5.36;
use ntheory qw(is_square is_semiprime factor vecall);
sub comma { reverse((reverse shift) =~ s/.{3}\K/,/gr) =~ s/^,//r }
sub table ($c, @V) { my $t = $c * (my $w = 5); (sprintf(('%' . $w . 'd') x @V, @V)) =~ s/.{1,$t}\K/\n/gr }
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] == @blum;
}
$C{$_ % 10}++ for @blum;
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>;
- Output:
The first fifty 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 26,828th Blum integer: 524,273 The 100,000th Blum integer: 2,075,217 The 200,000th Blum integer: 4,275,533 The 300,000th Blum integer: 6,521,629 The 400,000th Blum integer: 8,802,377 1: 100005 (25.001%) 3: 100067 (25.017%) 7: 99989 (24.997%) 9: 99939 (24.985%)
Phix
You can run this online here.
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
p4n3 &= 4*j-1
end if
end for
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
- Output:
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 26,828th Blum integer is: 524,273 The 100,000th Blum integer is: 2,075,217 The 200,000th Blum integer is: 4,275,533 The 300,000th Blum integer is: 6,521,629 The 400,000th Blum integer is: 8,802,377 Percentage 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
Python
''' python example for task rosettacode.org/wiki/Blum_integer '''
from sympy import factorint
def generate_blum():
''' Generate the Blum integers in series '''
for candidate in range(1, 10_000_000_000):
fexp = factorint(candidate).items()
if len(fexp) == 2 and sum(p[1] == 1 and p[0] % 4 == 3 for p in fexp) == 2:
yield candidate
print('First 50 Blum integers:')
lastdigitsums = [0, 0, 0, 0]
for idx, blum in enumerate(generate_blum()):
if idx < 50:
print(f'{blum: 4}', end='\n' if (idx + 1) % 10 == 0 else '')
elif idx + 1 in [26_828, 100_000, 200_000, 300_000, 400_000]:
print(f'\nThe {idx+1:,}th Blum number is {blum:,}.')
j = blum % 10
lastdigitsums[0] += j == 1
lastdigitsums[1] += j == 3
lastdigitsums[2] += j == 7
lastdigitsums[3] += j == 9
if idx + 1 == 400_000:
print('\n% distribution of the first 400,000 Blum integers is:')
for k, dig in enumerate([1, 3, 7, 9]):
print(f'{lastdigitsums[k]/4000:>8.5}% end in {dig}')
break
- Output:
Same as Wren example.
Quackery
sqrt
, isprime
, and eratosthenes
are defined at Sieve of Eratosthenes#Quackery.
from
, index
, incr
, and end
are defined at Loops/Increment loop index within loop body#Quackery.
In the line 600000 eratosthenes
, 600000
suggests foreknowledge of the value of the 26828th Blum integer, and while I was not the first person to complete this task, so knew that 524273
would suffice, I reasoned that if I had determined the generally linear relationship between n
and Blum(n)
from computing the first fifty Blum integers with a slower but unbounded primality test I would have felt confident that 600000
would suffice. See the graphs at https://oeis.org/A016105/graph for confirmation.
Precomputing the primes takes a hot minute, but the overhead is worth it for the overall speed gain.
blum
returns 0
(i.e. false
) if the number passed to it, and the value of the smaller divisor (non-zero is taken as true
) rather than false
and true
(i.e. 1
) as the divisor is available without extra calculation. So as well as listing the Blum integers, their factors are shown in the output.
600000 eratosthenes
[ dup sqrt
tuck dup * = ] is exactsqrt ( n --> n b )
[ dup isprime iff
[ drop false ] done
dup 4 mod 1 != iff
[ drop false ] done
dup exactsqrt iff
[ 2drop false ] done
temp put
3 from
[ 4 incr
index temp share > iff
[ drop false end ]
done
index isprime not if done
dup index /mod 0 != iff
drop done
isprime not if done
drop index end ]
temp release ] is blum ( n --> n )
[ dup blum
over echo
say " = "
dup echo
say " * "
/ echo cr ] is echoblum ( n --> )
say "The First 50 Blum integers:"
cr cr
[] 1 from
[ 4 incr
index blum if
[ index join ]
dup size 50 = if end ]
witheach echoblum
cr
say "The 26828th Blum integer:"
cr cr
0 1 from
[ 4 incr
index blum if 1+
dup 26828 = if
[ drop index end ] ]
echoblum
- Output:
The First 50 Blum integers: 21 = 3 * 7 33 = 3 * 11 57 = 3 * 19 69 = 3 * 23 77 = 7 * 11 93 = 3 * 31 129 = 3 * 43 133 = 7 * 19 141 = 3 * 47 161 = 7 * 23 177 = 3 * 59 201 = 3 * 67 209 = 11 * 19 213 = 3 * 71 217 = 7 * 31 237 = 3 * 79 249 = 3 * 83 253 = 11 * 23 301 = 7 * 43 309 = 3 * 103 321 = 3 * 107 329 = 7 * 47 341 = 11 * 31 381 = 3 * 127 393 = 3 * 131 413 = 7 * 59 417 = 3 * 139 437 = 19 * 23 453 = 3 * 151 469 = 7 * 67 473 = 11 * 43 489 = 3 * 163 497 = 7 * 71 501 = 3 * 167 517 = 11 * 47 537 = 3 * 179 553 = 7 * 79 573 = 3 * 191 581 = 7 * 83 589 = 19 * 31 597 = 3 * 199 633 = 3 * 211 649 = 11 * 59 669 = 3 * 223 681 = 3 * 227 713 = 23 * 31 717 = 3 * 239 721 = 7 * 103 737 = 11 * 67 749 = 7 * 107 The 26828th Blum integer: 524273 = 223 * 2351
Raku
use List::Divvy;
use Lingua::EN::Numbers;
sub is-blum(Int $n ) {
return False if $n.is-prime;
my $factor = find-factor($n);
return True if ($factor.is-prime && ( my $div = $n div $factor ).is-prime && ($div != $factor)
&& ($factor % 4 == 3) && ($div % 4 == 3));
False;
}
sub find-factor ( Int $n, $constant = 1 ) {
my $x = 2;
my $rho = 1;
my $factor = 1;
while $factor == 1 {
$rho *= 2;
my $fixed = $x;
for ^$rho {
$x = ( $x * $x + $constant ) % $n;
$factor = ( $x - $fixed ) gcd $n;
last if 1 < $factor;
}
}
$factor = find-factor( $n, $constant + 1 ) if $n == $factor;
$factor;
}
my @blum = lazy (2..Inf).hyper(:1000batch).grep: &is-blum;
say "The first fifty Blum integers:\n" ~
@blum[^50].batch(10)».fmt("%3d").join: "\n";
say '';
printf "The %9s Blum integer: %9s\n", .&ordinal-digit(:c), comma @blum[$_-1] for 26828, 1e5, 2e5, 3e5, 4e5;
say "\nBreakdown by last digit:";
printf "%d => %%%7.4f:\n", .key, +.value / 4e3 for sort @blum[^4e5].categorize: {.substr(*-1)}
- Output:
The first fifty 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 26,828th Blum integer: 524,273 The 100,000th Blum integer: 2,075,217 The 200,000th Blum integer: 4,275,533 The 300,000th Blum integer: 6,521,629 The 400,000th Blum integer: 8,802,377 Breakdown by last digit: 1 => %25.0013: 3 => %25.0168: 7 => %24.9973: 9 => %24.9848:
RPL
Blum integers are necessarily of the form 4k + 21, which allows to speed up the quest.
≪ FACTORS @ n FACTORS returns { p1 n1 p2 n2 .. } for n = p1n1 * p2n2 * .. CASE DUP SIZE 4 ≠ THEN DROP 0 END LIST→ DROP ROT * 1 ≠ THEN DROP2 0 END 4 MOD SWAP 4 MOD * 9 == END ≫ 'BLUM?' STO ≪ { } 17 DO 4 + IF DUP BLUM? THEN SWAP OVER + SWAP END UNTIL OVER SIZE 50 == END 50 SWAP DO 4 + IF DUP BLUM? THEN SWAP 1 + SWAP END UNTIL OVER 26828 == END NIP ≫ 'TASK' STO
- Output:
2: { 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 } 1: 524273
Runs in 10 minutes 55 on the iHP48 emulator.
Ruby
require 'prime'
BLUM_EXP = [1, 1]
blums = (1..).step(2).lazy.select do |c|
next if c % 5 == 0
primes, exps = c.prime_division.transpose
exps == BLUM_EXP && primes.all?{|pr| (pr-3) % 4 == 0}
end
n, m = 50, 26828
res = blums.first(m)
puts "First #{n} Blum numbers:"
res.first(n).each_slice(10){|s| puts "%4d"*s.size % s}
puts "\n#{m}th Blum number: #{res.last}"
- Output:
First 50 Blum numbers: 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 26828th Blum number: 524273
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
}
}
- Output:
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
Sidef
Takes about 30 seconds:
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}%)"
}
}
}
- Output:
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%)
Wren
import "./math" for Int
import "./fmt" for Fmt
var inc = [4, 2, 4, 2, 4, 6, 2, 6]
// Assumes n is odd.
var firstPrimeFactor = Fn.new { |n|
if (n == 1) return 1
if (n%3 == 0) return 3
if (n%5 == 0) return 5
var k = 7
var i = 0
while (k * k <= n) {
if (n%k == 0) {
return k
} else {
k = k + inc[i]
i = (i + 1) % 8
}
}
return n
}
var blum = List.filled(50, 0)
var bc = 0
var counts = { 1: 0, 3: 0, 7: 0, 9: 0 }
var i = 1
while (true) {
var p = firstPrimeFactor.call(i)
if (p % 4 == 3) {
var q = i / p
if (q != p && q % 4 == 3 && Int.isPrime(q)) {
if (bc < 50) blum[bc] = i
counts[i % 10] = counts[i % 10] + 1
bc = bc + 1
if (bc == 50) {
System.print("First 50 Blum integers:")
Fmt.tprint("$3d ", blum, 10)
System.print()
} else if (bc == 26828 || bc % 1e5 == 0) {
Fmt.print("The $,9r Blum integer is: $,9d", bc, i)
if (bc == 400000) {
System.print("\n\% distribution of the first 400,000 Blum integers:")
for (i in [1, 3, 7, 9]) {
Fmt.print(" $6.3f\% end in $d", counts[i]/4000, i)
}
return
}
}
}
}
i = (i % 5 == 3) ? i + 4 : i + 2
}
- Output:
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 26,828th Blum integer is: 524,273 The 100,000th Blum integer is: 2,075,217 The 200,000th Blum integer is: 4,275,533 The 300,000th Blum integer is: 6,521,629 The 400,000th Blum integer is: 8,802,377 % 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
XPL0
Simple minded brute force takes 93 seconds on Pi4.
int Prime1;
func Semiprime(N); \Return 'true' if N is semiprime
int N, F, C;
[C:= 0; F:= 2;
repeat if rem(N/F) = 0 then
[C:= C+1;
Prime1:= N;
N:= N/F;
]
else F:= F+1;
until F > N;
return C = 2;
];
int N, C, Prime2;
[Format(4,0);
N:= 3; C:= 0;
loop [if Semiprime(N) then
[if rem(Prime1/4) = 3 then
[Prime2:= N/Prime1;
if Prime2 # Prime1 and rem(Prime2/4) = 3 then
[C:= C+1;
if C <= 50 then
[RlOut(0, float(N));
if rem(C/10) = 0 then CrLf(0);
];
if rem(C/1000)=0 then ChOut(0, ^.);
if C >= 26828 then
[Text(0, "^m^jThe 26828th Blum integer is: ");
IntOut(0, N); CrLf(0);
quit;
];
];
];
];
N:= N+2;
];
]
- Output:
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
Faster factoring. Does stretch. Takes 9.2 seconds.
include xpllib; \for Print
int Prime1;
func Semiprime(N); \Returns 'true' if odd N >= 3 is semiprime
int N, D, C;
[C:= 0; D:= 3;
while D*D <= N do
[while rem(N/D) = 0 do
[if C = 2 then return false;
C:= C+1;
N:= N/D;
];
D:= D+2;
];
Prime1:= N;
return C = 1;
];
int N, C, I, Goal, Prime2;
int FD, DC(10); \final digit and digit counters
[Text(0, "First 50 Blum integers:^m^j");
N:= 3; C:= 0; Goal:= 100_000;
for I:= 0 to 9 do DC(I):= 0;
loop [if Semiprime(N) then
[if rem(Prime1/4) = 3 then
[Prime2:= N/Prime1;
if rem(Prime2/4) = 3 and Prime2 # Prime1 then
[C:= C+1;
if C <= 50 then
[Print("%5d", N);
if rem(C/10) = 0 then Print("\n");
];
if C = 26_828 then
Print("\nThe 26,828th Blum integer is: %7,d\n", N);
if C = Goal then
[Print("The %6,dth Blum integer is: %7,d\n", Goal, N);
if Goal = 400_000 then quit;
Goal:= Goal + 100_000;
];
FD:= rem(N/10);
DC(FD):= DC(FD)+1;
];
];
];
N:= N+2;
];
Print("\n% distribution of the first 400,000 Blum integers:\n");
for I:= 0 to 9 do
if DC(I) > 0 then
Print("%4.3f\% end in %d\n", float(DC(I)) / 4_000., I);
]
- Output:
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 26,828th Blum integer is: 524,273 The 100,000th Blum integer is: 2,075,217 The 200,000th Blum integer is: 4,275,533 The 300,000th Blum integer is: 6,521,629 The 400,000th Blum integer is: 8,802,377 % 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
- Programming Tasks
- Solutions by Programming Task
- ALGOL 68
- BASIC
- BASIC256
- FreeBASIC
- Gambas
- OxygenBasic
- PureBasic
- QB64
- C
- C sharp
- C++
- Common Lisp
- EasyLang
- Fortran
- FutureBasic
- Go
- Go-rcu
- Haskell
- J
- Java
- Jq
- Julia
- Mathematica
- Wolfram Language
- Maxima
- Nim
- Pascal
- Free Pascal
- Perl
- Ntheory
- Phix
- Python
- Quackery
- Raku
- RPL
- Ruby
- Scala
- Sidef
- Wren
- Wren-math
- Wren-fmt
- XPL0