
=encoding utf8

=head1 NAME

Sidef::Types::Number::Number

=head1 DESCRIPTION

The C<Number> class implements support for numerical operations, supporting integers, rationals, floating-points and complex numbers at arbitrary precision.

This class also implements many useful mathematical methods, from basic arithmetical operations, to advanced number-theoretic functions, including primality testing and prime factorization methods.

=head1 SYNOPSIS

    var a = Num(string)
    var b = Number(string, base)

=head1 INHERITS

Inherits methods from:

       * Sidef::Object::Object

=head1 METHODS

=head2 !

    n!

Factorial of C<n>. (C<1*2*3*...*n>)

Aliases: I<fac>, I<factorial>

=cut

=head2 !!

    n!!

Double-factorial of C<n>.

Aliases: I<dfac>, I<dfactorial>, I<double_factorial>

=cut

=head2 %

    n % k

Remainder of C<n/k>.

Aliases: I<mod>

=cut

=head2 %%

    n %% k

Returns true if C<n> is divisible by C<k>. False otherwise.

Aliases: I<is_div>

=cut

=head2 &

    a & b

Bitwise AND operation.

Aliases: I<and>

=cut

=head2 *

    a * b

Multiplication of C<a> and C<b>.

Aliases: I<mul>

=cut

=head2 **

    a**b

Exponentiation: C<a> to power C<b>.

Aliases: I<pow>

=cut

=head2 +

    a + b

Addition of C<a> and C<b>.

Aliases: I<add>

=cut

=head2 ++

    n++
    ++n
    n.inc

Increment C<n> by C<1> and return the result.

Aliases: I<inc>

=cut

=head2 -

    a - b

Subtraction of C<a> and C<b>.

Aliases: I<sub>

=cut

=head2 --

    n--
    --n
    n.dec

Decrement C<n> by C<1> and return the result.

Aliases: I<dec>

=cut

=head2 ..

    a .. b

Create an inclusive-inclusive C<RangeNum> object, from C<a> to C<b>.

Equivalent with:

    RangeNum(a, b)

Aliases: I<to>, I<upto>

=cut

=head2 ..^

    a ..^ b

Create an inclusive-exclusive C<RangeNum> object, from C<a> to C<b-1>.

Equivalent with:

    RangeNum(a, b-1)

Aliases: I<xto>, I<xupto>

=cut

=head2 /

    a / b

Division of C<a> by C<b>.

Aliases: I<÷>, I<div>

=cut

=head2 //

    a // b

Integer floor-division of C<a> and C<b>.

When C<a> and C<b> are integers, this is equivalent with:

    floor(a/b)

Aliases: I<idiv>

=cut

=head2 :

    a : b

Create a new complex number.

Equivalent with:

    Complex(a, b)

Aliases: I<pair>

=cut

=head2 <

    a < b

Returns true if C<a> is less than C<b>.

Aliases: I<lt>

=cut

=head2 <<

    a << b

Left shift C<a> by C<b> bits, which is equivalent with (assuming C<a> and C<b> are integers):

    floor(a * 2**b)

Aliases: I<lsft>, I<shift_left>

=cut

=head2 <=>

    a <=> b

Comparison of C<a> with C<b>. Returns C<-1> if C<a> is less than C<b>, C<0> if C<a> and C<b> are equal and C<+1> if C<a> is greater than C<b>.

Aliases: I<cmp>

=cut

=head2 approx_cmp

    a <~> b
    approx_cmp(a, b)
    approx_cmp(a, b, k)

Approximate comparison of C<a> and C<b>.

Equivalent with:

    a.round(k) <=> b.round(k)

When C<k> is omitted, it uses the default floating-point precision to deduce C<k>.

=cut

=head2 ==

    a == b

Equality check. Returns true if C<a> and C<b> are equal.

Aliases: I<eq>

=cut

=head2 >

    a > b

Returns true if C<a> is greater than C<b>. False otherwise.

Aliases: I<gt>

=cut

=head2 >>

    a >> b

Right shift C<a> by C<b> bits, which is equivalent with (assuming C<a> and C<b> are integers):

    floor(a / 2**b)

Aliases: I<rsft>, I<shift_right>

=cut

=head2 ^

    a ^ b

Bitwise XOR operation.

Aliases: I<xor>

=cut

=head2 ^..

    a ^.. b

Creates a reversed exclusive-inclusive C<RangeNum> object, from C<a-1> down to C<b>.

Equivalent with:

    RangeNum(a-1, b, -1)

Aliases: I<xdownto>

=cut

=head2 C

    Num.C
    Num.catalan_G

Returns the Catalan constant: 0.915965594177...

Aliases: I<catalan_G>

=cut

=head2 Y

    Num.Y()

Returns the Euler–Mascheroni constant: 0.5772156649...

Aliases: I<γ>, I<euler_gamma>

=cut

=head2 |

    a | b

Returns the bitwise OR or C<a> and C<b>.

Aliases: I<or>

=cut

=head2 ~

    ~a
    a.not

Returns the bitwise NOT of C<a>.

Aliases: I<not>

=cut

=head2 Γ

    Γ(n)
    gamma(n)

The gamma function, where for any positive integer C<n>:

    gamma(n) = (n-1)!

Aliases: I<gamma>

=cut

=head2 δ

    δ(a,b)

The Kronecker delta function, which returns 1 iff C<a==b> and 0 otherwise.

Aliases: I<kronecker_delta>

=cut

=head2 ζ

    ζ(s)
    zeta(s)

The Riemann Zeta function.

Aliases: I<zeta>

=cut

=head2 η

    η(s)
    eta(s)

The Dirichlet eta function, defined as:

    eta(s) = (1 - 2^(1-s)) * zeta(s)

Aliases: I<eta>

=cut

=head2 μ

    μ(n)
    mu(n)
    moebius(n)

The Möbius function: μ(n).

Aliases: I<mu>, I<mobius>, I<möbius>, I<moebius>

=cut

=head2 Π

    Π(...)
    prod(...)

Returns the product of a given list of numbers.

Aliases: I<prod>

=cut

=head2 π

    π(n)
    π(a,b)
    Num.π

Returns the PI numerical value:

    say Num.pi      #=> 3.1415...

When applied on a Number object (as C<n.pi> or C<pi(n)>), it returns the number of primes <= n:

    say 100.pi      #=> number of primes <= 100
    say pi(100)     #=> 25

When an additional argument is given, it returns the number of primes in the range C<a..b>:

    say pi(50, 100)     # number of primes in the range 50..100

Aliases: I<pi>

=cut

=head2 Σ

    Σ(...)
    sum(...)

Returns the sum of a given list of numbers.

Aliases: I<sum>

=cut

=head2 σ

    σ(n,k=1)
    sigma(n,k=1)

Returns the sum of the positive divisors of n, each raised to the power k.

Aliases: I<sigma>

=cut

=head2 τ

    Num.τ
    Num.tau

    τ(n)
    tau(n)

Returns the TAU constant (2*PI), or the number of positive divisors of n.

Example:

    say Num.tau     #=> 6.283185307179...
    say tau(120)    #=> 16

Aliases: I<tau>

=cut

=head2 φ

    Num.φ
    Num.phi

    φ(n)
    phi(n)

Returns the golden ratio constant PHI, or the Euler totient of n.

Example:

    say Num.phi       #=> 1.618033988749...
    say phi(180)      #=> 48

Aliases: I<phi>

=cut

=head2 Ψ

    Ψ(x)
    digamma(x)

The digamma function, defined as:

     Γ'(x)
    -------
     Γ(x)

Aliases: I<digamma>

=cut

=head2 Ω

    Ω(n,k=0)
    bigomega(n,k=0)

For k = 0, it returns the number of prime factors of n (counted with multiplicity).

In general, is equivalent with:

    n.prime_power_divisors.sum {|d| n**k / d**k }

Aliases: I<bigomega>, I<prime_power_sigma0>

=cut

=head2 ω

    ω(n,k=0)
    omega(n,k=0)

For k = 0, it returns the number of distinct prime factors of n.

In general, is equivalent with:

    n.prime_divisors.sum {|d| n**k / d**k }

Aliases: I<omega>, I<prime_sigma0>, I<prime_power_usigma0>

=cut

=head2 ≅

    a ≅ b

Returns true if C<a> and C<b> are approximately equal to each other.

Aliases: I<=~=>, I<approx_eq>

=cut

=head2 ≠

    a ≠ b

Returns true if C<a> and C<b> are different from each other.

Aliases: I<!=>, I<ne>

=cut

=head2 ≤

    a ≤ b

Returns true if C<a> is less than or equal to C<b>.

Aliases: I<E<lt>=>, I<le>

=cut

=head2 ≥

    a ≥ b

Returns true if C<a> is greather than or equal to C<b>.

Aliases: I<E<gt>=>, I<ge>

=cut

=head2 abs

    n.abs

The absolute value of C<n>.

=cut

=head2 abundancy

    abundancy(n)

Returns the abundancy index of n, defined as:

    sigma(n)/n

Aliases: I<abundancy_index>

=cut

=head2 acmp

    acmp(a,b)

Absolute comparison of C<a> and C<b>, defined as:

    abs(a) <=> abs(b)

=cut

=head2 acos

    n.acos

Inverse cosine of C<n> in radians.

=cut

=head2 acosh

    n.acosh

Inverse hyperbolic cosine of C<n>.

=cut

=head2 acot

    n.acot

Inverse cotangent of C<n> in radians.

=cut

=head2 acoth

    n.acoth

Inverse hyperbolic cotangent of C<n>.

=cut

=head2 acsc

    n.acsc

Inverse cosecant of C<n> in radians.

=cut

=head2 acsch

    n.acsch

Inverse hyperbolic cosecant of C<n>.

=cut

=head2 addmod

    addmod(a, b, m)

Modular integer addition: C<(a+b) % m>.

Example:

    say addmod(43, 97, 127)     # == (43+97)%127

=cut

=head2 agm

    agm(a, b)

Arithmetic-geometric mean of C<a> and C<b>.

=cut

=head2 ai

    x.ai
    Ai(x)

Airy function of the first kind: C<Ai(x)>.

Aliases: I<airy>

=cut

=head2 all_composite

    all_composite(...)

Returns true if all the given values are composite positive integers, by checking first for small factors, then running a B-PSW primality test.

=cut

=head2 all_prime

    all_prime(...)

Returns true if all the given values are prime numbers.

This is done by first running a primality pre-test on all values and returning early if one of the values is composite with a small factor. Otherwise, we run a B-PSW primality test on each value.

=cut

=head2 almost_prime_divisors

    n.almost_prime_divisors
    n.almost_prime_divisors(k)

Returns the k-almost prime divisors of n.

Example:

    say 5040.almost_prime_divisors(7)   #=> [720, 1008, 1680, 2520]

When C<k> is omitted, an array of arrays with the k-almost prime divisors of n, for each k in the range C<0..bigomega(n)>, is returned:

    say 120.almost_prime_divisors       #=> [[1], [2, 3, 5], [4, 6, 10, 15], [8, 12, 20, 30], [24, 40, 60], [120]]

=cut

=head2 almost_primes

    k.almost_primes(n)
    k.almost_primes(a,b)

Return an array with the k-almost primes <= n, or in the range C<a..b>.

Example:

    5.almost_primes(1e6)        # array of 5-almost primes <= 1e6
    5.almost_primes(1e5, 1e6)   # array of 5-almost primes in the range [1e5, 1e6]

=cut

=head2 approx_ge

    approx_ge(a, b)
    approx_ge(a, b, k)

True if C<a> is approximately greater than or equal to C<b>.

Equivalent with:

    a.round(k) >= b.round(k)

=cut

=head2 approx_gt

    approx_gt(a, b)
    approx_gt(a, b, k)

True if C<a> is approximately greater than C<b>.

Equivalent with:

    a.round(k) > b.round(k)

=cut

=head2 approx_le

    approx_le(a, b)
    approx_le(a, b, k)

True if C<a> is approximately less than or equal to C<b>.

Equivalent with:

    a.round(k) <= b.round(k)

=cut

=head2 approx_lt

    approx_lt(a, b)
    approx_lt(a, b, k)

True if C<a> is approximately less than C<b>.

Equivalent with:

    a.round(k) < b.round(k)

=cut

=head2 approx_ne

    approx_ne(a, b)
    approx_ne(a, b, k)

True if C<a> is approximately different than C<b>.

Equivalent with:

    a.round(k) != b.round(k)

=cut

=head2 as_bin

    n.as_bin

Returns a C<String> with the binary representation of C<n>.

Example:

    say 42.as_bin     # "101011"

=cut

=head2 as_dec

    n.as_dec
    n.as_dec(k)

Given a rational number C<n>, it returns its decimal expansion as a C<String> object, expanded at C<k> decimal places.

Example:

    say (1/17 -> as_dec(10))      # 0.05882352941
    say (1/17 -> as_dec(30))      # 0.0588235294117647058823529411765

Aliases: I<as_float>

=cut

=head2 asec

    n.asec

Inverse secant of C<n> in radians.

=cut

=head2 asech

    n.asech

Inverse hyperbolic secant of C<n>.

=cut

=head2 as_frac

    n.as_frac
    n.as_frac(base)

String-representation of C<n> as fraction.

Example:

    say 24.as_frac                 # 24/1
    say bernoulli(10).as_frac      # 5/66
    say bernoulli(12).as_frac(36)  # -j7/23u

If C<n> is an integer, it uses C<1> for the denominator.

=cut

=head2 as_hex

    n.as_hex

Returns a String representing the integer part of C<n> in hexadecimal (base 16).

Example:

    say 42.as_hex       # "2a"

Returns C<nil> when C<n> cannot be converted to an integer.

=cut

=head2 asin

    n.asin

Inverse sine of C<n> in radians.

=cut

=head2 asinh

    n.asinh

Inverse hyperbolic sine of C<n>.

=cut

=head2 as_int

    n.as_int
    n.as_int(base)

Returns a String representing the integer part of C<n> in a given base, where the base must be between 2 and 62.

When the base is omitted, it defaults to base 10.

Example:

    say 255.as_int      # "255"
    say 255.as_int(16)  # "ff"

Returns C<nil> when C<n> cannot be converted to an integer.

=cut

=head2 as_oct

    n.as_oct

Returns a String representing the integer part of C<n> in octal (base 8).

Example:

    say 42.as_oct   # 52

Returns C<nil> when C<n> cannot be converted to an integer.

=cut

=head2 as_rat

    n.as_rat
    n.as_rat(base)

Returns a rational string-representation of C<n> in a given base, where the base must be between 2 and 62.

When the base is omitted, it defaults to base 10.

Example:

    say as_rat(42)          # "42"
    say as_rat(2/4)         # "1/2"
    say as_rat(255, 16)     # "ff"

Returns C<nil> when C<n> cannot be converted to a rational number.

=cut

=head2 atan

    n.atan

Inverse tangent of C<n> in radians.

=cut

=head2 atan2

    atan2(a, b)

Four-quadrant inverse tangent of C<a> and C<b>.

=cut

=head2 atanh

    n.atanh

Inverse hyperbolic tangent of C<n>.

=cut

=head2 base

    n.base(b)

Returns a String-representation of C<n> in a given base C<b>, which must be between 2 and 62.

Aliases: I<in_base>

=cut

=head2 bell

    n.bell

Returns the n-th Bell number.

Aliases: I<bell_number>

=cut

=head2 bern

    n.bern
    bernoulli(n)
    bernoulli(n, x)

Returns the n-th Bernoulli number. When an additional argument is provided, it returns the n-th Bernoulli polynomial evaluated at C<x>.

Example:

    say bernoulli(10).as_rat       # B_10    = 5/66
    say bernoulli(10, 2).as_rat    # B_10(2) = 665/66

Aliases: I<bernfrac>, I<bernoulli>, I<bernoulli_number>

=cut

=head2 bernoulli_polynomial

    bernoulli_polynomial(n)
    bernoulli_polynomial(n, x)

Returns the n-th Bernoulli polynomial: C<B_n(x)>.

When C<x> is ommited, a Polynomial object is returned.

=cut

=head2 bernreal

    n.bernreal

Return an approximation to the n-th Bernoulli number as a floating-point number.

=cut

=head2 bessel_j

    bessel_j(x, n)

First order Bessel function: C<J_n(x)>.

=cut

=head2 bessel_y

    bessel_y(x, n)

Second order Bessel function: C<Y_n(x)>.

=cut

