### ----------------------------------------------------
### -- Math::NumberTheory
### -- Licenses: Artistic-2.0
### -- Authors: Anton Antonov
### -- File: lib/Math/NumberTheory.rakumod
### ----------------------------------------------------
### -- Chunk 2 of 2
 number of positive integers
# less than or equal to a number x which are not divisible
# by any of the first a primes.
# https://mathworld.wolfram.com/LegendresFormula.html

# Speeding up computations
state %legendre-formula-cache;
# Using the recurrence relation
sub legendre-formula(Numeric:D $x) {
    sub phi($x, $a) {
        return $x if $a == 0;
        return %legendre-formula-cache{$x}{$a} //= phi($x, $a - 1) - phi(($x div prime($a)), $a - 1);
    }

    sub pi($n) {
        return 0 if $n < 2;
        my $a = pi(($n.sqrt).floor);
        return phi($n, $a) + $a - 1;
    }

    return pi($x.floor);
}

#==========================================================
# Prime omega
#==========================================================
#| Give the number of prime factors counting multiplicities Ω(n) in the argument.
multi sub prime-omega($x) is export {*}
multi sub prime-omega(@x) {
    return @x.map({ prime-omega($_) }).List;
}
multi sub prime-omega(Int:D $x) {
    return $x == 1 ?? 0 !! factor-integer($x)».tail.sum;
}

multi sub prime-omega(Complex:D $x) {
    return $x == 1 ?? 0 !! factor-gaussian-integer($x, :!include-unit)».tail.sum;
}

#==========================================================
# Prime nu
#==========================================================
#| Give the number of distinct primes in the argument.
multi sub prime-nu($x) is export {*}
multi sub prime-nu(@x) {
    return @x.map({ prime-nu($_) }).List;
}
multi sub prime-nu(Int:D $x) {
    return $x == 1 ?? 0 !! factor-integer($x).elems;
}

multi sub prime-nu(Complex:D $x) {
    return $x == 1 ?? 0 !! factor-gaussian-integer($x, :!include-unit).elems;
}

#==========================================================
# MoebiusMu function
#==========================================================
# http://reference.wolfram.com/language/ref/MoebiusMu.html

# Give the Möbius function µ(n).
proto sub moebius-mu($n, Bool:D :gaussian(:$gaussian-integers) = False) is export {*}

multi sub moebius-mu(@n, Bool:D :gaussian(:$gaussian-integers) = False) {
    return @n.map({ moebius-mu($_, :$gaussian-integers) }).List;
}

multi sub moebius-mu(Int:D $n, Bool:D :gaussian(:$gaussian-integers) = False) {
    return moebius-mu($n + 0i) if $gaussian-integers;
    return 1 if $n == 1;
    my @exps = factor-integer($n)».tail;
    return @exps.max == 1 ?? (-1) ** @exps.elems !! 0;
}

multi sub moebius-mu(Complex:D $n, Bool:D :gaussian(:$gaussian-integers) = False) {
    my @exps =  factor-gaussian-integer($n, :!include-unit)».tail;
    return @exps.max == 1 ?? (-1) ** @exps.elems !! 0;
}

#==========================================================
# SquareFreeQ function
#==========================================================
# http://reference.wolfram.com/language/ref/SquareFreeQ.html

# Give True if the first argument is a square-free number, and False otherwise.
proto sub is-square-free($n, Bool:D :gaussian(:$gaussian-integers) = False) is export {*}

multi sub is-square-free(@n, Bool:D :gaussian(:$gaussian-integers) = False) {
    return @n.map({ is-square-free($_, :$gaussian-integers) }).List;
}

multi sub is-square-free(Int:D $n, Bool:D :gaussian(:$gaussian-integers) = False) {
    return moebius-mu($n, :$gaussian-integers) != 0 ?? True !! False;
}

multi sub is-square-free(Complex:D $n, Bool:D :gaussian(:$gaussian-integers) = False) {
    return moebius-mu($n) != 0 ?? True !! False;
}

