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

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

• Semiprime

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

Translation of: FreeBASIC

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

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

Gambas

Translation of: FreeBASIC

Public Prime1 As Integer

Public Sub Main() 
  
  Dim n As Integer = 3, c As Integer = 0, Prime2 As Integer
  
  Print "The first 50 Blum integers:" 
  Do 
    If isSemiprime(n) Then 
      If Prime1 Mod 4 = 3 Then 
        Prime2 = n / Prime1 
        If (Prime2 <> Prime1) And (Prime2 Mod 4 = 3) Then 
          c += 1 
          If c <= 50 Then 
            Print Format$(n, "####"); 
            If c Mod 10 = 0 Then Print 
          End If 
          If c >= 26828 Then  
            Print "\nThe 26828th Blum integer is: "; n 
            Break
          End If 
        End If 
      End If 
    End If 
    n += 2 
  Loop

End

Function isSemiprime(n As Integer) As Boolean 

  Dim d As Integer = 3, c As Integer = 0 
  While d * d <= n 
    While n Mod d = 0 
      If c = 2 Then Return False 
      n /= d 
      c += 1 
    Wend 
    d += 2 
  Wend 
  Prime1 = n 
  Return c = 1 

End Function

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#

Translation of: Java

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

EasyLang

Translation of: FreeBASIC

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

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.

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

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

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

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.

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.

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 
≫ '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

Translation of: Java

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

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