Blum integer: Difference between revisions
(Created Nim solution.) |
(Replace the "blum" sequence with an array.) |
||
Line 355: | Line 355: | ||
var blum |
var blum: array[50, int] |
||
var bc = 0 |
var bc = 0 |
||
var counts: CountTable[int] |
var counts: CountTable[int] |
Revision as of 14:35, 23 May 2023
- Definition
A positive integer n is a Blum integer if n = p x q is a semi-prime for which p and q are distinct primes congruent to 3 mod 4. In other words, p and q must be of the form 4t + 3 where t is some non-negative integer.
- Example
21 is a Blum integer because it has two prime factors: 3 (= 4 x 0 + 3) and 7 (= 4 x 1 + 3).
- Task
Find and show on this page the first 50 Blum integers.
Also show the 26,828th.
- Stretch
Find and show the 100,000th, 200,000th, 300,000th and 400,000th Blum integers.
For the first 400,000 Blum integers, show the percentage distribution by final decimal digit (to 3 decimal places). Clearly, such integers can only end in 1, 3, 7 or 9.
- Related task
- References
- Wikipedia article Blum integer
- OEIS sequence A016105: Blum integers
ALGOL 68
If running this with ALGOL 68G, you will need to specify a larger heap size with -heap 256M
on the command line.
Builds tables of the unique prime factor counts and the last prime factor to handle the stretch goal, uses prime factorisation to find the first 50 Blum integers.
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
#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
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.
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
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;//8802377;//
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;
var
BlumField : array [0..LIMIT] of boolean;
BlumPrimes : array[0..332398] of Uint32;
sieve : array[0..((LIMIT-3) DIV 4)] of boolean;
EndDigit : array[0..9] of Uint32;
n,idx,j,k,BlPrCnt : Uint64;
begin
//sieve primes 4*i+3
BlPrCnt:= 0;
IDX := 0;
repeat
if NOT(sieve[idx]) then
begin
n := idx*4+3;
BlumPrimes[BlPrCnt] := n;
inc(BlPrCnt);
j := idx+n;
if j > High(sieve) then
break;
while j <= High(sieve) do
begin
sieve[j] := true;
inc(j,n);
end;
end;
inc(idx);
until idx>High(sieve);
//collect the rest
for idx := idx+1 to High(sieve) do
if NOT(sieve[idx]) then
Begin
BlumPrimes[BlPrCnt] := idx*4+3;
inc(BlPrCnt);
end;
writeln('There ',BlPrCnt,' primes 4*n+3 to Limit ',LIMIT);
dec(BlPrCnt);
//generate Blum-Integers
For idx := 0 to BlPrCnt do
Begin
n := BlumPrimes[idx];
For j := idx+1 to BlPrCnt 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 Not(BlumField[idx]) do
inc(idx);
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
idx :=0;
j := 0 ;
k := 26828;
repeat
while Not(BlumField[idx]) do
inc(idx);
//count last decimal digit
inc(EndDigit[idx MOD 10]);
inc(j);
if j = k then
begin
writeln(CommaUint(j):8,'.th Blum integer is: ',CommaUint(idx):10);
if k < 100000 then
k := 100000
else
k += 100000;
end;
inc(idx);
until j >= 400000;
Writeln;
idx := 0;
For j := 0 to 9 do
inc(idx,EndDigit[j]);
For j := 0 to 9 do
if EndDigit[j] <> 0 then
writeln('Digit ',j:2,' rel. ',EndDigit[j]/idx*100:7:3,'%');
end.
- @TIO.RUN:
There 332398 primes 4*n+3 to Limit 10000000 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,828.th Blum integer is: 524,273 100,000.th Blum integer is: 2,075,217 200,000.th Blum integer is: 4,275,533 300,000.th Blum integer is: 6,521,629 400,000.th Blum integer is: 8,802,377 Digit 1 rel. 25.001% Digit 3 rel. 25.017% Digit 7 rel. 24.997% Digit 9 rel. 24.985% Real time: 0.108 s User time: 0.076 s Sys. time: 0.032 s CPU share: 98.90 % C-Version //-O3 -marchive=native Real time: 1.658 s User time: 1.612 s Sys. time: 0.033 s CPU share: 99.18 %
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
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