#==========================================================
# LiouvilleLambda function
#==========================================================
# http://reference.wolfram.com/language/ref/LiouvilleLambda.html

# Gives the Liouville lambda function λ(n).
proto sub liouville-lambda($n, Bool:D :gaussian(:$gaussian-integers) = False) is export {*}

multi sub liouville-lambda(@n, Bool:D :gaussian(:$gaussian-integers) = False) {
    return @n.map({ liouville-lambda($_, :$gaussian-integers) }).List;
}

multi sub liouville-lambda(Int:D $n, Bool:D :gaussian(:$gaussian-integers) = False) {
    return liouville-lambda($n + 0i) if $gaussian-integers;
    return (-1) ** prime-omega($n);
}

multi sub liouville-lambda(Complex:D $n, Bool:D :gaussian(:$gaussian-integers) = False) {
    return (-1) ** prime-omega($n);
}

#==========================================================
# CarmichaelLambda function
#==========================================================
# http://reference.wolfram.com/language/ref/CarmichaelLambda.html
# https://en.wikipedia.org/wiki/Carmichael_function#Recurrence_for_λ(n)

#| Give the Carmichael function λ(n).
#| C($n) -- An integer.
proto sub carmichael-lambda($n) is export {*}

multi sub carmichael-lambda(@n) {
    return @n.map({ carmichael-lambda($_) }).List;
}
multi sub carmichael-lambda(Int:D $n is copy) {
    return 0 if $n == 0;

    $n = abs($n);

    return 1 if $n == 1;

    my @factors = factor-integer(abs($n));

    my @lambdas = do for @factors -> $p {
        my $r = ($p.head - 1) * $p.head ** ($p.tail - 1);
        if $p.head == 2 && $p.tail ≥ 3 {
           $r /= 2
        }
        $r
    }
    return [lcm] @lambdas;
}

#==========================================================
# Kronecker delta
#==========================================================
#| Gives the Kronecker delta value, which is equal to 1 if all the arguments are equal, and 0 otherwise.
proto sub kronecker-delta(**@n) is export {*}

multi sub kronecker-delta() {
    return 1;
}
multi sub kronecker-delta(Numeric:D $n) {
    return $n == 0 ?? 1 !! 0;
}

multi sub kronecker-delta(+@n) {
    return kronecker-delta() if @n.elems == 0;
    return kronecker-delta(@n.head) if @n.elems == 1;
    return ([&&] @n.tail(*-1).map({ $_ == @n.head })) ?? 1 !! 0;
}

#==========================================================
# Simple extensions
#==========================================================
#| Gives true if the number is not prime.
proto sub is-composite($n) is export {*}
multi sub is-composite(Complex:D $n) { !is-prime-gaussian($n) }
multi sub is-composite(Int:D $n) { !is-prime($n) }

#| Gives true the number is a power of a prime.
sub is-prime-power(Int:D $n) is export {
    factor-integer($n).elems == 1;
}

#| Give Mangold lambda for an integer
sub mangold-lambda(Int:D $n) is export {
    my @res = factor-integer($n);
    return @res.elems == 1 ?? @res.head.head.log !! 0;
}

#==========================================================
# Coprime numbers
#==========================================================
#| Give true if the arguments are coprime.
proto sub are-coprime(|) is export {*}

multi sub are-coprime($n1, $n2) {
    ($n1 gcd $n2) == 1
}

multi sub are-coprime(*@ns) {
    die 'All elements are expected to be numbers.'
    unless @ns.all ~~ Numeric:D;

    for @ns.kv -> $i, $n1 {
        for @ns.kv -> $j, $n2 {
            if $i != $j {
                return False if $n1 gcd $n2 > 1
            }
        }
    }
    return True;
}

#==========================================================
# Twin, cousin, and sexy primes
#=========================================================