=head2 beta

    beta(a, b)

The beta function (also called the Euler integral of the first kind).

Defined as:

    beta(a, b) = gamma(a)*gamma(b) / gamma(a+b)

=cut

=head2 bit

    n.bit(k)
    n.getbit(k)

Returns 1 if bit C<k> of C<n> is set, and 0 if it is not set.

Return C<nil> when C<n> cannot be truncated to an integer or when C<k> is negative.

Example:

    say getbit(0b1001, 0)   # 1
    say getbit(0b1000, 0)   # 0

Aliases: I<getbit>, I<testbit>

=cut

=head2 bits

    n.bits()

Returns the binary digits of n.

The bits are ordered from the most significant bit to the least significant bit.

Example:

    say 1234.bits       #=> [1, 0, 0, 1, 1, 0, 1, 0, 0, 1, 0]

Equilvalent with:

    n.digits(2).flip

=cut

=head2 bit_scan0

    n.bit_scan0
    n.bit_scan0(k)

Scan C<n>, starting from bit index C<k>, towards more significant bits, until a 0-bit is found.

When C<k> is omitted, C<k=0> is assumed.

Returns C<nil> if C<n> cannot be truncated to an integer or if C<k> is negative.

=cut

=head2 bit_scan1

    n.bit_scan1
    n.bit_scan1(k)

Scan C<n>, starting from bit index C<k>, towards more significant bits, until a 1-bit is found.

When C<k> is omitted, C<k=0> is assumed.

Returns C<nil> if C<n> cannot be truncated to an integer or if C<k> is negative.

=cut

=head2 bsearch

    bsearch(n, {...})
    bsearch(a, b, {...})

Binary search from to C<0> to C<n>, or from C<a> to C<b>, which can be any arbitrary large integers.

The last argument is a block which does the comparisons.

This function finds a value C<k> such that f(k) = 0. Returns C<nil> otherwise.

    say bsearch(20,      {|k| k*k  <=> 49   })   #=> 7   (7*7  = 49)
    say bsearch(3, 1000, {|k| k**k <=> 3125 })   #=> 5   (5**5 = 3125)

=cut

=head2 bsearch_ge

    bsearch_ge(n, {...})
    bsearch_ge(a, b, {...})

Binary search from to C<0> to C<n>, or from C<a> to C<b>, which can be any arbitrary large integers.

The last argument is a block which does the comparisons.

This function finds a value C<k> such that f(k-1) < 0 and f(k) >= 0. Returns C<nil> otherwise.

    bsearch_ge(1e6,       { .exp <=> 1e+9 })  #  21   (exp( 21) >= 1e+9)
    bsearch_ge(-1e6, 1e6, { .exp <=> 1e-9 })  # -20   (exp(-20) >= 1e-9)

=cut

=head2 bsearch_le

    bsearch_le(n, {...})
    bsearch_le(a, b, {...})

Binary search from to C<0> to C<n>, or from C<a> to C<b>, which can be any arbitrary large integers.

The last argument is a block which does the comparisons.

This function finds a value C<k> such that f(k) <= 0 and f(k+1) > 0. Returns C<nil> otherwise.

    bsearch_le(1e6,       { .exp <=> 1e+9 })  #  20   (exp( 20) <= 1e+9)
    bsearch_le(-1e6, 1e6, { .exp <=> 1e-9 })  # -21   (exp(-21) <= 1e-9)

=cut

=head2 bsearch_max

    bsearch_max(n, {...})
    bsearch_max(a,b, {...})

Binary search, returning the largest integer value in the range C<a..b> that satisifes the given comparison function.

Example:

    say bsearch_max(1, 1e6, {|k| pi(k) <=> 100 })   #=> 546

where:

    n = 546 is the  largest value satisfying pi(n) <= 100

=cut

=head2 bsearch_min

    bsearch_min(n, {...})
    bsearch_min(a,b, {...})

Binary search, returning the smallest integer value in the range C<a..b> that satisifes the given comparison function.

Example:

    say bsearch_min(1, 1e6, {|k| pi(k) <=> 100 })   #=> 541

where:

    n = 541 is the smallest value satisfying pi(n) >= 100

=cut

=head2 bsearch_solve

    bsearch_solve(n, {...})
    bsearch_solve(a,b, {...})

It computes the inverse of any continuous function, given the range that includes the inverse value.

For floating-point values, the C<approx_cmp(a,b)> method (or the `<~>` operator) is recommended to be used for comparisons.

Example:

    say bsearch_inverse(100, {|x| exp(x) <~> 2 })          # solution to x for: exp(x) =   2
    say bsearch_inverse(200, {|x| x**2  <=> 43 })          # solution to x for:    x^2 =  43
    say bsearch_inverse(-10, 10, {|x| x**3 <~> -43 })      # solution to x for:    x^3 = -43
    say bsearch_inverse(300, 500, {|x| Li(x) <~> 100 })    # solution to x for:  Li(x) = 100

This method can also be used in computing approximations to some integer-specific functions:

    var n = 100000
    var v = 2*int(n*log(n) / log(log(n)))

    say nth_semiprime(n)                                        #=> 459577
    say bsearch_inverse(v, {|x| semiprime_count(x) <=> n })     #=> 459577.93302154541015625

Aliases: I<bsearch_inverse>

=cut

=head2 by

    n.by { ... }

Returns an array with C<n> elements >= 0 that satisfy the provided block of code.

Example:

    say 10.by { .is_prime }      # first 10 primes
    say 10.by { .is_square }     # first 10 squares

Aliases: I<first>

=cut

=head2 cadd

    cadd(a,b,x,y)

Complex arithmetic addition, defined as:

    cadd(a,b,x,y)   #=> (a+x, b+y)

Aliases: I<complex_add>

=cut

=head2 catalan

    n.catalan
    n.catalan(k)

Returns the n-th Catalan number.

If two arguments are provided, it returns the C<C(n,k)> entry in Catalan's triangle.

=cut

=head2 cbrt

    n.cbrt

Cube root of C<n>, as a floating-point value.

=cut

=head2 cdiv

    cdiv(a,b,x,y)

Complex arithmetic division, defined as:

    cdiv(a,b,x,y)   #=> ((a*x + b*y)/(x*x + y*y), (b*x - a*y)/(x*x + y*y))

Aliases: I<complex_div>

=cut

=head2 ceil

    n.ceil

Round C<n> towards positive Infinity.

Aliases: I<ceiling>

=cut

=head2 cfrac

    n.cfrac
    n.cfrac(k)

Compute C<k> terms of the simple continued-fraction expansion of C<n>.

Example:

    say sqrt(12).cfrac(6)    # [3, 2, 6, 2, 6, 2, 6]

Can also be used to compute very good rational approximations to a given real number:

    say Num.pi.cfrac(10).flip.reduce{|a,b| b + 1/a }.as_rat      # 4272943/1360120

When C<k> is omitted, it uses the default floating-point precision to deduce C<k>.

Aliases: I<as_cfrac>

=cut

=head2 chebyshevT

    chebyshevT(n)
    chebyshevT(n, x)

Compute the Chebyshev polynomials of the first kind: C<T_n(x)>, where C<n> must be a native integer.

Defined as:

    T(0, x) = 1
    T(1, x) = x
    T(n, x) = 2*x*T(n-1, x) - T(n-2, x)

When C<x> is ommited, a Polynomial object is returned.

Aliases: I<chebyshevt>

=cut

=head2 chebyshevTmod

    chebyshevTmod(n, x, m)

Compute the modular Chebyshev polynomials of the first kind: C<T_n(x) mod m>, where C<n> and C<m> must be integers (arbitrarily large).

=cut

=head2 chebyshevU

    chebyshevU(n)
    chebyshevU(n, x)

Compute the Chebyshev polynomials of the second kind: C<U_n(x)>, where C<n> must be a native integer.

Defined as:

    U(0, x) = 1
    U(1, x) = 2*x
    U(n, x) = 2*x*U(n-1, x) - U(n-2, x)

When C<x> is ommited, a Polynomial object is returned.

Aliases: I<chebyshevu>

=cut

=head2 chebyshevUmod

    chebyshevUmod(n, x, m)

Compute the modular Chebyshev polynomials of the second kind: C<U_n(x) mod m>, where C<n> and C<m> must be integers (arbitrarily large).

=cut

=head2 chr

    n.chr

Convert the integer C<n> into a character.

Example:

    say 97.chr      # "a"
    say 9786.chr    # "☺"

=cut

=head2 cinv

    cinv(a,b)

Complex arithmetic inversion, defined as:

    cinv(a,b)   #=> (a/(a*a + b*b), (-b)/(a*a + b*b))

Aliases: I<complex_inv>

=cut

=head2 cinvmod

    cinvmod(a,b,m)

Complex modular inversion modulo C<m>: returns a pair of integers (x,y) such that:

    cmod(cmul(a,b,x,y), m) == (1, 0)

Aliases: I<complex_invmod>

=cut

=head2 circular_permutations

    n.circular_permutations
    n.circular_permutations { ... }

Returns an array of arrays with the circular permutations of the integers in the range C<0..n-1>, or iterates over the circular permutations when a block is given.

Example:

    5.circular_permutations {|*a| say a }

=cut

=head2 cis

    cis(x)

Euler's formula applied on C<x>, defined as:

    cis(x) = cos(x) + sin(x)*i

=cut

=head2 clearbit

    n.clearbit(k)

Set the k-th bit of integer n to 0.

Example:

    say clearbit(0b1001, 0).as_bin  #=> 1000
    say clearbit(0b1100, 2).as_bin  #=> 1000

=cut

=head2 cmod

    cmod(a,b,m)

Complex arithmetic modular operation, defined as:

    cmod(a,b,m)     #=> (a%m, b%m)

Aliases: I<complex_mod>

=cut

=head2 cmul

    cmul(a,b,x,y)

Complex arithmetic multiplication, defined as:

    cmul(a,b,x,y)   #=> (a*x - b*y, a*y + b*x)

Aliases: I<complex_mul>

=cut

=head2 combinations

    n.combinations(k)
    n.combinations(k, { ... })

Returns an array with the k-combinations of the integers in the range C<0..n-1>, or iterates over the k-combinations when a block is given.

Example:

    5.combinations(2, {|*a| say a })

=cut

=head2 combinations_with_repetition

    n.combinations_with_repetition(k)
    n.combinations_with_repetition(k, { ... })

Returns an array with the k-combinations with repetition of the integers in the range C<0..n-1>, or iterates over the k-combinations with repetition when a block is given.

Example:

    5.combinations_with_repetition(2, {|*a| say a })

=cut

=head2 commify

    n.commify

Returns a string with thousands separators added after each group of 3 decimal digits.

Exmaple:

    say 1000.commify       #=> 1,000
    say 1e10.commify       #=> 10,000,000,000

=cut

=head2 complex

    n.complex
    complex(a,b)

Converts C<n> to a complex number, or creates a complex number from C<a> and C<b>.

Example:

    say complex(3, 4)       #=> 3+4i

This is equivalent with:

    say Complex(3, 4)       #=> 3+4i

=cut

=head2 complex_cmp

    complex_cmp(a,b,x,y)

Complex number comparison, defined as:

    (a <=> x) || (b <=> y)

=cut

=head2 complex_ipow

    complex_ipow(a,b,n)

Complex integer exponentiation: returns C<(x,y)> such that C<x+y*i = (a+b*i)^n>.

Example:

    say [complex_ipow(3,4,5)]       #=> [-237, -3116]

=cut

=head2 composite

    n.composite

Returns the n-th composite number (OEIS: A002808).

Example:

    say composite(10**9)        #=> 1053422339

Aliases: I<nth_composite>

=cut

=head2 composite_count

    composite_count(n)
    composite_count(a,b)

Returns the number of composite numbers <= n, or in the range C<a..b>.

Example:

    say composite_count(100)          # number of composites <= 100
    say composite_count(50, 100)      # number of composites in the range 50..100

=cut

=head2 composites

    composites(n)
    composites(a,b)

Returns an array with the composite numbers <= n, or in the range C<a..b>.

=cut

=head2 conj

    conj(x)

Complex conjugate of C<x>. For real integers, this is a fixed-point function.

=cut

=head2 consecutive_lcm

    consecutive_lcm(n)

Returns the least common multiple (LCM) of all the integers in the range C<1..n>.

Aliases: I<consecutive_integer_lcm>

=cut

=head2 convergents

    n.convergents(k)

Returns an array with the continued fraction convergents for a given real number C<n>, where C<k> is the number of convergenents to be computed and returned.

Example:

    say Num.pi.convergents(5)   #=> [3, 22/7, 333/106, 355/113, 103993/33102]

=cut

=head2 core

    core(n)

Squarefree part of n.

Example:

    say 30.of { .core }     #=> OEIS: A007913

Equivalent to PARI/GP C<core(n)> function.

=cut

=head2 cos

    cos(x)

Trigonometric cosine function.

=cut

=head2 cosh

    cosh(x)

Hyperbolic cosine function.

=cut

=head2 cot

    cot(x)

Trigonometric cotangent function.

=cut

=head2 coth

    coth(x)

Hyperbolic cotangent function.

=cut

=head2 cpow

    cpow(a,b,n)

Computes C<(a + b*i)^n>, where a,b are real numbers and n is an integer. Returns the real and imaginary part as a list.

Example:

    say [cpow(3, 4, 10)]    #=> [-9653287, 1476984]

Aliases: I<complex_pow>

=cut

=head2 cpowmod

    cpowmod(a,b,n,m)

Efficiently computes C<(a + b*i)^n mod m>, where a,b,n,m are all integers. Returns the real and imaginary part as a list.

Example:

    say [complex_powmod(3, 4, 1000, 1e6)]   #=> [585313, 426784]

Aliases: I<complex_powmod>

=cut

=head2 csc

    csc(x)

Trigonometric cosecant function.

=cut

=head2 csch

    csch(x)

Hyperbolic cosecant function.

=cut

=head2 csub

    csub(a,b,x,y)

Complex arithmetic subtration, defined as:

    csub(a,b,x,y)   #=> (a-x, b-y)

Aliases: I<complex_sub>

=cut

=head2 cube

    cube(x)

Returns the cube of C<x>. Equivalent with C<x**3>.

=cut

=head2 cyclotomic

    cyclotomic(n)
    cyclotomic(n,x)

Returns the n-th Cyclotomic polynomial evaluated at C<x>.

    say cyclotomic(12, 10)      #=> 9901

When C<x> is ommited, a Polynomial object is returned:

    say cyclotomic(12)          #=> x^4 - x^2 + 1

Aliases: I<cyclotomic_polynomial>

=cut

=head2 dconv

    n.dconv(f,g)

Returns the Dirichlet convolution of f and g.

Example:

    say 20.of { .dirichlet_convolution({.moebius}, {_}) }

Aliases: I<dirichlet_convolution>

=cut

=head2 de

    r.de
    r.denominator

Returns the denominator for a rational number C<r>.

Example:

    say denominator(43/97)      #=> 97

Aliases: I<denominator>

=cut

=head2 defs

    n.defs { ... }

Returns an array with the first n defined values returned by the given block. The block is called with k = 0,1,...

Example:

    10.defs {|k| k.is_prime ? k+1 : nil }      # array of p+1 for the first 10 primes p

=cut

=head2 deg2rad

    deg2rad(x)

Convert degrees to radians.

Example:

    say deg2rad(180)    #=> 3.14159...

=cut

=head2 derangements

    n.derangements

Returns an array of arrays with the derangements of the integers in the range C<0..n-1>, or iterates over the derangements when a block is given.

Example:

    5.derangements {|*a| say a }

Aliases: I<complete_permutations>

=cut

=head2 derivative

    derivative(n)

Arithmetic derivative of C<n>, defined for rationals and integers (positive and negative).

Example:

    say derivative(5040)                 #=> 15168
    say derivative(-5040/4323).as_frac   #=> -6240176/2076481

Aliases: I<arithmetic_derivative>

=cut

=head2 digit

    n.digit(k, b=10)

Returns the k-th digit of C<n> in a base C<b>.

Exmaples:

    say 1119.digit(0)         #=> 9
    say 1181.digit(1)         #=> 8
    say 1711.digit(2)         #=> 7
    say 6111.digit(3)         #=> 6
    say 1234.digit(4)         #=> 0

It also supports negative indices:

    say 9111.digit(-1)        #=> 9
    say 1811.digit(-2)        #=> 8
    say 1234.digit(-42)       #=> nil

