factorial_from_primorials.pl

#!/usr/bin/perl

# Daniel "Trizen" Șuteu
# Date: 21 March 2019
# https://github.com/trizen

# Given a positive integer n, find the unique integer k such that the product of the primorial of all its prime factors, is equal to n!.

# Example for n = 10:
#   a(10) = 5040 = 2^4 * 3^2 * 5 * 7

# By mapping each prime factor `p` to `primorial(p)`, we get:
#
#   primorial(2)^4 * primorial(3)^2 * primorial(5) * primorial(7) = 10!
#
# where `primorial(p)` is the product of primes <= p.

# OEIS sequence by Allan C. Wechsler (Mar 20 2019):
#   https://oeis.org/A307035

# Efficient formula:
#   a(n) = a(n-1) * (A319626(n) / A319627(n))

use 5.020;
use strict;
use warnings;

use Memoize qw(memoize);
use experimental qw(signatures);
use ntheory qw(factor_exp prev_prime);
use Math::AnyNum qw(:overload ipow factorial primorial prod);

memoize('f');

sub g($n) {

    my $prod = 1;
    foreach my $pp (factor_exp($n)) {
        my ($p, $e) = @$pp;

        if ($p > 2) {
            $prod *= (Math::AnyNum->new($p) / prev_prime($p))**$e;
        }
        else {
            $prod *= ipow($p, $e);
        }
    }

    return $prod;
}

sub f($n) {
    return 1 if ($n <= 1);
    f($n - 1) * g($n);
}

sub isok ($n, $v) {
    prod(map { ipow(primorial($_->[0]), $_->[1]) } factor_exp($v)) == factorial($n);
}

foreach my $n (1 .. 100) {
    my $v = f($n);
    say "a($n) = $v = ", join(' * ', map { "$_->[0]^$_->[1]" } factor_exp($v));
    isok($n, $v) or die "error for $n";
}

