# Blum integer

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).

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.

References

## 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
```

## C

Translation of: Wren
```#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++

```#include <algorithm>
#include <cstdint>
#include <iomanip>
#include <iostream>
#include <vector>

bool is_prime_type_3(const int32_t number) {
if ( number < 2 ) return false;
if ( number % 2 == 0 ) return false;
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;
}

int32_t least_prime_factor(const int32_t number) {
if ( number == 1 ) { return 1; }
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;
}

int main() {
int32_t blums[50];
int32_t blum_count = 0;
int32_t last_digit_counts[10] = {};
int32_t number = 1;

while ( blum_count < 400'000 ) {
const int32_t prime = least_prime_factor(number);
if ( prime % 4 == 3 ) {
const int32_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 ( int32_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 << "\nPercent 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
```

## FreeBASIC

```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
```
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

Translation of: Wren
Library: Go-rcu
```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.
```

## 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
}}
```

```   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

Works with: jq
```### 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
[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);

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

## 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

Translation of: Wren
```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

Translation of: Raku
Library: ntheory
```use v5.36;
use enum <false true>;
use ntheory <is_prime gcd>;

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) {
return false if \$n < 2 or is_prime \$n;
my \$factor = find_factor(\$n);
my \$div = int(\$n / \$factor);
return true if is_prime(\$factor) && is_prime(\$div) && (\$div != \$factor) && (\$factor%4 == 3) && (\$div%4 == 3);
false;
}

sub find_factor (\$n, \$constant = 1) {
my(\$x, \$rho, \$factor) = (2, 1, 1);
while (\$factor == 1) {
\$rho *= 2;
my \$fixed = \$x;
for (0..\$rho) {
\$x = ( \$x * \$x + \$constant ) % \$n;
\$factor = gcd((\$x-\$fixed), \$n);
last if 1 < \$factor;
}
}
\$factor = find_factor(\$n, \$constant+1) if \$n == \$factor;
\$factor;
}

my(\$i, @blum, %C);
my @nth = (26828, 1e5, 2e5, 3e5, 4e5);

while (++\$i) {
push @blum, \$i if is_blum \$i;
last if \$nth[-1] == scalar @blum;
}
\$C{substr \$_, -1, 1}++ 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

Translation of: Pascal

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.

## 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.

Works with: HP version 49
```≪ 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
```
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
```

## Wren

Library: Wren-math
Library: Wren-fmt
```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
```