=cut

=head2 digits

    n.digits(b=10)

Returns the digits of C<n> in base C<b>, ordered from the least significant digit to the most significant digit.

Example:

    say 1234.digits      #=> [4,3,2,1]
    say 1234.digits(20)  #=> [14, 1, 3]

The reverse operation is:

    b.digits2num(n.digits(b)) == n

=cut

=head2 digits2num

    b.digits2num(digits)

Convert an array of digits to a number in the base C<b>.

The array of digits are ordered from the least significant digit to the most significant digit, as returned by C<n.digits(b)>.

Example:

    say 10.digits2num([4,3,2,1])  #=> 1234

Aliases: I<from_digits>

=cut

=head2 digits_sum

    sumdigits(n, b=10)

Sum of base C<b> digits of C<n>.

Example:

    say sumdigits(1234)         #=> 10
    say sumdigits(1e5!, 100)    #=> 10658934

This is equivalent to:

    n.digits(b).sum

Aliases: I<sum_digits>, I<sumdigits>

=cut

=head2 dirichlet_sum

    n.dirichlet_sum(f, g, F, G)

The method computes the following sum in O(sqrt(n)) steps:

    Sum_{k=1..n} Sum_{d|k} f(d) * g(k/d)

Where C<F> and C<G> are the partial sums of C<f> and C<g>, respectively.

Example:

    say dirichlet_hyperbola(1e3,
            { .is_squarefree ? 1 : 0 },
            { _*_ },
            { .squarefree_count },
            { .faulhaber(2) }
    )

Aliases: I<dirichlet_hyperbola>

=cut

=head2 divides

    a.divides(b)
    a `divides` b

Returns true if C<a> divides C<b>.

Example:

    say 3.divides(15)   #=> true

=cut

=head2 divisor_map

    n.divisor_map {|d| ... }

Maps the divisors of n to the given block.

Example:

    say 24.divisor_map {|d| 1/d }  #=> [1, 1/2, 1/3, 1/4, 1/6, 1/8, 1/12, 1/24]

=cut

=head2 divisor_prod

    n.divisor_prod {|d| ... }

Product of the mapping of the positive divisors of n.

Equivalent with:

    n.divisor_map {|d| ... }.prod

Aliases: I<divisors_prod>

=cut

=head2 divisors

    n.divisors(k=n)

Return an array with the positive divisors of n <= k.

Example:

    say 120.divisors        #=> [1, 2, 3, 4, 5, 6, 8, 10, 12, 15, 20, 24, 30, 40, 60, 120]
    say 120.divisors(13)    #=> [1, 2, 3, 4, 5, 6, 8, 10, 12]

=cut

=head2 divisor_sum

    n.divisor_sum {|d| ... }

Sum of the mapping of the positive divisors of n.

Example:

    say 5040.divisor_sum {|d| euler_phi(d)**2 } #=> 2217854

Equivalent with:

    n.divisor_map {|d| ... }.sum

Aliases: I<divisors_sum>

=cut

=head2 divmod

    divmod(a, b)
    divmod(a, b, m)

When only two arguments are provided, it returns C<(a//b, a%b)>.

Example:

    say [divmod(23, 10)]    #=> [2,3]

When three arguments are given, it does integer modular division: C<(a/b) % m>.

Example:

    say divmod(43, 97, 127)     # == (43 * invmod(97, 127))%127

=cut

=head2 downto

    a.downto(b, step=1)

Returns a reverse range from C<a> down to C<b>, with an optional stepping value.

Example:

    say 10.downto(1)        #=> RangeNum(10, 1, -1)
    say 10.downto(1, 2)     #=> RangeNum(10, 1, -2)

=cut

=head2 dump

    n.dump

Returns a stringification version of C<n>.

Example:

    say dump(42)    #=> "42"
    say dump(3/4)   #=> "3/4"

=cut

=head2 e

    Num.e

Returns the I<e> mathematical constant: 2.71828...

=cut

=head2 each_almost_prime

    k.each_almost_prime(n, {...})
    k.each_almost_prime(a,b, {...})

Iterates over the k-almost prime numbers <= n, or in the range C<a..b>.

Example:

    11.almost_primes_each(1e7, {|n| say n })        # iterate over 11-almost primes <= 1e6
    11.almost_primes_each(1e6, 1e7, {|n| say n })   # iterate over 11-almost primes in the range [1e6, 1e7]

Aliases: I<almost_primes_each>

=cut

=head2 each_composite

    n.each_composite { ... }
    each_composite(a,b, { ... })

Iterate over the composite numbers <= n, or in the given range C<a..b>.

Example:

    # Iterate over the composite integers between 100 and 200
    each_composite(100, 200, {|c|
        say c
    })

    # Iterate over the composite integers <= 100
    100.each_composite {|c|
        say c
    }

Aliases: I<composites_each>

=cut

=head2 each_omega_prime

    k.omega_primes_each(n, { ... })
    k.omega_primes_each(a, b, { ... })

Iterates over the k-omega primes <= n, or in the range C<a..b>.

k-omega primes are numbers C<n> such that C<omega(n) == k>.

Example:

    1.omega_primes_each(100, {|n| say n })          # iterate over prime powers <= 100
    2.omega_primes_each(50, 100, {|n| say n })      # iterate over 2-omega primes in range 50..100

Aliases: I<omega_primes_each>

=cut

=head2 each_powerful

    k.powerful_each(n, { ... })
    k.powerful_each(a, b, { ... })

Iterates over the k-powerful numbers <= n, or in the range C<a..b>.

Example:

    2.powerful_each(100, {|n| say n })          # iterate over 2-powerful numbers <= 100
    2.powerful_each(50, 100, {|n| say n })      # iterate over 2-powerful in the range 50..100

Aliases: I<powerful_each>

=cut

=head2 each_prime

    n.each_prime {...}
    each_prime(a,b,{...})

It efficiently iterates over the primes in the given range, using a segmented prime sieve.

Example:

    # Iterate over the primes between 100 and 200
    each_prime(100, 200, {|p|
        say p
    })

    # Iterate over the primes <= 100
    100.each_prime {|p|
        say p
    }

It's consirably faster than using the C<next_prime(n)> method.

Aliases: I<primes_each>

=cut

=head2 each_semiprime

    n.each_semiprime { ... }
    each_semiprime(a,b, { ... })

Iterate over semiprimes <= n, or in the range C<a..b>:

    100.each_semiprime {|k| say k }         # iterate over semiprimes <= 100
    each_semiprime(50, 100, {|k| say k })   # iterate over semiprimes in the range [50, 100]

Aliases: I<semiprimes_each>

=cut

=head2 each_squarefree

    n.each_squarefree {...}
    each_squarefree(a,b,{...})

Iterates over the squarefree numbers in a given range.

Example:

    # Iterate over the squarefree numbers in the interval [100, 200]
    each_squarefree(100, 200, {|n|
        say n
    })

    # Iterate over the squarefree numbers <= 100
    100.each_squarefree {|n|
        say n
    }

Aliases: I<squarefree_each>

=cut

=head2 each_squarefree_almost_prime

    k.each_squarefree_almost_prime(n, {...})
    k.each_squarefree_almost_prime(a,b,{...})

Iterates over the squarefree k-almost primes <= n, or in the range C<a..b>.

Example:

    # Iterate over squarefree 3-almost primes <= 100
    3.squarefree_almost_primes_each(100, { .say })

    # Iterate over squarefree 3-almost primes in the range 50..100
    3.squarefree_almost_primes_each(50, 100, { . say })

Aliases: I<squarefree_almost_primes_each>

=cut

=head2 ecm_factor

    n.ecm_factor()
    n.ecm_factor(B)
    n.ecm_factor(B,reps)

Hendrik Lenstra's elliptic curve factorization method (ECM).

=cut

=head2 ei

    ei(x)

Exponential integral function.

Aliases: I<eint>

=cut

=head2 erf

    erf(x)

The Gauss error function.

=cut

=head2 erfc

    erfc(x)

The complementary error function.

=cut

=head2 euler

    euler(n)
    euler(n,x)

Returns the n-th Euler number:

    say 10.of {|n| euler_number(n) }   #=> [1, 0, -1, 0, 5, 0, -61, 0, 1385, 0]

Returns the n-th Euler polynomial evaluated at C<x>, when C<x> is given:

    say euler(10, 5)    #=> 1981100

Aliases: I<euler_number>

=cut

=head2 euler_polynomial

    euler_polynomial(n)
    euler_polynomial(n, x)

Returns the n-th Euler polynomial evaluated at C<x>.

When C<x> is ommited, a Polynomial object is returned.

=cut

=head2 eval

    eval(x,v)

Returns back C<x>.

=cut

=head2 exp

    exp(x)

Exponential function: C<e^x>.

=cut

=head2 exp10

    exp10(x)

Exponential function: C<10^x>.

=cut

=head2 exp2

    exp2(x)

Exponential function: C<2^x>.

=cut

=head2 exp_mangoldt

    exp_mangoldt(n)

Returns exp(Λ(n)); the exponential of the Mangoldt function.

=cut

=head2 expmod

    powmod(b, n, m)

Modular exponentiation: C<b^n mod m>, where C<b> is an integer or a rational number and C<n> and C<m> are both integers.

Example:

    say powmod(2, 42, 43)           #=> 1
    say powmod(3/4, 1234, 4171)     #=> 2138

Aliases: I<powmod>

=cut

=head2 expnorm

    expnorm(n, b=10)

Returns C<exp(n)> normalized in the range C<[0,1)>.

Example:

    say expnorm(log(2) * 20996011)          #=> 0.125976895450330105020494309574...

    say exp(log(10!))               #=> 3628800
    say expnorm(log(10!))           #=> 0.36288
    say expnorm(log(10!), 2)        #=> 0.86517333984375

    say 3628800.base(2)             #=> 1101110101111100000000
    say 0.86517333984375.base(2)    #=> 11011101011111/100000000000000

Complex numbers are also supported:

    say expnorm(log(-1))            #=> 0.06682015101903....+0.074398033695749....i
    say expnorm(log(-1234)).abs     #=> 0.1234

=cut

=head2 factor

    factor(n)

Returns an array with the prime factors of C<n>, where C<n> is a positive integer, such that C<n.factor.prod == n>.

Example:

    say 180.factor  #=> [2, 2, 3, 3, 5]

Aliases: I<factors>

=cut

=head2 factor_exp

    factor_exp(n)

Returns an array of pairs C<[p,k]> factors C<p^k> of n.

Example:

    say 180.factor_exp  #=> [[2, 2], [3, 2], [5, 1]]

Aliases: I<factors_exp>

=cut

=head2 factorialmod

    factorialmod(n,m)

Factorial of n modulo m (faster than C<n! % m>).

=cut

=head2 factorial_power

    factorial_power(n, p)

Returns the number of times C<p> divides C<n!>, where C<p> is a prime number.

Equivalent with (but more efficient):

    valuation(n!, p)

Aliases: I<factorial_valuation>

=cut

=head2 factorial_sum

    factorial_sum(n)

Left factorial of C<n>, defined as:

    Sum_{k=0..n-1} k!

Example:

    say 20.of { .factorial_sum }    # OEIS: A003422

Aliases: I<left_factorial>

=cut

=head2 factor_map

    n.factor_map {|p,k| ... }

Maps the prime-power factorization C<(p,k)> of n to the given block.

Example:

    say 5040.factor_map  {|p,k| p**k }  #=> [16, 9, 5, 7]

=cut

=head2 factor_prod

    n.factor_prod {|p,k| ... }

Product of the mapping of the prime-power factorization of n.

Example:

    say 5040.factor_prod {|p,k| (p-1) * p**(k-1) }  #=> 1152

Equivalent with:

    n.factor_map {|p,k| ... }.prod

Aliases: I<factors_prod>

=cut

=head2 factor_sum

    n.factor_sum {|p,k| ... }

Sum of the mapping of the prime-power factorization of n.

Equivalent with:

    n.factor_map {|p,k| ... }.sum

Aliases: I<factors_sum>

=cut

=head2 falling_factorial

    falling_factorial(n,k)

Falling factorial: C<(n)_k = n * (n - 1) * ... * (n - k + 1)>, defined as:

    binomial(n, k) * k!

For negative values of C<k>, falling factorial is defined as:

    falling_factorial(n, -k) = 1/falling_factorial(n + k, k)

When the denominator is zero, NaN is returned.

=cut

=head2 faulhaber

    faulhaber(n,k)

Sum of powers: C<1^k + 2^k + ... + n^k>, using Faulhaber's summation formula.

The value for C<k> must be a non-negative native integer.

Example:

    faulhaber(5, 2)   # 1^2 + 2^2 + 3^2 + 4^2 + 5^2 = 55

Aliases: I<faulhaber_sum>

=cut

=head2 faulhaber_polynomial

    faulhaber_polynomial(n)
    faulhaber_polynomial(n, x)

Computes the n-th Faulhaber polynomials evaluated at C<x>.

Defined in terms of the Bernoulli polynomials, as:

    faulhaber_polynomial(n,x) = (bernoulli_polynomial(n+1,x+1) - bernoulli_polynomial(n+1, 1))/(n+1)

When C<x> is ommited, a Polynomial object is returned.

=cut

=head2 faulhaber_range

    faulhaber_range(a, b, k)

Sum of powers: C<a^k + (a+1)^k + (a+2)^k + ... + b^k>, using Faulhaber's summation formula.

The value for C<k> must be a non-negative native integer.

Example:

    faulhaber_range(50, 100, 2)   # 50^2 + 51^2 + ... + 100^2 = 297925

=cut

=head2 fermat_factor

    n.fermat_factor(k=1e4)

Tries to factorize a given number using Fermat's factorization method (using at most C<k> iterations).

Works for odd composite non-power numbers C<n> that have two divisors close to C<sqrt(n)>.

=cut

=head2 fib

    fib(n)
    fib(n,k)

Returns the n-th Fibonacci number.

    say 10.of { .fib }      #=> [0, 1, 1, 2, 3, 5, 8, 13, 21, 34]

When an additional is provided, it returns the n-th Fibonacci number of k-th order.

    say 20.of { .fib(3) }   #=> Tribonacci numbers
    say 20.of { .fib(4) }   #=> Tetranacci numbers
    say 20.of { .fib(5) }   #=> Pentanacci numbers

Aliases: I<fibonacci>

=cut

=head2 fibmod

    fibmod(n,m)

Efficiently compute the n-th Fibonacci number modulo m.

Aliases: I<fibonacci_mod>, I<fibonaccimod>

=cut

=head2 flip

    n.flip(base=10)

Returns the reversal of C<n> in base C<b>. When C<b> is not given, it defaults to base 10.

Example:

    say 20.of { .flip }         # A004086
    say 20.of { .flip(2) }      # A030101

Aliases: I<reverse>

=cut

=head2 flipbit

    flipbit(n,k)

Flip the value of the k-th bit of integer C<n>.

Example:

    say flipbit(0b1000, 0).as_bin  #=> 1001
    say flipbit(0b1001, 0).as_bin  #=> 1000

=cut

=head2 floor

    floor(x)

Round C<x> towards -Infinity.

Example:

    say floor( 2.5)     #  2
    say floor(-2.5)     # -3

=cut

=head2 flt_factor

    n.flt_factor(base=2, tries=1e4)

Tries to find a factor of C<n>, using a new factorization method, inspired by Fermat's Little Theorem (FLT).

The method is particularily effective for numbers that have factors closed to each other, or have a factor C<k> for which C<znorder(2,k)> is small.

Example (try with base 3 and give up after 10^6 iterations):

    say flt_factor(2**64 + 1, 3, 1e6)   #=> [274177, 67280421310721]

=cut

=head2 fusc

    fusc(n)

Returns the n-th term in Stern's diatomic series (or Stern-Brocot sequence).

Example:

    say 30.of { .fusc }     #=> OEIS: A002487

=cut

=head2 gcd

    gcd(...)

The greatest common divisor of a list of integers.

=cut

=head2 gcdext

    gcdext(a,b)

The extended greatest common divisor of C<a> and C<b>, returning C<(u, v, d)>, where C<d = gcd(a,b)>, while C<u> and C<v> are the coefficients satisfying:

    u*a + v*b = d.

The value of C<d> is always non-negative.

=cut

=head2 geometric_sum

    geometric_sum(n,r)

Geometric sum: C<r^0 + r^1 + ... + r^n>, using the following formula:

    geometric_sum(n, r) = (r^(n+1) - 1) / (r - 1)

Example:

    say geometric_sum(5, 8)     # 8^0 + 8^1 + 8^2 + 8^3 + 8^4 + 8^5 = 37449

=cut

=head2 gpf

    gpf(n)

Returns the greatest prime factor of n.

Defined with the base-cases:

    gpf(0) = 0
    gpf(1) = 1

Example:

    say gpf(2**128 + 1)     #=> 5704689200685129054721

=cut

=head2 hamdist

    hamdist(a,b)

Returns the Hamming distance (number of bit-positions where the bits differ) between integers C<a> and C<b>.

=cut

=head2 harm

    harmonic(n)
    harmonic(n, k)

Returns the n-th Harmonic number C<H_n>. The harmonic numbers are the sum of reciprocals of the first C<n> natural numbers: C<1 + 1/2 + 1/3 + ... + 1/n>.

When an additional argument is given, it returns the n-th Harmonic number of the k-th order.

Aliases: I<harmfrac>, I<harmonic>, I<harmonic_number>

=cut

=head2 harmreal

    harmreal(n)
    harmreal(n, k)

Returns the n-th Harmonic number C<H_n> as a floating-point value, defined as:

    harmreal(n) = digamma(n+1) + γ

where C<γ> is the Euler-Mascheroni constant.

When an additional argument is given, it returns the n-th Harmonic number of the k-th order.

=cut

=head2 hclassno

    hclassno(n)

Returns the Hurwitz-Kronecker class number.

Example:

    say 30.of { .hclassno.nu }      # OEIS: A058305
    say 30.of { .hclassno.de }      # OEIS: A058306
    say 30.of { 12 * .hclassno }    # OEIS: A259825

=cut

=head2 hermiteH

    hermiteH(n)
    hermiteH(n, x)

Physicists' Hermite polynomials: C<H_n(x)>.

When C<x> is ommited, a Polynomial object is returned.

Aliases: I<hermite_polynomialH>

=cut

=head2 hermiteHe

    hermiteHe(n)
    hermiteHe(n, x)

Probabilists' Hermite polynomials: C<He_n(x)>.

When C<x> is ommited, a Polynomial object is returned.

Aliases: I<hermite_polynomialHe>

=cut

=head2 holf_factor

    n.holf_factor(reps)

Hart's OLF method (variant of Fermat's method).