sub related-primes(UInt:D $n, UInt:D :$step = 2) is export {
    # This is slow, but it would be nice if it is made to run faster.
    #my @ps = (1..(2 * $n * log($n)).floor)».&prime;

    my @ps = primes-up-to(floor(2 * $n * log($n)));
    my @res = (@ps (&) (@ps >>+>> $step)).keys.sort;
    @res = @res[^(min($n, @res.elems))];
    return @res.map({ ($_ - $step, $_) }).List;
}

#| Get the first n-pairs of primes that differ by 2.
sub twin-primes(UInt:D $n) is export {
    related-primes($n, :2step);
}

#| Get the first n-pairs of primes that differ by 4.
sub cousin-primes(UInt:D $n) is export {
    related-primes($n, :4step);
}

#| Get the first n-pairs of primes that differ by 6.
sub sexy-primes(UInt:D $n) is export {
    related-primes($n, :6step);
}

#==========================================================
# Integer partitions
#==========================================================
sub accel-asc(Int:D $n) {
    my @a = 0 xx ($n + 1);
    my $k = 1;
    my $y = $n - 1;

    gather {
        while $k != 0 {
            my $x = @a[$k - 1] + 1;
            $k--;

            while 2 * $x <= $y {
                @a[$k] = $x;
                $y -= $x;
                $k++;
            }

            my $l = $k + 1;
            while $x <= $y {
                @a[$k] = $x;
                @a[$l] = $y;
                take @a[0..($k + 1)];
                $x++;
                $y--;
            }

            @a[$k] = $x + $y;
            $y = $x + $y - 1;
            take @a[0..$k];
        }
    }
}

#| Give a list of all possible ways to partition the integer argument into smaller integers.
proto sub integer-paritions(Int:D $n, |) is export {*}

multi sub integer-paritions(Int:D $n, UInt:D $k) {
    integer-paritions($n, k-max => $k);
}

multi sub integer-paritions(Int:D $n, (Numeric:D $k-min, Numeric:D $k-max)) {
    integer-paritions($n, :$k-min, :$k-max);
}

multi sub integer-paritions(Int:D $n, UInt:D :$k-min = 1, Numeric:D :$k-max = Inf) {
    my @res = accel-asc($n);
    @res .= grep({ $k-min ≤ $_.elems ≤ $k-max });
    return @res».reverse».List.reverse.List;
}

#==========================================================
# Random prime
#==========================================================
# http://reference.wolfram.com/language/ref/RandomPrime.html
#| Give pseudo random prime numbers in specified ranges.
#| C<:$range> -- Integer or Range.
#| C<:$n> -- Number of primes to return.
proto sub random-prime(|) is export {*}
multi sub random-prime(Int:D $max, $n = Whatever) {
    die 'If the first argument is a number then it is expected to be an integer greater than 1.'
    unless $max > 1;
    return random-prime(2 .. $max, $n);
}

multi sub random-prime(Range:D $range, $n is copy = Whatever) {
    my ($min, $max) = $range.head, $range.tail;

    die 'Correct range argument is expected.'
    unless $min.defined && $max.defined;

    die 'The start of the range argument is expected to start with an integer greater than 1.'
    unless $min > 1;

    die 'The end of the range argument is expected to be greater than the range start.'
    unless $min < $max;

    die 'The second argument is expected to be a positive integer or Whatever.'
    unless $n ~~ Int:D && $n > 0 || $n.isa(Whatever);

    my $i = 0;
    my @res;
    my $n2 = $n.isa(Whatever) ?? 1 !! $n;
    for $min .. $max -> $c {
        if $c.is-prime {
            @res.push($c);
            $i++
        }
        last if $i ≥ $n2
    }

    return do if $n.isa(Whatever) {
        @res.head
    } elsif @res.elems == $n2 {
        @res.List
    } else {
        @res.roll($n2).List
    }
}

