Brilliant numbers: Difference between revisions
(→{{header|Phix}}: replaced with C++ translation, much faster even/and with 1e15 limit, even on 32 bit.) |
m (→{{header|Phix}}: see talk page) |
||
Line 531: | Line 531: | ||
<span style="color: #008080;">for</span> <span style="color: #000000;">primes</span> <span style="color: #008080;">in</span> <span style="color: #000000;">primes_by_digits</span> <span style="color: #008080;">do</span> |
<span style="color: #008080;">for</span> <span style="color: #000000;">primes</span> <span style="color: #008080;">in</span> <span style="color: #000000;">primes_by_digits</span> <span style="color: #008080;">do</span> |
||
<span style="color: #008080;">for</span> <span style="color: #000000;">i</span><span style="color: #0000FF;">=</span><span style="color: #000000;">1</span> <span style="color: #008080;">to</span> <span style="color: #7060A8;">length</span><span style="color: #0000FF;">(</span><span style="color: #000000;">primes</span><span style="color: #0000FF;">)</span> <span style="color: #008080;">do</span> |
<span style="color: #008080;">for</span> <span style="color: #000000;">i</span><span style="color: #0000FF;">=</span><span style="color: #000000;">1</span> <span style="color: #008080;">to</span> <span style="color: #7060A8;">length</span><span style="color: #0000FF;">(</span><span style="color: #000000;">primes</span><span style="color: #0000FF;">)</span> <span style="color: #008080;">do</span> |
||
<span style="color: #000080;font-style:italic;">--see talk page |
|||
⚫ | <span style="color: #008080;">for</span> <span style="color: #000000;">j</span><span style="color: #0000FF;">=</span><span style="color: #000000;"> |
||
-- for j=i to length(primes) do </span> |
|||
⚫ | |||
<span style="color: #000000;">brilliant_numbers</span> <span style="color: #0000FF;">&=</span> <span style="color: #000000;">primes</span><span style="color: #0000FF;">[</span><span style="color: #000000;">i</span><span style="color: #0000FF;">]*</span><span style="color: #000000;">primes</span><span style="color: #0000FF;">[</span><span style="color: #000000;">j</span><span style="color: #0000FF;">]</span> |
<span style="color: #000000;">brilliant_numbers</span> <span style="color: #0000FF;">&=</span> <span style="color: #000000;">primes</span><span style="color: #0000FF;">[</span><span style="color: #000000;">i</span><span style="color: #0000FF;">]*</span><span style="color: #000000;">primes</span><span style="color: #0000FF;">[</span><span style="color: #000000;">j</span><span style="color: #0000FF;">]</span> |
||
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span> |
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span> |
Revision as of 22:17, 22 April 2022
You are encouraged to solve this task according to the task description, using any language you may know.
Brilliant numbers are a subset of semiprime numbers. Specifically, they are numbers that are the product of exactly two prime numbers that both have the same number of digits when expressed in base 10.
Brilliant numbers are useful in cryptography and when testing prime factoring algorithms.
- E.G.
- 3 × 3 (9) is a brilliant number.
- 2 × 7 (14) is a brilliant number.
- 113 × 691 (78083) is a brilliant number.
- 2 × 31 (62) is semiprime, but is not a brilliant number (different number of digits in the two factors).
- Task
- Find and display the first 100 brilliant numbers.
- For the orders of magnitude 1 through 6, find and show the first brilliant number greater than or equal to the order of magnitude, and, its position in the series (or the count of brilliant numbers up to that point).
- Stretch
- Continue for larger orders of magnitude.
- See also
C++
<lang cpp>#include <algorithm>
- include <chrono>
- include <iomanip>
- include <iostream>
- include <locale>
- include <vector>
- include <primesieve.hpp>
auto get_primes_by_digits(uint64_t limit) {
primesieve::iterator pi; std::vector<std::vector<uint64_t>> primes_by_digits; std::vector<uint64_t> primes; for (uint64_t p = 10; p < limit;) { uint64_t prime = pi.next_prime(); if (prime > p) { primes_by_digits.push_back(std::move(primes)); p *= 10; } primes.push_back(prime); } return primes_by_digits;
}
int main() {
std::cout.imbue(std::locale(""));
auto start = std::chrono::high_resolution_clock::now();
auto primes_by_digits = get_primes_by_digits(1000000000);
std::cout << "First 100 brilliant numbers:\n"; std::vector<uint64_t> brilliant_numbers; for (const auto& primes : primes_by_digits) { for (auto i = primes.begin(); i != primes.end(); ++i) for (auto j = i; j != primes.end(); ++j) brilliant_numbers.push_back(*i * *j); if (brilliant_numbers.size() >= 100) break; } std::sort(brilliant_numbers.begin(), brilliant_numbers.end()); for (size_t i = 0; i < 100; ++i) { std::cout << std::setw(5) << brilliant_numbers[i] << ((i + 1) % 10 == 0 ? '\n' : ' '); }
std::cout << '\n'; uint64_t power = 10; size_t count = 0; for (size_t p = 1; p < 2 * primes_by_digits.size(); ++p) { const auto& primes = primes_by_digits[p / 2]; size_t position = count + 1; uint64_t min_product = 0; for (auto i = primes.begin(); i != primes.end(); ++i) { uint64_t p1 = *i; auto j = std::lower_bound(i, primes.end(), (power + p1 - 1) / p1); if (j != primes.end()) { uint64_t p2 = *j; uint64_t product = p1 * p2; if (min_product == 0 || product < min_product) min_product = product; position += std::distance(i, j); if (p1 >= p2) break; } } std::cout << "First brilliant number >= 10^" << p << " is " << min_product << " at position " << position << '\n'; power *= 10; if (p % 2 == 1) { size_t size = primes.size(); count += size * (size + 1) / 2; } }
auto end = std::chrono::high_resolution_clock::now(); std::chrono::duration<double> duration(end - start); std::cout << "\nElapsed time: " << duration.count() << " seconds\n";
}</lang>
- Output:
First 100 brilliant numbers: 4 6 9 10 14 15 21 25 35 49 121 143 169 187 209 221 247 253 289 299 319 323 341 361 377 391 403 407 437 451 473 481 493 517 527 529 533 551 559 583 589 611 629 649 667 671 689 697 703 713 731 737 767 779 781 793 799 803 817 841 851 869 871 893 899 901 913 923 943 949 961 979 989 1,003 1,007 1,027 1,037 1,067 1,073 1,079 1,081 1,121 1,139 1,147 1,157 1,159 1,189 1,207 1,219 1,241 1,247 1,261 1,271 1,273 1,333 1,343 1,349 1,357 1,363 1,369 First brilliant number >= 10^1 is 10 at position 4 First brilliant number >= 10^2 is 121 at position 11 First brilliant number >= 10^3 is 1,003 at position 74 First brilliant number >= 10^4 is 10,201 at position 242 First brilliant number >= 10^5 is 100,013 at position 2,505 First brilliant number >= 10^6 is 1,018,081 at position 10,538 First brilliant number >= 10^7 is 10,000,043 at position 124,364 First brilliant number >= 10^8 is 100,140,049 at position 573,929 First brilliant number >= 10^9 is 1,000,000,081 at position 7,407,841 First brilliant number >= 10^10 is 10,000,600,009 at position 35,547,995 First brilliant number >= 10^11 is 100,000,000,147 at position 491,316,167 First brilliant number >= 10^12 is 1,000,006,000,009 at position 2,409,600,866 First brilliant number >= 10^13 is 10,000,000,000,073 at position 34,896,253,010 First brilliant number >= 10^14 is 100,000,380,000,361 at position 174,155,363,187 First brilliant number >= 10^15 is 1,000,000,000,000,003 at position 2,601,913,448,897 Elapsed time: 0.167788 seconds
Go
<lang go>package main
import (
"fmt" "math" "rcu" "sort"
)
var primes = rcu.Primes(1e8 - 1)
type res struct {
bc interface{} next int
}
func getBrilliant(digits, limit int, countOnly bool) res {
var brilliant []int count := 0 pow := 1 next := math.MaxInt for k := 1; k <= digits; k++ { var s []int for _, p := range primes { if p >= pow*10 { break } if p > pow { s = append(s, p) } } for i := 0; i < len(s); i++ { for j := i; j < len(s); j++ { prod := s[i] * s[j] if prod < limit { if countOnly { count++ } else { brilliant = append(brilliant, prod) } } else { if next > prod { next = prod } break } } } pow *= 10 } if countOnly { return res{count, next} } return res{brilliant, next}
}
func main() {
fmt.Println("First 100 brilliant numbers:") brilliant := getBrilliant(2, 10000, false).bc.([]int) sort.Ints(brilliant) brilliant = brilliant[0:100] for i := 0; i < len(brilliant); i++ { fmt.Printf("%4d ", brilliant[i]) if (i+1)%10 == 0 { fmt.Println() } } fmt.Println() for k := 1; k <= 13; k++ { limit := int(math.Pow(10, float64(k))) r := getBrilliant(k, limit, true) total := r.bc.(int) next := r.next climit := rcu.Commatize(limit) ctotal := rcu.Commatize(total + 1) cnext := rcu.Commatize(next) fmt.Printf("First >= %18s is %14s in the series: %18s\n", climit, ctotal, cnext) }
}</lang>
- Output:
First 100 brilliant numbers: 4 6 9 10 14 15 21 25 35 49 121 143 169 187 209 221 247 253 289 299 319 323 341 361 377 391 403 407 437 451 473 481 493 517 527 529 533 551 559 583 589 611 629 649 667 671 689 697 703 713 731 737 767 779 781 793 799 803 817 841 851 869 871 893 899 901 913 923 943 949 961 979 989 1003 1007 1027 1037 1067 1073 1079 1081 1121 1139 1147 1157 1159 1189 1207 1219 1241 1247 1261 1271 1273 1333 1343 1349 1357 1363 1369 First >= 10 is 4 in the series: 10 First >= 100 is 11 in the series: 121 First >= 1,000 is 74 in the series: 1,003 First >= 10,000 is 242 in the series: 10,201 First >= 100,000 is 2,505 in the series: 100,013 First >= 1,000,000 is 10,538 in the series: 1,018,081 First >= 10,000,000 is 124,364 in the series: 10,000,043 First >= 100,000,000 is 573,929 in the series: 100,140,049 First >= 1,000,000,000 is 7,407,841 in the series: 1,000,000,081 First >= 10,000,000,000 is 35,547,995 in the series: 10,000,600,009 First >= 100,000,000,000 is 491,316,167 in the series: 100,000,000,147 First >= 1,000,000,000,000 is 2,409,600,866 in the series: 1,000,006,000,009 First >= 10,000,000,000,000 is 34,896,253,010 in the series: 10,000,000,000,073
J
<lang J>oprimes=: {{ NB. all primes of order y
p:(+i.)/-/\ p:inv +/\1 9*10^y
}}
obrill=: {{ NB. all brilliant numbers of order y primes
~.,*/~oprimes y
}}
brillseq=: {{ NB. sequences of brilliant numbers up through order y-1 primes
/:~;obrill each i.y
}}</lang>
Task examples:
<lang J> 10 10 $brillseq 2
4 6 9 10 14 15 21 25 35 49 121 143 169 187 209 221 247 253 289 299 319 323 341 361 377 391 403 407 437 451 473 481 493 517 527 529 533 551 559 583 589 611 629 649 667 671 689 697 703 713 731 737 767 779 781 793 799 803 817 841 851 869 871 893 899 901 913 923 943 949 961 979 989 1003 1007 1027 1037 1067 1073 1079
1081 1121 1139 1147 1157 1159 1189 1207 1219 1241 1247 1261 1271 1273 1333 1343 1349 1357 1363 1369 NB. order, index, value
(brillseq 4) (],.(I. 10^]) ([,.{) [) 1 2 3 4 5 6
1 3 10 2 10 121 3 73 1003 4 241 10201 5 2504 100013 6 10537 1018081</lang>
Stretch goal (results are order, index, value): <lang J> (brillseq 4) (],.(I. 10^]) ([,.{) [) ,7 7 124363 10000043
(brillseq 5) (],.(I. 10^]) ([,.{) [) 8 9
8 573928 100140049 9 7407840 1000000081</lang>
Julia
<lang julia> using Primes
function isbrilliant(n)
p = factor(n).pe return (length(p) == 1 && p[1][2] == 2) || length(p) == 2 && ndigits(p[1][1]) == ndigits(p[2][1]) && p[1][2] == p[2][2] == 1
end
function testbrilliants()
println("First 100 brilliant numbers:") foreach(p -> print(lpad(p[2], 5), p[1] % 20 == 0 ? "\n" : ""), enumerate(filter(isbrilliant, 1:1370))) bcount, results, positions = 0, zeros(Int, 9), zeros(Int, 9) for n in 1:10^10 if isbrilliant(n) bcount += 1 for i in 1:9 if n >= 10^i && results[i] == 0 results[i] = n positions[i] = bcount println("First >=", lpad(10^i, 12), " is", lpad(bcount, 8), " in the series: $n") end end end end return results, positions
end
testbrilliants()
</lang>
- Output:
First 100 brilliant numbers: 4 6 9 10 14 15 21 25 35 49 121 143 169 187 209 221 247 253 289 299 319 323 341 361 377 391 403 407 437 451 473 481 493 517 527 529 533 551 559 583 589 611 629 649 667 671 689 697 703 713 731 737 767 779 781 793 799 803 817 841 851 869 871 893 899 901 913 923 943 949 961 979 989 1003 1007 1027 1037 1067 1073 1079 1081 1121 1139 1147 1157 1159 1189 1207 1219 1241 1247 1261 1271 1273 1333 1343 1349 1357 1363 1369 First >= 10 is 4 in the series: 10 First >= 100 is 11 in the series: 121 First >= 1000 is 74 in the series: 1003 First >= 10000 is 242 in the series: 10201 First >= 100000 is 2505 in the series: 100013 First >= 1000000 is 10538 in the series: 1018081 First >= 10000000 is 124364 in the series: 10000043 First >= 100000000 is 573929 in the series: 100140049 First >= 1000000000 is 7407841 in the series: 1000000081
Perl
<lang perl>use strict; use warnings; use feature 'say'; use List::AllUtils <max head firstidx uniqint>; use ntheory <primes is_semiprime forsetproduct>;
sub table { my $t = shift() * (my $c = 1 + length max @_); ( sprintf( ('%'.$c.'d')x@_, @_) ) =~ s/.{1,$t}\K/\n/gr } sub comma { reverse ((reverse shift) =~ s/(.{3})/$1,/gr) =~ s/^,//r }
my(@B,@Br); for my $oom (1..5) {
my @P = grep { $oom == length } @{primes(10**$oom)}; forsetproduct { is_semiprime($_[0] * $_[1]) and push @B, $_[0] * $_[1] } \@P, \@P; @Br = uniqint sort { $a <=> $b } @Br, @B;
}
say "First 100 brilliant numbers:\n" . table 10, head 100, @Br;
for my $oom (1..9) {
my $key = firstidx { $_ > 10**$oom } @Br; printf "First >= %13s is position %9s in the series: %13s\n", comma(10**$oom), comma($key), comma $Br[$key];
}</lang>
- Output:
First 100 brilliant numbers: 4 6 9 10 14 15 21 25 35 49 121 143 169 187 209 221 247 253 289 299 319 323 341 361 377 391 403 407 437 451 473 481 493 517 527 529 533 551 559 583 589 611 629 649 667 671 689 697 703 713 731 737 767 779 781 793 799 803 817 841 851 869 871 893 899 901 913 923 943 949 961 979 989 1003 1007 1027 1037 1067 1073 1079 1081 1121 1139 1147 1157 1159 1189 1207 1219 1241 1247 1261 1271 1273 1333 1343 1349 1357 1363 1369 First >= 10 is position 4 in the series: 14 First >= 100 is position 10 in the series: 121 First >= 1,000 is position 73 in the series: 1,003 First >= 10,000 is position 241 in the series: 10,201 First >= 100,000 is position 2,504 in the series: 100,013 First >= 1,000,000 is position 10,537 in the series: 1,018,081 First >= 10,000,000 is position 124364 in the series: 10,000,043 First >= 100,000,000 is position 573929 in the series: 100,140,049 First >= 1,000,000,000 is position 7407841 in the series: 1,000,000,081
Faster approach
<lang perl>use 5.020; use strict; use warnings;
use ntheory qw(:all); use experimental qw(signatures);
sub is_briliant_number ($n) {
is_semiprime($n) || return; my @f = factor($n); length($f[0]) == length($f[1]);
}
sub next_brilliant_number ($n) {
++$n while not is_briliant_number($n); $n;
}
sub brilliant_numbers_count ($n) {
use integer;
my $count = 0; my $len = length(sqrtint($n));
foreach my $k (1 .. $len - 1) { my $pi = prime_count(10**($k - 1), 10**$k - 1); $count += binomial($pi, 2) + $pi; }
my $min = 10**($len - 1); my $max = 10**$len - 1;
my $pi_min = prime_count($min); my $pi_max = prime_count($max);
my $j = -1;
forprimes { if ($_*$_ <= $n) { $count += (($max <= $n/$_) ? $pi_max : prime_count($n/$_)) - $pi_min - ++$j; } else { lastfor; } } $min, $max;
return $count;
}
say "First 100 brilliant numbers:";
my @nums;
for (my $k = 1 ; scalar(@nums) < 100 ; ++$k) {
push(@nums, $k) if is_briliant_number($k);
}
while (@nums) {
my @slice = splice(@nums, 0, 10); say join ' ', map { sprintf("%4s", $_) } @slice;
}
say ;
foreach my $n (1 .. 13) {
my $v = next_brilliant_number(vecprod((10) x $n)); printf("First brilliant number >= 10^%d is %s", $n, $v); printf(" at position %s\n", brilliant_numbers_count($v));
}</lang>
- Output:
First 100 brilliant numbers: 4 6 9 10 14 15 21 25 35 49 121 143 169 187 209 221 247 253 289 299 319 323 341 361 377 391 403 407 437 451 473 481 493 517 527 529 533 551 559 583 589 611 629 649 667 671 689 697 703 713 731 737 767 779 781 793 799 803 817 841 851 869 871 893 899 901 913 923 943 949 961 979 989 1003 1007 1027 1037 1067 1073 1079 1081 1121 1139 1147 1157 1159 1189 1207 1219 1241 1247 1261 1271 1273 1333 1343 1349 1357 1363 1369 First brilliant number >= 10^1 is 10 at position 4 First brilliant number >= 10^2 is 121 at position 11 First brilliant number >= 10^3 is 1003 at position 74 First brilliant number >= 10^4 is 10201 at position 242 First brilliant number >= 10^5 is 100013 at position 2505 First brilliant number >= 10^6 is 1018081 at position 10538 First brilliant number >= 10^7 is 10000043 at position 124364 First brilliant number >= 10^8 is 100140049 at position 573929 First brilliant number >= 10^9 is 1000000081 at position 7407841 First brilliant number >= 10^10 is 10000600009 at position 35547995 First brilliant number >= 10^11 is 100000000147 at position 491316167 First brilliant number >= 10^12 is 1000006000009 at position 2409600866 First brilliant number >= 10^13 is 10000000000073 at position 34896253010
Phix
Replaced with C++ translation; much faster and now goes comfortably to 1e15 even on 32 bit. You can run this online here.
-- -- demo\rosetta\BrilliantNumbers.exw -- ================================= -- with javascript_semantics requires("1.0.2") -- (for in) atom t0 = time() function get_primes_by_digits(integer limit) sequence primes = get_primes_le(power(10,limit)), primes_by_digits = {} integer p = 10 while length(primes) do integer pi = abs(binary_search(p,primes))-1 primes_by_digits &= {primes[1..pi]} primes = primes[pi+1..$] p*= 10 end while return primes_by_digits end function sequence primes_by_digits = get_primes_by_digits(8) procedure first100() sequence brilliant_numbers = {} for primes in primes_by_digits do for i=1 to length(primes) do --see talk page -- for j=i to length(primes) do for j=1 to i do brilliant_numbers &= primes[i]*primes[j] end for end for if length(brilliant_numbers)>=100 then exit end if end for brilliant_numbers = sort(brilliant_numbers)[1..100] sequence j100 = join_by(brilliant_numbers,1,10," ","\n","%,5d") printf(1,"First 100 brilliant numbers:\n%s\n\n",{j100}) end procedure first100() atom pwr = 10, count = 0 for p=1 to 2*length(primes_by_digits)-1 do sequence primes = primes_by_digits[floor(p/2)+1] atom pos = count+1, min_product = 0 for i=1 to length(primes) do integer p1 = primes[i], j = abs(binary_search(floor((pwr+p1-1)/p1),primes,i)) if j<=length(primes) then -- (always is, I think) integer p2 = primes[j] atom prod = p1*p2 if min_product=0 or prod<min_product then min_product = prod end if pos += j-i if p1>=p2 then exit end if end if end for printf(1,"First brilliant number >= 10^%d is %,d at position %,d\n", {p, min_product, pos}) pwr *= 10; if odd(p) then integer size = length(primes) count += size * (size + 1) / 2; end if end for ?elapsed(time()-t0) {} = wait_key()
- Output:
First 100 brilliant numbers: 4 6 9 10 14 15 21 25 35 49 121 143 169 187 209 221 247 253 289 299 319 323 341 361 377 391 403 407 437 451 473 481 493 517 527 529 533 551 559 583 589 611 629 649 667 671 689 697 703 713 731 737 767 779 781 793 799 803 817 841 851 869 871 893 899 901 913 923 943 949 961 979 989 1,003 1,007 1,027 1,037 1,067 1,073 1,079 1,081 1,121 1,139 1,147 1,157 1,159 1,189 1,207 1,219 1,241 1,247 1,261 1,271 1,273 1,333 1,343 1,349 1,357 1,363 1,369 First brilliant number >= 10^1 is 10 at position 4 First brilliant number >= 10^2 is 121 at position 11 First brilliant number >= 10^3 is 1,003 at position 74 First brilliant number >= 10^4 is 10,201 at position 242 First brilliant number >= 10^5 is 100,013 at position 2,505 First brilliant number >= 10^6 is 1,018,081 at position 10,538 First brilliant number >= 10^7 is 10,000,043 at position 124,364 First brilliant number >= 10^8 is 100,140,049 at position 573,929 First brilliant number >= 10^9 is 1,000,000,081 at position 7,407,841 First brilliant number >= 10^10 is 10,000,600,009 at position 35,547,995 First brilliant number >= 10^11 is 100,000,000,147 at position 491,316,167 First brilliant number >= 10^12 is 1,000,006,000,009 at position 2,409,600,866 First brilliant number >= 10^13 is 10,000,000,000,073 at position 34,896,253,010 First brilliant number >= 10^14 is 100,000,380,000,361 at position 174,155,363,187 First brilliant number >= 10^15 is 1,000,000,000,000,003 at position 2,601,913,448,897 "3.3s"
Raku
1 through 7 are fast. 8 and 9 take a bit longer. <lang perl6>use Lingua::EN::Numbers;
- Find an abundance of primes to use to generate brilliants
my %primes = (2..100000).grep( &is-prime ).categorize: { .chars };
- Generate brilliant numbers
my @brilliant = lazy flat (1..*).map: -> $digits {
sort flat (^%primes{$digits}).race.map: { %primes{$digits}[$_] X× (flat %primes{$digits}[$_ .. *]) }
};
- The task
put "First 100 brilliant numbers:\n" ~ @brilliant[^100].batch(10)».fmt("%4d").join("\n") ~ "\n" ;
for 1 .. 7 -> $oom {
my $threshold = exp $oom, 10; my $key = @brilliant.first: :k, * >= $threshold; printf "First >= %13s is %9s in the series: %13s\n", comma($threshold), ordinal-digit(1 + $key, :u), comma @brilliant[$key];
}</lang>
- Output:
First 100 brilliant numbers: 4 6 9 10 14 15 21 25 35 49 121 143 169 187 209 221 247 253 289 299 319 323 341 361 377 391 403 407 437 451 473 481 493 517 527 529 533 551 559 583 589 611 629 649 667 671 689 697 703 713 731 737 767 779 781 793 799 803 817 841 851 869 871 893 899 901 913 923 943 949 961 979 989 1003 1007 1027 1037 1067 1073 1079 1081 1121 1139 1147 1157 1159 1189 1207 1219 1241 1247 1261 1271 1273 1333 1343 1349 1357 1363 1369 First >= 10 is 4ᵗʰ in the series: 10 First >= 100 is 11ᵗʰ in the series: 121 First >= 1,000 is 74ᵗʰ in the series: 1,003 First >= 10,000 is 242ⁿᵈ in the series: 10,201 First >= 100,000 is 2505ᵗʰ in the series: 100,013 First >= 1,000,000 is 10538ᵗʰ in the series: 1,018,081 First >= 10,000,000 is 124364ᵗʰ in the series: 10,000,043 First >= 100,000,000 is 573929ᵗʰ in the series: 100,140,049 First >= 1,000,000,000 is 7407841ˢᵗ in the series: 1,000,000,081
Rust
<lang rust>// [dependencies] // primal = "0.3" // indexing = "0.4.1"
fn get_primes_by_digits(limit: usize) -> Vec<Vec<usize>> {
let mut primes_by_digits = Vec::new(); let mut power = 10; let mut primes = Vec::new(); for prime in primal::Primes::all().take_while(|p| *p < limit) { if prime > power { primes_by_digits.push(primes); primes = Vec::new(); power *= 10; } primes.push(prime); } primes_by_digits
}
fn main() {
use indexing::algorithms::lower_bound; use std::time::Instant;
let start = Instant::now();
let primes_by_digits = get_primes_by_digits(1000000000);
println!("First 100 brilliant numbers:"); let mut brilliant_numbers = Vec::new(); for primes in &primes_by_digits { for i in 0..primes.len() { let p1 = primes[i]; for j in i..primes.len() { let p2 = primes[j]; brilliant_numbers.push(p1 * p2); } } if brilliant_numbers.len() >= 100 { break; } } brilliant_numbers.sort(); for i in 0..100 { let n = brilliant_numbers[i]; print!("{:4}{}", n, if (i + 1) % 10 == 0 { '\n' } else { ' ' }); }
println!(); let mut power = 10; let mut count = 0; for p in 1..2 * primes_by_digits.len() { let primes = &primes_by_digits[p / 2]; let mut position = count + 1; let mut min_product = 0; for i in 0..primes.len() { let p1 = primes[i]; let n = (power + p1 - 1) / p1; let j = lower_bound(&primes[i..], &n); let p2 = primes[i + j]; let product = p1 * p2; if min_product == 0 || product < min_product { min_product = product; } position += j; if p1 >= p2 { break; } } println!("First brilliant number >= 10^{p} is {min_product} at position {position}"); power *= 10; if p % 2 == 1 { let size = primes.len(); count += size * (size + 1) / 2; } }
let time = start.elapsed(); println!("\nElapsed time: {} milliseconds", time.as_millis());
}</lang>
- Output:
First 100 brilliant numbers: 4 6 9 10 14 15 21 25 35 49 121 143 169 187 209 221 247 253 289 299 319 323 341 361 377 391 403 407 437 451 473 481 493 517 527 529 533 551 559 583 589 611 629 649 667 671 689 697 703 713 731 737 767 779 781 793 799 803 817 841 851 869 871 893 899 901 913 923 943 949 961 979 989 1003 1007 1027 1037 1067 1073 1079 1081 1121 1139 1147 1157 1159 1189 1207 1219 1241 1247 1261 1271 1273 1333 1343 1349 1357 1363 1369 First brilliant number >= 10^1 is 10 at position 4 First brilliant number >= 10^2 is 121 at position 11 First brilliant number >= 10^3 is 1003 at position 74 First brilliant number >= 10^4 is 10201 at position 242 First brilliant number >= 10^5 is 100013 at position 2505 First brilliant number >= 10^6 is 1018081 at position 10538 First brilliant number >= 10^7 is 10000043 at position 124364 First brilliant number >= 10^8 is 100140049 at position 573929 First brilliant number >= 10^9 is 1000000081 at position 7407841 First brilliant number >= 10^10 is 10000600009 at position 35547995 First brilliant number >= 10^11 is 100000000147 at position 491316167 First brilliant number >= 10^12 is 1000006000009 at position 2409600866 First brilliant number >= 10^13 is 10000000000073 at position 34896253010 First brilliant number >= 10^14 is 100000380000361 at position 174155363187 First brilliant number >= 10^15 is 1000000000000003 at position 2601913448897 Elapsed time: 676 milliseconds
Sidef
<lang ruby>func is_briliant_number(n) {
n.is_semiprime && (n.factor.map{.len}.uniq.len == 1)
}
func brilliant_numbers_count(n) {
var count = 0 var len = n.isqrt.len
for k in (1 .. len-1) { var pi = prime_count(10**(k-1), 10**k - 1) count += binomial(pi, 2)+pi }
var min = (10**(len - 1)) var max = (10**len - 1)
each_prime(min, max, {|p| count += prime_count(p, max `min` idiv(n, p)) })
return count
}
say "First 100 brilliant numbers:"
100.by(is_briliant_number).each_slice(10, {|*a|
say a.map { '%4s' % _}.join(' ')
})
say
for n in (1 .. 12) {
var v = (10**n .. Inf -> first_by(is_briliant_number)) printf("First brilliant number >= 10^%d is %s", n, v) printf(" at position %s\n", brilliant_numbers_count(v))
}</lang>
- Output:
First 100 brilliant numbers: 4 6 9 10 14 15 21 25 35 49 121 143 169 187 209 221 247 253 289 299 319 323 341 361 377 391 403 407 437 451 473 481 493 517 527 529 533 551 559 583 589 611 629 649 667 671 689 697 703 713 731 737 767 779 781 793 799 803 817 841 851 869 871 893 899 901 913 923 943 949 961 979 989 1003 1007 1027 1037 1067 1073 1079 1081 1121 1139 1147 1157 1159 1189 1207 1219 1241 1247 1261 1271 1273 1333 1343 1349 1357 1363 1369 First brilliant number >= 10^1 is 10 at position 4 First brilliant number >= 10^2 is 121 at position 11 First brilliant number >= 10^3 is 1003 at position 74 First brilliant number >= 10^4 is 10201 at position 242 First brilliant number >= 10^5 is 100013 at position 2505 First brilliant number >= 10^6 is 1018081 at position 10538 First brilliant number >= 10^7 is 10000043 at position 124364 First brilliant number >= 10^8 is 100140049 at position 573929 First brilliant number >= 10^9 is 1000000081 at position 7407841 First brilliant number >= 10^10 is 10000600009 at position 35547995 First brilliant number >= 10^11 is 100000000147 at position 491316167 First brilliant number >= 10^12 is 1000006000009 at position 2409600866
Wren
<lang ecmascript>import "./math" for Int import "./seq" for Lst import "./fmt" for Fmt
var primes = Int.primeSieve(1e7-1)
var getBrilliant = Fn.new { |digits, limit, countOnly|
var brilliant = [] var count = 0 var pow = 1 var next = Num.maxSafeInteger for (k in 1..digits) { var s = primes.where { |p| p > pow && p < pow * 10 }.toList for (i in 0...s.count) { for (j in i...s.count) { var prod = s[i] * s[j] if (prod < limit) { if (countOnly) { count = count + 1 } else { brilliant.add(prod) } } else { next = next.min(prod) break } } } pow = pow * 10 } return countOnly ? [count, next] : [brilliant, next]
}
System.print("First 100 brilliant numbers:") var brilliant = getBrilliant.call(2, 10000, false)[0] brilliant.sort() brilliant = brilliant[0..99] for (chunk in Lst.chunks(brilliant, 10)) Fmt.print("$4d", chunk) System.print() for (k in 1..12) {
var limit = 10.pow(k) var res = getBrilliant.call(k, limit, true) var total = res[0] var next = res[1] Fmt.print("First >= $,17d is $,15r in the series: $,17d", limit, total + 1, next)
}</lang>
- Output:
First 100 brilliant numbers: 4 6 9 10 14 15 21 25 35 49 121 143 169 187 209 221 247 253 289 299 319 323 341 361 377 391 403 407 437 451 473 481 493 517 527 529 533 551 559 583 589 611 629 649 667 671 689 697 703 713 731 737 767 779 781 793 799 803 817 841 851 869 871 893 899 901 913 923 943 949 961 979 989 1003 1007 1027 1037 1067 1073 1079 1081 1121 1139 1147 1157 1159 1189 1207 1219 1241 1247 1261 1271 1273 1333 1343 1349 1357 1363 1369 First >= 10 is 4th in the series: 10 First >= 100 is 11th in the series: 121 First >= 1,000 is 74th in the series: 1,003 First >= 10,000 is 242nd in the series: 10,201 First >= 100,000 is 2,505th in the series: 100,013 First >= 1,000,000 is 10,538th in the series: 1,018,081 First >= 10,000,000 is 124,364th in the series: 10,000,043 First >= 100,000,000 is 573,929th in the series: 100,140,049 First >= 1,000,000,000 is 7,407,841st in the series: 1,000,000,081 First >= 10,000,000,000 is 35,547,995th in the series: 10,000,600,009 First >= 100,000,000,000 is 491,316,167th in the series: 100,000,000,147 First >= 1,000,000,000,000 is 2,409,600,866th in the series: 1,000,006,000,009
XPL0
<lang XPL0> func NumDigits(N); \Return number of digits in N int N, Cnt; [Cnt:= 0; repeat N:= N/10;
Cnt:= Cnt+1;
until N = 0; return Cnt; ];
func Brilliant(N); \Return 'true' if N is a brilliant number int N, Limit, Cnt, F; int A(3); [Limit:= sqrt(N); Cnt:= 0; F:= 2; loop [if rem(N/F) = 0 then
[A(Cnt):= F; Cnt:= Cnt+1; if Cnt > 2 then quit; N:= N/F; ] else F:= F+1; if F > N then quit; if F > Limit then [A(Cnt):= N; Cnt:= Cnt+1; quit; ]; ];
if Cnt # 2 then return false; return NumDigits(A(0)) = NumDigits(A(1)); ];
int Cnt, N, Mag; [Format(5, 0); Cnt:= 0; N:= 4; loop [if Brilliant(N) then
[RlOut(0, float(N)); Cnt:= Cnt+1; if Cnt >= 100 then quit; if rem(Cnt/10) = 0 then CrLf(0); ]; N:= N+1; ];
CrLf(0); CrLf(0); Format(7, 0); Cnt:= 0; N:= 4; Mag:= 10; loop [if Brilliant(N) then
[Cnt:= Cnt+1; if N >= Mag then [Text(0, "First >= "); RlOut(0, float(Mag)); Text(0, " is "); RlOut(0, float(Cnt)); Text(0, " in series: "); RlOut(0, float(N)); CrLf(0); if Mag >= 1_000_000 then quit; Mag:= Mag*10; ]; ]; N:= N+1; ];
]</lang>
- Output:
4 6 9 10 14 15 21 25 35 49 121 143 169 187 209 221 247 253 289 299 319 323 341 361 377 391 403 407 437 451 473 481 493 517 527 529 533 551 559 583 589 611 629 649 667 671 689 697 703 713 731 737 767 779 781 793 799 803 817 841 851 869 871 893 899 901 913 923 943 949 961 979 989 1003 1007 1027 1037 1067 1073 1079 1081 1121 1139 1147 1157 1159 1189 1207 1219 1241 1247 1261 1271 1273 1333 1343 1349 1357 1363 1369 First >= 10 is 4 in series: 10 First >= 100 is 11 in series: 121 First >= 1000 is 74 in series: 1003 First >= 10000 is 242 in series: 10201 First >= 100000 is 2505 in series: 100013 First >= 1000000 is 10538 in series: 1018081