=cut

=head2 hyperfactorial

    Number.hyperfactorial()

Hyperfactorial of C<n>, defined as C<Prod_{k=1..n} k^k>.

=cut

=head2 hyperfactorial_ln

    hyperfactorial_ln(n)

Natural logarithm of C<hyperfactorial(n)>, where C<n> is a non-negative integer.

Aliases: I<lnhyperfactorial>, I<hyperfactorial_log>

=cut

=head2 hypot

    hypot(x,y)

The value of the hypotenuse for catheti C<x> and C<y>, defined as:

    sqrt(x**2 + y**2)

Also defined for complex numbers.

=cut

=head2 i

    x.i

Multiplies C<x> by the imaginary unit C<i>.

Example:

    say 42.i         #=> 42i
    say 42i.i        #=> -42

=cut

=head2 iadd

    iadd(a,b)

Integer addition: C<a+b>.

=cut

=head2 icbrt

    icbrt(n)

Integer cube root of n.

=cut

=head2 idiv_ceil

    idiv_ceil(a,b)

Integer division of integers C<a> and C<b>, rounded towards +Infinity.

When C<a> and C<b> are integers, this is equivalent with:

    ceil(a/b)

=cut

=head2 idiv_round

    idiv_round(a,b)

Integer division of integers C<a> and C<b>, rounded towards nearest integer.

When C<a> and C<b> are integers, this is equivalent with:

    floor(a/b + 1/2)

=cut

=head2 idiv_trunc

    idiv_trunc(a,b)

Integer division of integers C<a> and C<b>, rounded towards 0.

When C<a> and C<b> are integers, this is equivalent with:

    trunc(a/b)

=cut

=head2 ilog

    ilog(n,b)

Integer logarithm of C<n> to base C<b>, satisfying:

    b**ilog(n,b) <= n < b**(ilog(n,b)+1)

=cut

=head2 ilog10

    ilog10(n)

Integer logarithm of C<n> to base C<10>, equivalent to:

    ilog(n,10)

=cut

=head2 ilog2

    ilog2(n)

Integer logarithm of C<n> to base C<2>, equivalent to:

    ilog(n,2)

=cut

=head2 im

    x.im

The imaginary part of complex number C<x>. Return 0 when C<x> is a real number.

Aliases: I<imag>, I<imaginary>

=cut

=head2 imod

    imod(a,m)

Integer remainder of C<a> when divided by C<m>: C<a % m>.

=cut

=head2 imul

    imul(a,b)

Integer multiplication: C<a*b>.

=cut

=head2 inf

    Num.inf

Returns the positive Infinity special floating-point value.

=cut

=head2 int

    int(x)
    trunc(x)

Truncate C<x> to an integer.

Aliases: I<to_i>, I<to_int>, I<trunc>

=cut

=head2 inv

    inv(x)

Multiplicative inverse of C<x>: C<1/x>.

=cut

=head2 inverse_phi

    n.inverse_phi

Returns all the solutions C<x> to the Euler totient function: phi(x) = n.

Aliases: I<inverse_totient>

=cut

=head2 inverse_phi_len

    n.inverse_phi_len

Returns the number of solutions to the Euler totient function: phi(x) = n.

Equivalent with:

    n.inverse_phi.len

Aliases: I<inverse_totient_len>

=cut

=head2 inverse_phi_max

    n.inverse_phi_max

Returns the largest solution x to the Euler totient function: phi(x) = n.

Equivalent with:

    n.inverse_phi.max

Returns C<nil> if there are no solutions.

Aliases: I<inverse_euler_phi_max>

=cut

=head2 inverse_phi_min

    n.inverse_phi_min

Returns the smallest solution x to the Euler totient function: phi(x) = n.

Equivalent with:

    n.inverse_phi.min

Returns C<nil> if there are no solutions.

Aliases: I<inverse_euler_phi_min>

=cut

=head2 inverse_polygonal

    polygonal_inverse(n)

Returns an array of pairs C<[r,k]> such that C<polygonal(r,k) = n>.

Example:

    say polygonal_inverse(4012)   #=> [[2, 4012], [4, 670], [8, 145], [4012, 2]]

Aliases: I<polygonal_inverse>

=cut

=head2 inverse_psi

    n.inverse_psi

Returns all the solutions C<x> to the Dedekin psi function: psi(x) = n.

Example:

    say inverse_psi(120)    #=> [75, 76, 87, 95]

Aliases: I<inverse_dedekind_psi>

=cut

=head2 inverse_psi_len

    n.inverse_psi_len

Returns the number of solutions to Dedekind's psi function: psi(x) = n.

Equivalent to C<n.inverse_psi.len>, but much faster.

Aliases: I<inverse_dedekind_psi_len>

=cut

=head2 inverse_psi_max

    n.inverse_psi_max

Returns the largest solution x to the Dedekind psi function: psi(x) = n.

Equivalent with:

    n.inverse_psi.max

Returns C<nil> if there are no solutions.

Aliases: I<inverse_dedekind_psi_max>

=cut

=head2 inverse_psi_min

    n.inverse_psi_min

Returns the smallest solution x to the Dedekind psi function: psi(x) = n.

Equivalent with:

    n.inverse_psi.min

Returns C<nil> if there are no solutions.

Aliases: I<inverse_dedekind_psi_min>

=cut

=head2 inverse_sigma

    n.inverse_sigma(k=1)

The method returns all the numbers C<m> for which C<sigma(m,k) = n>, where C<n> is given.

Example:

    say inverse_sigma(42)           # [20, 26, 41]
    say inverse_sigma(22100, 2)     # [120, 130, 141]

Also works with arbitrary large integers:

    say inverse_sigma(9325257382230393314439814176)

=cut

=head2 inverse_sigma_len

    n.inverse_sigma_len(k=1)

Returns the number of solutions to the sigma sum of divisors function: sigma_k(x) = n.

Equivalent to C<n.inverse_sigma(k).len>, but much faster.

=cut

=head2 inverse_sigma_max

    n.inverse_sigma_max

Returns the largest solution x to the sum of divisors function: sigma(x) = n.

Equivalent with:

    n.inverse_sigma.max

Returns C<nil> if there are no solutions.

=cut

=head2 inverse_sigma_min

    n.inverse_sigma_min

Returns the smallest solution x to the sum of divisors function: sigma(x) = n.

Equivalent with:

    n.inverse_sigma.min

Returns C<nil> if there are no solutions.

=cut

=head2 inverse_uphi

    n.inverse_uphi

Returns an array with all the solutions C<x> to C<uphi(x) = n>.

Example:

    say inverse_uphi(120)       #=> [121, 143, 144, 155, 164, 183, 220, 231, 240, 242, 286, 310, 366, 462]

=cut

=head2 inverse_usigma

    n.inverse_usigma

Returns an array with all the possible solutions for C<x> in usigma(x) = n.

Example:

    say inverse_usigma(120)                         #=> [60, 87, 92, 95, 99]
    say inverse_usigma(5040).len                    #=> 38
    say inverse_usigma(5040).map { .usigma }.uniq   #=> [5040]

=cut

=head2 invmod

    n.invmod(m)

Returns the modular inverse of n modulo m.

=cut

=head2 ipolygonal_root

    ipolygonal_root(n,k)

First integer k-gonal root of n.

Example:

    say ipolygonal_root(n, 5)                   # integer pentagonal root
    say ipolygonal_root(polygonal(10, 5), 5)    # prints: "10"

=cut

=head2 ipolygonal_root2

    ipolygonal_root2(n,k)

Second integer k-gonal root of C<n>.

Example:

    say ipolygonal_root2(n, 5)                   # second integer pentagonal root
    say ipolygonal_root2(polygonal(-10, 5), 5)   # prints: "-10"

=cut

=head2 ipow

    ipow(b,n)

Integer exponentiation: C<b^n>.

=cut

=head2 ipow10

    ipow10(n)

Integer exponentiation: C<10^n>.

=cut

=head2 ipow2

    ipow2(n)

Integer exponentiation: C<2^n>.

=cut

=head2 iquadratic_formula

    iquadratic_formula(a,b,c)

Returns a list of integer solutions C<(x_1, x_2)> to the quadratic equation: C<a*x^2 + b*x + c = 0>, defined as:

    floor((-b ± isqrt(b^2 - 4ac)) / (2a))

Example:

    say [iquadratic_formula(13, -42, -34)]  #=> [-3, -1]

Aliases: I<integer_quadratic_formula>

=cut

=head2 irand

    irand(n)
    irand(a,b)

Returns a pseudorandom integer in range C<0..n> (all inclusive), or in the range C<a..b> (all inclusive).

If C<a> is greater than C<b>, the returned integer will be in the range C<b..a>.

Example:

    irand(10)        # a pseudorandom integer in the interval [0, 10]
    irand(10, 20)    # a pseudorandom integer in the interval [10, 20]

=cut

=head2 iroot

    n.iroot(k)

Integer k-th root of n.

=cut

=head2 irootrem

    n.irootrem(k)

Returns a list with the integer k-th root of n and k-th root remainder of n.

Equivalent with:

    (n.iroot(k), n - n.iroot(k)**k)

=cut

=head2 is_abundant

    n.is_abundant()

Returns true when sigma(n) > 2*n.

=cut

=head2 is_aks_prime

    n.is_aks_prime()

Return true if C<n> passes the Agrawal-Kayal-Saxena (AKS) primality test.

=cut

=head2 is_almost_prime

    n.is_almost_prime(k=2)

Return true if C<n> is a k-almost prime (i.e.: true iff C<n> is the product of C<k> not necessarily distinct primes).

Equivalently, k-almost primes are numbers n such that C<bigomega(n) == k>.

Example:

    say 20.by { .is_almost_prime(1) }   # primes
    say 20.by { .is_almost_prime(2) }   # semiprimes
    say 20.by { .is_almost_prime(3) }   # 3-almost primes

=cut

=head2 is_between

    n.is_between(min, max)

Returns a true value when `n >= min` and `n <= max`.

=cut

=head2 is_bpsw_prime

    n.is_bpsw_prime()

Returns true if C<n> passes the B-PSW primality test (extra-strong variant).

=cut

=head2 is_carmichael

    n.is_carmichael()

Returns true if is a Carmichael number.

=cut

=head2 is_chebyshev

    n.is_chebyshev()

Returns true if C<n> is an odd composite Chebyshev pseudoprime, as defined by OEIS A175530.

Aliases: I<is_chebyshev_psp>, I<is_chebyshev_pseudoprime>

=cut

=head2 is_complex

    x.is_complex

Returns true if C<x> is a complex number.

Example:

    say is_complex(complex(4))        # false (is real)
    say is_complex(complex(4i))       # false (is imaginary)
    say is_complex(complex(3+4i))     # true

=cut

=head2 is_composite

    n.is_composite

Returns true if C<n> is a positive > 1 composite number.

=cut

=head2 is_congruent

    n.is_congruent(a,m)

Returns true when n is congruent to a modulo m.

Example:

    say 99923.is_congruent(-2, 5)   #=> true
    say 99923.is_congruent(3, 5)    #=> true

Also defined for rationals, floats and complex numbers:

    say is_congruent(124, 1/4, 3/4) #=> true

=cut

=head2 is_coprime

    is_coprime(a,b)

Returns true if C<gcd(a,b) = 1>.

=cut

=head2 is_cube

    n.is_cube

Return true if n is a cube number (if it can be written as b^3, for some integer b).

=cut

=head2 is_cyclic

    n.is_cyclic

Returns true when C<gcd(phi(n), n) = 1>, where C<phi(n)> is the Euler totient function.

Example:

    say 30.by { .is_cyclic }    # OEIS: A003277

=cut

=head2 is_ecpp_prime

    n.is_ecpp_prime

Return true if C<n> can be proved prime using the Elliptic Curve Primality Proving algorithm.

=cut

=head2 iseed

    iseed(n)

Re-seed the C<irand()> function. The value for C<n> can be any arbitrary large integer.

=cut

=head2 is_euler_psp

    n.is_euler_psp(bases...)

Return true if C<n> is an Euler-Jacobi pseudoprime, given a list of bases.

Aliases: I<is_euler_pseudoprime>

=cut

=head2 is_even

    n.is_even

Returns true if C<n> is an integer divisible by 2.

=cut

=head2 is_fib

    n.is_fib

Returns true if n is a Fibonacci number. False otherwise.

Aliases: I<is_fibonacci>

=cut

=head2 is_fib_psp

    n.is_fib_psp(P=1, Q=-1)

Returns true if C<n> passes the Lucas test to the C<U>, using the parameters C<P> and C<Q>.

Example:

    say 10.by { .is_composite && .is_lucasU_psp }               # Fibonacci pseudoprimes
    say 10.by { .is_odd_composite && .is_lucasU_psp(2, -1) }    #=> OEIS: A327651

Aliases: I<is_lucasu_psp>, I<is_lucasU_psp>, I<is_fibonacci_psp>, I<is_lucasU_pseudoprime>, I<is_fibonacci_pseudoprime>

=cut

=head2 is_frobenius_psp

    n.is_frobenius_psp(a,b)

Return true if C<n> is a Frobenius probable prime with respect to the polynomial C<x^2 - ax + b>.

Aliases: I<is_frobenius_pseudoprime>

=cut

=head2 is_fundamental

    n.is_fundamental

Returns true if C<n> is a fundamental discriminant.

=cut

=head2 is_gaussian_prime

    is_gaussian_prime(a,b)

Returns true if a+b*i is a Gaussian prime.

Example:

    say is_gaussian_prime(3, 4)       #=> false
    say is_gaussian_prime(13, 42)     #=> true

Provided by L<Math::Prime::Util::GMP> >= 0.52.

=cut

=head2 is_imag

    x.is_imag

Returns true if C<x> is an imaginary number.

Example:

    say is_imag(complex(4))           # false (is real)
    say is_imag(complex(4i))          # true
    say is_imag(complex(3+4i))        # false (is complex)

=cut

=head2 is_imprimitive_carmichael

    n.is_imprimitive_carmichael

Returns true if C<n> is an imprimitive Carmichael numbers, as defined by OEIS: A328935.

