### ----------------------------------------------------
### -- Math::Sequences
### -- Licenses: Artistic-2.0
### -- Authors: Aaron Sherman, JJ Merelo, Anton Antonov
### -- File: lib/Math/Sequences/Integer.rakumod
### ----------------------------------------------------
### -- Chunk 1 of 3
# Integer sequences
#
# This module is absolutely not meant to be exhausive.
# It should contain only those sequences that are frequently needed.

unit module Math::Sequences::Integer;

class Integers is Range is export {
    multi method new(
            $min = -Inf,
            $max = Inf,
            :$excludes-min = True,
            :$excludes-max = True) {
        #Seems as if this doesn't get passed on correctly...
        #nextwith $min, $max, :$excludes-min, :$excludes-max
        self.Range::new($min, $max, :$excludes-min, :$excludes-max);
    }

    method is-default {
        self.min == -Inf and self.max == Inf and
            self.excludes-min and self.excludes-max;
    }

    method symbol { "ℤ" }
    method gist {
        if self.is-default {
            self.symbol;
        } else {
            my $emin = self.excludes-min ?? '^' !! '';
            my $emax = self.excludes-max ?? '^' !! '';
            "{self.^name}({self.min}{$emin}..{$emax}{self.max})"
        }
    }

    method !params {
        self.min.raku, self.max.raku,
            |(<excludes-min excludes-max>.map: -> $param {
                ":{$param}(" ~ self."$param"() ~ ')'
            });
    }

    method raku {
        "Integers.new(" ~ self!params.join(",") ~ ")"
    }

    method !min-countable {
        self.min.succ > self.min;
    }

    # An iterator based on the current range (full of fail)
    method iterator {
        class :: does Iterator {
            has $.min;
            has $.max;
            has $.offset;
            has $.countable;
            has $.endcmp;
            has $.cur = 0;

            method is-lazy()  { True  }

            method pull-one() {
                my $count = $!cur++;
                my $next = $!min + $count + $!offset;
                my $err = "Cannot count from non-finite starting point";
                my &endcmp = $!endcmp;

                if      endcmp($next, $!max)          { IterationEnd }
                elsif   $count == 0 and $!offset == 0 { $!min }
                elsif   ! $!countable                 { fail $err }
                else                                  { $next }
            }
        }.new(
            :min(self.min), :max(self.max),
            :countable(self!min-countable),
            :offset(self.excludes-min ?? 1 !! 0),
            :endcmp(self.excludes-max ?? {$^cur > $^end} !! {$^cur >= $^end}));
    }

    # All of the integers >= $n
    method from(Int $n) {
        my $min = $n;
        my $max = self.max;
        my $excludes-min = self.excludes-min;
        my $excludes-max = self.excludes-max;

        self.WHAT.new($min, $max, :$excludes-min, :$excludes-max);
    }

    method of { ::Int }
    method Numeric { Inf }
    method is-int { True }
    method infinite { True }
    method elems { Inf }
    method Str { self.gist }
}

# Naturals can mean 1..Inf or 0..Inf.
# For ease of naming, we define 𝕀 using Naturals and ℕ as the same
# sequence, counting from 1.
class Naturals is Integers is export {
    method new(
            $min = 0,
            $max = Inf,
            :$excludes-min = False,
            :$excludes-max = True) {
        nextwith $min, $max, :$excludes-min, :$excludes-max
    }

    method symbol { '𝕀' }
    method is-default {
        self.min == 0 and self.max == Inf and
            !self.excludes-min and self.excludes-max;
    }

    method from($n where * >= 0) { callsame }

    # Triangular counting
    method triangle {
        lazy gather for |self -> $n {
            for self.min .. $n -> $k {
                take ($n, $k);
            }
        }
    }
}

our constant \ℤ is export = Integers.new;
our constant \𝕀 is export = Naturals.new;
our constant \ℕ is export = 𝕀.from(1);

our constant \Z is export = Integers.new;
our constant \I is export = Naturals.new;
our constant \N is export = 𝕀.from(1);

#####
# Utilities for the OEIS entires:

# via:
# https://github.com/perl6/perl6-examples/blob/master/categories/best-of-rosettacode/binomial-coefficient.pl
sub infix:<choose>($n, $p) is export(:support) { [*] ($n ... 0) Z/ 1 .. $p }
sub infix:<ichoose>($n,$p --> Int) is export(:support) { ($n choose $p).floor }

# Per OEIS A000111:
# 2*a(n+1) = Sum_{k=0..n} binomial(n, k)*a(k)*a(n-k).
sub euler-up-down($i is copy) is export(:support) {
    $i++;
    if $i < 2 {
        1
    } else {
        my $sum = [+] (0..$i).map: -> $k {
            ($i ichoose $k) * euler-up-down($k) * euler-up-down($i-$k) };
        $sum div 2;
    }
}

# Per OEIS A000111:
# From Wikipedia:
#  The Euler zigzag numbers are related to Entringer numbers, from which the zigzag numbers may be computed. The Entringer numbers can be defined recursively as follows:[3]

#    E ( 0 , 0 ) = 1 {\displaystyle E(0,0)=1} E(0,0)=1
#    E ( n , 0 ) = 0 for  n > 0 {\displaystyle E(n,0)=0\qquad {\mbox{for }}n>0} E(n,0)=0\qquad {\mbox{for }}n>0
#    E ( n , k ) = E ( n , k − 1 ) + E ( n − 1 , n − k ) {\displaystyle E(n,k)=E(n,k-1)+E(n-1,n-k)} E(n,k)=E(n,k-1)+E(n-1,n-k).
#
#  The nth zigzag number is equal to the Entringer number E(n, n).

my %Entringer;
multi sub Entringer(0, 0 --> 1) { };

multi sub Entringer($ where * > 0, 0 --> 0) { };