__END__
a(1) = 1 =
a(2) = 2 = 2^1
a(3) = 3 = 3^1
a(4) = 12 = 2^2 * 3^1
a(5) = 20 = 2^2 * 5^1
a(6) = 60 = 2^2 * 3^1 * 5^1
a(7) = 84 = 2^2 * 3^1 * 7^1
a(8) = 672 = 2^5 * 3^1 * 7^1
a(9) = 1512 = 2^3 * 3^3 * 7^1
a(10) = 5040 = 2^4 * 3^2 * 5^1 * 7^1
a(11) = 7920 = 2^4 * 3^2 * 5^1 * 11^1
a(12) = 47520 = 2^5 * 3^3 * 5^1 * 11^1
a(13) = 56160 = 2^5 * 3^3 * 5^1 * 13^1
a(14) = 157248 = 2^6 * 3^3 * 7^1 * 13^1
a(15) = 393120 = 2^5 * 3^3 * 5^1 * 7^1 * 13^1
a(16) = 6289920 = 2^9 * 3^3 * 5^1 * 7^1 * 13^1
a(17) = 8225280 = 2^9 * 3^3 * 5^1 * 7^1 * 17^1
a(18) = 37013760 = 2^8 * 3^5 * 5^1 * 7^1 * 17^1
a(19) = 41368320 = 2^8 * 3^5 * 5^1 * 7^1 * 19^1
a(20) = 275788800 = 2^10 * 3^4 * 5^2 * 7^1 * 19^1
a(21) = 579156480 = 2^9 * 3^5 * 5^1 * 7^2 * 19^1
a(22) = 1820206080 = 2^10 * 3^5 * 5^1 * 7^1 * 11^1 * 19^1
a(23) = 2203407360 = 2^10 * 3^5 * 5^1 * 7^1 * 11^1 * 23^1
a(24) = 26440888320 = 2^12 * 3^6 * 5^1 * 7^1 * 11^1 * 23^1
a(25) = 73446912000 = 2^12 * 3^4 * 5^3 * 7^1 * 11^1 * 23^1
a(26) = 173601792000 = 2^13 * 3^4 * 5^3 * 7^1 * 13^1 * 23^1
a(27) = 585906048000 = 2^10 * 3^7 * 5^3 * 7^1 * 13^1 * 23^1
a(28) = 3281073868800 = 2^12 * 3^7 * 5^2 * 7^2 * 13^1 * 23^1
a(29) = 4137006182400 = 2^12 * 3^7 * 5^2 * 7^2 * 13^1 * 29^1
a(30) = 20685030912000 = 2^12 * 3^7 * 5^3 * 7^2 * 13^1 * 29^1
a(31) = 22111584768000 = 2^12 * 3^7 * 5^3 * 7^2 * 13^1 * 31^1
a(32) = 707570712576000 = 2^17 * 3^7 * 5^3 * 7^2 * 13^1 * 31^1
a(33) = 1667845251072000 = 2^16 * 3^8 * 5^3 * 7^1 * 11^1 * 13^1 * 31^1
a(34) = 4362056810496000 = 2^17 * 3^8 * 5^3 * 7^1 * 11^1 * 17^1 * 31^1
a(35) = 10178132557824000 = 2^17 * 3^7 * 5^3 * 7^2 * 11^1 * 17^1 * 31^1
a(36) = 91603193020416000 = 2^17 * 3^9 * 5^3 * 7^2 * 11^1 * 17^1 * 31^1
a(37) = 109332843282432000 = 2^17 * 3^9 * 5^3 * 7^2 * 11^1 * 17^1 * 37^1
a(38) = 244391061454848000 = 2^18 * 3^9 * 5^3 * 7^2 * 11^1 * 19^1 * 37^1
a(39) = 433238699851776000 = 2^17 * 3^10 * 5^3 * 7^2 * 13^1 * 19^1 * 37^1
a(40) = 5776515998023680000 = 2^20 * 3^9 * 5^4 * 7^2 * 13^1 * 19^1 * 37^1
a(41) = 6401004214026240000 = 2^20 * 3^9 * 5^4 * 7^2 * 13^1 * 19^1 * 41^1
a(42) = 26884217698910208000 = 2^20 * 3^10 * 5^3 * 7^3 * 13^1 * 19^1 * 41^1
a(43) = 28195642952515584000 = 2^20 * 3^10 * 5^3 * 7^3 * 13^1 * 19^1 * 43^1
a(44) = 177229755701526528000 = 2^22 * 3^10 * 5^3 * 7^2 * 11^1 * 13^1 * 19^1 * 43^1
a(45) = 664611583880724480000 = 2^20 * 3^11 * 5^4 * 7^2 * 11^1 * 13^1 * 19^1 * 43^1
a(46) = 1609059624132280320000 = 2^21 * 3^11 * 5^4 * 7^2 * 11^1 * 13^1 * 23^1 * 43^1
a(47) = 1758739589167841280000 = 2^21 * 3^11 * 5^4 * 7^2 * 11^1 * 13^1 * 23^1 * 47^1
a(48) = 42209750140028190720000 = 2^24 * 3^12 * 5^4 * 7^2 * 11^1 * 13^1 * 23^1 * 47^1
a(49) = 82731110274455253811200 = 2^24 * 3^12 * 5^2 * 7^4 * 11^1 * 13^1 * 23^1 * 47^1
a(50) = 459617279302529187840000 = 2^25 * 3^10 * 5^4 * 7^4 * 11^1 * 13^1 * 23^1 * 47^1
a(51) = 901556970939576483840000 = 2^24 * 3^11 * 5^4 * 7^4 * 11^1 * 17^1 * 23^1 * 47^1
a(52) = 4261905680805270650880000 = 2^26 * 3^11 * 5^4 * 7^4 * 13^1 * 17^1 * 23^1 * 47^1
a(53) = 4805978746439986053120000 = 2^26 * 3^11 * 5^4 * 7^4 * 13^1 * 17^1 * 23^1 * 53^1
a(54) = 32440356538469905858560000 = 2^24 * 3^14 * 5^4 * 7^4 * 13^1 * 17^1 * 23^1 * 53^1
a(55) = 84962838553135467724800000 = 2^24 * 3^13 * 5^5 * 7^3 * 11^1 * 13^1 * 17^1 * 23^1 * 53^1
a(56) = 951583791795117238517760000 = 2^27 * 3^13 * 5^4 * 7^4 * 11^1 * 13^1 * 17^1 * 23^1 * 53^1
a(57) = 1595302239185931841044480000 = 2^26 * 3^14 * 5^4 * 7^4 * 11^1 * 13^1 * 19^1 * 23^1 * 53^1
a(58) = 4022936081425393338286080000 = 2^27 * 3^14 * 5^4 * 7^4 * 11^1 * 13^1 * 19^1 * 29^1 * 53^1
a(59) = 4478362807624494470922240000 = 2^27 * 3^14 * 5^4 * 7^4 * 11^1 * 13^1 * 19^1 * 29^1 * 59^1
a(60) = 44783628076244944709222400000 = 2^28 * 3^14 * 5^5 * 7^4 * 11^1 * 13^1 * 19^1 * 29^1 * 59^1
a(61) = 46301717163575281818009600000 = 2^28 * 3^14 * 5^5 * 7^4 * 11^1 * 13^1 * 19^1 * 29^1 * 61^1
a(62) = 98989878073850602507468800000 = 2^29 * 3^14 * 5^5 * 7^4 * 11^1 * 13^1 * 19^1 * 31^1 * 61^1
a(63) = 311818115932629397898526720000 = 2^27 * 3^16 * 5^4 * 7^5 * 11^1 * 13^1 * 19^1 * 31^1 * 61^1
a(64) = 19956359419688281465505710080000 = 2^33 * 3^16 * 5^4 * 7^5 * 11^1 * 13^1 * 19^1 * 31^1 * 61^1
a(65) = 39307980675143584704783974400000 = 2^33 * 3^15 * 5^5 * 7^5 * 13^2 * 19^1 * 31^1 * 61^1
a(66) = 185309051754248327893981593600000 = 2^33 * 3^16 * 5^5 * 7^4 * 11^1 * 13^2 * 19^1 * 31^1 * 61^1
a(67) = 203536171598928491293389619200000 = 2^33 * 3^16 * 5^5 * 7^4 * 11^1 * 13^2 * 19^1 * 31^1 * 67^1
a(68) = 1064650743748241339073114931200000 = 2^35 * 3^16 * 5^5 * 7^4 * 11^1 * 13^1 * 17^1 * 19^1 * 31^1 * 67^1
a(69) = 1933181613648122431474866585600000 = 2^34 * 3^17 * 5^5 * 7^4 * 11^1 * 13^1 * 17^1 * 23^1 * 31^1 * 67^1
a(70) = 9021514197024571346882710732800000 = 2^35 * 3^16 * 5^5 * 7^5 * 11^1 * 13^1 * 17^1 * 23^1 * 31^1 * 67^1
a(71) = 9560112059533500979532424806400000 = 2^35 * 3^16 * 5^5 * 7^5 * 11^1 * 13^1 * 17^1 * 23^1 * 31^1 * 71^1
a(72) = 172082017071603017631583646515200000 = 2^36 * 3^18 * 5^5 * 7^5 * 11^1 * 13^1 * 17^1 * 23^1 * 31^1 * 71^1
a(73) = 176929397834183384325431073177600000 = 2^36 * 3^18 * 5^5 * 7^5 * 11^1 * 13^1 * 17^1 * 23^1 * 31^1 * 73^1
a(74) = 422347594829986143228448368230400000 = 2^37 * 3^18 * 5^5 * 7^5 * 11^1 * 13^1 * 17^1 * 23^1 * 37^1 * 73^1
a(75) = 1759781645124942263451868200960000000 = 2^36 * 3^17 * 5^7 * 7^5 * 11^1 * 13^1 * 17^1 * 23^1 * 37^1 * 73^1
a(76) = 7867259119382094824843646074880000000 = 2^38 * 3^17 * 5^7 * 7^5 * 11^1 * 13^1 * 19^1 * 23^1 * 37^1 * 73^1
a(77) = 17307970062640608614656021364736000000 = 2^38 * 3^17 * 5^6 * 7^5 * 11^2 * 13^1 * 19^1 * 23^1 * 37^1 * 73^1
a(78) = 61364621131180339633780439384064000000 = 2^38 * 3^18 * 5^6 * 7^5 * 11^1 * 13^2 * 19^1 * 23^1 * 37^1 * 73^1
a(79) = 66408288621414340151625407004672000000 = 2^38 * 3^18 * 5^6 * 7^5 * 11^1 * 13^2 * 19^1 * 23^1 * 37^1 * 79^1
a(80) = 1770887696571049070710010853457920000000 = 2^42 * 3^17 * 5^7 * 7^5 * 11^1 * 13^2 * 19^1 * 23^1 * 37^1 * 79^1
a(81) = 8965118963890935920469429945630720000000 = 2^38 * 3^21 * 5^7 * 7^5 * 11^1 * 13^2 * 19^1 * 23^1 * 37^1 * 79^1
a(82) = 19868642028082614742661979879505920000000 = 2^39 * 3^21 * 5^7 * 7^5 * 11^1 * 13^2 * 19^1 * 23^1 * 41^1 * 79^1
a(83) = 20874649219377937008113219367075840000000 = 2^39 * 3^21 * 5^7 * 7^5 * 11^1 * 13^2 * 19^1 * 23^1 * 41^1 * 83^1
a(84) = 175347053442774670868151042683437056000000 = 2^40 * 3^22 * 5^6 * 7^6 * 11^1 * 13^2 * 19^1 * 23^1 * 41^1 * 83^1
a(85) = 382166654939380692917765093028003840000000 = 2^40 * 3^21 * 5^7 * 7^6 * 11^1 * 13^1 * 17^1 * 19^1 * 23^1 * 41^1 * 83^1
a(86) = 801617861580164380266531658546544640000000 = 2^41 * 3^21 * 5^7 * 7^6 * 11^1 * 13^1 * 17^1 * 19^1 * 23^1 * 43^1 * 83^1
a(87) = 1516103346901615240938875093338030080000000 = 2^40 * 3^22 * 5^7 * 7^6 * 11^1 * 13^1 * 17^1 * 19^1 * 29^1 * 43^1 * 83^1
a(88) = 19059584932477448743231572601963806720000000 = 2^43 * 3^22 * 5^7 * 7^5 * 11^2 * 13^1 * 17^1 * 19^1 * 29^1 * 43^1 * 83^1
a(89) = 20437386252897505278886867006925045760000000 = 2^43 * 3^22 * 5^7 * 7^5 * 11^2 * 13^1 * 17^1 * 19^1 * 29^1 * 43^1 * 89^1
a(90) = 153280396896731289591651502551937843200000000 = 2^42 * 3^23 * 5^8 * 7^5 * 11^2 * 13^1 * 17^1 * 19^1 * 29^1 * 43^1 * 89^1
a(91) = 253609383956409951869823395131388067840000000 = 2^42 * 3^23 * 5^7 * 7^6 * 11^1 * 13^2 * 17^1 * 19^1 * 29^1 * 43^1 * 89^1
a(92) = 1228003332841563977474934334320405381120000000 = 2^44 * 3^23 * 5^7 * 7^6 * 11^1 * 13^2 * 17^1 * 23^1 * 29^1 * 43^1 * 89^1
a(93) = 1969039826797680170778774018824098283520000000 = 2^43 * 3^24 * 5^7 * 7^6 * 11^1 * 13^2 * 17^1 * 23^1 * 31^1 * 43^1 * 89^1
a(94) = 4304412644627486884958250180685238108160000000 = 2^44 * 3^24 * 5^7 * 7^6 * 11^1 * 13^2 * 17^1 * 23^1 * 31^1 * 47^1 * 89^1
a(95) = 8018023553717867726883015042452894515200000000 = 2^44 * 3^23 * 5^8 * 7^6 * 11^1 * 13^2 * 19^1 * 23^1 * 31^1 * 47^1 * 89^1
a(96) = 384865130578457650890384722037738936729600000000 = 2^48 * 3^24 * 5^8 * 7^6 * 11^1 * 13^2 * 19^1 * 23^1 * 31^1 * 47^1 * 89^1
a(97) = 419459749057420136363677730760232324300800000000 = 2^48 * 3^24 * 5^8 * 7^6 * 11^1 * 13^2 * 19^1 * 23^1 * 31^1 * 47^1 * 97^1
a(98) = 1644282216305086934545616704580110711259136000000 = 2^49 * 3^24 * 5^6 * 7^8 * 11^1 * 13^2 * 19^1 * 23^1 * 31^1 * 47^1 * 97^1
a(99) = 5813712121935843090000573348336820014809088000000 = 2^47 * 3^26 * 5^6 * 7^7 * 11^2 * 13^2 * 19^1 * 23^1 * 31^1 * 47^1 * 97^1
a(100) = 64596801354842701000006370537075777942323200000000 = 2^49 * 3^24 * 5^8 * 7^7 * 11^2 * 13^2 * 19^1 * 23^1 * 31^1 * 47^1 * 97^1