Example:

    say 325533792014488126487416882038879701391121.is_imprimitive_carmichael   # true

The method efficiently tries to factorize large Carmichael numbers, using the C<miller_factor(n)> method.

=cut

=head2 is_inf

    x.is_inf

Returns true if C<x> equals positive infinity (C<Inf>).

=cut

=head2 is_int

    x.is_int

Returns true if C<x> is an integer.

=cut

=head2 is_khashin_psp

    n.is_khashin_psp

Return true if C<n> passes the Frobenius test of Sergey Khashin.

Aliases: I<is_khashin_pseudoprime>, I<is_frobenius_khashin_psp>, I<is_frobenius_khashin_pseudoprime>

=cut

=head2 is_lucas

    n.is_lucas

Returns true if n is a Lucas number. False otherwise.

=cut

=head2 is_lucas_carmichael

    n.is_lucas_carmichael

Returns true if C<n> is a Lucas-Carmichael number.

Example:

    say 10.by(:is_lucas_carmichael)                               # OEIS: A006972
    say is_lucas_carmichael(58735331016965175152455996272482303)  # true

=cut

=head2 is_lucas_psp

    n.is_lucas_psp

Return true if C<n> is a Lucas pseudoprime.

Aliases: I<is_lucas_pseudoprime>

=cut

=head2 is_lucasv_psp

    n.is_lucasv_psp(P=1, Q=-1)

Returns true if C<n> passes the Lucas test to the C<V> sequence, using the paramters C<P> and C<Q>.

Example:

    say 10.by { .is_composite && .is_lucasV_psp }               # Bruckman-Lucas pseudoprimes
    say 10.by { .is_odd_composite && .is_lucasV_psp(2, -1) }    #=> OEIS: A330276

Aliases: I<is_lucasV_psp>, I<is_bruckman_lucas_psp>, I<is_lucasV_pseudoprime>, I<is_bruckman_lucas_pseudoprime>

=cut

=head2 is_mersenne_prime

    is_mersenne_prime(p)

Returns true if C<2^p - 1> is a Mersenne prime.

Example:

    say 607.is_mersenne_prime       # prints true; (2**607 - 1) is a Mersenne prime.

=cut

=head2 is_mone

    x.is_mone

Returns true if C<x> equals -1.

=cut

=head2 is_nan

    x.is_nan

Returns true if C<x> holds the Not-a-Number special value (C<NaN>).

=cut

=head2 is_neg

    x.is_neg

Returns true if C<x> is negative.

Aliases: I<is_negative>

=cut

=head2 is_ninf

    x.is_ninf

Returns true if C<x> equals negative infinity (C<-Inf>).

=cut

=head2 is_nm1_prime

    n.is_nm1_prime

Return true if C<n> can be proved prime using the factorization of C<n-1>.

Aliases: I<is_pm1_prime>, I<is_nminus1_prime>

=cut

=head2 is_np1_prime

    n.is_np1_prime

Return true if C<n> can be proved prime using the factorization of C<n+1>.

Aliases: I<is_pp1_prime>, I<is_nplus1_prime>

=cut

=head2 is_odd

    n.is_odd

Returns true when C<n> is an integer not divisible by 2.

=cut

=head2 is_odd_composite

    n.is_odd_composite

Returns true when C<n> is an odd composite integer.

=cut

=head2 is_omega_prime

    n.is_omega_prime(k=2)

Return true if C<n> is a k-omega prime (i.e.: true if C<n> is divisible by exactly C<k> different primes).

Equivalently, k-omega primes are numbers n such that C<omega(n) == k>.

Example:

    say 20.by { .is_omega_prime(1) }   # prime powers
    say 20.by { .is_omega_prime(2) }   # numbers n such that omega(n) == 2

=cut

=head2 is_one

    n.is_one

Returns true when C<n> equals 1.

=cut

=head2 is_over_psp

    n.is_over_psp

Returns true if C<n> is an overpseudoprime to base C<b>. Multiple bases can also be provided.

An overpseudoprime to base C<b> is a also a strong Fermat pseudoprime to base C<b> and a super-pseudoprime to base C<b>, where znorder(b,n) == znorder(b,p) for every p|n.

Example:

    say 10.by { .is_composite && .is_over_psp }         # overpseudoprimes to base 2
    say 10.by { .is_composite && .is_over_psp(3) }      # overpseudoprimes to base 3

Aliases: I<is_over_pseudoprime>, I<is_overpseudoprime>

=cut

=head2 is_palindrome

    n.is_palindrome(b=10)

Returns true if the given number C<n> is palindromic in the given base C<b>. When no base is given, it defaults to 10.

Example:

    # Numbers that are palindromic in bases 2 and 10 (OEIS: A007632)
    say 1e6.range.grep{ .is_palindrome(2) && .is_palindrome(10) }

Aliases: I<is_palindromic>

=cut

=head2 is_pell_lucas_psp

    n.is_pell_lucas_psp

It returns true if C<V_n(2, -1) = 2 (mod n)>.

Example:

    say 20.by { .is_pell_lucas_pseudoprime }                               #=> OEIS: A270342
    say 20.by { .is_pell_lucas_pseudoprime && .is_composite }              #=> OEIS: A335668
    say 20.by { .is_pell_lucas_pseudoprime && .is_composite && .is_odd }   #=> OEIS: A330276

Aliases: I<is_pell_lucas_pseudoprime>

=cut

=head2 is_pell_psp

    n.is_pell_psp

These are odd numbers that satisfy:

    U_n(2, -1) = (2|n) (mod n)

Example:

    say 10.by { .is_pell_pseudoprime && .is_composite }  # OEIS: A099011

Aliases: I<is_pell_pseudoprime>

=cut

=head2 is_perrin_psp

    n.is_perrin_psp

Returns true if C<n> passes the Perrin primality test.

Aliases: I<is_perrin_pseudoprime>

=cut

=head2 is_plumb_psp

    n.is_plumb_psp

Return true if C<n> passes Colin Plumb's Euler Criterion primality test.

Aliases: I<is_euler_plumb_psp>, I<is_plumb_pseudoprime>, I<is_euler_plumb_pseudoprime>

=cut

=head2 is_polygonal

    is_polygonal(n,k)

Returns true if C<n> is a first k-gonal number.

Example:

    say is_polygonal(145, 5)      #=> 1 ("145" is a pentagonal number)
    say is_polygonal(155, 5)      #=> 0

=cut

=head2 is_polygonal2

    is_polygonal2(n,k)

Returns true when C<n> is a second k-gonal number.

Example:

    say is_polygonal2(145, 5)      #=> 0
    say is_polygonal2(155, 5)      #=> 1 ("155" is a second-pentagonal number)

=cut

=head2 is_pos

    x.is_pos

Returns true when C<x> is a positive integer.

Aliases: I<is_positive>

=cut

=head2 is_pow

    n.is_power
    n.is_power(k)

When k is omitted, it true if n is a perfect power.

When k is given, it returns true if n can be expressed as C<n = b^k> for some b >= 1.

Example:

    say 225.is_power          # true: 225 == 15**2
    say 100.is_power(2)       # true: 100 is square (10**2)
    say 125.is_power(3)       # true: 125 is a cube ( 5**3)

Aliases: I<is_power>, I<is_perfect_power>

=cut

=head2 is_powerfree

    n.is_powerfree(k=2)

Returns true when all the exponents in the prime-power factorization of n are < k.

=cut

=head2 is_powerful

    n.is_powerful(k=2)

Returns true when all the exponents in the prime-power factorization of n are >= k.

=cut

=head2 is_power_of

    n.is_power_of(b)

Return true if n is a power of b, such that n = b^k for some k >= 0:

    n.is_power_of(b)    # true if n == b^k for some k >= 0

Example:

    say 1000.range.grep { .is_power_of(2) }     # powers of 2
    say 1000.range.grep { .is_power_of(5) }     # powers of 5

=cut

=head2 is_practical

    n.is_practical()

Returns true if n is a practical number (OEIS: A005153).

Example:

    say 20.by { .is_practical }

=cut

=head2 is_prime

    n.is_prime

Returns true if C<n> is a prime number.

=cut

=head2 is_prime_power

    n.is_prime_power

Returns true if n is a power of some prime.

=cut

=head2 is_primitive_root

    n.is_primitive_root(m)

Returns true if n is a primitive root modulo m.

=cut

=head2 is_prob_prime

    is_prob_prime(n)

Returns true if n is probably prime. The method does some trial division, then it performs a B-PSW test.

=cut

=head2 is_prob_squarefree

    n.is_prob_squarefree(k)

Returns true iff C<n> is not a power and is not divisible by a square <= k.

=cut

=head2 is_prov_prime

    is_prov_prime(n)

Returns true if n is provable prime.

Aliases: I<is_provable_prime>

=cut

=head2 is_psp

    n.is_psp(bases...)

Returns true if C<n> is a Fermat pseudoprime to the provided bases.

Aliases: I<is_fermat_psp>, I<is_pseudoprime>, I<is_fermat_pseudoprime>

=cut

=head2 isqrt

    n.isqrt

Integer square root of n.

=cut

=head2 isqrtrem

    n.isqrtrem

Returns a list with the integer square root of n and the square root remainder of n.

Equivalent with:

    (n.isqrt, n - n.isqrt**2)

=cut

=head2 is_rat

    x.is_rat

Returns true if C<x> is a rational number.

=cut

=head2 is_real

    x.is_real

Returns true if C<x> is a real number.

    say is_real(complex(4))           # true
    say is_real(complex(4i))          # false (is imaginary)
    say is_real(complex(3+4i))        # false (is complex)

=cut

=head2 is_rough

    n.is_rough(k)

Returns true if all prime factors of n are >= k.

Example:

    say 30.by { .is_rough(3) }  #=> OEIS: A005408
    say 30.by { .is_rough(5) }  #=> OEIS: A007310
    say 30.by { .is_rough(7) }  #=> OEIS: A007775
    # ...
    say 30.by { .is_rough(23) } #=> OEIS: A166063

=cut

=head2 is_safe_prime

    n.is_safe_prime

It returns true if both C<n> and C<(n-1)/2> are prime.

Example:

    say 30.by { .is_safe_prime }    #=> OEIS: A005385

=cut

=head2 is_semiprime

    n.is_semiprime

Returns true if n has exactly two prime factors (not necessarily distinct).

=cut

=head2 is_smooth

    n.is_smooth(k)

Returns true if all the prime factors of C<n> are <= k. False otherwise.

=cut

=head2 is_smooth_over_prod

    n.is_smooth_over_prod(k)

Returns true when C<n> is smooth over the prime factors of C<k>.

=cut

=head2 is_sqr

    n.is_sqr
    is_square(n)

Returns true if C<n> is a perfect square integer.

Aliases: I<is_square>, I<is_perfect_square>

=cut

=head2 is_square_free

    n.is_square_free

Returns true if the prime factorization of n does not include duplicated factors (i.e.: n is not divisible by a square).

Aliases: I<is_squarefree>

=cut

=head2 is_stronger_lucas_psp

    n.is_stronger_lucas_psp

Return true if C<n> is an extra-strong Lucas pseudoprime.

Aliases: I<is_extra_strong_lucas_psp>, I<is_stronger_lucas_pseudoprime>, I<is_extra_strong_lucas_pseudoprime>

=cut

=head2 is_strong_fib

    n.is_strong_fib

Teturns true if C<n> is a strong Fibonacci pseudoprime, satisfying:

    V_n(P,Q) = P (mod)

for Q = -1 and all P.

Odd composite integer n is a strong Fibonacci pseudoprime iff:

    1) n is a Carmichael number: p-1 | n-1
    2) 2(p + 1) | (n − 1) or 2(p + 1) | (n − p)

for each prime p|n.

Example:

    say is_strong_fibonacci_pseudoprime(443372888629441)    #=> true
    say is_strong_fibonacci_pseudoprime(39671149333495681)  #=> true

Aliases: I<is_strong_fib_psp>, I<is_strong_fibonacci>, I<is_strong_fibonacci_psp>, I<is_strong_fibonacci_pseudoprime>

=cut

=head2 is_strongish_lucas_psp

    n.is_strongish_lucas_psp

Return true if C<n> is almost an extra-strong Lucas pseudoprime.

Aliases: I<is_strongish_lucas_pseudoprime>

=cut

=head2 is_strong_lucas_psp

    n.is_strong_lucas_psp

Return true if C<n> is a strong Lucas pseudoprime.

Aliases: I<is_strong_lucas_pseudoprime>

=cut

=head2 is_strong_psp

    n.is_strong_psp(bases...)

Return true if C<n> is a strong pseudoprime.

Aliases: I<miller_rabin>, I<is_strong_fermat_psp>, I<is_strong_pseudoprime>, I<is_strong_fermat_pseudoprime>

=cut

=head2 is_super_psp

    n.is_super_psp(bases...)

It returns true if the given value of n is a super-pseudoprime to the given bases. When no base is given, the base 2 is used (which represents the Super-Poulet numbers: A050217)

Example:

    # Super-Poulet numbers (OEIS: A050217)
    say 1e4.range.grep { .is_super_pseudoprime }.grep{ .is_composite }

    # Super-Poulet numbers to base 3 (OEIS: A328662)
    say 1e4.range.grep { .is_super_pseudoprime(3) }.grep{ .is_composite }

    # Super-Poulet numbers to base 2 and 3
    say 1e5.range.grep { .is_super_pseudoprime(2, 3) }.grep{ .is_composite }

Aliases: I<is_super_pseudoprime>, I<is_superpseudoprime>

=cut

=head2 is_totient

    n.is_totient

Given an integer C<n>, returns true if there exists an integer C<x> such that C<euler_phi(x) == n>.

=cut

=head2 isub

    isub(a,b)

Integer subtraction: C<a-b>.

=cut

=head2 is_underwood_psp

    n.is_underwood_psp

Return true if C<n> passes the efficient Frobenius test of Paul Underwood.

Aliases: I<is_underwood_pseudoprime>, I<is_frobenius_underwood_psp>, I<is_frobenius_underwood_pseudoprime>

=cut

=head2 is_zero

    x.is_zero

Returns true when C<x> equals 0.

=cut

=head2 jacobi

    jacobi(a,n)

Returns the Jacobi symbol: C<(a|n)>.

=cut

=head2 jordan_totient

    jordan_totient(n,k)

Jordan's totient C<J_k(n)>, which is a generalization of Euler's totient function.

=cut

=head2 kronecker

    kronecker(a,n)

Returns the Kronecker symbol: C<(a|n)>.

=cut

=head2 laguerre

    laguerre(n)
    laguerre(n, x)

Laguerre polynomials: C<L_n(x)>.

When C<x> is ommited, a Polynomial object is returned.

Aliases: I<laguerreL>, I<laguerre_polynomial>

=cut

=head2 lambda

    lambda(n)

Carmichael lambda function: C<λ(n)>, defined as the smallest positive integer C<m> such that:

    powmod(a, m, n) == 1

for every integer C<a> between C<1> and C<n> that is coprime to C<n>.

Alias: I<carmichael_lambda>.

=cut

=head2 lambert_w

    lambert_w(x)

The Lambert-W function. When the value of C<x> is less than C<-1/e>, it returns a complex number.

    It also accepts a complex number as input.

Identities (assuming x>0):

    LambertW(exp(x)*x) = x
    LambertW(log(x)*x) = log(x)

=cut

=head2 lcm

    lcm(...)

Least common multiple of a list of integers.

=cut

=head2 legendre

    legendre(a,p)

Returns the Legendre symbol: C<(a|p)>.

=cut

=head2 legendreP

    legendreP(n)
    legendreP(n,x)

Legendre polynomials: C<P_n(x)>.

When C<x> is ommited, a Polynomial object is returned.

Aliases: I<legendrep>, I<legendre_polynomial>

=cut

=head2 legendre_phi

    n.legendre_phi(k)

Returns the count of numbers <= C<n> that are not divisible by the first C<k> primes.

Equivalent with:

    prime(k+1).rough_count(n)

=cut

=head2 len

    n.len(b=10)

Returns the number of digits of the integer part of C<n> in a given base..

    say 5040.len        #=> 4
    say 5040.len(2)     #=> 13

Aliases: I<size>, I<length>

=cut

=head2 lgamma

    lgamma(x)

Natural logarithm of C<abs(Γ(x))>.

Aliases: I<gamma_abs_log>

=cut

=head2 lgrt

    lgrt(x)

