Blum integer: Difference between revisions

From Rosetta Code
Content added Content deleted
(J draft)
Line 294: Line 294:
Same as Wren example.
Same as Wren example.
</pre>
</pre>

=={{header|J}}==

Implementation: <syntaxhighlight lang=J>blumseq=: {{
b=. 1e3
p=. (#~ 1 p: ])<.3 4 p.i.b
r=. /:~ (</~ #&, */~) p
L=. done=. 0
while. -.done do.
L0=. L
P=. (#~ 1 p: ])<.3 4 p. b+1+i.1e3
R1=. p (</ #&, */) P
R2=. (</~ #&, */~) P
b=. b+1e3
p=. p,P
r=. /:~ r, R1, R2
L=. r I. {.R2
if. y<L do.
L=. y{r
done=. (L=L0)*0=#$L0
r=. (y+1){.r
end.
end.
y{.r
}}</syntaxhighlight>

In other words, generate some blum numbers and keep adding to them until the smallest number being added to the sequence exceeds the largest value within the desired sequence length.

Task examples:

<syntaxhighlight lang=J> 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
524621</syntaxhighlight>


=={{header|Julia}}==
=={{header|Julia}}==

Revision as of 22:16, 23 May 2023

Blum integer is a draft programming task. It is not yet considered ready to be promoted as a complete task, for reasons that should be found in its talk page.
Definition

A positive integer n is a Blum integer if n = p x q is a semi-prime for which p and q are distinct primes congruent to 3 mod 4. In other words, p and q must be of the form 4t + 3 where t is some non-negative integer.

Example

21 is a Blum integer because it has two prime factors: 3 (= 4 x 0 + 3) and 7 (= 4 x 1 + 3).

Task

Find and show on this page the first 50 Blum integers.

Also show the 26,828th.

Stretch

Find and show the 100,000th, 200,000th, 300,000th and 400,000th Blum integers.

For the first 400,000 Blum integers, show the percentage distribution by final decimal digit (to 3 decimal places). Clearly, such integers can only end in 1, 3, 7 or 9.

Related task
References


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. 

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
        FOR j FROM i + i BY i TO UPB upfc DO
            upfc[ j ] +:= 1;
            lpf[  j ]  := i
        OD
    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 WHILE b count < 400 000 DO
        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 unique prime factors - could be a Blum integer            #
            IF lpf[ i ] MOD 4 = 3 THEN
                # the last prime factor mod 4 is three                      #
                IF ( i OVER lpf[ i ] ) MOD 4 = 3 THEN
                    # and so is the other one                               #
                    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.

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: 

blumseq=: {{
  b=. 1e3
  p=. (#~ 1 p: ])<.3 4 p.i.b
  r=. /:~ (</~ #&, */~) p
  L=. done=. 0
  while. -.done do.
    L0=. L
    P=. (#~ 1 p: ])<.3 4 p. b+1+i.1e3
    R1=. p (</ #&, */) P
    R2=. (</~ #&, */~) P    
    b=. b+1e3
    p=. p,P
    r=. /:~ r, R1, R2
    L=. r I. {.R2
    if. y<L do.
      L=. y{r
      done=. (L=L0)*0=#$L0
      r=. (y+1){.r
    end.
  end.
  y{.r
}}

In other words, generate some blum numbers and keep adding to them until the smallest number being added to the sequence exceeds the largest value within the desired sequence length.

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
524621

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

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.5} % 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 %

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:

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
Retrieved from "https://rosettacode.org/wiki/Blum_integer?oldid=343020"