multi sub Entringer($n, $k) { %Entringer{"$n,$k"} //= Entringer($n, $k - 1) + Entringer($n - 1, $n - $k) };


# Per OEIS A000123
sub binpart($n) { $n ?? binpart($n-1) + binpart($n div 2) !! 1 }

our %BROKEN = %(A003418 => "dispatch error in core");

sub factorial($n) is export(:support) { ([*] 1..$n) or 1 }

sub factors($n is copy, :%map) is export(:support) {
    gather do {
        if %map{$n}:exists {
            take %map{$n};
        } elsif $n < 4 {
            take $n;
        } else {
            .take for prime-factors($n);

            ### Inline factoring code from Raku module Prime::Factor ########
            ### https://raku.land/zef:thundergnat/Prime::Factor
            ### Used with permission.

            sub prime-factors ( Int $n where * > 0 ) {
                return $n if $n.is-prime;
                return [] if $n == 1;
                my $factor = find-factor( $n );
                sort flat prime-factors( $factor ), prime-factors( $n div $factor );
            }

            sub find-factor ( Int $n, $constant = 1 ) {
                return 2 unless $n +& 1;
                # magic number below: product of primes 3 through 43
                if (my $gcd = $n gcd 6541380665835015) > 1 {
                    return $gcd if $gcd != $n
                }
                my $x      = 2;
                my $rho    = 1;
                my $factor = 1;
                while $factor == 1 {
                    $rho = $rho +< 1;
                    my $fixed = $x;
                    my int $i = 0;
                    while $i < $rho {
                        $x = ( $x * $x + $constant ) % $n;
                        $factor = ( $x - $fixed ) gcd $n;
                        last if 1 < $factor;
                        $i = $i + 1;
                    }
                }
                $factor = find-factor( $n, $constant + 1 ) if $n == $factor;
                $factor
            }

            ### End inlined code ###############################################
        }
    }
}

sub divisors($n) is export(:support) {
    gather do {
        take 1;
        if $n > 2 {
            for 2 .. ($n div 2) -> $i {
                take $i if $n %% $i;
            }
        }
        take $n if $n != 1;
    }
}

# Helper which fixes factors(1) to be empty
sub prime-factors($n) is export(:support) {
    factors($n).grep: {$^d ≥ 2}
}

# The prime signature of a number is the Bag of the positive
# exponents that appear in its prime factorization.
# https://en.wikipedia.org/wiki/Prime_signature
sub prime-signature($n --> Bag:D) is export(:support) {
    prime-factors($n).Bag\ # prime factorization bag
        .values.Bag        # exponent bag
}

sub sigma($n, $exponent=1) is export(:support) {
    [+] divisors($n).map: -> $j { $j ** $exponent };
}

sub planar-partitions($n) is export(:support) {
    if $n <= 1 {
        1;
    } else {
        ([+] (1..$n).map: -> $k {
            planar-partitions($n-$k) * sigma($k, 2)
        }) div $n;
    }
}

#| The number of sequences that sum to $target where elements
#| <= $n and elements in each sum are distinct.
sub strict-partitions(Int:D $n, Int :$target=$n) is export(:support) {
    return strict-partitions($n.abs, :$target).map({.map: {-$_}}) if $n < 0;
    return [$n] if $n==0 and not $target;
    $target > 0 or die "\$target($target) cannot be < 0";
    gather loop (my $i = min($n,$target); $i > 0; $i--) {
        given ($target)-$i -> $remain {
            when  $remain < 0 { next }
            when $remain == 0 { take [$i] }
            when      $i <= 1 { next } # only if $remain != 0 per above
            default {
                for strict-partitions($i-1, :target($remain)) -> $rest {
                    # Merge and add to results
                    take ($i.Array, $rest.Array).map: |*;
                }
            }
        }
    }
}

sub totient ($n) is export(:support) {
    +(^$n).grep: * gcd $n == 1
}

sub moebius ($n) is export(:support) {
    given $n.&prime-signature {
        when *.elems    == 0 { 1 #`{ constant one   } }
        when *.keys.max == 1 { .{1} %% 2 ?? 1 !! -1   }
        default              { 0 #`{ non-squarefree } }
    }
}

# The number of k-ary necklaces of length n.
sub necklaces ($n, :ary($k) = 2) {
    return 1 if $n == 0;
    $n R/ sum divisors($n).map: -> $d { totient($d) * $k**($n / $d) }
}

# https://mail.python.org/pipermail/edu-sig/2012-December/010721.html
sub Pi-digits is export(:support) {
    my ($q, $r, $t) = 1, 180, 60;
    gather for 2..* -> $j {
        my ($u,$y) =
            3*(3*$j+1)*(3*$j+2),
            ($q*(27*$j-12)+5*$r) div (5*$t);
        take $y;
        ($q, $r, $t) =
            10*$q*$j*(2*$j-1),
            10*$u*($q*(5*$j-2)+$r-$y*$t),
            $t*$u;
    }
}

# May as well, while we're here...
sub FatPi($digits=100) is export {
  FatRat.new(+([~] Pi-digits[^($digits)]), 10**($digits-1));
}

sub Eulers-number ( Int $terms = 500 ) is export(:support) {
    # Generates decimal digits of e accurate at least up to term.
    # Returns first 500 decimal digits by default as a trade-off
    # between completeness and run time.
    (sum map { FatRat.new(1,factorial($_)) }, ^(ceiling($terms * .66) max 100))\
    .substr(0, $terms+2).FatRat
}

# Stirling numbers of the first kind
multi Sterling1 (0, 0) is export(:support) { 1 }
multi Sterling1 (Int \n where * > 0, 0) is export(:support) { 0 }
multi Sterling1 ( 0, Int \k where * > 0) is export(:support) { 0 }
multi Sterling1 (Int \n, Int \k) is export(:support) {
    state %seen;
    (%seen{"{n - 1}|{k - 1}"} //= Sterling1(n - 1, k - 1)) -
    (n - 1) * (%seen{"{n - 1}|{k}"} //= Sterling1(n - 1, k))
}

# Stirling numbers of the second kind
multi Sterling2 (0, 0) is export(:support) { 1 }
multi Sterling2 (Int \n where * > 0, 0) is export(:support) { 0 }
multi Sterling2 (Int \n, Int \k where * == n) is export(:support) { 1 }
multi Sterling2 (Int \n, Int \k) is export(:support) {
    1/factorial(k) * sum (0 .. k).map: -> \j {
        (-1)**j * (k choose j) * (k - j)**n
    }
}

sub Horadam( Int $p, Int $q, Int $r, Int $s ) {
  [ $p, $q, {$^n1 × $r + $^n2 × $s} … ∞]
}


# If we don't yet have a formula for a given sequence, we use &NOSEQ in a
# range to define where our canned data ends. Because we use "fail", the
# failure only happens if you try to use a value past the end of our known
# entries.
sub NOSEQ { fail "This sequence has not yet been defined" }

# Needed for other sequences
our @A010051 is export = ℕ.map: { .is-prime ?? 1 !! 0 };
our @A080040 is export = 2, 2, -> $a, $b {2*$a + 2*$b} ... *;

# These are the "core" OEIS sequences as defined here:
# http://oeis.org/wiki/Index_to_OEIS:_Section_Cor
#

# groups
our @A000001 is export = 0, 1, 1, 1, 2, 1, 2, 1, 5, 2, 2, 1, &NOSEQ ... *;
# Kolakoski
our @A000002 is export = 1, 2, 2, -> $i {
    (state @a).push(|(((@a ?? @a[*-1] !! $i)%2+1) xx $i)); @a.shift} ... *;
# A000004 / 0's
our @A000004 is export = 0 xx *;
# A000005 / divisors
our @A000005 is export = ℕ.map: { divisors($^n).elems };
# A000007 / 0^n
our @A000007 is export = 𝕀.map: -> $n { 0 ** $n };
# A000009 / distinct partitions
our @A000009 is export = 𝕀.map: { strict-partitions($^i).elems };
# A000010 / totient
our @A000010 is export = ℕ.map: &totient;
# A000012 / 1's
our @A000012 is export = 1 xx *;
# A000014 / series-reduced trees
our @A000014 is export = 0, 1, 1, 0, 1, 1, 2, 2, 4, 5, 10, 14, &NOSEQ ... *;
# A000019 / prim. perm. groups
our @A000019 is export = 1, 1, 2, 2, 5, 4, 7, 7, 11, 9, 8, 6, &NOSEQ ... *;
# A000027 / natural numbers
our @A000027 is export = |ℕ; # We chose 𝕀[0]=0, OEIS chose 𝕀[0]=1
# A000029 / bracelets
our @A000029 is export = 𝕀.map: anon sub ($n) {
    return 1 if $n == 0;
    ½*necklaces($n) + ($n %% 2 ??
        3 * 2**($n/2 - 2) !!
        2**(($n - 1) / 2)
    )
}
# A000031 / necklaces
our @A000031 is export = 𝕀.map: &necklaces;
# A000032 / Lucas
our @A000032 is export = 2, 1, * + * ... *;
our @sequence-Lucas is export =  @A000032;
# A000035 / 0101...
our @A000035 is export = |(0,1) xx *;
# A000040 / primes
our @A000040 is export = lazy 𝕀.grep: {.is-prime};
# A000041 / partitions
our @A000041 is export = 1, 1, 2, 3, 5, 7, 11, 15, 22, 30, 42, &NOSEQ ... *;
# A000043 / Mersenne
our @A000043 is export = lazy 𝕀.grep: { .is-prime and (2**$_-1).is-prime };
# A000045 / Fibonacci
our @A000045 is export = 0, 1, * + * ... *;
# A000048 / necklaces-two-color-interchangable
our @A000048 is export = 𝕀.map: anon sub ($n) {
    return 1 if $n == 0;
    (2*$n) R/ sum divisors($n).grep(* !%% 2).map: -> $d { moebius($d) * 2**($n / $d) }
}
# A000055 / trees
our @A000055 is export = 1, 1, 1, 1, 2, 3, 6, 11, 23, 47, 106, &NOSEQ ... *;
# A000058 / Sylvester
our @A000058 is export = 2, { $_**2 - $_ + 1 } ... *;
# A000069 / odious
our @A000069 is export = lazy 𝕀.grep: -> $n { ([+] $n.base(2).comb) !%% 2 };
# A000079 / 2^n
our @A000079 is export = 1, * * 2 ... *;
# A000081 / rooted trees
our @A000081 is export = 0, 1, 1, 2, 4, 9, 20, 48, 115, 286, &NOSEQ ... *;
# A000085 / self-inverse perms.
our @A000085 is export = 1, 1, -> $a,$b { state $n = 1; $b+($n++)*$a } ... *;
# A000088 / graphs
our @A000088 is export = 1, 1, 2, 4, 11, 34, 156, 1044, 12346, &NOSEQ ... *;
# A000105 / polyominoes
our @A000105 is export = 1, 1, 1, 2, 5, 12, 35, 108, 369, 1285, &NOSEQ ... *;
# A000108 / Catalan
our @A000108 is export = lazy 𝕀.map: {(2*$^n choose $^n)/($^n+1)};
# A000109 / polyhedra
our @A000109 is export = 1, 1, 1, 2, 5, 14, 50, 233, 1249, 7595, 49566, 339722,
                         2406841, 17490241, 129664753, 977526957, 7475907149,
                         57896349553, 453382272049, 3585853662949, 28615703421545,
                         &NOSEQ ... *;
# A000110 / Bell
our @A000110 is export = lazy gather {
    my @bells = lazy [1], -> @b {
        my @c = @b.tail;
        @c.push: @b[$_] + @c[$_] for ^@b;
        @c
    } ... *;
    @bells.map: { take .head };
};
# A000111 / Euler
our @A000111 is export = lazy 𝕀.map: -> $n {Entringer($n, $n)};
# A000112 / posets
our @A000112 is export = 1, 1, 2, 5, 16, 63, 318, 2045, 16999, 183231, 2567284,
                         46749427, 1104891746, 33823827452, 1338193159771,
                         68275077901156, 4483130665195087, &NOSEQ ... *;
# A000120 / 1's in n
our @A000120 is export = lazy 𝕀.map: -> $n {$n.base(2).comb.grep({+$_}).elems};
# A000123 / binary partitions
our @A000123 is export = lazy 𝕀.map: &binpart;
# A000124 / Lazy Caterer
our @A000124 is export = lazy 𝕀.map: -> $n {($n * ($n+1)) / 2 + 1};
# A000129 / Pell
our @A000129 is export = 0, 1, * + 2 * * ... *;
our @Pell-sequence is export = @A000129;
# A000140 / Kendall-Mann
our @A000140 is export = 1, &NOSEQ ... *;
# A000142 / n!
our @A000142 is export = 𝕀.map: -> $n { factorial($n) };
# A000161 / partitions into 2 squares
our @A000161 is export = flat 1, {
    state $n++;
    my $k = (0 .. *).map({.²}).first: * >= $n, :k;
    my @sq = (0 ..^ $k).map({.²});
    my $cnt = ($n == $k²) ?? 1 !! 0;
    my %seen;
    for @sq {
        next if %seen{$_};
        if $n - $_ ∈ @sq {
            $cnt++;
            %seen{$n - $_} = True ;
        }
    }
    $cnt
}...*;
# A000166 / derangements
our @A000166 is export = lazy 1, -> $a {state $n++; $n*$a + (-1)**$n } ... *;
# A000169 / labeled rooted trees
our @A000169 is export = {state $n++; $n**($n - 1) } ... *;
# A000182 / tangent
our @A000182 is export = lazy (1, 3 ... *).map: -> $n {Entringer($n, $n)}
# A000203 / sigma
our @A000203 is export = ℕ.map: -> $n { sigma($n) };
# A000204 / Lucas
our @A000204 is export = 1, 3, *+* ... *;
# A000217 / triangular
our @A000217 is export = 𝕀.map: -> $n {($n*($n+1)) div 2};
# A000219 / planar partitions
our @A000219 is export = 𝕀.map: -> $n { planar-partitions($n) };
# A000225 / 2^n-1
our @A000225 is export = 𝕀.map: -> $n {2**$n-1};
# A000244 / 3^n
our @A000244 is export = 𝕀.map: -> $n {3**$n};
# A000262 / sets of lists
our @A000262 is export = 1, 1, -> $a, $b {state $n++; (2*($n+1)-1)*$b - $n*($n - 1) * $a } ... *;
# A000272 / n^(n-2)
our @A000272 is export = 𝕀.map: -> $n {$n ?? $n**($n-2) !! 1};
# A000273 / directed graphs
our @A000273 is export = 1, &NOSEQ ... *;
# A000290 / n^2
our @A000290 is export = 𝕀.map: -> $n {$n**2};
# A000292 / tetrahedral
our @A000292 is export = 𝕀.map: -> $n { ($n*($n+1)*($n+2)) div 6 };
# A000302 / 4^n
our @A000302 is export = 𝕀.map: -> $n {4**$n}
# A000311 / Schroeder's fourth
our @A000311 is export = 0, 1, -> $a {
    state $n = 1;
    ++$n;
    -($n-1) * $a + sum (1 .. ($n - 1)).map: -> $k {
        ($n choose $k) * @A000311[$k] *@A000311[$n - $k]
    }
} ... *;
# A000312 / mappings
our @A000312 is export = lazy 1, {state $n++; $n ** $n } ... *;
# A000326 / pentagonal
our @A000326 is export = lazy 0, {state $n++; $n*(3*$n-1)/2 } ... *;
# A000330 / square pyramidal
our @A000330 is export = lazy 0, {state $n++; $n*($n+1)*(2*$n+1)/6 } ... *;
# A000364 / Euler or secant
our @A000364 is export = lazy (0, 2 ... *).map: -> $n {Entringer($n, $n)}
# A000396 / perfect
our @A000396 is export = lazy gather for @A000040 #`{primes} {
                             my $n = 2**$_ - 1;
                             take $n * 2**($_ - 1) if $n.is-prime;
                         }
# A000521 / j
our @A000521 is export = 1, &NOSEQ ... *;
# A000578 / n^3
our @A000578 is export = 𝕀.map: -> $n {$n ** 3}
# A000583 / n^4
our @A000583 is export = 𝕀.map: -> $n {$n ** 4}
# A000593 / sum odd divisors
our @A000593 is export = {state $n++; sum $n.&divisors.grep: * % 2 } ... *;
# A000594 / Ramanujan tau
our @A000594 is export = 1, {
    (state $i = 1)++;
    my $s = 0;
    my $t = 1;
    my $u = 0;
    for 1..$i -> $j {
        $t += 9 * $j;
        $u += $j;
        last if $i <= $u;
        $s += (-1) ** ($j % 2 + 1) * (2 * $j + 1) * ($i - $t) * @A000594[$i-$u-1];
    }
    $s / ($i - 1);
} ... *;
# A000602 / hydrocarbons
our @A000602 is export = 1, &NOSEQ ... *;
# A000609 / threshold functions
our @A000609 is export = 2, 4, 14, 104, 1882, 94572, 15028134,
                         8378070864, 17561539552946, 144130531453121108,
                         &NOSEQ ... *;
# A000670 / preferential arrangements
our @A000670 is export = 1, &NOSEQ ... *;
# A000688 / abelian groups
our @A000688 is export = (1..*).map: {
    # @A000041 NYI. Hardcoded list, will fail at term 2 ** 173526
    state @a = <1 1 2 3 5 7 11 15 22 30 42 56 77 101 135 176 231 297 385 490 627
          792 1002 1255 1575 1958 2436 3010 3718 4565 5604 6842 8349 10143 12310
          14883 17977 21637 26015 31185 37338 44583 53174 63261 75175 89134
          105558 124754 147273 173525>».Int;
    my @f = .&factors.Bag.grep( *.value > 1 )».value || 0;
    [*] @f.map: { @a[$_] }
};
# A000720 / pi(n)
our @A000720 is export = [\+] @A010051;
# A000793 / Landau
our @A000793 is export = 1, { strict-partitions(++$).map({[lcm] $_}).max } ... *;
# A000796 / Pi
our