Returns the "logarithm-root" of C<x>, such that C<lgrt(x) ** lgrt(x) =~= x>.

Example:

    say lgrt(100)       #=> 3.59728502354041750549765225178229
    say lgrt(-100)      #=> 3.70202936660214594290193962952737+1.34823128471151901327831464969872i

=cut

=head2 li

    li(x)

Returns the logarithmic integral of x.

Example:

    say 100.li  # prints: 30.12614158...

=cut

=head2 li2

    li2(x)

Dilogarithm function, defined as the integral of C<-log(1-t)/t> from 0 to x.

=cut

=head2 liouville

    liouville(n)

The Liouville function.

Equivalent to:

    (-1)**omega(n)

=cut

=head2 liouville_sum

    n.liouville_sum

Computes partial sums of the Liouville lambda function.

Example:

    say liouville_sum(10**9)    #=> -25216
    say liouville_sum(10**10)   #=> -116026

=cut

=head2 ln

    x.ln

Natural logarithm of C<x>.

=cut

=head2 ln2

    Num.ln2

Returns the natural logarithm of 2 constant.

=cut

=head2 lnbern

    lnbern(n)

Returns the natural logarithm of the n-th Bernoulli number.

Aliases: I<bern_log>, I<lnbernreal>, I<bernoulli_log>

=cut

=head2 lngamma

    lngamma(x)

Natural logarithm of C<Γ(x)>.

Aliases: I<gamma_log>

=cut

=head2 lnsuperfactorial

    lnsuperfactorial(n)

Natural logarithm of C<superfactorial(n)>.

Aliases: I<superfactorial_ln>, I<superfactorial_log>

=cut

=head2 log

    log(x)
    log(x, b)

Natural logarithm of C<x> to base I<e>, or to a given base C<b>.

=cut

=head2 log10

    x.log10

Logarithm of C<x> to base 10.

=cut

=head2 log2

    x.log2

Logarithm of C<x> to base 2.

=cut

=head2 logarithmic_derivative

    logarithmic_derivative(n)

Return the logarithmic derivative of C<n>, defined as:

    derivative(n)/n

=cut

=head2 lpf

    lpf(n)

Returns the least prime factor of n.

Defined with the base-cases:

    lpf(0) = 0
    lpf(1) = 1

Example:

    say lpf(fibonacci(1234))    #=> 234461

=cut

=head2 lsb

    lsb(n)

Returns the least significant bit of n.

Example:

    say 0b110010101111000000.lsb    # 6

=cut

=head2 lucas

    lucas(n)

Returns the n-th Lucas number.

=cut

=head2 lucas_factor

    n.lucas_factor(j=1, tries=100)

Effective in factoring Carmichael numbers, Fermat pseudoprimes, Lucas pseudoprimes and Lucas-Carmichael numbers.

Example:

    say lucas_factor(2425361208749736840354501506901183117777758034612345610725789878400467)

Aliases: I<lucas_miller_factor>

=cut

=head2 lucas_mod

    lucasmod(n,m)

Efficiently compute the n-th Lucas number modulo m.

Aliases: I<lucasmod>

=cut

=head2 lucasU

    lucasU(P, Q, n)

The Lucas C<U_n(P, Q)> function.

    say 20.of{|n| lucasU(1, -1, n) }    # the Fibonacci numbers
    say 20.of{|n| lucasU(2, -1, n) }    # the Pell numbers
    say 20.of{|n| lucasU(1, -2, n) }    # the Jacobsthal numbers

Aliases: I<lucasu>

=cut

=head2 lucasUmod

    lucasUmod(P,Q,n,m)

Efficiently compute the Lucas C<U_n(P, Q)> function modulo m.

Aliases: I<lucasumod>

=cut

=head2 lucasUVmod

    lucasUVmod(P,Q,n,m)

Efficiently compute the Lucas C<U_n(P, Q)> and C<V_n(P,Q)> functions modulo m.

Equivalent with:

    (lucasUmod(P,Q,n,m), lucasVmod(P,Q,n,m))

Aliases: I<lucasuvmod>

=cut

=head2 lucasV

    lucasV(P, Q, n)

The Lucas C<V_n(P, Q)> function.

    say 20.of{|n| lucasV(1, -1, n) }    # the Lucas numbers
    say 20.of{|n| lucasV(2, -1, n) }    # the Pell-Lucas numbers
    say 20.of{|n| lucasV(1, -2, n) }    # the Jacobsthal-Lucas numbers

Aliases: I<lucasv>

=cut

=head2 lucasvmod

    lucasVmod(P,Q,n,m)

Efficiently compute the Lucas C<V_n(P, Q)> function modulo m.

Aliases: I<lucasVmod>

=cut

=head2 make_coprime

    n.make_coprime(k)

Returns the largest divisor of n that is coprime to k.

=cut

=head2 mangoldt

    mangoldt(n)

The Mangoldt function. For the exponential values, see C<exp_mangoldt>.

=cut

=head2 max

    max(...)

Returns the maximum value from a list of numbers.

=cut

=head2 mertens

    mertens(n)
    mertens(a,b)

Returns the Mertens functions, which is defined as the partial sums of the Moebius function:

    Sum_{k=1..n} moebius(k)

When an additional argument is given, the returned result is:

    Sum_{k=a..b} moebius(k)

Example:

    say mertens(100000)     # equilvanet with: (1..100000 -> sum { .moebius })
    say mertens(21, 123)    # equivalent with: (21..123 -> sum { .moebius })

=cut

=head2 mfac

    mfac(n,k)
    mfactorial(n,k)

The generalized multi-factorial of n.

Example:

    say 15.of { .mfac(2) }    # double-factorials (OEIS: A006882)
    say 15.of { .mfac(3) }    # triple-factorials (OEIS: A007661)

Aliases: I<mfactorial>, I<multi_factorial>

=cut

=head2 miller_factor

    n.miller_factor(tries=100)

Effective in factoring Carmichael numbers and Fermat pseudoprimes.

It returns an array with the factors of C<n>. However, sometimes, not all factors are prime.

Example:

    say miller_factor(58571442634534443082821160508299574798027946748324125518533225605795841)

Aliases: I<miller_rabin_factor>

=cut

=head2 miller_rabin_random

    n.miller_rabin_random(k)

Return true if C<n> passes the Miller-Rabin primality test with C<k> random bases.

=cut

=head2 min

    min(...)

Returns the smallest value from a list of numbers.

=cut

=head2 mone

    Num.mone

Returns the -1 value.

=cut

=head2 motzkin

    motzkin(n)

Returns the n-th Motzkin number. (OEIS: A001006)

Example:

    say 10.of { .motzkin }   #=> [1, 1, 2, 4, 9, 21, 51, 127, 323, 835]

=cut

=head2 msb

    msb(n)

Returns the most significant bit of n.

Example:

    say 0b110010101111000000.msb    # 17

=cut

=head2 mulmod

    mulmod(a,b,m)

Modular integer multiplication: C<(a*b) % m>.

Example:

    say mulmod(43, 97, 127)     # == (43*97)%127

=cut

=head2 multinomial

    multinomial(...)

The multinomial coefficient, given a list of native integers.

Example:

    say multinomial(1, 4, 4, 2)     #=> 34650

=cut

=head2 nan

    Num.nan

Returns the Not-a-Number special value (C<NaN>).

=cut

=head2 n_composites

    n_composites(n, start=4)

Returns n consecutive composite numbers starting from C<start>.

Example:

    say n_composites(10)        #=> [4, 6, 8, 9, 10, 12, 14, 15, 16, 18]
    say n_composites(5, 50)     #=> [50, 51, 52, 54, 55]

Aliases: I<ncomposites>, I<next_composites>

=cut

=head2 nd

    n.st({...})
    n.nd({...})
    n.rd({...})
    n.th({...})

It returns the n-th value for which the provided block evaluates to a true value, starting couting from 0.

Example:

    say 100.th { .is_prime }    # 100-th prime

Also aliased as C<.st>, C<.nd> and C<.rd>:

    say 1.st { .is_prime }      # first prime
    say 2.nd { .is_prime }      # second prime
    say 3.rd { .is_prime }      # third prime

Aliases: I<rd>, I<st>, I<th>

=cut

=head2 neg

    x.neg

Negates the sign of <Cx> (equivalent with: C<-x>).

=cut

=head2 new

    Number(string, base=10)
    Num.new(string, base=10)

Create a new Number object, given a string and a base.

Aliases: I<call>

=cut

=head2 next_composite

    n.next_composite

Given a non-negative integer C<n>, it returns the next composite number after C<n>.

=cut

=head2 next_palindrome

    n.next_palindrome(b=10)

Efficiently returns the next palindrome in base-b greater than n.

Example:

    # Iterate over the base-10 palindromic numbers < 10^6
    for (var n = 0; n < 1e6; n = n.next_palindrome) {
        say n
    }

=cut

=head2 next_pow

    n.next_pow(b)

Returns the perfect power after C<n>, with base C<b>.

Equivalent with:

    b**(1+ilog(n,b))

Aliases: I<next_power>

=cut

=head2 next_pow2

    n.next_pow2

Returns the next power of 2 larger than C<n>.

Equivalent with:

    n.next_pow(2)

Aliases: I<next_power2>

=cut

=head2 next_prime

    n.next_prime

Returns the next prime larger than C<n>.

=cut

=head2 next_twin_prime

    n.next_twin_prime

Returns next twin prime number larger than C<n>.

Provided by L<Math::Prime::Util::GMP> >= 0.52.

=cut

=head2 ninf

    Num.ninf

Returns the negative infinity special value (C<-Inf>).

=cut

=head2 nok

    nok(n,k)
    binomial(n,k)

Returns the binomial coefficient C<n> over C<k>, also called the "choose" function.

Equivalent to:

                        n!
    binomial(n, k) = --------
                     k!(n-k)!

Aliases: I<binomial>

=cut

=head2 norm

    norm(x)

Returns the normalized value of C<x>: C<abs(x)^2>.

=cut

=head2 n_primes

    n_primes(n, start=2)

Returns an array containing C<n> consecutive primes that are `>= start` (if omitted, then start = 2).

Example:

    say n_primes(10)        #=> [2, 3, 5, 7, 11, 13, 17, 19, 23, 29]
    say n_primes(10, 1000)  #=> [1009, 1013, 1019, 1021, 1031, 1033, 1039, 1049, 1051, 1061]

The value for C<start> can be any arbitary large integer:

    say n_primes(10, 2**128 + 1)    # first 10 primes > 2^128

Aliases: I<nprimes>, I<next_primes>

=cut

=head2 nu

    r.nu
    r.numerator

Returns the numerator of rational number C<r>.

Example:

    say numerator(43/97)     #=> 43

Aliases: I<numerator>

=cut

=head2 nude

    nude(r)

Returns a list with the numerator and the denominator of a rational number C<r>.

Example:

    say [nude(43/97)]       #=> [43, 97]

=cut

=head2 num2perm

    num2perm(n,k)

Given a non-negative integer C<n> and integer C<k>, return the rank C<k> lexicographic permutation of C<n> elements.

C<k> will be interpreted as mod C<n!>.

The inverse value for C<k> is given by C<Array#perm2num>:

    say num2perm(5, 43)             #=> [1, 4, 0, 3, 2]
    say num2perm(5, 43).perm2num    #=> 43

=cut

=head2 numify

    n.numify

Returns a raw native number representation for the self-number.

Can be used for assigning values to C<Num!PREC> variable.

Example:

    local Num!PREC = 42.numify  # set precision to 42 bits
    say sqrt(2)                 # 1.414213562

Native numbers can also be used when indexing an array:

    var arr = [42, 43, 44]
    var idx = 2.numify
    say arr[idx]            #=> 44

Although this may or may not be actually faster.

=cut

=head2 of

    n.of {|k| ... }

Returns an array with n elements mapped to the given block. The block is called with k = 0,1,2...,n-1.

Example:

    say 10.of { _*_ }       #=> first 10 squares
    say 10.of { .fib }      #=> first 10 Fibonacci numbers

=cut

=head2 omega_prime_divisors

    n.omega_prime_divisors
    n.omega_prime_divisors(k)

Returns the k-omega prime divisors of n.

Example:

    say 5040.omega_prime_divisors(4)    #=> [210, 420, 630, 840, 1260, 1680, 2520, 5040]

When C<k> is omitted, an array of arrays with the k-omega prime divisors of n, for each k in the range C<0..omega(n)>, is returned:

    say 120.omega_prime_divisors        #=> [[1], [2, 3, 4, 5, 8], [6, 10, 12, 15, 20, 24, 40], [30, 60, 120]]

=cut

=head2 omega_prime_count

    k.omega_prime_count(n)
    k.omega_prime_count(a,b)

Returns the number of k-omega primes <= n, or in the range C<a..b>.

Example:

    say 1.omega_prime_count(100)            # number prime powers <= 100
    say 2.omega_prime_count(50, 100)        # number of 2-omega primes in range 50..100

Aliases: I<omega_primepi>

=cut

=head2 omega_primes

    k.omega_primes(n)
    k.omega_primes(a,b)

Returns an array with k-omega primes <= n, or in the range C<a..b>.

k-omega primes are numbers C<n> such that C<omega(n) == k>.

Example:

    say 1.omega_primes(100)         # prime powers <= 100
    say 2.omega_primes(50, 100)     # 2-omega primes in range 50..100

=cut

=head2 one

    Num.one

Returns the 1 value.

=cut

=head2 partitions

    partitions(n)

Returns the number of partitions of n.

Aliases: I<partition_number>

=cut

=head2 parts

    x.parts

Returns an array with the real and imaginary parts of C<x>.

Example:

    say parts(5)        #=> [5, 0]
    say parts(3+4i)     #=> [3, 4]

=cut

=head2 pbrent_factor

    n.pbrent_factor(reps)

Pollard-Brent rho factorization method.

=cut

=head2 perfect_power

    n.perfect_power

Returns the largest power C<k> of C<n> for which there exists an integer C<r>, such that: C<n = r^k>.

Example:

    say perfect_power(15**5)    #=> 5

=cut

=head2 perfect_root

    n.perfect_root

Returns the smallest root C<r> of C<n> for which there exists an integer C<k>, such that: C<n = r^k>.

Example:

    say perfect_root(15**5)    #=> 15

=cut

=head2 permutations

    n.permutations

Returns an array of arrays with the permutations of the integers in the range C<0..n-1>, or iterates over the permutations when a block is given.

Example:

    5.permutations {|*a| say a }

=cut

=head2 pi_k

    k.almost_prime_count(n)
    k.almost_prime_count(a,b)

Returns the number of k-almost primes <= n, or in the range C<a..b>.

Example:

    say 1.almost_prime_count(100)           # count of primes <= 100
    say 2.almost_prime_count(50, 100)       # count of semiprimes in range 50..100

Aliases: I<almost_primepi>, I<almost_prime_count>

=cut

=head2 pm1_factor

    n.pm1_factor(B)

Pollard p-1 factorization method.

Aliases: I<pminus1_factor>

=cut

=head2 pn_primes

    pn_primes(n)
    pn_primes(a,b)

Returns the first n prime numbers, or the primes in the range prime(a)..prime(b).

Example:

    say pn_primes(25)        # the first 25 primes
    say pn_primes(100, 110)  # the primes from 100-th prime to 110-th prime (inclusive)

=cut

=head2 pn_primorial

    pn_primorial(n)

Returns the product of the first n primes.

=cut

=head2 polygonal

    polygonal(n,k)

Returns the n-th k-gonal number. When C<n> is negative, it returns the second k-gonal number.

Example:

   say 10.of {|n| polygonal( n, 3) }  # triangular numbers
   say 10.of {|n| polygonal( n, 5) }  # pentagonal numbers
   say 10.of {|n| polygonal(-n, 5) }  # second pentagonal numbers

=cut

=head2 polygonal_root

    polygonal_root(n,k)

Returns the k-gonal root of C<n>. Also defined for complex numbers.

Example:

    say polygonal_root(n, 3)      # triangular root
    say polygonal_root(n, 5)      # pentagonal root

=cut

=head2 polygonal_root2

    polygonal_root2(n,k)

Returns the second k-gonal root of C<n>. Also defined for complex numbers.

Example:

    say polygonal_root2(n, 5)      # second pentagonal root

=cut

=head2 polymod

    n.polymod(...)

Returns a list of mod results corresponding to the divisors in given list.

Example:

    say [120.polymod(10)]    # (0, 12)
    say [120.polymod(10,10)] # [0, 2, 1]

