1234567891011121314151617181920212223242526272829303132333435363738394041424344454647484950515253545556575859606162636465666768697071727374757677787980818283848586878889 |
- #!/usr/bin/perl
- # Author: Trizen
- # Date: 28 February 2021
- # Edit: 11 April 2024
- # https://github.com/trizen
- # Fast recursive algorithm for generating all the k-powerful numbers in a given range [A,B].
- # A positive integer n is considered k-powerful, if for every prime p that divides n, so does p^k.
- # Example:
- # 2-powerful = a^2 * b^3, for a,b >= 1
- # 3-powerful = a^3 * b^4 * c^5, for a,b,c >= 1
- # 4-powerful = a^4 * b^5 * c^6 * d^7, for a,b,c,d >= 1
- # OEIS:
- # https://oeis.org/A001694 -- 2-powerful numbers
- # https://oeis.org/A036966 -- 3-powerful numbers
- # https://oeis.org/A036967 -- 4-powerful numbers
- # https://oeis.org/A069492 -- 5-powerful numbers
- # https://oeis.org/A069493 -- 6-powerful numbers
- use 5.036;
- use ntheory qw(:all);
- sub divceil ($x, $y) { # ceil(x/y)
- (($x % $y == 0) ? 0 : 1) + divint($x, $y);
- }
- sub powerful_numbers ($A, $B, $k = 2) {
- my @powerful;
- sub ($m, $r) {
- my $from = 1;
- my $upto = rootint(divint($B, $m), $r);
- if ($r <= $k) {
- if ($A > $m) {
- # Optimization by Dana Jacobsen (from Math::Prime::Util::PP)
- my $l = divceil($A, $m);
- if (($l >> $r) == 0) {
- $from = 2;
- }
- else {
- $from = rootint($l, $r);
- $from++ if (powint($from, $r) != $l);
- }
- }
- foreach my $j ($from .. $upto) {
- push @powerful, mulint($m, powint($j, $r));
- }
- return;
- }
- foreach my $v ($from .. $upto) {
- gcd($m, $v) == 1 or next;
- is_square_free($v) or next;
- __SUB__->(mulint($m, powint($v, $r)), $r - 1);
- }
- }
- ->(1, 2 * $k - 1);
- sort { $a <=> $b } @powerful;
- }
- require Math::Sidef;
- my $A = int rand 1e5;
- my $B = int rand 1e7;
- foreach my $k (2 .. 5) {
- say "Testing: k = $k";
- my @a1 = powerful_numbers($A, $B, $k);
- my @a2 = Math::Sidef::powerful($k, $A, $B);
- my @a3 = grep { $_ >= $A } powerful_numbers(1, $B, $k);
- "@a1" eq "@a2" or die "error for: powerful_numbers($A, $B, $k)";
- "@a1" eq "@a3" or die "error for: powerful_numbers($A, $B, $k)";
- }
- say join(', ', powerful_numbers(1e6 - 1e4, 1e6, 2)); #=> 990025, 990125, 990584, 991232, 992016, 994009, 995328, 996004, 996872, 998001, 998784, 1000000
|