#==========================================================
# Real digits
#==========================================================
# http://reference.wolfram.com/language/ref/RealDigits.html
#| Real digits
#| C<$x> -- Number to convert.
#| C<:$base> -- Conversion base.
#| C<:$n> -- Digit exponent to start with.
#| C<:$tol> -- Tolerance to stop the conversion with.
#| C<:$length> -- Max number of digits.
proto sub real-digits(Numeric:D $x, *@args, *%args) is export {*}

multi sub real-digits($x, Numeric:D $b = 10, $n = Whatever, *%args) {
    return real-digits($x ~~ Int:D ?? $x.FatRat !! $x, :$b, :$n, |%args);
}

multi sub real-digits(Numeric:D $x is copy,
                      Numeric:D :base(:$b) = 10,
                      :start(:first-digit-exponent(:$n)) is copy = Whatever,
                      Numeric:D :tolerance(:$tol) = 10e-14,
                      :l(:len(:$length)) is copy = Inf) {
    if $x == 0 { return 0, 0 }
    $x = abs($x);
    my @digits;
    my $current = $x;

    if $length.isa(Whatever) { $length = Inf }
    die 'The argument $length is expected to be a positive integer or Whatever.'
    unless ($length === Inf || $length ~~ Int:D) && $length > 0;

    if $n.isa(Whatever) { $n = $x.log($b).floor }
    die 'The third argument is expected to be a number or Whatever.' unless $n ~~ Numeric:D;

    if $x ~~ FatRat:D {
        my $exp = $n.FatRat;
        my $bf = $b.FatRat;

        my $bfe = [*] ($bf xx $exp);
        while $current / $x > $tol && @digits.elems < $length {
            my $r = $current / $bfe;
            my $digit = do if $r.round(10 ** -100) == 1 { 1 } else { ($current / $bfe).floor };
            @digits .= push($digit);
            $current -= $digit.FatRat * $bfe;
            $exp--;
            $bfe = $bfe / $bf;
        }
    } else {
        my $exp = $n;
        my $bf = $b;

        my $bfe = [*] ($bf xx $exp);
        while $current / $x > $tol && @digits.elems < $length {
            my $digit = floor($current / $bfe);
            @digits .= push($digit);
            $current -= $digit * $bfe;
            $exp--;
            $bfe = $bfe / $bf;
        }
    }

    return @digits, $n + 1;
}

#==========================================================
# Phi number system
#==========================================================
# https://resources.wolframcloud.com/FunctionRepository/resources/PhiNumberSystem/
# https://mathworld.wolfram.com/PhiNumberSystem.html

# Using N[Sqrt[5], 100] in Wolfram Language
constant $sqrt5 = 2.236067977499789696409173668731276235440618359611525724270897245410520925637804899414414408378782275.FatRat;

# Fibonacci 401 and 400
constant $fibonacci401 = 284812298108489611757988937681460995615380088782304890986477195645969271404032323901.FatRat;
constant $fibonacci400 = 176023680645013966468226945392411250770384383304492191886725992896575345044216019675.FatRat;

#| Golden ratio (phi)
our constant \phi is export = $fibonacci401 / $fibonacci400;
our constant \ϕ is export = $fibonacci401 / $fibonacci400;
#our constant \ϕ is export = (1.FatRat + $sqrt5) / 2.FatRat;
sub golden-ratio(Bool:D :pre(:$pre-computed) = True) is export {
    return $pre-computed ?? \phi !! (1 + sqrt(5)) / 2;
}

#| Phi number system.
#| C<$n> -- An integer number to convert.
#| C<:$tol> -- Tolerance to stop the conversion with.
#| C<:$length> -- Max number of digits.
sub phi-number-system(Int:D $n, Numeric:D :tolerance(:$tol) = 10e-16, :l(:len(:$length)) is copy = Whatever) is export {
    if $length.isa(Whatever) { $length = 2 * ($sqrt5 * abs($n)).log(ϕ).floor + 1; }
    my ($digits, $exp) = real-digits($n.FatRat, ϕ, :$tol, :$length);
    return $exp <<->> ($digits.grep(*== 1, :k) >>+>> 1);
}