Particularly useful for:

    var (sec, min, hours, days) = seconds.polymod(60, 60, 24)

=cut

=head2 popcount

    n.popcount

Number of 1's in binary representation of n.

This value is also known as the Hamming weight value.

Aliases: I<hammingweight>

=cut

=head2 power_count

    n.power_count

Returns the number of perfect powers <= n.

Example:

    say power_count(10**6)      #=> 1111
    say power_count(10**20)     #=> 10004650118

Aliases: I<perfect_power_count>

=cut

=head2 power_divisors

    power_divisors(n)
    power_divisors(n, k)

Returns the power divisors of n.

Example:

    say 5040.power_divisors    #=> [1, 4, 8, 9, 16, 36, 144]

Equivalent with:

    n.divisors.grep { .is_power }

When an additional argument is given, it returns the k-th power divisors of n.

Example:

    say 10!.power_divisors(3)  #=> [1, 8, 27, 64, 216, 1728]
    say 10!.power_divisors(4)  #=> [1, 16, 81, 256, 1296, 20736]

Equivalent with:

    n.divisors.grep { .is_power(k) }

=cut

=head2 powerfree_count

    k.powerfree_count(n)
    k.powerfree_count(a,b)

It efficiently counts the number of k-powerfree numbers <= n, or in the range C<a..b>:

    say 3.powerfree_count(100)       #=> count of cubefree numbers <= 100
    say 3.powerfree_count(50, 100)   #=> count of cubefree numbers in the range 50..100

=cut

=head2 powerfree_part

    k.powerfree_part(n)

Returns the k-powerfree part of C<n>.

    say 30.of { 2.powerfree_part(_) }       # squarefree part (OEIS: A007913)
    say 30.of { 3.powerfree_part(_) }       # cubefree part (OEIS: A050985)

=cut

=head2 powerfree_part_sum

    k.powerfree_part_sum(n)
    k.powerfree_part_sum(a,b)

Returns the sum of the k-powerfree parts of integers <= n or in the range C<a..b>:

    say 2.powerfree_part_sum(100)            # sum of squarefree part of n for n <= 100
    say 3.powerfree_part_sum(50, 100)        # sum of cubefree part of n for n in the range 50..100

=cut

=head2 powerfree_sum

    k.powerfree_sum(n)
    k.powerfree_sum(a,b)

Returns the sum of k-powerfree numbers <= n or in the range C<a..b>:

    say 2.powerfree_sum(100)            # sum of squarefree numbers <= 100
    say 3.powerfree_sum(50, 100)        # sum of cubefree numbers in the range 50..100

=cut

=head2 powerful

    k.powerful(n)
    k.powerful(a,b)

It efficiently generates all the k-powerful numbers <= n, or in the range C<a..b>.

Example:

    say 2.powerful(100)          #=> 2-powerful numbers <= 100
    say 2.powerful(50, 100)      #=> 2-powerful numbers in the range 50..100

=cut

=head2 powerful_count

    k.powerful_count(n)
    k.powerful_count(a,b)

It efficiently counts the number of k-powerful numbers <= n, or in the range C<a..b>;

    say 2.powerful_count(100)       #=> count of 2-powerful numbers <= 100
    say 2.powerful_count(50, 100)   #=> count of 2-powerful numbers in the range 50..100

=cut

=head2 power_udivisors

    power_udivisors(n)
    power_udivisors(n,k)

Returns the unitary power divisors of C<n>, such that each divisors C<d> is a perfect power and C<gcd(n/d, d) = 1>.

Example:

    say 15!.power_udivisors         #=> [1, 49, 125, 729, 2048, 35721, 91125]

Equivalent with:

    n.udivisors.grep { .is_power }

When an additional argument is given, it returns the unitary k-th power divisors of C<n>.

Example:

    say 15!.power_udivisors(2)      #=> [1, 49, 729, 35721]
    say 15!.power_udivisors(3)      #=> [1, 125, 729, 91125]

Equivalent with:

    n.udivisors.grep { .is_power(k) }

Aliases: I<power_unitary_divisors>, I<unitary_power_divisors>

=cut

=head2 pp1_factor

    n.pp1_factor(B)

Williams' p+1 factorization method.

Aliases: I<pplus1_factor>

=cut

=head2 prev_prime

    n.prev_prime

Returns the previous prime smaller than n.

=cut

=head2 prho_factor

    n.prho_factor(reps)

Pollard rho factorization method.

=cut

=head2 primality_pretest

    n.primality_pretest()

The method returns true when C<n> passes the internal primality pretest (when n is large enough and it has no small factors).

=cut

=head2 prime

    prime(n)

Returns the n-th prime number.

Example:

    say prime(100)      #  100th prime: 541
    say prime(1000)     # 1000th prime: 7919

Aliases: I<nth_prime>

=cut

=head2 prime_divisors

    prime_divisors(n)

Returns the unique prime factors of C<n>.

=cut

=head2 prime_lower

    n.prime_lower

Lower bound for the n-th prime.

Aliases: I<nth_prime_lower>

=cut

=head2 primepi

    primepi(n)
    primepi(a,b)

Returns the number of primes <= n, or in the range C<a..b>.

Aliases: I<prime_count>, I<count_primes>

=cut

=head2 primepi_lower

    n.primepi_lower

Lower bound for prime_count(n).

Aliases: I<prime_count_lower>

=cut

=head2 primepi_upper

    n.primepi_upper

Upper bound for prime_count(n).

Aliases: I<prime_count_upper>

=cut

=head2 prime_power

    n.prime_power

Returns the exponent C<k> if C<n> is a power of the form C<n = p^k> for some prime C<p>. Returns 1 otherwise.

Example:

    say prime_power(15)     #=> 1
    say prime_power(43**5)  #=> 5

=cut

=head2 prime_power_count

    prime_power_count(n)
    prime_power_count(a,b)

Returns the number of prime powers <= n, or in the range C<a..b>.

Example:

    say prime_power_count(1e15)     # number of prime powers <= 10^15
    say prime_power_count(1e6, 1e8) # number of prime powers in [10^6, 10^8]

=cut

=head2 prime_power_divisors

    n.prime_power_divisors()

Returns the prime power divisors of n.

Example:

    say prime_power_divisors(5040)  #=> [2, 3, 4, 5, 7, 8, 9, 16]

Equivalent with:

    n.divisors.grep{.is_prime_power}

=cut

=head2 prime_power_sigma

    n.prime_power_sigma(k=1)

Sum of the prime power divisors of n, each divisor raised to power k.

Example:

    say prime_power_sigma(10!)      #=> 667
    say prime_power_sigma(10!, 2)   #=> 95459

Equivalent with:

    n.prime_power_divisors.sum {|d| d**k }

=cut

=head2 prime_power_udivisors

    prime_power_udivisors(n)

Returns the unitary prime power divisors of n.

Example:

    say prime_power_udivisors(10!)  #=> [7, 25, 81, 256]

The product of the unitary prime power divisors of a number, is the number itself.

This method is equivalent with:

    n.factor_map {|p,k| p**k }.sort

Aliases: I<prime_power_unitary_divisors>, I<unitary_prime_power_divisors>

=cut

=head2 prime_power_usigma

    Number.prime_power_usigma(k=1)

Sum of the unitary prime power divisors of n, each divisor raised to power C<k>.

Example:

    say prime_power_usigma(10!)     #=> 369     (== 7 + 25 + 81 + 256)
    say prime_power_usigma(10!, 2)  #=> 72771   (== 7^2 + 25^2 + 81^2 + 256^2)

Equivalent with:

    n.factor_map {|p,e| p**(e * k) }.sum

=cut

=head2 prime_root

    n.prime_root

Returns the prime C<p> if n can be expressed as C<n = p^k> for some prime number C<p> and integer C<k>. Returns C<n> otherwise.

Example:

    say prime_root(15)      #=> 15
    say prime_root(43**5)   #=> 43

=cut

=head2 primes

    primes(n)
    primes(a,b)

Returns an array with the prime numbers <= n, or in the range C<a..b>.

=cut

=head2 prime_sigma

    n.prime_sigma(k=1)

Sum of the unique prime divisors of n, each divisor raised to power k.

Example:

    say prime_sigma(100!)       #=> 1060
    say prime_sigma(100!, 2)    #=> 65796

=cut

=head2 prime_sum

    prime_sum(n)
    prime_sum(a,b)

Sum of the prime numbers <= n, or in the range C<a..b>.

Aliases: I<primes_sum>, I<sum_primes>

=cut

=head2 prime_udivisors

    n.prime_udivisors

Returns the unique unitary prime factors of C<n>.

Aliases: I<prime_unitary_divisors>, I<unitary_prime_divisors>

=cut

=head2 prime_upper

    n.prime_upper

Upper bound for the n-th prime.

Aliases: I<nth_prime_upper>

=cut

=head2 prime_usigma

    n.prime_usigma(k=1)

Sum of the unique unitary prime divisors of n, each divisor raised to power k.

Example:

    say prime_usigma(100!)      #=> 732
    say prime_usigma(100!, 2)   #=> 55330

=cut

=head2 prime_usigma0

    prime_usigma0(n)

The number of unique unitary prime divisors of n.

=cut

=head2 primitive_part

    n.primitive_part(f)

Returns the primitive part of C<f(n)>, for n > 0, such that:

    a(n) = primitive part of f(n)
    f(n) = Product_{d|n} a(d)

Example:

    func f(n) { n.fib }
    func a(n) { n.primitive_part(f) }

    say 20.of { a(_) }
    say 20.of { f(_) }
    say 20.of { .divisors.prod {|d| a(d) } }

=cut

=head2 primorial

    Number.primorial()

Returns the

=cut

=head2 primorial_deflation

    n.primorial_deflation()

Primorial deflation of n, satisfying:

    primorial_inflation(primorial_deflation(n)) = n

Defined as:

    primorial_deflation(n) = A319626(n)/A319627(n)

=cut

=head2 primorial_inflation

    n.primorial_inflation()

Primorial inflation of n: fully multiplicative with a(p) = primorial(p).

Defined as:

    primorial_inflation(n) = A108951(n)

The method also accepts C<n> to be a fraction, which is computed as:

    primorial_inflation(a/b) = primorial_inflation(a)/primorial_inflation(b)

=cut

=head2 psi

    n.psi(k=1)
    dedekind_psi(n,k=1)

Dedekind psi function.

Example:

    say 10.of { .dedekind_psi    }   #=> [0, 1, 3, 4, 6, 6, 12, 8, 12, 12]
    say 10.of { .dedekind_psi(2) }   #=> [0, 1, 5, 10, 20, 26, 50, 50, 80, 90]

Aliases: I<dedekind_psi>

=cut

=head2 qnr

    qnr(n)

Returns the least quadratic non-residue of n.

Example:

    say qnr(17676352761153241)      #=> 37
    say qnr(172138573277896681)     #=> 41

Aliases: I<quadratic_nonresidue>

=cut

=head2 qs_factor

    qs_factor(n)

Pomerance's Quadratic Sieve factorization method (SIMPQS).

=cut

=head2 quadratic_formula

    quadratic_formula(a,b,c)

Returns a list of solutions C<(x_1, x_2)> to the quadratic equation: C<a*x^2 + b*x + c = 0>, defined as:

    (-b ± sqrt(b^2 - 4ac)) / (2a)

Example:

    say [quadratic_formula(13, -42, -34)]     #=> [3.9011...., -0.6704...]

=cut

=head2 rad

    rad(n)

Returns the radical of C<n>, which is the largest squarefree divisor of C<n>.

Example:

    say rad(2**5 * 3**9)    #=> 6

Equivalent to:

    n.factor.uniq.prod

=cut

=head2 rad2deg

    rad2deg(x)

Convert radians to degrees.

Example:

    say rad2deg(Num.pi)     #=> 180

=cut

=head2 ramanujan_sum

    ramanujan_sum(n,q)

The Ramanujan sum function, defined as:

    c_q(n) = μ(q/gcd(n,q)) * φ(q) / φ(q/gcd(n,q))

Example:

    say 20.of {|q| ramanujan_sum(2, q) }    #=> OEIS: A086831

=cut

=head2 ramanujan_tau

    ramanujan_tau(n)

Ramanujan's tau function.

=cut

=head2 rand

    rand(n)
    rand(a,b)

Returns a pseudorandom floating-point in the interval C<[0,n)>, or in the interval C<[a,b)>.

=cut

=head2 random_bytes

    n.random_bytes

Returns an array with C<n> random values between 0 and 255.

=cut

=head2 random_maurer_nbit_prime

    random_maurer_nbit_prime(n)

Generate a random n-bit Marurer prime.

Aliases: I<random_nbit_maurer_prime>

=cut

=head2 random_nbit_prime

    random_nbit_prime(n)

Returns a random n-bit prime.

Example:

    say random_nbit_prime(100)      # 100-bit random prime

=cut

=head2 random_nbit_strong_prime

    random_nbit_strong_prime(n)

Generate a random n-bit strong prime number.

Aliases: I<random_strong_nbit_prime>

=cut

=head2 random_ndigit_prime

    random_ndigit_prime(n)

Returns a random n-digit prime.

Example:

    say random_ndigit_prime(100)    # 100-digit random prime

=cut

=head2 random_prime

    random_prime(n)
    random_prime(a,b)

Returns a random prime <= n, or in the range C<a..b>.

Example:

    say random_prime(100)         # random prime in range [2, 100]
    say random_prime(100, 1000)   # random prime in range [100, 1000]

=cut

=head2 random_safe_prime

    random_safe_prime(n)

Returns a random n-bit safe-prime.

Example:

    say random_safe_prime(512)   #=> 512-bit safe prime

Provided by L<Math::Prime::Util::GMP> >= 0.52.

=cut

=head2 random_string

    n.random_string

Returns a string with C<n> random characters (bytes).

=cut

=head2 range

    range(n)
    range(a,b)
    range(a,b,step)

Creates a new RangeNumber object.

Example:

    say range(10)           #=> RangeNum(0, 9, 1)
    say range(1, 10)        #=> RangeNum(1, 10, 1)
    say range(1, 10, 2)     #=> RangeNum(1, 10, 2)

=cut

=head2 rat

    x.rat

Convert C<x> to a rational number. Returns C<NaN> when this conversion is not possible.

Aliases: I<to_r>, I<to_rat>

=cut

=head2 rat_approx

    x.rat_approx

Given a real number C<x>, it returns a very good (sometimes exact) rational approximation to C<n>, computed with continued fractions.

Example:

    say rat_approx(3.14).as_frac        #=> 22/7
    say rat_approx(zeta(-5)).as_frac    #=> -1/252

Returns NaN when C<x> is not a real number.

=cut

=head2 re

    re(z)

Returns the real part of a complex number C<z>.

Aliases: I<real>

=cut

=head2 reals

    reals(z)

Returns a list with the real and imaginary parts of C<z>.

=cut

=head2 remdiv

    n.remdiv(k)

Removes all occurrences of the divisor C<k> from integer C<n>.

Equivalent with:

    n / k**valuation(n,k)

Aliases: I<remove>

=cut

=head2 rising_factorial

    rising_factorial(n,k)

Rising factorial: C<n^(k) = n * (n + 1) * ... * (n + k - 1)>, defined as:

    binomial(n + k - 1, k) * k!

For negative values of C<k>, rising factorial is defined as:

    rising_factorial(n, -k) = 1/rising_factorial(n - k, k)

=cut

=head2 root

    root(n,k)

The k-th root of n, defined as:

    n**(1/k)

=cut

=head2 rotate

    rotate(n, k, b=10)

Rotate the digits of C<n> to the left when C<k> is positive or to the right otherwise, in a given base C<b>, or base 10 when no base is given.

Example:

    say rotate(12345, 2)         #=> 34512
    say rotate(12345, -1)        #=> 51234

=cut

=head2 rough_count

    k.rough_count(n)
    k.rough_count(a,b)

Returns the number of k-rough numbers <= n, or in the range C<a..b>.

Example:

    say 97.rough_count(1e6)    #=> 122005
    say 23.rough_count(1e9)    #=> 171024023

Also works with arbitarily large n:

    say 43.rough_count(1e34)   #=> 1450936704022016442012254601096989

=cut

=head2 rough_part

    n.rough_part(k)

It returns the k-rough part of n that contains all the prime factors p|n such that p >= k.

Example:

    say 15.of {|n| rough_part(n!, 3) }      #=> OEIS: A049606

=cut

=head2 round

    round(x,k=0)

Rounds C<x> to the k-th decimal place.

A negative argument rounds that many digits after the decimal point, while a positive argument rounds that many digits before the decimal point.

Example:

    say round(1234.567)             #=> 1235
    say round(1234.567, 2)          #=> 1200
    say round(3.123+4.567i, -2)     #=> 3.12+4.57*i

Aliases: I<roundf>

=cut

=head2 run

    run(a,b)

Given two arguments, it returns the second one.

=cut

=head2 sec

    sec(x)

Trigonometric secant function.

=cut

=head2 secant_number

    secant_number(n)

Return the n-th secant/zig number. (OEIS: A000364)

=cut

=head2 sech

    sech(x)

Hyperbolic secant function.

=cut

=head2 seed

    seed(n)

Re-seed the C<rand()> function. The value for C<n> can be any arbitrary large integer.

=cut

=head2 semiprime

    semiprime(n)

Returns the n-th semiprime number.

Example:

    say 10.of {|k| semiprime(10**k) }    #=> OEIS: A114125

Aliases: I<nth_semiprime>

=cut

=head2 semiprime_count

    semiprime_count(n)
    semiprime_count(a,b)

Counts the number of semiprimes <= n (OEIS: A072000), or in the range C<a..b>.

Example:

    say 20.of {|k| semiprime_count(2**k) }

=cut

=head2 semiprimes

    semiprimes(n)
    semiprimes(a,b)

Returns an array with the semiprimes <= n, or in the range C<a..b>.

Example:

    say semiprimes(100)             # semiprimes <= 100
    say semiprimes(50, 100)         # semiprimes in the range [50, 100]

=cut

=head2 setbit

    setbit(n,k)

Set the k-th bit of integer C<n> to 1.

Example:

    say setbit(0b1000, 0).as_bin    #=> 1001
    say setbit(0b1000, 2).as_bin    #=> 1100

=cut

=head2 sgn

    sgn(x)

Returns the sign of x.

Defined as:

    x / abs(x)

Aliases: I<sign>

=cut

=head2 sigma0

    tau(n)
    sigma0(n)

Returns the number of positive divisors of n.

=cut

=head2 sin

    sin(x)

Trigonometric sine function.

=cut

=head2 sin_cos

    sin_cos(x)

Returns a list with the values C<(sin(x), cos(x))>.

=cut

=head2 sinh

    sinh(x)

Hyperbolic sine function.

=cut

=head2 smooth_count

    k.smooth_count(n)
    k.smooth_count(a,b)

Returns the number of k-smooth numbers <= n, or in the range C<a..b>.

Example:

    say 30.of {|n| 13.smooth_count(10**n) }   #=> OEIS: A106629

=cut

=head2 smooth_part

    n.smooth_part(k)

It efficiently returns the largest divisor of n that is k-smooth.

Example:

    say 20.of {|n| smooth_part(n!, 3) }     #=> OEIS: A118381

=cut

=head2 solve_pell

    solve_pell(d, k=1)

Find the smallest pair of integers C<(x,y)> that satisfy the Pell equation: C<x^2 - d*y^2 = k>.

Example:

    say [solve_pell(863)]       #=> [18524026608, 630565199]
    say [solve_pell(953, -1)]   #=> [2746864744, 88979677]

Return C<(nil,nil)> when a solution does not exist.

=cut

=head2 sopf

    sopf(n)

Sum of the distinct primes dividing n.

Example:

    say 30.of { .sopf }     #=> OEIS: A008472

=cut

=head2 sopfr

    sopfr(n)

Sum of primes dividing n (with repetition).

Example:

    say 30.of { .sopfr }    #=> OEIS: A001414

=cut

=head2 sqr

    sqr(x)

Returns the square of C<x>. Equivalent with C<x*x>.

Aliases: I<square>

=cut

=head2 sqrt

    sqrt(x)

Returns the square root of C<x>. Equivalent with C<x**(1/2)>.

=cut

=head2 sqrt_cfrac

    n.sqrt_cfrac
    n.sqrt_cfrac(k)

Returns the expansion of the continued fraction for square root of n.

When an additional argument C<k> is specified, it includes only the first C<k> terms from the expansion period.

Example:

    say 28.sqrt_cfrac               #=> [5, 3, 2, 3, 10]
    say 28.sqrt_cfrac(2)            #=> [5, 3, 2]
    say 12345678910.sqrt_cfrac(10)  #=> [111111, 9, 26, 1, 2, 3, 1, 1, 2, 8, 1]

=cut

=head2 sqrt_cfrac_period

    n.sqrt_cfrac_period

Returns the period of the expansion of the continued fraction for square root of n.

Example:

    say 28.sqrt_cfrac_period        #=> [3, 2, 3, 10]

=cut

=head2 sqrt_cfrac_period_each

    n.sqrt_cfrac_period_each { ... }
    n.sqrt_cfrac_period_each({ ... }, max_iter)

Iterate over the period of the expansion of the continued fraction for square root of n.

Example:

    28.sqrt_cfrac_period_each {|r| say r }

=cut

=head2 sqrt_cfrac_period_len

    n.sqrt_cfrac_period_len

Returns the length of the period of continued fraction for square root of n. (OEIS: A003285)

=cut

=head2 sqrtmod

    sqrtmod(a,m)

Modular square root of C<a>, returning a solution C<r>, such that C<r^2 = a (mod m)>.

Example:

    say sqrtmod(544, 800)         #=> 512
    say sqrtmod(436, 1752)        #=> 1010

=cut

=head2 sqrtmod_all

    sqrtmod_all(a, n)

Given two integers C<a> and C<n>, return a sorted array of all modular square roots of C<a> mod C<|n|>. If no square root exists, an empty array is returned.

Example:

    say sqrtmod_all(4095, 8469)     #=> [1110, 1713, 3933, 4536, 6756, 7359]

=cut

=head2 square_divisors

    square_divisors(n)

Returns the square divisors of n.

Example:

    say 5040.square_divisors    #=> [1, 4, 9, 16, 36, 144]

Equivalent with:

    n.divisors.grep { .is_square }

=cut

=head2 squarefree

    squarefree(n)
    squarefree(a,b)

Returns an array with the squarefree numbers <= n, or in the range C<a..b>.

Example:

    # Squarefree numbers in the interval [100, 200]
    say squarefree(100, 200)

    # Squarefree numbers <= 100
    say squarefree(100)

=cut

=head2 squarefree_almost_primes

    k.squarefree_almost_primes(n)
    k.squarefree_almost_primes(a,b)

Returns an array with the squarefree k-almost primes <= n, or in the range C<a..b>.

Example:

    say 2.squarefree_almost_primes(100)      #=> squarefree semiprimes <= 100
    say 2.squarefree_almost_primes(50, 100)  #=> squarefree semiprimes in the range 50..100

    say 3.squarefree_almost_primes(100)      #=> squarefree 3-almost primes <= 100
    say 3.squarefree_almost_primes(50, 100)  #=> squarefree 3-almost primes in the range 50..100

=cut

=head2 square_free_count

    square_free_count(n)
    square_free_count(a,b)

Returns the number of squarefree integers <= n, or in the range C<a..b>.

Example:

    # The number of square-free integers between 10^5 and 10^6
    say square_free_count(10**5, 10**6) # 547132

    # The number of square-free integers <= 2^40
    say square_free_count(2**40)        # 668422917419

Aliases: I<squarefree_count>

=cut

=head2 squarefree_divisors

    squarefree_divisors(n)

Returns the square-free divisors of n.

Example:

    say squarefree_divisors(120)  #=> [1, 2, 3, 5, 6, 10, 15, 30]

=cut

=head2 squarefree_pi_k

    k.squarefree_almost_prime_count(n)
    k.squarefree_almost_prime_count(a,b)

Returns the number of squarefree k-almost primes <= n, or in the range C<a..b>.

Example:

    say 1.squarefree_almost_prime_count(1000)          # count of primes <= 1000
    say 2.squarefree_almost_prime_count(50, 1000)      # count of squarefree semiprimes in range 50..1000

Aliases: I<squarefree_almost_primepi>, I<squarefree_almost_prime_count>

=cut

=head2 squarefree_sigma

    squarefree_sigma(n,k=1)

Sum of the square free divisors of C<n>, each divisor raised to power C<k>.

Example:

    say squarefree_sigma(5040)      #=> 576
    say squarefree_sigma(5040, 2)   #=> 65000

=cut

=head2 squarefree_sum

    squarefree_sum(n)
    squarefree_sum(a,b)

Returns the sum of squarefree numbers <= n or in the range C<a..b>. (OEIS: A066779)

Example:

    say squarefree_sum(1e12)    #=> 303963551353876732927386
    say squarefree_sum(1e13)    #=> 30396355090144154315002969

=cut

=head2 squarefree_udivisors

    squarefree_udivisors(n)

Returns the unitary squarefree divisors of C<n>.

Example:

    say squarefree_udivisors(5040)   #=> [1, 5, 7, 35]

Aliases: I<squarefree_unitary_divisors>, I<unitary_squarefree_divisors>

=cut

=head2 squarefree_usigma

    n.squarefree_usigma(k=1)

Sum of the unitary squarefree divisors of n, each divisor raised to power k.

Example:

    say 5040.squarefree_usigma     #=> 48
    say 5040.squarefree_usigma(2)  #=> 1300

=cut

=head2 squarefree_usigma0

    squarefree_usigma0(n)

Returns the number of unitary squarefree divisors of n.

Example:

    say squarefree_usigma0(5040)    #=> 4

=cut

=head2 square_sigma

    n.square_sigma(k=1)

Sum of the square divisors of C<n>, each divisors raised to the power C<k>.

Example:

    say 5040.square_sigma       #=> 210
    say 5040.square_sigma(2)    #=> 22386

Equivalent with:

    n.square_divisors.sum {|d| d**k }

=cut

=head2 square_sigma0

    square_sigma0(n)

Returns the number of square divisors of C<n>.

=cut

=head2 squares_r

    squares_r(n, k=2)

Sum of squares function C<r_k(n)>: returns the number of ways or representing C<n> as a sum of C<k> squares.

Example:

    say 30.of { .squares_r(2) }     # OEIS: A004018
    say 30.of { .squares_r(3) }     # OEIS: A005875
    say 30.of { .squares_r(4) }     # OEIS: A000118

Aliases: I<sum_of_squares_count>

=cut

=head2 square_udivisors

    square_udivisors(n)

Returns the unitary square divisors of C<n>, such that each divisors C<d> is a square and C<gcd(n/d, d) = 1>.

Example:

    say 5040.square_udivisors   #=> [1, 9, 16, 144]

Aliases: I<square_unitary_divisors>, I<unitary_square_divisors>

=cut

=head2 square_usigma

    n.square_usigma(k=1)

Sum of the unitary square divisors of C<n>, each divisors raised to power C<k>.

Example:

    say square_usigma(5040)     #=> 170
    say square_usigma(5040, 2)  #=> 21074

Equivalent with:

    n.square_udivisors.sum {|d| d**k }

=cut

=head2 square_usigma0

    square_usigma0(n)

Returns the number of unitary square divisors of C<n>.

=cut

=head2 squfof_factor

    n.squfof_factor(reps)

Shanks' SQUFOF method.

=cut

=head2 stirling

    stirling(n,k)

Stirling numbers of the first kind.

Aliases: I<stirling1>

=cut

=head2 stirling2

    stirling2(n,k)

Stirling numbers of the second kind.

=cut

=head2 stirling3

    stirling3(n,k)

Stirling numbers of the third kind (also known as Lah numbers).

=cut

=head2 subfactorial

    subfactorial(n,k=0)

Number of derangements of a set with n elements having exactly k fixed points. Also known as the rencontres numbers.

Example:

    say 20.of { .subfactorial }        #=> OEIS: A000166
    say 20.of { .subfactorial(2) }     #=> OEIS: A000387

=cut

=head2 submod

    submod(a,b,m)

Modular intger subtraction: C<(a-b)%m>.

Example:

    say submod(43, 97, 127)     #=> (43-97)%127

=cut

=head2 subsets

    n.subsets
    n.subsets(k)
    n.subsets { ... }
    n.subsets(k, { ... })

Returns an array with the subsets of the integers in the range C<0..n-1>, or iterates over the k-subsets when a block is given.

Example:

    say 5.subsets           # all k-subsets of 0..4
    say 5.subsets(2)        # all 2-subsets of 0..4

    5.subsets {|*a| say a }         # iterate over all k-subsets of 0..4
    5.subsets(2, {|*a| say a })     # iterate over all 2-subsets of 0..4

=cut

=head2 sum_remainders

    sum_remainders(n, v)

Returns the following sum: C<Sum_{k=1..n} v % k>, computed in C<O(sqrt(v))> steps.

Example:

    say 20.of {|n| sum_remainders(n, n) }        #=> OEIS: A004125
    say 20.of {|n| sum_remainders(n, n.prime) }  #=> OEIS: A099726

=cut

=head2 superfactorial

    superfactorial(n)

Returns the product of first n factorials.

Example:

    say 10.of { .superfactorial }   #=> OEIS: A000178

=cut

=head2 superprimorial

    superprimorial(n)

Returns the product of first n primorials.

Example:

    say 10.of { .superprimorial }   #=> OEIS: A006939

=cut

=head2 tan

    tan(x)

Trigonometric tangent function.

=cut

=head2 tangent_number

    tangent_number(n)

Returns the n-th tangent/zag number. (OEIS: A000182).

=cut

=head2 tanh

    tanh(x)

Hyperbolic tangent function.

=cut

=head2 times

    n.times {|k| ... }

Call a given block of code C<n> times with k = 0,1,2,...,n-1.

Example:

    5.times {|k| say "k = #{k}" }

=cut

=head2 to_f

    x.to_f

Convert C<x> to a floating-point value number.

Aliases: I<float>, I<to_float>

=cut

=head2 to_n

    x.to_n

Fixed-point function: returns C<x>.

Aliases: I<to_num>

=cut

=head2 to_poly

    x.to_poly

Converts C<x> to a Polynomial object.

=cut

=head2 to_s

    x.to_s

Converts C<x> to a String object.

Aliases: I<to_str>

=cut

=head2 totient

    phi(n)
    totient(n)

Euler's totient function.

Aliases: I<euler_phi>, I<eulerphi>, I<euler_totient>

=cut

=head2 trial_factor

    n.trial_factor(limit)

Trial division.

=cut

=head2 tuples

    n.tuples(k)
    n.tuples(k, { ... })

Returns an array with the k-tuples of the integers in the range C<0..n-1>, or iterates over the k-tuples when a block is given.

Example:

    5.tuples(2, {|*a| say a })

Aliases: I<variations>

=cut

=head2 tuples_with_repetition

    n.tuples_with_repetition(k)
    n.tuples_with_repetition(k, { ... })

Returns an array with the k-tuples with repetition of the integers in the range C<0..n-1>, or iterates over the k-tuples with repetition when a block is given.

Example:

    5.tuples_with_repetition(2, {|*a| say a })

Aliases: I<variations_with_repetition>

=cut

=head2 udivisors

    udivisors(n)

Returns the unitary divisors of C<n>.

These are the divisors C<d> of C<n> such that C<gcd(n/d, d) = 1>.

Example:

    say 120.udivisors   #=> [1, 3, 5, 8, 15, 24, 40, 120]

Aliases: I<unitary_divisors>

=cut

=head2 uphi

    uphi(n)

The unitary totient function. (OEIS: A047994)

=cut

=head2 usigma

    n.usigma(k=1)

Returns the sum of the unitary divisors of C<n>, each divisor raised to the power C<k>.

=cut

=head2 usigma0

    usigma0(n)

Returns the number of unitary divisor of n.

Equivalently, the number of squarefree divisors of n.

Example:

    say usigma0(5040)   #=> 16

Aliases: I<squarefree_sigma0>

=cut

=head2 valuation

    n.valuation(k)

Returns the number of times n is divisible by k.

Example:

    say valuation(2**32, 4)     # prints: 16

=cut

=head2 zero

    Num.zero

Returns the number 0.

=cut

=head2 znorder

    znorder(a,m)

The smallest positive integer k such that C<powmod(a, k, m) == 1>.

Return C<NaN> if C<a> is not coprime to C<m>.

=cut

=head2 znprimroot

    znprimroot(n)

The smallest primitive root of (Z/nZ)^*.

=cut
