NAME

Sidef::Types::Number::Number - Arbitrary-precision numeric type providing extensive mathematical and number-theoretic functionality.

DESCRIPTION

The 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.

The following class-variables can be changed during runtime:

Num!PREC            = 192        # precision for floating-point numbers
Num!ROUND           = 0          # rounding mode for floating-point numbers
Num!VERBOSE         = false      # true to enable verbose/debug mode
Num!USE_YAFU        = false      # true to use YAFU for factoring large integers
Num!USE_PFGW        = false      # true to use PFGW64 as a primality pretest for large enough n
Num!USE_PARI_GP     = false      # true to use PARI/GP in several methods
Num!USE_FACTORDB    = false      # true to use factordb.com for factoring large integers
Num!USE_PRIMESUM    = false      # true to use Kim Walisch's primesum in prime_sum(n)
Num!USE_PRIMECOUNT  = false      # true to use Kim Walisch's primecount in prime_count(n)
Num!USE_CONJECTURES = false      # true to use conjectured methods for better performance
Num!SPECIAL_FACTORS = true       # true to try to find factors of special form in factor(n)

The supported rounding modes for floating-point numbers, are:

Num!ROUND = 0   # Round to nearest.
Num!ROUND = 1   # Round towards zero.
Num!ROUND = 2   # Round towards +infinity.
Num!ROUND = 3   # Round towards -infinity.
Num!ROUND = 4   # Round away from zero. (with MPFR >= 3.0.0)
Num!ROUND = 5   # Faithful rounding. (with MPFR >= 4.0.0)

The values can also be modified only in a local scope, by using the local keyword:

func f(n) {
    local Num!PREC = 1024
    # do some work
}

In the above example, the f(n) function will use 1024 bits of precision for floating-point numbers, while outside the function, the default precision will be used.

NOTE: the local scope extends to any other functions or methods that the are called from the local scope.

SYNOPSIS

# Basic number creation
var a = Num("3.14159")
var b = Number("1010", 2)      # Binary: 10
var c = 42                     # Integer literal

# Arithmetic operations
say (2**100 + 1)               # Large integer arithmetic
say as_frac(1/3 + 1/6)         # Rational arithmetic: 1/2

# Number-theoretic functions
say 100.prime                  # 100th prime: 541
say factor(2**64 + 1)          # Prime factorization
say euler_phi(100)             # Euler's totient: 40

# Floating-point with custom precision
local Num!PREC = 1024         # set precision to 1024 bits
say sqrt(2)                   # High-precision square root

INHERITS

Inherits methods from Sidef::Object::Object.

OPERATORS & COMPARISONS

!=

a != b
ne(a, b)

Returns true if a and b are different.

say (42 != 13)        #=> true
say (5 != 5)          #=> false

Aliases: ≠, ne

*

a * b
mul(a, b)

Multiplies a and b.

say (6 * 7)           #=> 42
say (2/3 * 3/4)       #=> 1/2

Aliases: mul

**

a ** b
pow(a, b)

Raises a to the power b.

say 2**10           #=> 1024
say 5**3            #=> 125
say pow(2, 100)     #=> 1267650600228229401496703205376

Aliases: pow

+

a + b
add(a, b)

Adds two numbers a and b.

say (42 + 13)       #=> 55
say (1/3 + 1/6)     #=> 1/2
say add(2.5, 3.7)   #=> 6.2

Aliases: add

++

n++
++n
n.inc

Increments n by 1 and returns the result.

var x = 5
say ++x             #=> 6
say x++             #=> 6 (returns 6, x is now 7)
say x.inc           #=> 8

The inc method-call does not modify the variable in-place.

Aliases: inc

-

a - b
sub(a, b)

Subtracts b from a.

say (50 - 8)          #=> 42
say (1/2 - 1/3)       #=> 1/6

Aliases: sub

--

n--
--n
n.dec

Decrements n by 1 and returns the result.

var x = 5
say --x             #=> 4
say x--             #=> 4 (returns 4, x is now 3)
say x.dec           #=> 2

The dec method-call does not modify the variable in-place.

Aliases: dec

..

a .. b

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

Equivalent with:

RangeNum(a, b)

Aliases: to, upto

..^

a ..^ b

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

Equivalent with:

RangeNum(a, b-1)

Aliases: xto, xupto

/

a / b
div(a, b)

Divides a by b.

say (10 / 3)          #=> 10/3
say (15 / 5)          #=> 3
say (7.5 / 2.5)       #=> 3

Aliases: ÷, div

//

a // b
idiv(a, b)

Integer floor division. Equivalent to floor(a/b) for integers.

say (17 // 5)         #=> 3
say (-17 // 5)        #=> -4
say idiv(100, 7)      #=> 14

Aliases: divint, fld, idiv, idiv_floor

<

a < b
lt(a, b)

Returns true if a is less than b.

say (3 < 5 )          #=> true
say (10 < 7)          #=> false

Aliases: lt

<=

a <= b
le(a, b)

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

say (3 <= 5)          #=> true
say (5 <= 5)          #=> true

Aliases: ≤, le

<=>

a <=> b
cmp(a, b)

Three-way comparison: returns -1 if a < b, 0 if a == b, +1 if a > b.

say (5 <=> 10)        #=> -1
say (10 <=> 10)       #=> 0
say (15 <=> 10)       #=> 1

Aliases: cmp

<~>

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

Approximate comparison of a and b.

Equivalent with:

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

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

say (0.1f + 0.2f <~> 0.3)   #=> 0 (equal once rounded)

Aliases: approx_cmp

==

a == b
eq(a, b)

Returns true if a and b are equal.

say (42 == 42)        #=> true
say (1/2 == 0.5)      #=> true

Aliases: eq

=~=

a =~= b

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

Aliases: ≅, approx_eq

>

a > b
gt(a, b)

Returns true if a is greater than b.

say (10 > 7)          #=> true
say (3 > 5)           #=> false

Aliases: gt

>=

a >= b
ge(a, b)

Returns true if a is greater than or equal to b.

say (10 >= 7)         #=> true
say (7 >= 7)          #=> true

Aliases: ≥, ge

^..

a ^.. b

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

Equivalent with:

RangeNum(a-1, b, -1)

Aliases: xdownto

abs

n.abs

The absolute value of n. Works for integers, rationals, floating-points and complex numbers.

say abs(-42)         #=> 42
say abs(-1/3)        #=> 1/3

acmp

acmp(a,b)

Absolute comparison of a and b, defined as:

abs(a) <=> abs(b)

say acmp(-10, 3)      #=> 1
say acmp(-3, 10)      #=> -1

approx_ge

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

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

Equivalent with:

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

approx_gt

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

True if a is approximately greater than b.

Equivalent with:

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

approx_le

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

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

Equivalent with:

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

approx_lt

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

True if a is approximately less than b.

Equivalent with:

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

approx_ne

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

True if a is approximately different than b.

Equivalent with:

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

divides

a.divides(b)
a `divides` b

Returns true if a divides b.

say 3.divides(15)   #=> true

divmod

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

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

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

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

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

geometric_sum

geometric_sum(n,r)

Geometric sum: 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

hypot

hypot(x,y)

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

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

Also defined for complex numbers.

say hypot(3, 4)       #=> 5

iabs

iabs(n)

Integer absolute value, by first truncating n to an integer.

say iabs(-4.7)        #=> 4

Aliases: absint

iadd

iadd(a,b)

Integer addition, by first truncating a and b to integers: int(a) + int(b).

say iadd(2.9, 3.1)    #=> 5

Aliases: addint

icmp

x.icmp(y)

Integer comparison, by first truncating x and y to integers.

say 5.5.icmp(5.1)     #=> 0 (both truncate to 5)

Aliases: cmpint

imod

imod(a,m)

Integer remainder of a when divided by m, by first truncating a and m to integers: int(a) % int(m).

say imod(17.9, 5)     #=> 2

Aliases: modint

imul

imul(a,b)

Integer multiplication, by first truncating a and b to integers: int(a) * int(b).

say imul(2.9, 3.9)    #=> 6

Aliases: mulint

ineg

ineg(n)

Integer negation, by first truncating n to an integer.

say ineg(4.7)         #=> -4

Aliases: negint

inv

inv(n)

Multiplicative inverse of n: 1/n.

say inv(4)            #=> 1/4
say inv(1/3)          #=> 3

is_divisible

n %% k
n.is_divisible(k)

Returns true if n is divisible by k. False otherwise.

say 15.is_divisible(5)   #=> true
say (17 %% 5)            #=> false

Aliases: %%, is_divisible_by, is_div

isub

isub(a,b)

Integer subtraction, by first truncating a and b to integers: int(a) - int(b).

say isub(10.9, 3.1)   #=> 7

Aliases: subint

max

max(...)

Returns the maximum value from a list of numbers.

say max(3, 7, 2, 9, 1)    #=> 9

Aliases: vecmax

min

min(...)

Returns the smallest value from a list of numbers.

say min(3, 7, 2, 9, 1)    #=> 1

Aliases: vecmin

mod

n % k
mod(n, k)

Returns the remainder of n divided by k.

say (17 % 5)          #=> 2
say (-17 % 5)         #=> 3
say 100.mod(7)        #=> 2

Aliases: %

neg

x.neg

Negates the sign of x (equivalent with: -x).

say 5.neg             #=> -5
say (-5 -> neg)       #=> 5

prod

prod(...)

Returns the product of a given list of numbers.

say prod(1, 2, 3, 4, 5)   #=> 120

Aliases: Π, vecprod

remdiv

n.remdiv(k)

Removes all occurrences of the divisor k from integer n.

Equivalent with:

n / k**valuation(n,k)

Aliases: remove

sgn

sgn(x)

Returns the sign of x, as -1, 0 or 1.

Defined as:

x / abs(x)      # for x != 0

say sgn(-42)          #=> -1
say sgn(0)            #=> 0
say sgn(1/3)          #=> 1

Aliases: sign

smod

x.smod(m)

Returns the residual mod m such that it is within half of the modulus.

say smod(1, 6)      #=> 1
say smod(4, 6)      #=> -2

sum

sum(...)

Returns the sum of a given list of numbers.

say sum(1, 2, 3, 4, 5)    #=> 15

Aliases: Σ, vecsum

sum_remainders

sum_remainders(n, v)

Returns the following sum: Sum_{k=1..n} v % k, computed in O(sqrt(v)) steps.

say 20.of {|n| sum_remainders(n, n) }        #=> OEIS: A004125
say 20.of {|n| sum_remainders(n, n.prime) }  #=> OEIS: A099726

Negative values of v are also supported.

BITWISE OPERATIONS

&

a & b
and(a, b)

Performs bitwise AND operation on integers a and b.

say (0b1100 & 0b1010)     #=> 8 (0b1000)
say 12.and(10)            #=> 8

Aliases: and

<<

a << b
shift_left(a, b)

Left shifts a by b bits. Equivalent to a × 2^b for integers.

say (5 << 3)            #=> 40 (5 × 8)
say shift_left(1, 10)   #=> 1024

Aliases: lsft, shift_left

>>

a >> b
shift_right(a, b)

Right shifts a by b bits. Equivalent to floor(a / 2^b) for integers.

say (40 >> 3)              #=> 5 (40 / 8)
say shift_right(1024, 10)  #=> 1

Aliases: rsft, shift_right

^

a ^ b
xor(a, b)

Performs bitwise XOR operation on integers a and b.

say (0b1100 ^ 0b1010)   #=> 6 (0b0110)
say 12.xor(10)          #=> 6

Aliases: xor

bit

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

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

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

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

Aliases: getbit, testbit

bit_scan0

n.bit_scan0
n.bit_scan0(k)

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

When k is omitted, k=0 is assumed.

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

say 0b1011.bit_scan0     #=> 2

bit_scan1

n.bit_scan1
n.bit_scan1(k)

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

When k is omitted, k=0 is assumed.

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

say 0b1100.bit_scan1     #=> 2

bits

n.bits

Returns the binary digits of n.

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

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

Equivalent with:

n.digits(2).flip

clearbit

n.clearbit(k)

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

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

flipbit

flipbit(n,k)

Flip the value of the k-th bit of integer n.

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

hamdist

hamdist(a,b)

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

say hamdist(0b1010, 0b0011)   #=> 2

lsb

lsb(n)

Returns the index of the least significant bit of n that is nonzero.

say 0b110010101111000000.lsb    # 6

msb

msb(n)

Returns the index of the most significant bit of n.

say 0b110010101111000000.msb    # 17

popcount

n.popcount

Number of 1's in binary representation of n.

This value is also known as the Hamming weight value.

say 0b1011.popcount   #=> 3

Aliases: hammingweight

rotate

rotate(n, k, b=10)

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

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

setbit

setbit(n,k)

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

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

|

a | b
or(a, b)

Performs bitwise OR operation on integers a and b.

say (0b1100 | 0b1010)   #=> 14 (0b1110)
say 12.or(10)           #=> 14

Aliases: or

~

~a
not(a)

Returns the bitwise complement of a.

say as_bin(~0b1010)      #=> -1011
say 10.not               #=> -11

Aliases: not

TYPE CONVERSION & FORMATTING

as_bin

n.as_bin

Returns a String with the binary representation of n.

say 42.as_bin     # "101011"

as_dec

n.as_dec
n.as_dec(k)

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

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

Aliases: as_float

as_frac

n.as_frac
n.as_frac(base)

String-representation of n as fraction.

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

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

as_hex

n.as_hex

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

say 42.as_hex       # "2a"

Returns nil when n cannot be converted to an integer.

as_int

n.as_int
n.as_int(base)

Returns a String containing the integer representation of n in the given base.

The base must be between 2 and 62 and defaults to base 10. The sign is preserved for negative values.

This method first converts n to an integer. Returns nil when that conversion is not possible.

say 255.as_int         #=> 255
say 255.as_int(16)     #=> ff
say (-255).as_int(16)  #=> -ff

as_oct

n.as_oct

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

say 42.as_oct   # 52

Returns nil when n cannot be converted to an integer.

as_rat

n.as_rat
n.as_rat(base)

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

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

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

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

base

n.base(b)

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

say 255.base(16)      #=> ff
say 10.base(2)        #=> 1010

Aliases: in_base

cfrac

n.cfrac
n.cfrac(k)

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

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 k is omitted, it uses the default floating-point precision to deduce k.

Aliases: as_cfrac

chr

n.chr

Convert the integer n into a character.

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

commify

n.commify

Returns a string with thousands separators inserted in the decimal representation of n.

The sign and any fractional part are preserved.

Example:

say 1000.commify          #=> 1,000
say 1e10.commify          #=> 10,000,000,000
say (-12345.678).commify  #=> -12,345.678

convergents

n.convergents(k)

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

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

de

r.de
r.denominator

Returns the denominator for a rational number r.

say denominator(43/97)      #=> 97

Aliases: denominator

digit

n.digit(k, b=10)

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

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

digital_root

digital_root(n, b=10)

Returns the digital root of n, with respect to base b. (OEIS A010888)

say 30.of { .digital_root }               #=> OEIS: A010888
say 30.of { .digital_root(.isqrt+1) }     #=> OEIS: A122197

Also known as "repeated digital sum".

Both n and b can be arbitrary large, as long as b > 1.

digits

n.digits(b=10)

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

Negative numbers are treated by their absolute value. If b <= 1, it returns an empty array.

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

The reverse operation is:

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

digits2num

b.digits2num(digits)

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

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

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

Aliases: from_digits

digits_sum

sumdigits(n, b=10)

Sum of base b digits of n.

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

This is equivalent to:

n.digits(b).sum

Aliases: sum_digits, sumdigits

dump

n.dump

Returns a stringification version of n.

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

egypt_greedy

egypt_greedy(p/q)

Greedy algorithm for Egyptian fraction expansion (also called the Fibonacci-Sylvester algorithm): at each step, extract the largest unit fraction less than the target and replace the target with the remainder.

say egypt_greedy(9/10)      #=> [2, 3, 15]
say egypt_greedy(5/121)     #=> [25, 757, 763309, 873960180913, 1527612795642093418846225]

Returns the array of denominators, such that egypt_greedy(p/q).sum{|d| 1/d } equals p/q.

flip

n.flip(base=10)

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

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

Aliases: reverse

int

int(n)
trunc(n)

Truncate n to an integer (rounding towards zero).

say int(4.7)          #=> 4
say int(-4.7)         #=> -4
say int(7/2)          #=> 3

Aliases: to_i, to_int, trunc

len

n.len(b=10)

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

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

Aliases: size, length

lift

n.lift

Returns the self object.

new

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

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

say Number("42")        #=> 42
say Number("101", 2)    #=> 5
say Number("ff", 16)    #=> 255

Aliases: call

nu

r.nu
r.numerator

Returns the numerator of rational number r.

say numerator(43/97)     #=> 43

Aliases: numerator

nude

nude(r)

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

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

numify

n.numify

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

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

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.

rat

x.rat

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

say 0.5.rat           #=> 1/2
say 4.rat             #=> 4

Aliases: to_r, to_rat

rat_approx

x.rat_approx
x.rat_approx(max_den)

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

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

It also accepts and optional argument max_den, returning the best rational approximation to <x> with denominator q <= max_den.

say rat_approx(Num.pi, 1e4).as_frac    #=> 355/113
say rat_approx(Num.pi, 1e9).as_frac    #=> 2549491779/811528438

Returns NaN when x is not a real number.

Aliases: best_rational_approx

to_f

x.to_f

Convert x to a floating-point value number.

say to_f(1/3)         #=> 0.333333333333333333333333333333333333333333333333
say 42.to_f           #=> 42

Aliases: f, float, to_float

to_n

x.to_n

Fixed-point function: returns x.

Aliases: to_num

to_poly

x.to_poly

Converts x to a Polynomial object.

to_s

x.to_s

Converts x to a String object.

Aliases: to_str

SEARCH, RANGES & ITERATION

bsearch

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

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

The last argument is a block which does the comparisons.

This function finds a value k such that f(k) = 0. Returns 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)

bsearch_ge

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

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

The last argument is a block which does the comparisons.

This function finds a value k such that f(k-1) < 0 and f(k) >= 0. Returns 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)

bsearch_le

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

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

The last argument is a block which does the comparisons.

This function finds a value k such that f(k) <= 0 and f(k+1) > 0. Returns 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)

bsearch_max

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

Binary search, returning the largest integer value in the range a..b that satisfies the given comparison function.

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

where:

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

bsearch_min

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

Binary search, returning the smallest integer value in the range a..b that satisfies the given comparison function.

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

where:

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

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 approx_cmp(a,b) method (or the `<~>` operator) is recommended to be used for comparisons.

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: bsearch_inverse

by

n.by { ... }

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

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

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,2,...

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

downto

a.downto(b, step=1)

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

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

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 counting from 0.

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

Also aliased as .st, .nd and .rd:

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

Aliases: rd, st, th

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.

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

range

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

Creates a new RangeNumber object.

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

run

run(a,b)

Given two arguments, it returns the second one.

say run(1, 42)      #=> 42

times

n.times {|k| ... }

Call a given block of code n times with k = 0,1,2,...,n-1.

5.times {|k| say "k = #{k}" }

ROUNDING, ROOTS, POWERS & LOGARITHMS

cbrt

n.cbrt

Cube root of n, as a (complex) floating-point value.

say cbrt(27)          #=> 3
say cbrt(-8)          #=> 1 + 1.73205080756887729352744634150587236694280525381i

ceil

n.ceil

Round n towards positive Infinity.

say ceil(2.1)         #=> 3
say ceil(-2.1)        #=> -2

Aliases: ceiling

cube

cube(x)

Returns the cube of x. Equivalent with x**3.

say cube(5)           #=> 125

exp

exp(x)
exp(b, x)

Exponential function: e^x.

When two arguments are given, it does floating-point exponentiation: b^x.

say exp(1)            #=> 2.718281828459045
say exp(2, 10)        #=> 1024

exp10

exp10(x)

Exponential function: 10^x.

say exp10(3)          #=> 1000

exp2

exp2(x)

Exponential function: 2^x.

say exp2(10)          #=> 1024

expnorm

expnorm(n, b=10)

Returns exp(n) normalized in the range [0,1). The base b can be any value != {0,1}.

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

Equivalent with:

exp(n) / b**k        # for some positive integer value of k

floor

floor(x)

Round x towards -Infinity.

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

icbrt

icbrt(n)

Integer cube root of n.

say icbrt(27)         #=> 3
say icbrt(-8)         #=> -2

Aliases: cbrtint

idiv_ceil

idiv_ceil(a,b)

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

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

ceil(a/b)

say idiv_ceil(7, 2)   #=> 4
say idiv_ceil(-7, 2)  #=> -3

Aliases: cld

idiv_round

idiv_round(a,b)

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

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

floor(a/b + 1/2)

say idiv_round(7, 2)  #=> 4
say idiv_round(5, 2)  #=> 3

Aliases: rdd

idiv_trunc

idiv_trunc(a,b)

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

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

trunc(a/b)

say idiv_trunc(7, 2)   #=> 3
say idiv_trunc(-7, 2)  #=> -3

Aliases: trd

ilog

ilog(n,b)

Integer logarithm of n to base b, satisfying:

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

say ilog(100, 3)      #=> 4    (3**4 = 81 <= 100 < 3**5 = 243)

Aliases: logint

ilog10

ilog10(n)

Integer logarithm of n to base 10, equivalent to:

ilog(n,10)

say ilog10(999)       #=> 2
say ilog10(1000)      #=> 3

ilog2

ilog2(n)

Integer logarithm of n to base 2, equivalent to:

ilog(n,2)

say ilog2(1023)       #=> 9
say ilog2(1024)       #=> 10

ipow

ipow(b,n)

Integer exponentiation: b^n.

say ipow(2, 10)       #=> 1024

Aliases: powint

ipow10

ipow10(n)

Integer exponentiation: 10^n.

say ipow10(6)         #=> 1000000

ipow2

ipow2(n)

Integer exponentiation: 2^n.

say ipow2(10)         #=> 1024

iroot

n.iroot(k)

Integer k-th root of n: the largest integer r such that r**k <= n.

say iroot(100, 2)     #=> 10
say iroot(30, 3)      #=> 3

Aliases: rootint

irootrem

n.irootrem(k)

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

say [30.irootrem(3)]   #=> [3, 3]

Equivalent with:

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

isqrt

n.isqrt

Integer square root of n: the largest integer r such that r**2 <= n.

say isqrt(100)        #=> 10
say isqrt(99)         #=> 9

Aliases: sqrtint

isqrtrem

n.isqrtrem

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

say [99.isqrtrem]     #=> [9, 18]

Equivalent with:

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

lgrt

lgrt(x)

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

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

ln

x.ln

Natural logarithm of x.

say ln(exp(1))        #=> 1
say ln(1)             #=> 0

ln2

Num.ln2

Returns the natural logarithm of 2 constant.

say Num.ln2           #=> 0.6931471805599453

log

log(x)
log(x, b)

Natural logarithm of x to base e, or to a given base b.

say log(exp(1))       #=> 1
say log(100, 10)      #=> 2

log10

x.log10

Logarithm of x to base 10.

say log10(1000)       #=> 3

log2

x.log2

Logarithm of x to base 2.

say log2(1024)        #=> 10

perfect_root

n.perfect_root

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

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

root

root(n,k)

The k-th root of n, defined as n**(1/k).

say root(8, 3)        #=> 2
say root(16, 4)       #=> 2

round

round(x,k=0)

Rounds 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.

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

Aliases: roundf

sqr

sqr(x)

Returns the square of x. Equivalent with x*x.

say sqr(9)            #=> 81

Aliases: square

sqrt

sqrt(x)

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

say sqrt(16)          #=> 4
say sqrt(2)           #=> 1.414213562373095
say sqrt(-1)          #=> 1i

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 k is specified, it includes only the first k terms from the expansion period.

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]

sqrt_cfrac_period

n.sqrt_cfrac_period

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

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

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.

28.sqrt_cfrac_period_each {|r| say r }

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)

say sqrt_cfrac_period_len(28)   #=> 4

sqrtQ

sqrtQ(n)

Returns a Quadratic object, representing sqrt(n).

say sqrtQ(2)**2       #=> Quadratic(2, 0, 2, 0)

TRIGONOMETRIC & HYPERBOLIC FUNCTIONS

acos

n.acos

Inverse cosine of n in radians.

say acos(1)           #=> 0
say acos(0)           #=> 1.570796326794897

acosh

n.acosh

Inverse hyperbolic cosine of n.

say acosh(1)          #=> 0

acot

n.acot

Inverse cotangent of n in radians.

say acot(1)           #=> 0.7853981633974483

acoth

n.acoth

Inverse hyperbolic cotangent of n.

say acoth(2)          #=> 0.5493061443340548

acsc

n.acsc

Inverse cosecant of n in radians.

say acsc(1)           #=> 1.570796326794897

acsch

n.acsch

Inverse hyperbolic cosecant of n.

say acsch(1)          #=> 0.881373587019543

asec

n.asec

Inverse secant of n in radians.

say asec(1)           #=> 0

asech

n.asech

Inverse hyperbolic secant of n.

say asech(1)          #=> 0

asin

n.asin

Inverse sine of n in radians.

say asin(1)           #=> 1.570796326794897

asinh

n.asinh

Inverse hyperbolic sine of n.

say asinh(0)          #=> 0

atan

n.atan

Inverse tangent of n in radians.

say atan(1)           #=> 0.7853981633974483

atan2

atan2(a, b)

Four-quadrant inverse tangent of a and b (the angle, in radians, of the point (b,a)).

say atan2(1, 1)       #=> 0.7853981633974483
say atan2(1, -1)      #=> 2.356194490192345

atanh

n.atanh

Inverse hyperbolic tangent of n.

say atanh(0)          #=> 0

cos

cos(x)

Trigonometric cosine function.

say cos(0)            #=> 1
say cos(Num.pi)       #=> -1

cosh

cosh(x)

Hyperbolic cosine function.

say cosh(0)           #=> 1

cot

cot(x)

Trigonometric cotangent function.

say cot(Num.pi/4)     #=> 1

coth

coth(x)

Hyperbolic cotangent function.

say coth(1)           #=> 1.313035285499331

csc

csc(x)

Trigonometric cosecant function.

say csc(Num.pi/2)     #=> 1

csch

csch(x)

Hyperbolic cosecant function.

say csch(1)           #=> 0.8509181282393215

deg2rad

deg2rad(x)

Convert degrees to radians.

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

rad2deg

rad2deg(x)

Convert radians to degrees.

say rad2deg(Num.pi)     #=> 180

sec

sec(x)

Trigonometric secant function.

say sec(0)            #=> 1

sech

sech(x)

Hyperbolic secant function.

say sech(0)           #=> 1

sin

sin(x)

Trigonometric sine function.

say sin(0)            #=> 0
say sin(Num.pi/2)     #=> 1

sin_cos

sin_cos(x)

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

say [sin_cos(0)]      #=> [0, 1]

sinh

sinh(x)

Hyperbolic sine function.

say sinh(0)           #=> 0

tan

tan(x)

Trigonometric tangent function.

say tan(0)            #=> 0
say tan(Num.pi/4)     #=> 1

tanh

tanh(x)

Hyperbolic tangent function.

say tanh(0)           #=> 0

SPECIAL FUNCTIONS & CONSTANTS

agm

agm(a, b)

Arithmetic-geometric mean of a and b.

say agm(1, sqrt(2))   #=> 1.198140234735592

ai

x.ai
Ai(x)

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

Aliases: airy

bessel_j

bessel_j(x, n)

First order Bessel function: J_n(x).

say bessel_j(1, 0)    #=> 0.7651976865579666

bessel_y

bessel_y(x, n)

Second order Bessel function: Y_n(x).

say bessel_y(1, 0)    #=> 0.08825696421567696

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)

say beta(2, 3)        #=> 1/12

CatalanG

Num.CatalanG

Returns the Catalan constant:

say Num.CatalanG      #=> 0.915965594177219015054603514932384110774149374282

Aliases: C

digamma

digamma(x)

The digamma function, defined as: Gamma'(x)/Gamma(x).

say digamma(1)        #=> -0.5772156649015329

Aliases: Ψ

e

Num.e

Returns the e mathematical constant:

say Num.e             #=> 2.718281828459045

ei

ei(x)

Exponential integral function.

say ei(1)             #=> 1.895117816355937

Aliases: eint

erf

erf(x)

The Gauss error function.

say erf(1)            #=> 0.8427007929497149

erfc

erfc(x)

The complementary error function.

say erfc(1)           #=> 0.1572992070502851

eta

eta(s)

The Dirichlet eta function, defined as:

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

Aliases: η

EulerGamma

Num.EulerGamma

Returns the Euler–Mascheroni constant: 0.5772156649...

Aliases: γ, y

gamma

gamma(n)
Num.gamma

When called as Num.gamma, it returns the Euler-Mascheroni constant:

say Num.gamma   #=> 0.5772156...

When given a real numeric value, it computes the gamma function as a floating-point value. For positive integers, it satisfies the identity gamma(n) = (n-1)!.

say gamma(6)      #=> 120
say gamma(1/2)    #=> 1.7724538...

Aliases: Γ

harm

harmonic(n)
harmonic(n, k)

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

say harmonic(10).as_rat     # 7381/2520

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

say harmonic(10, 2).as_rat   # 55991/2520

Aliases: harmfrac, harmonic, harmonic_number

harmreal

harmreal(n)
harmreal(n, k)

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

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

where γ is the Euler-Mascheroni constant.

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

hclassno

hclassno(n)

Returns the Hurwitz-Kronecker class number.

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

ellcard

ellcard(A, B, p)

Returns the number of points on the elliptic curve:

E : y² = x³ + Ax + B

over the finite field F_p, including the point at infinity.

The modulus p must be a positive prime. The curve is assumed to be nonsingular, that is:

4A³ + 27B² ≠ 0 (mod p).

Example:

say ellcard(2, 3, 5)      # => 7
say ellcard(2, 3, 7)      # => 6
say ellcard(0, 1, 11)     # => 12

By Hasse's theorem, the returned value always satisfies:

|#E(F_p) - (p + 1)| ≤ 2√p

Aliases: elliptic_curve_cardinality

lambert_w

lambert_w(x)

The Lambert-W function. When the value of x is less than -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)

lgamma

lgamma(x)

Natural logarithm of abs(Γ(x)).

Aliases: gamma_abs_log

li

li(x)

Returns the logarithmic integral of x.

say 100.li  # prints: 30.12614158...

li2

li2(x)

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

lnbern

lnbern(n)
lnbern(n,x)

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

When a second argument x is given, it returns the natural logarithm of the n-th Bernoulli polynomial, evaluated at x:

log(B_n(x))     # as a floating-point number

This is usually faster than log(bernreal(n,x)) for large n, but it may be inaccurate in precision in the least significant digits for certain inputs.

Aliases: bern_log, lnbernreal, bernoulli_log

lngamma

lngamma(x)

Natural logarithm of Γ(x).

say lngamma(6)        #=> 4.787491742782046

Aliases: gamma_log

mone

Num.mone

Returns the -1 value.

say Num.mone          #=> -1

nan

Num.nan

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

say Num.nan                #=> NaN
say (Num.nan == Num.nan)   #=> false

inf

Num.inf

Returns the positive Infinity special floating-point value (Inf).

say Num.inf           #=> Inf
say (1 / Num.inf)     #=> 0

ninf

Num.ninf

Returns the negative infinity special value (-Inf).

say Num.ninf          #=> -Inf

one

Num.one

Returns the 1 value.

say Num.one           #=> 1

pi

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

Returns the PI numerical value:

say Num.pi      #=> 3.1415...

When applied on a Number object (as n.pi or 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 a..b:

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

Aliases: π

zero

Num.zero

Returns the number 0.

say Num.zero          #=> 0

zeta

zeta(s)

The Euler-Riemann zeta function. It is currently defined only for real arguments.

say zeta(2)    #=> 1.644934...
say zeta(4)    #=> 1.082323...

Aliases: ζ

COMPLEX NUMBERS

arg

z.arg
z.angle
z.phase

Returns the argument (phase angle) of the complex number in radians.

say arg(1+1i)        #=> 0.7853981633974483

cadd

cadd(a,b,x,y)

Complex arithmetic addition, defined as:

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

Aliases: complex_add

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: complex_div

cinv

cinv(a,b)

Complex arithmetic inversion, defined as:

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

Aliases: complex_inv

cinvmod

cinvmod(a,b,m)

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

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

Aliases: complex_invmod

cis

cis(x)

Euler's formula applied on x, defined as:

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

cmod

cmod(a,b,m)

Complex arithmetic modular operation, defined as:

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

Aliases: complex_mod

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: complex_mul

complex

n.complex
complex(a,b)

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

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

This is equivalent with:

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

complex_cmp

complex_cmp(a,b,x,y)

Complex number comparison, defined as:

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

complex_ipow

complex_ipow(a,b,n)

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

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

conj

conj(x)

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

say conj(3+4i)        #=> 3 - 4i
say 5.conj            #=> 5

cpow

cpow(a,b,n)

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

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

Aliases: complex_pow

cpowmod

cpowmod(a,b,n,m)

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

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

Aliases: complex_powmod

csub

csub(a,b,x,y)

Complex arithmetic subtraction, defined as:

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

Aliases: complex_sub

i

x.i

Multiplies x by the imaginary unit i.

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

im

x.im

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

say im(3+4i)          #=> 4
say 5.im              #=> 0

Aliases: imag, imaginary

is_complex

x.is_complex

Returns true if x is a complex number.

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

is_imag

x.is_imag

Returns true if x is an imaginary number.

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

is_real

x.is_real

Returns true if 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)

norm

norm(x)

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

say norm(3+4i)        #=> 25
say norm(-5)          #=> 25

pair

a.pair(b)

Create a new complex number.

Equivalent with:

Complex(a, b)

Aliases: :, pair

parts

x.parts

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

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

re

re(z)

Returns the real part of a complex number z.

say re(3+4i)          #=> 3
say 5.re              #=> 5

Aliases: real

reals

reals(z)

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

say [reals(3+4i)]     #=> [3, 4]

roots_of_unity

roots_of_unity(n)

Returns an array with the complex n-th roots of 1, as floating-point approximations.

say roots_of_unity(4).map { .round(-5) }   #=> [1, 1i, -1, -1i]

MODULAR ARITHMETIC & CONGRUENCES

addmod

addmod(a, b, m)

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

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

addmulmod

addmulmod(a, b, c, m)

Modular operation: (a + b*c) % m.

say addmulmod(43, 97, 2, 127)     # == (43 + 97*2) % 127

expmod

powmod(b, n, m)

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

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

Aliases: powmod

geometric_summod

geometric_summod(n, r, m)

Geometric sum modulo m: r^0 + r^1 + ... + r^n (mod m), using the following formula:

geometric_summod(n, r, m) = ((powmod(r, n+1, m)-1) * invmod(r-1, m)) % m

invmod

n.invmod(m)

Returns the modular inverse of n modulo m.

say 3.invmod(7)       #=> 5     (since 3*5 % 7 == 1)

is_congruent

n.is_congruent(a,m)

Returns true when n is congruent to a modulo m.

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

is_primitive_root

n.is_primitive_root(m)

Returns true if n is a primitive root modulo m.

say 3.is_primitive_root(7)    #=> true
say 2.is_primitive_root(7)    #=> false

dedekind_sum

dedekind_sum(h, k)

Returns the Dedekind sum s(h, k) as an exact rational number.

The Dedekind sum is defined by

s(h, k) = Σ ((n/k)) * ((hn/k))

where the sum is taken over n = 1, 2, ..., k-1, and the sawtooth function is:

((x)) = x - floor(x) - 1/2    if x is not an integer
      = 0                     otherwise

The second argument k must be a positive integer.

say dedekind_sum(1, 4)    # s(1,4) = 1/8
say dedekind_sum(1, 5)    # s(1,5) = 1/5
say dedekind_sum(2, 3)    # s(2,3) = -1/18

legendre

legendre(a,p)

Returns the Legendre symbol: (a|p).

say legendre(2, 7)    #=> 1
say legendre(3, 7)    #=> -1

jacobi

jacobi(a,n)

Returns the Jacobi symbol: (a|n).

say jacobi(2, 15)     #=> 1
say jacobi(7, 15)     #=> -1

kronecker

kronecker(a,n)

Returns the Kronecker symbol: (a|n).

say kronecker(2, 8)   #=> 0
say kronecker(5, 8)   #=> -1

hilbert

hilbert(a, b, p)

Returns the Hilbert symbol (a, b)_p, whose value is either 1 or -1.

For a prime p, the Hilbert symbol determines whether the quadratic equation

z² = ax² + by²

has a nontrivial solution over the field of p-adic numbers Q_p. It is a fundamental invariant in algebraic number theory, quadratic forms, and local class field theory.

The arguments a and b must be nonzero, and p must be a positive prime. If either a or b is zero, or if p <= 0, the function returns 0.

say hilbert(2, 3, 5)     # => 1
say hilbert(2, 3, 2)     # => -1
say hilbert(-1, -1, 2)   # => -1
say hilbert(-1, -1, 5)   # => 1

The Hilbert symbol is symmetric:

hilbert(a, b, p) == hilbert(b, a, p)

and satisfies the multiplicative identities

(ab, c)_p = (a, c)_p * (b, c)_p
(a, bc)_p = (a, b)_p * (a, c)_p

Aliases: hilbert_symbol

coredisc

coredisc(n)

Returns the fundamental discriminant associated with n. The result is obtained from the squarefree kernel of n, preserving the sign of n:

D = core(n)      if core(n) ≡ 1 (mod 4)
  = 4*core(n)    otherwise

The returned value is the discriminant of the quadratic field Q(sqrt(n)), making it useful in algebraic number theory and quadratic forms.

say coredisc(12)    #=> 12
say coredisc(20)    #=> 5
say coredisc(-5)    #=> -20
say coredisc(-12)   #=> -3

Aliases: fundamental_discriminant

reduce_quadratic_form

reduce_quadratic_form(A, B, C)

Returns the Gauss-reduced form equivalent to the positive definite binary quadratic form

Ax² + Bxy + Cy²

The input form must be positive definite, that is,

A > 0
B² - 4AC < 0

If these conditions are not satisfied, the function returns an empty array.

The returned form (a, b, c) is properly equivalent to the input form and satisfies the classical reduction conditions:

|b| ≤ a ≤ c,

with the additional tie-breaking rules

b ≥ 0    if |b| = a
b ≥ 0    if a = c

Thus, each equivalence class of positive definite binary quadratic forms has a unique reduced representative.

say reduce_quadratic_form(5, 6, 5)    #=> [4, 4, 5]
say reduce_quadratic_form(7, 10, 4)   #=> [1, 0, 3]

Reduced quadratic forms play a central role in the theory of binary quadratic forms, ideal class groups, and quadratic fields.

kronecker_delta

kronecker_delta(a,b)

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

say kronecker_delta(3, 3)   #=> 1
say kronecker_delta(3, 4)   #=> 0

Aliases: δ

muladdmod

muladdmod(a, b, c, m)

Modular operation: (a*b + c) % m.

say muladdmod(43, 97, 5, 127)     # == (43*97 + 5) % 127

muladdmulmod

muladdmulmod(a, b, c, d, m)

Modular operation: (a*b + c*d) % m.

say muladdmulmod(43, 97, 5, 3, 127)   # == (43*97 + 5*3) % 127

mulmod

mulmod(a,b,m)

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

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

mulsubmod

mulsubmod(a, b, c, m)

Modular operation: (a*b - c) % m.

say mulsubmod(43, 97, 5, 127)     # == (43*97 - 5) % 127

mulsubmulmod

mulsubmulmod(a, b, c, d, m)

Modular operation: (a*b - c*d) % m.

say mulsubmulmod(43, 97, 5, 3, 127)   # == (43*97 - 5*3) % 127

polymod

n.polymod(...)

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

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)

qnr

qnr(n)

Returns the least quadratic nonresidue of n.

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

Aliases: quadratic_nonresidue

rootmod

rootmod(a, k, n)

Returns a solution x to the following congruence:

x^k == a (mod n)

Return NaN if a solution does not exist.

say rootmod(2131968,5,10000000000)          #=> 7811648
say rootmod(15,3,1000000000000000000117)    #=> 72574612502199260377

rootmod_all

rootmod_all(a, k, n)

Returns an array with all the solutions x to the following congruence:

x^k == a (mod n)

If no solution exists, an empty array is returned.

say rootmod_all(243, 5, 1000)    #=> [3, 203, 403, 603, 803]

sqrtmod

sqrtmod(a,m)

Modular square root of a, returning a solution r, such that r^2 = a (mod m).

say sqrtmod(544, 800)         #=> 288
say sqrtmod(436, 1752)        #=> 134

sqrtmod_all

sqrtmod_all(a, n)

Returns an array with all the solutions x to the following congruence:

x^2 == a (mod n)

If no solution exists, an empty array is returned.

say sqrtmod_all(4095, 8469)     #=> [1110, 1713, 3933, 4536, 6756, 7359]

submod

submod(a,b,m)

Modular integer subtraction: (a-b) % m.

say submod(43, 97, 127)     #=> (43-97)%127

submulmod

submulmod(a, b, c, m)

Modular operation: (a - b*c) % m.

say submulmod(43, 97, 5, 127)     # == (43 - 97*5) % 127

znlog

znlog(a, g, m)

Returns the integer k that solves the congruence: a = g^k (mod m).

Returns NaN if a solution is not found.

say znlog(9, 2, 11)   #=> 6    (since 2**6 % 11 == 9)

znorder

znorder(a,m)

The smallest positive integer k such that powmod(a, k, m) == 1.

Return NaN if a is not coprime to m.

Aliases: multiplicative_order

znprimroot

znprimroot(n)

The smallest primitive root of (Z/nZ)^*.

COMBINATORICS & PARTITIONS

bell

n.bell

Returns the n-th Bell number.

say 10.bell           #=> 115975

Aliases: bell_number

bellmod

bellmod(n, m)

Returns the n-th Bell number modulo m.

say bellmod(10, 1000)   #=> 975

binomial

binomial(n,k)

Returns the binomial coefficient for choosing k objects from n objects, also called the "choose" function.

For the usual case of nonnegative integers, it counts the number of ways to choose k objects from n objects.

say binomial(5, 2)    #=> 10
say binomial(10, 0)   #=> 1
say binomial(5, 8)    #=> 0

Aliases: nok

binomialmod

binomialmod(n, k, m)

Returns the binomial coefficient binomial(n,k) modulo m.

say binomialmod(1e10, 1e5, 20!)     #=> 286953611424480000

Equivalent with (but much faster):

binomial(n,k) % m

eulerian

eulerian(n, k)

Returns the Eulerian number A(n, k), which counts the number of permutations of {1, 2, ..., n} having exactly k ascents (equivalently, exactly k descents).

say 30.of {|n| eulerian(n, 1) }                     #=> OEIS: A000295
say 30.of {|n| eulerian(n, 2) }                     #=> OEIS: A000460
say 10.of {|n| ^n->map {|k| eulerian(n,k) }.flat }  #=> OEIS: A008292

The first few rows of the Eulerian triangle are:

1
1
1 1
1 4 1
1 11 11 1
1 26 66 26 1

Aliases: eulerian_number

circular_permutations

n.circular_permutations
n.circular_permutations { ... }

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

5.circular_permutations {|*a| say a }

combinations

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

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

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

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 0..n-1, or iterates over the k-combinations with repetition when a block is given.

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

derangements

n.derangements

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

5.derangements {|*a| say a }

Aliases: complete_permutations

double_factorial

n!!
double_factorial(n)

Double-factorial of n. (OEIS A006882)

say 8.double_factorial    #=> 384   (8*6*4*2)
say 7.double_factorial    #=> 105   (7*5*3*1)

Aliases: !!, dfac, dfactorial

factorial

n!
n.factorial

Computes the factorial of n: 1*2*3*...*n.

Defined for nonnegative integers.

say 0.factorial     #=> 1
say 5.factorial     #=> 120
say 10.factorial    #=> 3628800
say 100.factorial   #=> 93326215443944152681699...

Aliases: !, fac

factorial_sum

factorial_sum(n)

Left factorial of n (OEIS A003422), defined as:

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

Example:

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

Aliases: left_factorial

factorial_valuation

factorial_valuation(n, k)

Returns the number of times n! is divisible by k, where k is either prime or composite.

Equivalent with (but more efficient):

valuation(n!, k)

say factorial_valuation(10, 2)   #=> 8    (10! has 2^8 as a factor)

Aliases: factorial_power

factorialmod

factorialmod(n,m)

Returns the factorial of n modulo m.

Equivalent with (but much faster):

factorial(n) % m

say factorialmod(20, 1000)   #=> 0

falling_factorial

falling_factorial(n,k)

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

binomial(n, k) * k!

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

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

When the denominator is zero, NaN is returned.

say falling_factorial(5, 3)    #=> 60    (5*4*3)

farey

farey(n)

Generates the Farey sequence of maximum denominator n.

say farey(5)    #=> [0, 1/5, 1/4, 1/3, 2/5, 1/2, 3/5, 2/3, 3/4, 4/5, 1]

farey_neighbors

farey_neighbors(n, p/q)

Returns the neighbors of p/q in the Farey sequence of max denominator n.

say [farey_neighbors(5, 2/5)]            #=> [1/3, 1/2]
say [farey_neighbors(2**32, 43/97)]      #=> [1903954555/4294967252, 1903954563/4294967270]

fubini

fubini(n)

Returns the n-th Fubini number. (OEIS A000670)

say 10.of { .fubini }    #=> OEIS: A000670

Aliases: fubini_number

fubini_numbers

fubini_numbers(n)

Returns an array containing the Fubini numbers with indices in the range 0..n.

say fubini_numbers(5)   #=> [1, 1, 3, 13, 75, 541]

hyperfactorial

hyperfactorial(n)

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

say hyperfactorial(4)     #=> 27648    (1^1 * 2^2 * 3^3 * 4^4)

hyperfactorial_ln

hyperfactorial_ln(n)

Natural logarithm of hyperfactorial(n), where n is a nonnegative integer.

say hyperfactorial_ln(4)  #=> 10.22730867160378

Aliases: lnhyperfactorial, hyperfactorial_log

kempner

kempner(n)

Returns the smallest positive integer m such that n divides factorial(m).

say 30.of { .kempner }      #=> OEIS: A002034

Aliases: smarandache

lnsuperfactorial

lnsuperfactorial(n)

Natural logarithm of superfactorial(n).

say lnsuperfactorial(4)  #=> 5.662960480135946

Aliases: superfactorial_ln, superfactorial_log

mfac

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

The generalized multi-factorial of n.

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

Aliases: mfactorial, multi_factorial

multinomial

multinomial(...)

The multinomial coefficient, given a list of native integers.

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

multiplicative_partitions

multiplicative_partitions(n)
multiplicative_partitions(n, max_value)
multiplicative_partitions(n, max_value, max_sum)

Returns an array of arrays of integers > 1 whose product is n.

say multiplicative_partitions(6)        #=> [[6], [2, 3]]
say multiplicative_partitions(8)        #=> [[8], [2, 4], [2, 2, 2]]
say multiplicative_partitions(30)       #=> [[30], [2, 15], [3, 10], [5, 6], [2, 3, 5]]

The number of such partitions is given by OEIS A001055.

When the max_value argument is provided, it generates only the partitions with values <= max_value.

say multiplicative_partitions(30, 10)        #=> [[3, 10], [5, 6], [2, 3, 5]]

When the max_sum argument is provided, it generates only the partitions whose sum is <= max_sum.

say multiplicative_partitions(30, nil, 11)   #=> [[5, 6], [2, 3, 5]]

Aliases: mpartitions

multisets

multisets(n, k)
multisets(n, k, max_sum)

Returns an array of arrays of integers of size n, with integers in the range 1..k.

say multisets(3, 2)
#=> [[1, 1, 1], [2, 1, 1], [2, 2, 1], [2, 2, 2]]

When the max_sum argument is provided, it generates only the multisets whose sum is <= max_sum.

say multisets(3, 4, 6)
#=> [[1, 1, 1], [2, 1, 1], [2, 2, 1], [2, 2, 2], [3, 1, 1], [3, 2, 1], [4, 1, 1]]

NOTE: these are also called "combinations with replacements".

necklaces

necklaces(n,k)

Number of n-bead necklaces with k colors when turning over is not allowed.

say 30.of { .necklaces(2) }     #=> OEIS: A000031

Defined as:

(1/n)*Sum_{d|n} phi(n/d)*k^d

necklaces_aperiodic

necklaces_aperiodic(n,k)

Number of n-bead necklaces with beads of k colors when turning over is not allowed and with primitive period n.

say 30.of { .necklaces_aperiodic(2) }   #=> OEIS: A001037

Also number of degree-n irreducible polynomials over GF(k).

Defined as:

(1/n) * Sum_{d|n} mu(n/d) * k^d

frobenius_number

var g = frobenius_number(a, b, ...)

Returns the Frobenius number of a set of positive integers.

The Frobenius number is the largest integer that cannot be expressed as a non-negative integer linear combination of the given values.

The input values must be positive and have greatest common divisor equal to 1. Otherwise, the Frobenius number is undefined and the function returns NaN.

say frobenius_number(4, 7)        # => 17
say frobenius_number(5, 8)        # => 27
say frobenius_number(6, 9, 20)    # => 43
say frobenius_number(1)           # => -1

num2perm

num2perm(n,k)

Given a nonnegative integer n and integer k, return the rank k lexicographic permutation of n elements.

k will be interpreted as mod n!.

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

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

partition_count

n.partition_count
n.partition_count(k)
n.partition_count(set)

Returns the number of integer partitions of n.

With a single argument, this function returns the partition number p(n), that is, the number of ways to write n as a sum of positive integers, where the order of the summands is ignored.

Examples:

say 0.partition_count      # => 1
say 5.partition_count      # => 7
say 10.partition_count     # => 42
say 100.partition_count    # => 190569292

Partitions into exactly k parts

When the second argument is an integer k, the function returns the number of partitions of n into exactly k positive parts.

say 7.partition_count(1)   # => 1
say 7.partition_count(2)   # => 3
say 7.partition_count(3)   # => 4
say 7.partition_count(7)   # => 1
say 7.partition_count(8)   # => 0

Restricted partitions

When the second argument is an array, it is interpreted as the set of allowed part sizes. The function then counts only those partitions whose parts belong to the specified set.

# Parts from {2, 3, 5}
say 10.partition_count([2, 3, 5])  # => 4
# 10 = 5+5
#    = 5+3+2
#    = 3+3+2+2
#    = 2+2+2+2+2

# Odd parts only
say 8.partition_count([1, 3, 5, 7])   # => 6

Aliases: partition_number

partitions

partitions(n)
partitions(n, max_value)

Return an array of arrays with the additive partitions of n, such that their sum equals n.

say partitions(4)       #=> [[4], [3, 1], [2, 2], [2, 1, 1], [1, 1, 1, 1]]

When the max_value argument is provided, it returns only the partitions with value <= max_value.

say partitions(4, 2)    #=> [[2, 2], [2, 1, 1], [1, 1, 1, 1]]

It uses the "Mathematica" ordering of the partitions (graded reverse lexicographic ordering).

permutations

n.permutations

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

5.permutations {|*a| say a }

rising_factorial

rising_factorial(n,k)

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

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

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

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

say rising_factorial(5, 3)   #=> 210    (5*6*7)

stirling

stirling(n,k)

Stirling numbers of the first kind.

say stirling(5, 2)    #=> -50

Aliases: stirling1

stirling2

stirling2(n,k)

Stirling numbers of the second kind.

say stirling2(5, 2)   #=> 15

stirling3

stirling3(n,k)

Stirling numbers of the third kind (also known as Lah numbers).

say stirling3(5, 2)   #=> 240

strict_partitions

strict_partitions(n)
strict_partitions(n, max_value)

Return an array of arrays with the additive strict partitions of n (partitions where all integers in a single partition are strictly unique/different), such that their sum equals n.

say strict_partitions(6)       #=> [[6], [5, 1], [4, 2], [3, 2, 1]]

When the max_value argument is provided, it returns only the strict partitions with value <= max_value.

say strict_partitions(12, 5)   #=> [[5, 4, 3], [5, 4, 2, 1]]

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.

say 20.of { .subfactorial }        #=> OEIS: A000166
say 20.of { .subfactorial(2) }     #=> OEIS: A000387

Aliases: rencontres

subsets

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

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

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

superfactorial

superfactorial(n)

Returns the product of first n factorials. (OEIS A000178)

say 10.of { .superfactorial }   #=> OEIS: A000178

tuples

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

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

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

Aliases: variations

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 0..n-1, or iterates over the k-tuples with repetition when a block is given.

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

Aliases: variations_with_repetition

SEQUENCES, POLYNOMIALS & FIGURATE NUMBERS

bern

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

Returns the n-th Bernoulli number of the first kind (with B_1 = -1/2). When an additional argument is provided, it returns the n-th Bernoulli polynomial evaluated at x.

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

For Bernoulli numbers of the second kind (with B_1 = 1/2) use:

bernoulli(n, 1)     # here B_1 = 1/2

Aliases: bernfrac, bernoulli, bernoulli_number

bernoulli_numbers

bernoulli_numbers(n)

Returns an array containing the Bernoulli numbers with indices in the range 0..n.

say bernoulli_numbers(5)    #=> [1, -1/2, 1/6, 0, -1/30, 0]
say bernoulli_numbers(6)    #=> [1, -1/2, 1/6, 0, -1/30, 0, 1/42]

bernoulli_polynomial

bernoulli_polynomial(n)
bernoulli_polynomial(n, x)

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

say bernoulli_polynomial(3, 2)   #=> 3

When x is omitted, a Polynomial object is returned.

bernreal

bernreal(n)
bernreal(n,x)

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

say bernreal(10)     #=> 0.0757575757575757575757575757575757575757575757576

When a second argument x is given, it return the n-th Bernoulli polynomial evaluated at x:

B_n(x)      # as a floating-point number

catalan

n.catalan
n.catalan(k)

Returns the n-th Catalan number.

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

say 10.catalan          #=> 16796
say catalan(4, 2)       #=> 9

centered_polygonal

n.centered_polygonal(k)

Returns the n-th centered k-gonal number.

say 15.of {|n| centered_polygonal(n, 3) }  # centered triangular numbers
say 15.of {|n| centered_polygonal(n, 6) }  # centered hexagonal numbers

centered_polygonal_root

n.centered_polygonal_root(k)

Returns the centered k-gonal root of n. Also defined for complex numbers.

say centered_polygonal_root(n, 3)      # centered triangular root
say centered_polygonal_root(n, 5)      # centered pentagonal root

centered_pyramidal

n.centered_pyramidal(k)

Returns the n-th centered k-gonal pyramidal number.

say centered_pyramidal(4, 6)     #=> 54

centered_pyramidal_root

n.centered_pyramidal_root(k)

Returns the centered k-gonal pyramidal root of n. Also defined for complex numbers.

say centered_pyramidal_root(54, 6)   #=> 4

chebyshevT

chebyshevT(n)
chebyshevT(n, x)

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

say chebyshevT(3, 2)      #=> 26

Defined as:

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

When x is omitted, a Polynomial object is returned.

Aliases: chebyshevt

chebyshevTmod

chebyshevTmod(n, x, m)

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

say chebyshevTmod(3, 2, 100)     #=> 26

chebyshevU

chebyshevU(n)
chebyshevU(n, x)

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

say chebyshevU(3, 2)      #=> 56

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 x is omitted, a Polynomial object is returned.

Aliases: chebyshevu

chebyshevUmod

chebyshevUmod(n, x, m)

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

say chebyshevUmod(3, 2, 100)     #=> 56

collatz

n.collatz

Returns the number of halving and tripling steps to reach 1 in the 3x+1 Collatz problem.

say 27.collatz        #=> 111

cyclotomic

cyclotomic(n)
cyclotomic(n,x)

Returns the n-th Cyclotomic polynomial evaluated at x.

say cyclotomic(12, 10)      #=> 9901

When x is omitted, a Polynomial object is returned:

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

Aliases: cyclotomic_polynomial

cyclotomicmod

cyclotomicmod(n, x, m)

Returns the n-th Cyclotomic polynomial evaluated at x, computed modulo m.

say cyclotomicmod(30!, 5040, 2**128 + 1)

It returns NaN when moebius(d) = -1 and gcd(n/d, m) != 1, for some d|n.

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 x, when x is given:

say euler(10, 5)    #=> 1981100

Aliases: euler_number

euler_numbers

euler_numbers(n)

Returns an array containing the Euler numbers with indices in the range 0..n.

say euler_numbers(9)    #=> [1, 0, -1, 0, 5, 0, -61, 0, 1385, 0]
say euler_numbers(10)   #=> [1, 0, -1, 0, 5, 0, -61, 0, 1385, 0, -50521]

euler_polynomial

euler_polynomial(n)
euler_polynomial(n, x)

Returns the n-th Euler polynomial evaluated at x.

say euler_polynomial(3, 2)   #=> 9/4

When x is omitted, a Polynomial object is returned.

faulhaber

faulhaber(n,k)

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

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

The value for k must be a nonnegative native integer.

Aliases: faulhaber_sum

faulhaber_polynomial

faulhaber_polynomial(n)
faulhaber_polynomial(n, x)

Computes the n-th Faulhaber polynomials evaluated at 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 x is omitted, a Polynomial object is returned.

faulhaber_range

faulhaber_range(a, b, k)

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

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

The value for k must be a nonnegative native integer.

fibbinary

fibbinary(n)

Returns the n-th Fibbinary number. (OEIS A003714)

say 10.of { .fibbinary }            # [0, 1, 2, 4, 5, 8, 9, 10, 16, 17]
say 30.of { .fibbinary.as_bin }     # Zeckendorf representation of n (OEIS: A014417)

zeckendorf

var terms = n.zeckendorf

Returns the Zeckendorf representation of n as an array of Fibonacci numbers.

By Zeckendorf's theorem, every positive integer can be written uniquely as a sum of distinct, non-consecutive Fibonacci numbers. This function returns the Fibonacci terms in decreasing order.

say 10.zeckendorf     # => [8, 2]
say 42.zeckendorf     # => [34, 8]
say 100.zeckendorf    # => [89, 8, 3]

The returned terms always satisfy:

n == terms.sum

and no two consecutive Fibonacci numbers appear in the representation.

fibmod

fibmod(n,m)
fibmod(n,k,m)

Efficiently computes the n-th Fibonacci number modulo m.

fibmod(n,m) == (fibonacci(n) % m)

When three arguments are given, it computes the n-th k-th order Fibonacci number modulo m.

fibmod(n,k,m) == (fibonacci(n,k) % m)

Aliases: fibonacci_mod, fibonaccimod

fibonacci

fibonacci(n)
fibonacci(n,k)

Returns the n-th Fibonacci number. (OEIS A000045)

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

When an additional argument k is provided, it returns the n-th k-th order Fibonacci number. For k > 2, the sequence starts with zeros up to index k-2.

say 20.of { .fibonacci(3) }   #=> Tribonacci numbers
say 20.of { .fibonacci(4) }   #=> Tetranacci numbers
say 20.of { .fibonacci(5) }   #=> Pentanacci numbers
say fibonacci(2, 5)           #=> 0

Aliases: fib

fibonacci_inverse

fibonacci_inverse(n)

Returns the index of the Fibonacci number n.

If n = F(k), where F(k) denotes the k-th Fibonacci number, the function returns k. If n is not a Fibonacci number, the function returns -1.

The Fibonacci indexing used is

F(0) = 0
F(1) = 1
F(2) = 1
F(3) = 2
F(4) = 3
...

Examples:

say 144.fibonacci_inverse    # => 12
say 145.fibonacci_inverse    # => -1

Aliases: fib_inverse

fusc

fusc(n)

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

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

genocchi

genocchi(n)

Returns the n-th Genocchi number (of first kind). (OEIS A001469)

Defined as:

gennochi(n) = 2*(1 - 4^n) * bernoulli(2*n)

hermiteH

hermiteH(n)
hermiteH(n, x)

Physicists' Hermite polynomials: H_n(x).

say hermiteH(3, 2)    #=> 40

When x is omitted, a Polynomial object is returned.

Aliases: hermite_polynomialH

hermiteHe

hermiteHe(n)
hermiteHe(n, x)

Probabilists' Hermite polynomials: He_n(x).

say hermiteHe(3, 2)   #=> 2

When x is omitted, a Polynomial object is returned.

Aliases: hermite_polynomialHe

ipolygonal_root

ipolygonal_root(n,k)

First integer k-gonal root of n.

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

ipolygonal_root2

ipolygonal_root2(n,k)

Second integer k-gonal root of n.

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

isolve_quadratic

isolve_quadratic(a,b,c)

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

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

Example:

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

Aliases: integer_quadratic_formula, iquadratic_formula

is_centered_polygonal

n.is_centered_polygonal(k)

Returns true if n is a centered k-gonal number.

say 15.by { .is_centered_polygonal(3) }  # centered triangular numbers
say 15.by { .is_centered_polygonal(6) }  # centered hexagonal numbers

is_chebyshev

n.is_chebyshev

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

Aliases: is_chebyshev_psp, is_chebyshev_pseudoprime

is_fib

n.is_fib

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

Aliases: is_fibonacci

is_lucas

n.is_lucas

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

is_polygonal

is_polygonal(n,k)

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

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

is_polygonal2

is_polygonal2(n,k)

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

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

is_pyramidal

n.is_pyramidal(k)

Returns true if n is a k-gonal pyramidal number.

say 15.by { .is_pyramidal(3) }  #=> tetrahedral numbers
say 15.by { .is_pyramidal(5) }  #=> pentagonal pyramidal numbers

jacobsthal

jacobsthal(n)

Returns the n-th Jacobsthal number. (OEIS A001045)

Equivalent to: lucasU(1, -2, n).

jacobsthal_lucas

jacobsthal_lucas(n)

Returns the n-th Jacobsthal-Lucas number. (OEIS A014551)

Equivalent to: lucasV(1, -2, n).

laguerre

laguerre(n)
laguerre(n, x)

Laguerre polynomials: L_n(x).

say laguerre(3, 2)    #=> -1/3

When x is omitted, a Polynomial object is returned.

Aliases: laguerreL, laguerre_polynomial

legendreP

legendreP(n)
legendreP(n,x)

Legendre polynomials: P_n(x).

say legendreP(3, 2)   #=> 17

When x is omitted, a Polynomial object is returned.

Aliases: legendrep, legendre_polynomial

lucas

lucas(n)

Returns the n-th Lucas number.

say 10.of { .lucas }  #=> [2, 1, 3, 4, 7, 11, 18, 29, 47, 76]

lucas_mod

lucasmod(n,m)

Efficiently compute the n-th Lucas number modulo m.

say lucas_mod(100, 1000)    #=> (lucas(100) % 1000)

Aliases: lucasmod

lucasU

lucasU(P, Q, n)

The Lucas 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: lucasu

lucasUmod

lucasUmod(P,Q,n,m)

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

say lucasUmod(1, -1, 100, 1000)   #=> (lucasU(1,-1,100) % 1000)

Aliases: lucasumod

lucasUVmod

lucasUVmod(P,Q,n,m)

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

Equivalent with:

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

Aliases: lucasuvmod

lucasV

lucasV(P, Q, n)

The Lucas 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: lucasv

lucasvmod

lucasVmod(P,Q,n,m)

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

say lucasvmod(1, -1, 100, 1000)   #=> (lucasV(1,-1,100) % 1000)

Aliases: lucasVmod

motzkin

motzkin(n)

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

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

padovan

padovan(n)

Returns the n-th Padovan number (OEIS A000931), defined as:

a(0) = 1
a(1) = 0
a(2) = 0
a(n) = a(n-2) + a(n-3)

Example:

say 10.of { .padovan }   #=> [1, 0, 0, 1, 0, 1, 1, 1, 2, 2]

padovanmod

padovanmod(n,m)

Returns the n-th Padovan number, modulo m.

say padovanmod(1000, 97)   #=> (padovan(1000) % 97)

pell

pell(n)

Returns the n-th Pell number. (OEIS A000129)

say 10.of { .pell }   #=> [0, 1, 2, 5, 12, 29, 70, 169, 408, 985]

Equivalent to: lucasU(2, -1, n).

pell_lucas

pell_lucas(n)

Returns the n-th Pell-Lucas number. (OEIS A002203)

say 10.of { .pell_lucas }   #=> [2, 2, 6, 14, 34, 82, 198, 478, 1154, 2786]

Equivalent to: lucasV(2, -1, n).

perrin

perrin(n)

Returns the n-th Perrin number (OEIS A001608), defined as:

a(0) = 3
a(1) = 0
a(2) = 2
a(n) = a(n-2) + a(n-3)

Example:

say 10.of { .perrin }   #=> [3, 0, 2, 3, 2, 5, 5, 7, 10, 12]

perrinmod

perrinmod(n,m)

Returns the n-th Perrin number, modulo m.

say perrinmod(1000, 97)   #=> (perrin(1000) % 97)

This function can be used as a pseudo-primality test for n > 1:

perrinmod(n, n) == 0    #=> n is a prime number or a Perrin pseudoprime

The list of Perrin pseudoprimes is given by OEIS A013998.

pisano_period

pisano_period(n)

Returns the n-th Pisano number: period of Fibonacci numbers mod n.

say pisano_period(10!)          #=> 86400
say pisano_period(2**128 + 1)   #=> 28356863910078205764000346543980814080

polygonal

polygonal(n,k)

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

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

polygonal_inverse

polygonal_inverse(n)

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

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

Aliases: inverse_polygonal

polygonal_root

polygonal_root(n,k)

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

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

polygonal_root2

polygonal_root2(n,k)

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

say polygonal_root2(n, 5)      # second pentagonal root

pyramidal

n.pyramidal(k)

Returns the n-th k-gonal pyramidal number.

say 15.of {|n| pyramidal(n, 3) }  # tetrahedral numbers
say 15.of {|n| pyramidal(n, 5) }  # pentagonal pyramidal numbers

pyramidal_root

n.pyramidal_root(k)

Returns the k-gonal pyramidal root of n. Also defined for complex numbers.

say pyramidal_root(pyramidal(1234, 3), 3)     #=> 1234

ramanujan_tau

ramanujan_tau(n)

Ramanujan's tau function.

say ramanujan_tau(1)    #=> 1
say ramanujan_tau(2)    #=> -24
say ramanujan_tau(3)    #=> 252

secant_number

secant_number(n)

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

say secant_number(0)    #=> 1
say secant_number(1)    #=> 1
say secant_number(2)    #=> 5
say secant_number(3)    #=> 61

solve_pell

solve_pell(D)
solve_pell(D, N)

Finds all fundamental integer solutions (x, y) of the generalized Pell equation:

x^2 - D*y^2 = N

where D is a positive non-square integer and N is an arbitrary integer (default: 1).

For the classical Pell equations (N = ±1), the method computes the fundamental solution using the continued fraction expansion of sqrt(D). When available and advantageous, PARI/GP is used for large inputs to accelerate the computation. (with Num!USE_PARI_GP = true)

For arbitrary values of N, the method implements the LMM (Lagrange–Matthews–Mollin) continued-fraction algorithm to enumerate all primitive solution classes. Each solution is then reduced to its fundamental representative under multiplication by the fundamental unit of Q(sqrt(D)).

Returns an array of [x, y] pairs sorted by increasing x. If the equation has no integer solutions, an empty array is returned.

say solve_pell(61)      #=> [[1766319049, 226153980]]
say solve_pell(61, -1)  #=> [[29718, 3805]]
say solve_pell(10, 9)   #=> [[3, 0], [7, 2], [7, -2]]

# Iterate over solutions
var (D, N) = (10, 9)
var U = Quadratic(solve_pell(D)[0]..., D)
var V = Quadratic(solve_pell(D, N)[0]..., D)
for n in (0..5) { say (U**n * V) }

solve_lcg

solve_lcg(n, r, m)

Return the smallest solution for x to the following linear congruence:

n*x == r (mod m)

Example:

say solve_lcg(143, 44, 231)     #=> 10

Return NaN when no solution exists (i.e.: when r is not divisible by gcd(n, m)).

Aliases: solve_linear_congruence

solve_lcg_all

solve_lcg_all(n, r, m)

Return an array with the values of x satisfying the following linear congruence:

n*x == r (mod m)

Example:

say solve_lcg_all(3, 12, 15)    #=> [4, 9, 14]

Aliases: linear_congruence

solve_quadratic

solve_quadratic(a,b,c)

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

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

Example:

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

Aliases: quadratic_formula

solve_cubic

solve_cubic(a, b, c, d)

Returns a list of (complex) solutions (x_1, x_2, x_3) to the cubic equation:

a*x^3 + b*x^2 + c*x + d = 0

Aliases: cubic_formula

solve_quadraticQ

solve_quadraticQ(a,b,c)

Returns a list of Quadratic objects, which are solutions (x_1, x_2) to the quadratic equation: a*x^2 + b*x + c = 0.

say [solve_quadraticQ(3,4,5)]     #=> [Quadratic(-2/3, 1/6, -44, 0), Quadratic(-2/3, -1/6, -44, 0)]

Aliases: quadratic_formulaQ

solve_quadratic_mod

solve_quadratic_mod(a, b, c, m)

Returns all solutions to the quadratic congruence:

a*x^2 + b*x + c ≡ 0 (mod m)

The result is returned as an array containing all distinct solutions modulo m, sorted in ascending order. If the congruence has no solutions, an empty array is returned.

If the leading coefficient a is zero, the congruence degenerates into the linear congruence:

b*x + c ≡ 0 (mod m)

which is solved automatically using solve_lcg_all.

When m is odd and the leading coefficient a is invertible modulo m, the function uses the modular analogue of the quadratic formula. After normalizing the equation by a, it completes the square to obtain:

(x + b/(2a))^2 ≡ (b/2a)^2 - c/a (mod m)

and computes all modular square roots using sqrtmod_all.

For the general case (including even moduli or noninvertible leading coefficients), the congruence is transformed into the equivalent square-root problem:

k^2 ≡ b^2 - 4ac (mod 4am)

Solutions are then reconstructed from the corresponding square roots modulo 4am.

Examples:

# x^2 - 1 ≡ 0 (mod 13)
say solve_quadratic_mod(1, 0, -1, 13)    #=> [1, 12]

# Linear congruence: 5x + 2 ≡ 0 (mod 11)
say solve_quadratic_mod(0, 5, 2, 11)     #=> [4]

# 3*x^2 + 4*x + 5 ≡ 0 (mod 124)
say solve_quadratic_mod(3, 4, 5, 124)    #=> [47, 55, 109, 117]

Aliases: modular_quadratic_formula, quadratic_congruence

solve_cubic_mod

solve_cubic_mod(A, B, C, D, m)

Returns all solutions to the cubic polynomial congruence:

A*x^3 + B*x^2 + C*x + D ≡ 0 (mod m)

The result is returned as an array containing all distinct solutions modulo m, sorted in ascending order. If no solutions exist, or if the equation cannot be transformed into a monic cubic modulo m, an empty array is returned.

The function requires that the modular inverses of 2, 3, 27, and the leading coefficient A exist modulo m. If any of these inverses do not exist, the function returns an empty array.

# x^3 - 1 ≡ 0 (mod 13)
say solve_cubic_mod(1, 0, 0, -1, 13)     # [1, 3, 9]

# x^3 + 2x + 5 ≡ 0 (mod 17)
say solve_cubic_mod(1, 0, 2, 5, 17)      # [2]

solve_quadratic_form

var solutions = solve_quadratic_form(d, n)

Given a positive integer n and a positive integer d, returns an array with [x,y] solutions to the equation:

x^2 + d*y^2 = n

Example:

say solve_quadratic_form(13, 97)      #=> []
say solve_quadratic_form(18, 43*97)   #=> [[11, 15], [61, 5]]

Returns an empty array if there are no solutions.

Aliases: cornacchia

solve_binary_quadratic_form

var solutions = solve_binary_quadratic_form(A, B, C, n)

Returns all integer solutions [x, y] of the binary quadratic equation:

Ax² + Bxy + Cy² = n.

The result is returned as an array of pairs [x, y], sorted lexicographically. If no integer solutions exist, an empty array is returned.

The current implementation applies to positive definite binary quadratic forms, namely those satisfying:

A > 0
4AC - B² > 0.

It transforms the equation into the equivalent form:

X² + DY² = 4An,

where

D = 4AC - B²,
X = 2Ax + By,
Y = y.

The transformed equation is solved using solve_quadratic_form, after which all valid sign combinations are examined to recover the corresponding integer solutions [x, y].

Examples:

# x² + y² = 25
say solve_binary_quadratic_form(1, 0, 1, 25)
# => [[-5,0], [-4,-3], [-4,3], [-3,-4], [-3,4],
#     [0,-5], [0,5], [3,-4], [3,4], [4,-3],
#     [4,3], [5,0]]

# x² + xy + y² = 13
say solve_binary_quadratic_form(1, 1, 1, 13)
# => [[-4, 1], [-4, 3], [-3, -1], [-3, 4],
#     [-1, -3], [-1, 4], [1, -4], [1, 3], [3, -4],
#     [3, 1], [4, -3], [4, -1]]

Every returned pair satisfies:

A*x² + B*x*y + C*y² == n

solve_ternary_quadratic_form

var solutions = solve_ternary_quadratic_form(a, b, c, n)

Returns all the positive solutions to the ternary quadratic form (up to rearrangement):

a*x^2 + b*y^2 + c*z^2 = n

Example:

say solve_ternary_quadratic_form(1, 1, 1, 18)   #=> [[0, 3, 3], [1, 1, 4]]
say solve_ternary_quadratic_form(1, 1, 2, 18)   #=> [[0, 0, 3], [0, 4, 1], [1, 3, 2], [3, 3, 0]]
say solve_ternary_quadratic_form(1, 2, 2, 18)   #=> [[0, 0, 3], [4, 0, 1]]
say solve_ternary_quadratic_form(1, 2, 3, 18)   #=> [[0, 3, 0], [2, 1, 2], [4, 1, 0]]

It does not support degenerate forms where a coefficient is zero or negative. It requires:

a > 0, b > 0, c > 0, n > 0

sum_of_polygonals

sum_of_polygonals(n, k)

Returns all representations as a sum of two k-gonal numbers: n = P(x,k) + P(y,k) with 0 ≤ x ≤ y, as an array of [x, y] arrays, sorted by x.

Example:

say sum_of_polygonals(18736, 3)      # sums of two triangular numbers
say sum_of_polygonals(99025, 4)      # sums of two squares numbers
say sum_of_polygonals(10622, 5)      # sums of two pentagonal numbers

sum_of_squares

sum_of_squares(n)

Returns all representations as a sum of two squares: n = x² + y² with 0 ≤ x ≤ y, as an array of [x, y] arrays, sorted by x.

Example:

say sum_of_squares(99025)

Output:

[[41, 312], [48, 311], [95, 300], [104, 297], [183, 256], [220, 225]]

sum_of_cubes

sum_of_cubes(n)

Returns all integer solutions (x, y) satisfying:

x^3 + y^3 = n

The result is returned as an array of pairs, where each pair is itself an array [x, y].

Returns an empty array when no integer solutions exist or when n = 0.

This implementation uses the factorization identity:

x^3 + y^3 = (x + y)(x^2 - xy + y^2)

to reduce the problem to testing the divisors of |n|. For each divisor A = x + y, the algorithm derives the corresponding quadratic equation, checks whether its discriminant is a perfect square, and reconstructs the integer solution when possible.

The algorithm naturally supports both positive and negative values of n.

Examples:

say sum_of_cubes(35)      # [[2, 3]]
say sum_of_cubes(1729)    # [[1, 12], [9, 10]]
say sum_of_cubes(-728)    # [[-12, 10], [-9, 1], [-8, -6]]
say sum_of_cubes(30)      # []

primitive_pythagorean_triples

primitive_pythagorean_triples(limit)

Returns all primitive Pythagorean triples whose perimeter does not exceed limit.

Each returned triple is an array [a, b, c], satisfying:

a² + b² = c²
gcd(a, b, c) = 1
a < b < c

A primitive Pythagorean triple is one whose three components are coprime. Every primitive Pythagorean triple appears exactly once in the returned array.

say primitive_pythagorean_triples(50)
#=> [[3, 4, 5], [5, 12, 13], [8, 15, 17]]

tangent_number

tangent_number(n)

Returns the n-th tangent/zag number. (OEIS A000182)

say 5.of { .tangent_number }   #=> [0, 1, 2, 16, 272]

tetranacci

tetranacci(n)

Returns the n-th Tetranacci number. (OEIS A000078)

Alias for: fibonacci(n,4).

say 10.of { .tetranacci }   #=> [0, 0, 0, 1, 1, 2, 4, 7, 13, 24]

tribonacci

tribonacci(n)

Returns the n-th Tribonacci number. (OEIS A000073)

Alias for: fibonacci(n, 3).

say 10.of { .tribonacci }   #=> [0, 0, 1, 1, 2, 4, 7, 13, 24, 44]

NUMBER THEORY: PRIMES & PRIMALITY

all_prime

all_prime(...)

Returns true if all the given values are prime numbers.

say all_prime(2, 3, 5, 7)   #=> true
say all_prime(2, 3, 4, 5)   #=> false

is_abs_euler_psp

n.is_abs_euler_psp

Returns true if n is an absolute Euler pseudoprime: p-1 | (n-1)/2 for every p|n, where n is a squarefree composite integer. (OEIS A033181)

say 10.by { .is_abs_euler_psp }     #=> OEIS: A033181

Aliases: is_absolute_euler_psp

is_aks_prime

n.is_aks_prime

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

is_balanced_prime

is_balanced_prime(n, order=1)

Returns true if n is a balanced prime of a given order.

say 10.by { .is_balanced_prime }       #=> OEIS: A006562
say 10.by { .is_balanced_prime(2) }    #=> OEIS: A082077
say 10.by { .is_balanced_prime(3) }    #=> OEIS: A082078

A prime number p_k is a balanced prime, of order n, iff the following equality holds:

p_k = 1/(2*n) * Sum_{i=1..n} (p_{k-i} + p_{k+i})

is_bfsw_psp

n.is_bfsw_psp

Returns true if n passes a slightly stronger and faster variant of the BFW primality test, based on the Lucas-V function, checking the following congruences:

V_{n+1} == 2*Q (mod n)
Q^(n+1) == Q^2 (mod n)

where Q = -2 and P is the first integer >= 2 such that D = P^2 - 4*Q gives kronecker(D, n) == -1.

There are no known composites that pass this test.

is_bpsw_prime

n.is_bpsw_prime

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

is_ecpp_prime

n.is_ecpp_prime

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

is_emirp

n.is_emirp(base = 10)

Returns true if n is prime and its reversal is a different prime (in the given base). Also known as emirps.

say 10.by { .is_emirp }        #=> decimal emirps (OEIS: A006567)
say 10.by { .is_emirp(2) }     #=> binary emirps (OEIS: A080790)

is_euler_psp

n.is_euler_psp(bases...)

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

Aliases: is_euler_pseudoprime

is_fib_psp

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

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

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

Aliases: is_lucasu_psp, is_lucasU_psp, is_fibonacci_psp, is_lucasU_pseudoprime, is_fibonacci_pseudoprime

is_frobenius_psp

n.is_frobenius_psp(a,b)

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

Aliases: is_frobenius_pseudoprime

is_gaussian_prime

is_gaussian_prime(a,b)

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

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

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

is_khashin_psp

n.is_khashin_psp

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

Aliases: is_khashin_pseudoprime, is_frobenius_khashin_psp, is_frobenius_khashin_pseudoprime

is_llr_prime

n.is_llr_prime(k)

Returns true if 2^n * k - 1 is prime, using the Lucas-Lehmer-Riesel (LLR) primality test.

say 15.by { .is_llr_prime(3) }   # numbers n such that 2^n * 3 - 1 is prime

is_lucas_psp

n.is_lucas_psp

Return true if n is a Lucas pseudoprime.

Aliases: is_lucas_pseudoprime

is_lucasv_psp

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

Returns true if n passes the Lucas test to the V sequence, using the parameters P and Q.

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: is_lucasV_psp, is_bruckman_lucas_psp, is_lucasV_pseudoprime, is_bruckman_lucas_pseudoprime

is_mersenne_prime

is_mersenne_prime(p)

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

say 31.is_mersenne_prime    #=> true (2^31 - 1 is prime)
say 11.is_mersenne_prime    #=> false

Uses the Lucas-Lehmer test for efficiency. =cut

is_nm1_prime

n.is_nm1_prime

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

Aliases: is_pm1_prime, is_nminus1_prime

is_np1_prime

n.is_np1_prime

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

Aliases: is_pp1_prime, is_nplus1_prime

is_over_psp

n.is_over_psp

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

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

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: is_over_pseudoprime, is_overpseudoprime

is_pell_lucas_psp

n.is_pell_lucas_psp

It returns true if V_n(2, -1) = 2 (mod n). (OEIS A270342)

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: is_pell_lucas_pseudoprime

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: is_pell_pseudoprime

is_perrin_psp

n.is_perrin_psp

Returns true if n passes the Perrin primality test.

Aliases: is_perrin_pseudoprime

is_plumb_psp

n.is_plumb_psp

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

Aliases: is_euler_plumb_psp, is_plumb_pseudoprime, is_euler_plumb_pseudoprime

is_prime

n.is_prime

Returns true if n is a prime number.

Values less than 2 are not prime.

say 17.is_prime           #=> true
say 18.is_prime           #=> false
say 1.is_prime            #=> false
say is_prime(2**127 - 1)  #=> true

For small native integers, this uses fast deterministic tests. For larger integers, it uses strong probable-prime tests together with additional checks from the underlying libraries.

is_prob_prime

n.is_prob_prime

Returns true if n is probably prime. Performs trial division followed by a B-PSW test.

say is_prob_prime(1003200318091)   #=> true

Faster than is_prime but may (extremely rarely) give false positives for large numbers.

is_proth_prime

n.is_proth_prime(k)

Returns true if 2^n * k + 1 is prime, using the Proth primality test.

say 15.by { .is_proth_prime(3) }   # numbers n such that 2^n * 3 + 1 is prime

is_prov_prime

n.is_prov_prime

Returns true if n is definitely a prime.

Aliases: is_provable_prime

is_psp

n.is_psp(bases...)

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

Each base must be coprime to n.

Aliases: is_fermat_psp, is_pseudoprime, is_fermat_pseudoprime

is_safe_prime

n.is_safe_prime

It returns true if both n and (n-1)/2 are prime. (OEIS A005385)

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

is_sexy_prime

is_sexy_prime(n)

Returns true if n and n+6 are both prime numbers.

say is_sexy_prime(5)     #=> true    (5 and 11 are both prime)
say is_sexy_prime(9)     #=> false

is_sophie_germain

is_sophie_germain(n)

Returns true if n and 2*n + 1 are both prime numbers.

say is_sophie_germain(11)    #=> true    (11 and 23 are both prime)
say is_sophie_germain(7)     #=> false

Aliases: is_sophie_germain_prime

is_strong_fib

n.is_strong_fib

Returns true if 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: is_strong_fib_psp, is_strong_fibonacci, is_strong_fibonacci_psp, is_strong_fibonacci_pseudoprime

is_strong_lucas_psp

n.is_strong_lucas_psp

Return true if n is a strong Lucas pseudoprime.

Aliases: is_strong_lucas_pseudoprime

is_strong_psp

n.is_strong_psp(bases...)

Return true if n is a strong Fermat pseudoprime to the given bases.

Each base must be coprime to n.

Aliases: miller_rabin, is_strong_fermat_psp, is_strong_pseudoprime, is_strong_fermat_pseudoprime

is_stronger_lucas_psp

n.is_stronger_lucas_psp

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

Aliases: is_extra_strong_lucas_psp, is_stronger_lucas_pseudoprime, is_extra_strong_lucas_pseudoprime

is_strongish_lucas_psp

n.is_strongish_lucas_psp

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

Aliases: is_strongish_lucas_pseudoprime

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)

# 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: is_super_pseudoprime, is_superpseudoprime

is_twin_prime

is_twin_prime(n)

Returns true if n and n+2 are both prime numbers.

say is_twin_prime(11)    #=> true    (11 and 13 are both prime)
say is_twin_prime(7)     #=> false

is_underwood_psp

n.is_underwood_psp

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

Aliases: is_underwood_pseudoprime, is_frobenius_underwood_psp, is_frobenius_underwood_pseudoprime

is_vpsp

n.is_vpsp

Returns true if n passes the Lucas-V BFW test: V_{n+1} == 2*Q (mod n) for carefully chosen values of P and Q.

The composite numbers that pass this test are extremely rare. (OEIS A365514)

When combined with a strong Fermat pseudoprime test to base 2, there are no known composites that pass both tests.

Reference:

Strengthening the Baillie-PSW primality test
https://arxiv.org/abs/2006.14425

Aliases: is_bfw_psp

next_prime

n.next_prime

Returns the next prime number greater than n.

say 10.next_prime           #=> 11
say 100.next_prime          #=> 101
say (10**20 -> next_prime)  #=> 100000000000000000039

next_primes

n.next_primes(start=2)

Returns an array of n consecutive primes starting from start (if omitted, then start = 2).

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

The value for start can be any arbitrarily large integer:

say 10.next_primes(2**128 + 1)    # first 10 primes >= 2^128

next_twin_prime

n.next_twin_prime

Returns next twin prime number larger than n.

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

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).

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

prev_prime

n.prev_prime

Returns the previous prime number smaller than n.

say 11.prev_prime       #=> 7
say 100.prev_prime      #=> 97

Returns nil if n <= 2.

prev_primes

n.prev_primes(start)

Returns an array with n consecutive decreasing prime numbers starting from start.

say 5.prev_primes(100)  #=> [97, 89, 83, 79, 73]

primality_pretest

n.primality_pretest

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

prime

prime(n)
nth_prime(n)

Returns the n-th prime number (1-indexed).

say prime(1)            #=> 2
say prime(100)          #=> 541
say prime(1000000)      #=> 15485863

For convenience, prime(0) returns 1, although 1 is not a prime.

say prime(0)            #=> 1

Aliases: nth_prime

prime_cluster

prime_cluster(lo, hi, *diffs)

Returns an array of primes p such that p+diffs[0], p+diffs[1], etc, are also primes.

say prime_cluster(1, 100, 2)    # twin primes <= 100

Aliases: sieve_prime_cluster

prime_lower

n.prime_lower

Lower bound for the n-th prime.

Aliases: nth_prime_lower

prime_upper

n.prime_upper

Upper bound for the n-th prime.

Aliases: nth_prime_upper

primepi

pi(n)
primepi(n)
primepi(a, b)

Counts the number of primes less than or equal to n, or in the range [a, b].

In the two-argument form, the lower bound is effectively clamped to 2, so values below 2 contribute nothing.

say primepi(1)         #=> 0
say primepi(100)       #=> 25
say primepi(50, 100)   #=> 10
say pi(10**12)         #=> 37607912018

Aliases: prime_count, count_primes

primepi_lower

n.primepi_lower

Lower bound for prime_count(n).

Aliases: prime_count_lower

primepi_upper

n.primepi_upper

Upper bound for prime_count(n).

Aliases: prime_count_upper

primes

primes(n)
primes(a, b)

Returns an array of prime numbers <= n, or in the range [a,b].

say primes(20)          #=> [2, 3, 5, 7, 11, 13, 17, 19]
say primes(10, 20)      #=> [11, 13, 17, 19]

primes_each

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

Iterates over the primes in the given range.

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

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

Aliases: each_prime

primes_sum

primes_sum(n)
primes_sum(a,b,k=1)

Sum of the prime numbers <= n, or in the range a..b, each prime raised to the k-th power.

Aliases: prime_sum, sum_primes

strong_fermat_psp

k.strong_fermat_psp(base, n)
k.strong_fermat_psp(base, a, b)

Returns an array with all the strong Fermat pseudoprimes to the given base that are <= n or in the range a..b and have k distinct prime factors.

say 3.strong_fermat_psp(2, 1e6)         # 3-omega strong Fermat psp to base 2
say 4.strong_fermat_psp(2, 1e6, 1e7)    # 4-omega strong Fermat psp to base 2 in range

strong_fermat_psp_each

k.strong_fermat_psp_each(base, n, { ... })
k.strong_fermat_psp_each(base, a, b, { ... })

Iterates over the k-omega strong Fermat pseudoprimes <= n or in the range a..b.

3.strong_fermat_psp_each(2, 1e6, {|n| say n })          # iterate over 3-omega strong Fermat psp to base 2 <= 1e6
4.strong_fermat_psp_each(2, 1e7, 1e8, {|n| say n })     # iterate over 4-omega strong Fermat psp to base 2 in range 1e7..1e8

Aliases: each_strong_fermat_psp

twin_primes

twin_primes(n)
twin_primes(a,b)

Returns an array of primes p such that p+2 is also prime, where p <= n or in range a..b.

say twin_primes(100)        #=> [3, 5, 11, 17, 29, 41, 59, 71]
say twin_primes(50, 200)    #=> [59, 71, 101, 107, 137, 149, 179, 191, 197]

NUMBER THEORY: PRIME-RELATED FAMILIES

almost_prime_count

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

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

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: almost_primepi, pi_k

almost_prime_divisors

n.almost_prime_divisors
n.almost_prime_divisors(k)

Returns the k-almost prime divisors of n.

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

When k is omitted, an array of arrays with the k-almost prime divisors of n, for each k in the range 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]]

almost_prime_sum

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

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

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

almost_primes

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

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

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]

almost_primes_each

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

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

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: each_almost_prime

carmichael

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

Returns an array with all the Carmichael numbers <= n or in the range a..b that have exactly k prime factors.

say 3.carmichael(1e4)           # 3-Carmichael numbers <= 1e4
say 4.carmichael(1e4, 1e6)      # 4-Carmichael numbers in range 1e4..1e6

carmichael_each

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

Iterates over the Carmichael <= n or in the range a..b that have exactly k prime factors.

3.carmichael_each(1e5, {|n| say n })          # iterate over 3-Carmichael numbers <= 1e5
4.carmichael_each(1e5, 1e6, {|n| say n })     # iterate over 4-Carmichael in range 1e5..1e6

Aliases: each_carmichael

carmichael_strong_fermat

k.carmichael_strong_fermat(base, n)
k.carmichael_strong_fermat(base, a, b)

Returns an array with all the Carmichael numbers <= n or in the range a..b that are also strong Fermat pseudoprimes to base and have exactly k prime factors.

say 3.carmichael_strong_fermat(2, 1e7)       # 3-Carmichael strong Fermat pseudoprimes to base 2 <= 1e7
say 4.carmichael_strong_fermat(2, 1e7, 1e9)  # 4-Carmichael strong Fermat pseudoprimes to base 2, in range 1e7..1e9

Aliases: strong_fermat_carmichael, carmichael_strong_psp

carmichael_strong_fermat_each

k.carmichael_strong_fermat_each(base, n, { ... })
k.carmichael_strong_fermat_each(base, a, b, { ... })

Iterates over the Carmichael <= n or in the range a..b that are also strong Fermat pseudoprimes to base and have exactly k prime factors.

3.carmichael_strong_fermat_each(2, 1e7, {|n| say n })          # iterate over 3-Carmichael strong Fermat to base 2 <= 1e7
4.carmichael_strong_fermat_each(2, 1e5, 1e9, {|n| say n })     # iterate over 4-Carmichael strong Fermat to base 2 in range 1e5..1e9

Aliases: carmichael_strong_psp_each, each_carmichael_strong_psp, each_carmichael_strong_fermat, each_strong_fermat_carmichael, strong_fermat_carmichael_each

fermat_psp

k.fermat_psp(base, n)
k.fermat_psp(base, a, b)

Returns an array with all the k-omega Fermat pseudoprimes to the given base that are <= n or in the range a..b.

say 3.fermat_psp(2, 10000)            # 3-omega Fermat psp to base 2
say 4.fermat_psp(2, 10000, 100000)    # 4-omega Fermat psp to base 2 in range

fermat_psp_each

k.fermat_psp_each(base, n, { ... })
k.fermat_psp_each(base, a, b, { ... })

Iterates over the k-omega Fermat pseudoprimes <= n or in the range a..b.

3.fermat_psp_each(2, 1e5, {|n| say n })          # iterate over 3-Fermat psp to base 2 <= 1e5
4.fermat_psp_each(2, 1e5, 1e6, {|n| say n })     # iterate over 4-Fermat psp to base 2 in range 1e4..1e6

Aliases: each_fermat_psp

is_almost_prime

n.is_almost_prime(k=2)

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

Equivalently, k-almost primes are numbers n that satisfy bigomega(n) == k.

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

In order for n to have at least k prime factors, without having any prime factors less than or equal to B, it implies that n must be greater than B^k, since all the prime factors of n are greater than B.

By setting Num!USE_CONJECTURES = true, the function uses a conjectured approach based on Pollard's rho method to find a larger bound for B, which requires O(sqrt(B)) steps to find a prime factor less than B. Therefore, if we take B = 10^12, after c*10^6 iterations (we use c=2) of the Pollard rho method without success in finding a prime factor of n, it's very probable that n has no prime factor less than 10^12.

By enabling the conjectured approach, the function becomes about 5x faster than the rigorous method, for large enough inputs.

is_carmichael

n.is_carmichael

Determine if n is a Carmichael number or not.

is_imprimitive_carmichael

n.is_imprimitive_carmichael

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

say 325533792014488126487416882038879701391121.is_imprimitive_carmichael   # true

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

is_lucas_carmichael

n.is_lucas_carmichael

Returns true if n is a Lucas-Carmichael number. (OEIS A006972)

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

is_omega_prime

n.is_omega_prime(k=2)

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

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

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

In order for n to have at least k prime factors, without having any prime factors less than or equal to B, it implies that n must be greater than B^k, since all the prime factors of n are greater than B.

By setting Num!USE_CONJECTURES = true, the function uses a conjectured approach based on Pollard's rho method to find a larger bound for B, which requires O(sqrt(B)) steps to find a prime factor less than B. Therefore, if we take B = 10^12, after c*10^6 iterations (we use c=2) of the Pollard rho method without success in finding a prime factor of n, it's very probable that n has no prime factor less than 10^12.

By enabling the conjectured approach, the function becomes about 5x faster than the rigorous method, for large enough inputs.

is_prime_power

n.is_prime_power

Returns true if n is a power of some prime.

is_semiprime

n.is_semiprime

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

linear_forms_primes

linear_forms_primes(A, B, [a_1, b_1], ..., [a_k, b_k])

Returns an array with values m in the range [A, B], such that m * a_1 + b_1, m * a_2 + b_2, ..., m * a_k + b_k are all prime.

say linear_forms_primes(1, 100, [1, 0])                     # primes in range 1..100
say linear_forms_primes(1, 500, [1, 0], [1, 2])             # lesser of twin primes in range 1..500
say linear_forms_primes(1, 200, [2,-1], [4,-1], [6, -1])    # OEIS: A124486

lnsuperprimorial

lnsuperprimorial(n)

Natural logarithm of superprimorial(n).

Aliases: superprimorial_ln, superprimorial_log

lucas_carmichael

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

Returns an array with all the Lucas-Carmichael numbers <= n or in the range a..b that have exactly k prime factors.

say 3.lucas_carmichael(10000)         # 3-Lucas-Carmichael numbers
say 4.lucas_carmichael(10000, 100000) # 4-Lucas-Carmichael numbers in range

lucas_carmichael_each

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

Iterates over the Lucas-Carmichael <= n or in the range a..b that have exactly k prime factors.

3.lucas_carmichael_each(1e5, {|n| say n })          # iterate over 3-Lucas-Carmichael numbers <= 1e5
4.lucas_carmichael_each(1e5, 1e6, {|n| say n })     # iterate over 4-Lucas-Carmichael in range 1e4..1e6

Aliases: each_lucas_carmichael

next_almost_prime

n.next_almost_prime(k=2)

Returns the next k-almost prime greater than n.

next_omega_prime

n.next_omega_prime(k=2)

Returns the next k-omega prime greater than n.

next_prime_power

n.next_prime_power

Given a nonnegative integer n, returns the next prime power (p^k with k >= 1) greater than n.

next_semiprime

n.next_semiprime

Returns the next semiprime greater than n.

next_sphenic

n.next_sphenic

Returns the next number that is the product of three distinct primes.

next_squarefree_almost_prime

n.next_squarefree_almost_prime(k=2)

Returns the next squarefree k-almost prime greater than n.

next_squarefree_semiprime

n.next_squarefree_semiprime

Returns the next squarefree semiprime greater than n.

nth_almost_prime

nth_almost_prime(n, k=2)

Returns the n-th k-almost prime.

say nth_almost_prime(1e7, 2)   #=> 56168169
say nth_almost_prime(1e7, 3)   #=> 41657362
say nth_almost_prime(1e7, 4)   #=> 47997635

nth_omega_prime

nth_omega_prime(n, k=2)

Returns the n-th k-omega prime.

say nth_omega_prime(1e7, 2)    #=> 42314023
say nth_omega_prime(1e7, 3)    #=> 28013887
say nth_omega_prime(1e7, 4)    #=> 39780102

nth_prime_power

nth_prime_power(n)

Returns the n-th prime power p^k with k >= 1.

say nth_prime_power(1e12)   #=> 29996212395727

nth_sphenic

nth_sphenic(n)

Returns the n-th number that is the product of three distinct primes. (OEIS A007304)

say nth_sphenic(1e7)     #=> 48108421

nth_squarefree_almost_prime

nth_squarefree_almost_prime(n, k=2)

Returns the n-th squarefree k-almost prime.

say nth_squarefree_almost_prime(1e7, 2)     #=> 56173891
say nth_squarefree_almost_prime(1e7, 3)     #=> 48108421
say nth_squarefree_almost_prime(1e7, 4)     #=> 81556446

omega_prime_count

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

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

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: omega_primepi

omega_prime_divisors

n.omega_prime_divisors
n.omega_prime_divisors(k)

Returns the k-omega prime divisors of n.

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

When k is omitted, an array of arrays with the k-omega prime divisors of n, for each k in the range 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]]

omega_prime_sum

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

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

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

omega_primes

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

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

k-omega primes are numbers n that satisfy omega(n) == k.

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

omega_primes_each

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

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

k-omega primes are numbers n that satisfy omega(n) == k.

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: each_omega_prime

pn_primorial

pn_primorial(n)

Returns the product of the first n primes.

prev_almost_prime

n.prev_almost_prime(k=2)

Returns the previous k-almost prime smaller than n.

prev_omega_prime

n.prev_omega_prime(k=2)

Returns the previous k-omega prime smaller than n.

prev_prime_power

n.prev_prime_power

Returns the previous prime power smaller than n.

prev_semiprime

n.prev_semiprime

Returns the previous semiprime smaller than n.

prev_sphenic

n.prev_sphenic

Returns the previous number that is the product of three distinct primes.

prev_squarefree_almost_prime

n.prev_squarefree_almost_prime(k=2)

Returns the previous squarefree k-almost prime smaller than n.

prev_squarefree_semiprime

n.prev_squarefree_semiprime

Returns the previous squarefree semiprime smaller than n.

prime_power

n.prime_power

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

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

prime_power_count

prime_power_count(n)
prime_power_count(a,b)

Returns the count of prime powers <= n, or in the range a..b.

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]

prime_power_count_lower

prime_power_count_lower(n)

Lower bound for prime_power_count(n).

prime_power_count_upper

prime_power_count_upper(n)

Upper bound for prime_power_count(n).

prime_power_lower

prime_power_lower(n)

Lower bound for the n-th prime power.

Aliases: nth_prime_power_lower

prime_power_sum

prime_power_sum(n)
prime_power_sum(a,b)

Returns the sum of prime powers <= n, or in the range a..b.

Aliases: prime_powers_sum

prime_power_upper

prime_power_upper(n)

Upper bound for the n-th prime power.

Aliases: nth_prime_power_upper

prime_powers

prime_powers(n)
prime_powers(a,b)

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

say prime_powers(100)          # prime powers <= 100
say prime_powers(50, 100)      # prime powers in the range [50, 100]

prime_powers_each

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

Iterate over prime powers <= n, or in the range a..b:

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

Aliases: each_prime_power

prime_root

n.prime_root

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

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

semiprime

semiprime(n)

Returns the n-th semiprime number. (OEIS A001358)

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

Aliases: nth_semiprime

semiprime_count

semiprime_count(n)
semiprime_count(a,b)

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

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

semiprime_sum

semiprime_sum(n)
semiprime_sum(a,b)

Returns the sum of semiprimes <= n, or in the range a..b.

Aliases: semiprimes_sum

semiprimes

semiprimes(n)
semiprimes(a,b)

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

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

semiprimes_each

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

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

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

Aliases: each_semiprime

sphenic

sphenic(n)
sphenic(a, b)

Returns an array with the numbers with exactly three distinct prime factors <= n or in range a..b. (OEIS A007304)

say sphenic(200)    #=> OEIS: A007304

Aliases: sphenic_numbers

sphenic_count

sphenic_count(n)
sphenic_count(a,b)

Count of the numbers with exactly three distinct prime factors <= n or in range a..b. (OEIS A215218)

say sphenic_count(1e10)      #=> 1997171674
say sphenic_count(1e11)      #=> 19522428788

sphenic_each

sphenic_each(n, {|k| ... })
sphenic_each(a, b, {|k| ... })

Iterate over the numbers with exactly three distinct prime factors <= n or in range a..b.

sphenic_each(500, 1000, {|k| say k })

squarefree_almost_prime_count

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

Returns the count of squarefree k-almost primes <= n, or in the range a..b.

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: squarefree_almost_primepi, squarefree_pi_k

squarefree_almost_prime_sum

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

Returns the sum of squarefree k-almost primes <= n, or in the range a..b.

say 1.squarefree_almost_prime_sum(1000)          # sum of primes <= 1000
say 2.squarefree_almost_prime_sum(50, 1000)      # sum of squarefree semiprimes in range 50..1000

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 a..b.

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

squarefree_almost_primes_each

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

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

# 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: each_squarefree_almost_prime

squarefree_fermat_psp

k.squarefree_fermat_psp(base, n)
k.squarefree_fermat_psp(base, a, b)

Returns an array with all the squarefree Fermat pseudoprimes to the given base that are <= n or in the range a..b and have exactly k prime factors.

say 3.squarefree_fermat_psp(2, 10000)           # squarefree 3-Fermat psp to base 2
say 4.squarefree_fermat_psp(2, 10000, 100000)   # squarefree 4-Fermat psp to base 2 in range

squarefree_fermat_psp_each

k.squarefree_fermat_psp_each(base, n, { ... })
k.squarefree_fermat_psp_each(base, a, b, { ... })

Iterates over the squarefree Fermat pseudoprimes <= n or in the range a..b that have exactly k prime factors.

3.squarefree_fermat_psp_each(2, 1e5, {|n| say n })          # iterate over squarefree 3-Fermat psp to base 2 <= 1e5
4.squarefree_fermat_psp_each(2, 1e5, 1e6, {|n| say n })     # iterate over squarefree 4-Fermat psp to base 2 in range 1e4..1e6

Aliases: each_squarefree_fermat_psp

squarefree_semiprime_count

squarefree_semiprime_count(n)
squarefree_semiprime_count(a,b)

Counts the number of squarefree semiprimes <= n (OEIS A072613), or in the range a..b.

say 10.of {|k| squarefree_semiprime_count(10**k) }  # OEIS: A036351

Aliases: squarefree_semiprimes_count

squarefree_semiprime_sum

squarefree_semiprime_sum(n)
squarefree_semiprime_sum(a,b)

Returns the sum of squarefree semiprimes <= n, or in the range a..b.

Aliases: squarefree_semiprimes_sum

squarefree_semiprimes

squarefree_semiprimes(n)
squarefree_semiprimes(a,b)

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

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

squarefree_semiprimes_each

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

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

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

Aliases: each_squarefree_semiprime

squarefree_strong_fermat_psp

k.squarefree_strong_fermat_psp(base, n)
k.squarefree_strong_fermat_psp(base, a, b)

Returns an array with all the squarefree strong Fermat pseudoprimes to the given base that are <= n or in the range a..b and have exactly k prime factors.

say 3.squarefree_strong_fermat_psp(2, 1e6)        # squarefree strong 3-Fermat psp to base 2 <= 1e6
say 4.squarefree_strong_fermat_psp(2, 1e7, 1e8)   # squarefree strong 4-Fermat psp to base 2, in range 1e7..1e8

squarefree_strong_fermat_psp_each

k.squarefree_strong_fermat_psp_each(base, n, { ... })
k.squarefree_strong_fermat_psp_each(base, a, b, { ... })

Iterates over the squarefree strong Fermat pseudoprimes <= n or in the range a..b that have exactly k prime factors.

3.squarefree_strong_fermat_psp_each(2, 1e6, {|n| say n })          # iterate over squarefree strong 3-Fermat psp to base 2 <= 1e6
4.squarefree_strong_fermat_psp_each(2, 1e7, 1e8, {|n| say n })     # iterate over squarefree strong 4-Fermat psp to base 2 in range 1e7..1e8

Aliases: each_squarefree_strong_fermat_psp

NUMBER THEORY: FACTORIZATION, GCD & LCM

chebyshev_factor

n.chebyshev_factor
n.chebyshev_factor(B)
n.chebyshev_factor(B,x)

Integer factorization method, based on Chebyshev polynomials, which have the following nesting property:

T_{m n}(x) = T_m(T_n(x))

which are efficiently computed using the Lucas V function V_n(P,Q):

T_n(x) = (1/2) * V_n(2x, 1)

This method is particularly effective for numbers that have a prime factor p such that p-1 or p+1 is B-smooth.

say chebyshev_factor(1124075136413 * 3556516507813, 4000, 3)

The product of the factors will give back n. However, some factors may be composite.

consecutive_lcm

consecutive_lcm(n)

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

Aliases: consecutive_integer_lcm

cop_factor

n.cop_factor(tries=n.ilog2)

Congruence of Powers (CoP) factorization method, trying to find algebraic factors of n.

say cop_factor((5**48 + 1)*(3**120 + 1))

An additional argument can be given to limit the number of iterations.

The product of the factors will give back n. However, some factors may be composite.

cyclotomic_factor

n.cyclotomic_factor
n.cyclotomic_factor(bases...)

Factorization method, based on cyclotomic polynomials, trying to find algebraic factors of n.

say cyclotomic_factor(2**120 + 1)
say cyclotomic_factor((10**258 - 1)/9 - 10**(258/2) - 1, 10)

An optional list of bases can be given to restrict the search only to the given bases.

The product of the factors will give back n. However, some factors may be composite.

dop_factor

n.dop_factor(tries=n.ilog2)

Difference of Powers (DoP) factorization method, trying to find algebraic factors of n.

say dop_factor(2**120 + 1)

An additional argument can be given to limit the number of iterations.

The product of the factors will give back n. However, some factors may be composite.

ecm_factor

n.ecm_factor
n.ecm_factor(B1)
n.ecm_factor(B1, curves)

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

The product of the factors will give back n. However, some factors may be composite.

factor

factor(n)
factor(n, { ... })

Returns the prime factors of n (with repetition), sorted in ascending order.

For n <= 0, it returns an empty array.

say 180.factor          #=> [2, 2, 3, 3, 5]
say 2023.factor         #=> [7, 17, 17]
say factor(2**67 - 1)   #=> [193707721, 761838257287]
say factor(1)           #=> []

An optional block can specify custom factorization methods:

say factor(10**120 - 10**40, {|k| k.ecm_factor })
say factor(10**120 - 10**40, {|k| k.fermat_factor })

The block is expected to return an array with zero or more factors of k. These factors do not need to be prime. The block is applied recursively to each new composite factor until no further splitting is possible.

Aliases: factors

factor_exp

factor_exp(n)

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

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

Aliases: factors_exp

factor_map

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

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

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

factor_prod

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

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

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

Equivalent with:

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

Aliases: factors_prod

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: factors_sum

factor_upto

n.factor_upto(limit)

Returns an array of factors of n, using trial division, the Pollard Rho method and the Elliptic Curve method, based on heuristic bounds to find all the prime factors of n that are less than limit.

The product of the factors will give back n. However, the last factor may be composite.

NOTE: some factors may be larger than limit. Also, it's not guaranteed that a prime factor less than limit is always found, even when it exists.

fermat_factor

n.fermat_factor(k=1e4)

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

Works for odd composite nonpower numbers n that have two divisors close to sqrt(n).

The product of the factors will give back n. However, some factors may be composite.

fibonacci_factor

n.fibonacci_factor(upto = 2*n.ilog2)

Tries to find special factors of a Fibonacci-like number.

say fibonacci_factor(480.fib)
say fibonacci_factor(480.lucas)

The product of the factors will give back n. However, some factors may be composite.

Aliases: fib_factor

flt_factor

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

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

This method is particularly effective for numbers that have factors close to each other, or have a factor k for which 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]

The product of the factors will give back n. However, some factors may be composite.

is_coprime

is_coprime(a,b)

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

gcd

gcd(...)

Returns the greatest common divisor of a list of integers.

With one argument, it returns the absolute value of that argument. By convention, gcd() is 0.

say gcd(54, 24)    #=> 6
say gcd(-54, 24)   #=> 6
say gcd(-42)       #=> 42
say gcd()          #=> 0

gcd_factors

n.gcd_factors([a, b, ...])

Given a positive integer and an array of integers, it tries to find nontrivial factors of n, checking each gcd(n, array[0]), gcd(n, array[1]), etc.

var n = 43*43*97*503
var a = [19*43*97, 1, 13*41*43*101]
say gcd_factors(n, a)                #=> [43, 43, 97, 503]

The product of the factors gives back n. However, some factors may be composite.

gcdext

gcdext(a,b)

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

u*a + v*b = d.

The value of d is always nonnegative.

gcud

gcud(...)

Returns the greatest common unitary divisor of a list of integers.

germain_factor

germain_factor(n)

Tries to factorize n, using the Sophie Germain identity:

x^4 + 4y^4 = (x^2 - 2xy + 2y^2) * (x^2 + 2xy + 2y^2)

The product of the factors will give back n. However, some factors may be composite.

Aliases: sophie_germain_factor

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

gpf_sum

gpf_sum(n)
gpf_sum(a,b)

Returns the sum of largest prime factors of numbers from 1 to n, or in the range a..b. (OEIS A088822)

say 30.of { .gpf_sum }   #=> OEIS: A088822

holf_factor

n.holf_factor(tries=1e4)

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

The product of the factors will give back n. However, some factors may be composite.

lcm

lcm(...)

Returns the least common multiple of a list of integers.

With one argument, it returns the absolute value of that argument. By convention, lcm() is 1.

say lcm(6, 10)     #=> 30
say lcm(-6, 10)    #=> 30
say lcm(21)        #=> 21
say lcm()          #=> 1

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

lpf_sum

lpf_sum(n)
lpf_sum(a,b)

Returns the sum of least prime factors of numbers from 1 to n, or in the range a..b. (OEIS A088821)

say 30.of { .lpf_sum }     #=> OEIS: A088821

lucas_factor

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

A factorization method, using the modular Lucas sequences, which is effective in factoring Carmichael numbers, Fermat pseudoprimes, Lucas pseudoprimes and Lucas-Carmichael numbers.

say lucas_factor(2425361208749736840354501506901183117777758034612345610725789878400467)

The product of the factors will give back n. However, some factors may be composite.

Aliases: lucas_miller_factor

mbe_factor

n.mbe_factor(tries=10)

Tries to find a factor of n, by using the "Modular Binary Exponentiation" factorization method (randomized version).

This method is particularly effective for numbers that contain a prime factor p such that p-1 is sufficiently smooth.

The running time of the method, is: O(tries * log(n)^2).

Example:

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

The product of the factors will give back n. However, some factors may be composite.

miller_factor

n.miller_factor(tries=100)

A factorization method, based on the Miller-Rabin primality test. Effective in factoring Carmichael numbers and Fermat pseudoprimes.

It returns an array with the factors of n. The product of the factors will give back n. However, some factors may be composite.

Aliases: miller_rabin_factor

pell_factor

n.pell_factor(tries=1e4)

Pell factorization method, trying to find a factor of n.

This method is particularly effective for numbers that have factors close to sqrt(n).

say pell_factor(10**120 - 10**40)

The product of the factors will give back n. However, some factors may be composite.

phi_finder_factor

n.phi_finder_factor(tries=1e4)

Tries to find a factor of n, using the Phi-finder factorization method due to Kyle Kloster (2010).

This method is particularly effective for semiprimes n = p*q such that p and q are relatively close to each other.

Example:

say phi_finder_factor(622882096110539)          #=> [23099599, 26965061]
say phi_finder_factor(132750061135361, 1e5)     #=> [9369673, 14168057]

The product of the factors will give back n. However, some factors may be composite.

pm1_factor

n.pm1_factor
n.pm1_factor(B)

Pollard p-1 factorization method.

The product of the factors will give back n. However, some factors may be composite.

Aliases: pminus1_factor

pp1_factor

n.pp1_factor(B)

Williams' p+1 factorization method.

The product of the factors will give back n. However, some factors may be composite.

Aliases: pplus1_factor

prime_divisors

prime_divisors(n)

Returns the unique prime factors of n.

prime_signature

n.prime_signature

Given a non-negative integer n, returns the prime signature of n.

This is the exponents of the prime-power factorization of n, sorted in descending order.

say prime_signature(324)                                      #=> [4, 2]  (324 = 2^2 * 3^4)
say 30.of { .prime_signature.prod_kv{|k,v| prime(k+1)**v } }  #=> OEIS: A046523

prime_signature_inverse

prime_signature_inverse(A, B, prime_signature)

Returns an array of numbers in range [A, B] that have a given prime signature.

say prime_signature_inverse(1, 1e3, [1])        #=> prime numbers
say prime_signature_inverse(1, 1e3, [2])        #=> squares of numbers
say prime_signature_inverse(1, 1e3, [2,1])      #=> squares of primes times another prime
say prime_signature_inverse(1, 1e3, [2,1,1])    #=> numbers of the form p^2*q*r

Aliases: prime_signature_numbers.

prime_signature_inverse_len

prime_signature_inverse_len(A, B, prime_signature)

Count numbers in range [A, B] that have a given prime signature.

say 10.of {|n| prime_signature_inverse_len(1, 10**n, [1]) }      # prime number count
say 10.of {|n| prime_signature_inverse_len(1, 10**n, [1,1]) }    # A036351
say 10.of {|n| prime_signature_inverse_len(1, 10**n, [1,1,1]) }  # A215218

Aliases: prime_signature_count.

qs_factor

qs_factor(n)

Pomerance's Quadratic Sieve factorization method (SIMPQS).

The product of the factors will give back n. However, some factors may be composite.

rho_brent_factor

n.rho_brent_factor
n.rho_brent_factor(tries)

Pollard-Brent rho factorization method.

The product of the factors will give back n. However, some factors may be composite.

rho_factor

n.rho_factor
n.rho_factor(tries)

Pollard rho factorization method.

The product of the factors will give back n. However, some factors may be composite.

special_factor

n.special_factor(tries=1)

Tries to find special factors of n, using various special-factorization methods, including trial division, p-1 method, p+1 method, HOLF method, Fermat method, Pell method, Miller-Rabin method and the Lucas-Miller method.

An additional optional argument can be given to increase the number of tries by the given factor. For example, tries = 2 will double the number of tries.

The method returns an array with the factors of n:

say special_factor((3**120 + 1) * (5**240 - 1))

The product of the factors will give back n. However, some factors may be composite.

Aliases: special_factors

squfof_factor

n.squfof_factor(tries=1e4)

Shanks' SQUFOF method.

The product of the factors will give back n. However, some factors may be composite.

trial_factor

n.trial_factor(limit)

Trial division.

The product of the factors will give back n. However, some factors may be composite.

valuation

n.valuation(k)

Returns the number of times n is divisible by k.

say valuation(2**32, 4)     # prints: 16

NUMBER THEORY: DIVISORS & DIVISOR SUMS

abundancy

abundancy(n)

Returns the abundancy index of n, defined as:

sigma(n)/n

Aliases: abundancy_index

aliquot

aliquot(n)

Returns the sum of divisors of n that are less than n.

Equivalent with:

sigma(n) - n

antidivisor_count

antidivisor_count(n)

Returns the number of antidivisors of n.

say antidivisor_count(13)   #=> 4
say antidivisor_count(27)   #=> 5

antidivisor_sum

antidivisor_sum(n, k=1)

Returns the sum of anti-divisors of n.

When k is specified, it computes the sum of d**k for all antidivisors d of n.

Aliases: antidivisor_sigma

antidivisors

antidivisors(n)

Returns an array with the antidivisors of n sorted from 1..n.

Antidivisors of n are numbers that do not divide n by the largest possible margin.

say antidivisors(24)    #=> [7, 16]
say antidivisors(128)   #=> [3, 5, 15, 17, 51, 85]

bdivisors

bdivisors(n)

Returns an array with the bi-unitary divisors of n.

say 48.bdivisors    #=> [1, 2, 3, 6, 8, 16, 24, 48]

These are divisors d of n that satisfy:

gcud(n/d, d) = 1

Aliases: biudivisors, bi_unitary_divisors

bsigma

bsigma(n, k=1)

Returns the sum of the bi-unitary divisors of n, each divisor raised to the k-th power. (OEIS A188999)

say 20.of { .bsigma }    #=> OEIS: A188999

Aliases: biusigma

bsigma0

bsigma0(n)

Returns the count of the bi-unitary divisors of n. (OEIS A286324)

say 20.of { .bsigma0 }   #=> OEIS: A286324

Aliases: biusigma0

bsigma0_sum

bsigma0_sum(n)

Sum of number of bi-unitary divisors of k for k=1..n.

say 30.of { .bsigma0_sum }    #=> OEIS: A306069

bsigma_sum

bsigma_sum(n,k=1)

Sum of the bsigma(x,k) function, for x = 1..n.

say 30.of { .bsigma_sum(0) }   # OEIS: A306069
say 30.of { .bsigma_sum    }   # OEIS: A307159

cube_divisors

n.cube_divisors

Returns the cube divisors of n.

say 10!.cube_divisors    #=> [1, 8, 27, 64, 216, 1728]

Equivalent with:

n.divisors.grep { .is_cube }

cube_sigma

n.cube_sigma(k=1)

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

say 10!.cube_sigma       #=> 2044
say 10!.cube_sigma(2)    #=> 3037530

Equivalent with:

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

cube_sigma0

n.cube_sigma0

Returns the count of cube divisors of n.

cube_udivisors

n.cube_udivisors

Returns the unitary cube divisors of n, such that each divisor d is a cube and gcd(n/d, d) = 1.

say 15!.cube_udivisors   #=> [1, 125, 729, 91125]

cube_usigma

n.cube_usigma(k=1)

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

say cube_usigma(15!)     #=> 91980
say cube_usigma(15!, 2)  #=> 8304312692

Equivalent with:

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

cube_usigma0

n.cube_usigma0

Returns the count of unitary cube divisors of n.

cubefree_divisors

n.cubefree_divisors

Returns an array with the cubefree divisors of n.

say cubefree_divisors(120)  #=> [1, 2, 3, 4, 5, 6, 10, 12, 15, 20, 30, 60]

cubefree_sigma

n.cubefree_sigma(k=1)

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

say cubefree_sigma(5040)      #=> 4368
say cubefree_sigma(5040, 2)   #=> 2484300

cubefree_sigma0

n.cubefree_sigma0

Returns the count of cubefree divisors of n.

cubefree_udivisors

n.cubefree_udivisors

Returns the unitary cubefree divisors of n.

say cubefree_udivisors(5040)   #=> [1, 5, 7, 9, 35, 45, 63, 315]

cubefree_usigma

n.cubefree_usigma(k=1)

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

say 5040.cubefree_usigma     #=> 480
say 5040.cubefree_usigma(2)  #=> 106600

cubefree_usigma0

n.cubefree_usigma0

Returns the number of unitary cubefree divisors of n.

say cubefree_usigma0(5040)    #=> 8

cubefull_divisors

cubefull_divisors(n)

Return an array with the cubefull (3-powerful) divisors of n.

cubefull_sigma

cubefull_sigma(n, j=1)

Sum of the cubefull (3-powerful) divisors of n, each divisor raised to the power j.

cubefull_sigma0

cubefull_sigma0(n)

Count of the cubefull (3-powerful) divisors of n.

cubefull_udivisors

cubefull_udivisors(n)

Return an array with the unitary cubefull (3-powerful) divisors of n.

cubefull_usigma

cubefull_usigma(n, j=1)

Sum of the unitary cubefull (3-powerful) divisors of n, each divisor raised to the power j.

cubefull_usigma0

cubefull_usigma0(n)

Count of the unitary cubefull (3-powerful) divisors of n.

divisors

n.divisors(k=n)

Returns an array with the positive divisors of n that are less than or equal to k.

When k is omitted, all positive divisors of n are returned in ascending order. Non-positive inputs produce an empty array.

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]
say 0.divisors          #=> []

divisors_each

n.divisors_each {|d| ... }

Lazily iterate over the divisors d of n. (not in sorted order)

24.divisors_each {|d| say d }

Aliases: divisor_each

divisors_map

n.divisors_map {|d| ... }

Maps the divisors of n to the given block.

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

Aliases: divisor_map

divisors_prod

n.divisors_prod {|d| ... }

Product of the mapping of the positive divisors of n.

Equivalent with:

n.divisors_map {|d| ... }.prod

Aliases: divisor_prod

divisors_sum

n.divisors_sum {|d| ... }

Sum of the mapping of the positive divisors of n.

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

Equivalent with:

n.divisors_map {|d| ... }.sum

Aliases: divisor_sum

edivisors

edivisors(n)

Returns an array with the exponential divisors (or e-divisors) of n.

say edivisors(5040)     #=> [210, 420, 630, 1260, 1680, 5040]

Aliases: exponential_divisors

esigma

esigma(n, k=1)

Returns the sum of the exponential divisors (or e-divisors) of n, each divisor raised to k-th power. (OEIS A051377)

say 20.of { .esigma }       #=> OEIS: A051377

esigma0

esigma0(n)

Returns the count of the exponential divisors (or e-divisors) of n. (OEIS A049419)

say 20.of { .esigma0 }      #=> OEIS: A049419

idivisors

n.idivisors

Returns an array with the infinitary divisors (or i-divisors) of n.

say 96.idivisors        #=> [1, 2, 3, 6, 16, 32, 48, 96]

Aliases: infinitary_divisors

inverse_count

inverse_count(target, count_func, bsearch_method='bsearch_min')

Finds the smallest integer m >= 0 such that count_func(m) >= target.

This method implements the following search strategy:

1. Exponential Search: Rapidly doubles the search range to find an upper bound hi such that count_func(hi) >= target.

2. Binary Search: Refines the range between lo and hi using the specified bsearch_method (defaults to 'bsearch_min').

The complexity is logarithmic relative to the value of the result, making it suitable for inverting monotonic functions like prime counts or smooth counts. Returns NaN if target is negative.

say inverse_count(25, { .prime_count })             #=> 97
say inverse_count(100, { 10.smooth_count(_) })      #=> 450

isigma

isigma(n, k=1)

Returns the sum of the infinitary divisors of n, each divisor raised to the k-th power. (OEIS A049417)

say 20.of { .isigma }       #=> OEIS: A049417

isigma0

isigma0(n)

Returns the count of the infinitary divisors of n. (OEIS A037445)

say 20.of { .isigma0 }      #=> OEIS: A037445

isigma0_sum

isigma0_sum(n)

Sum of number of infinitary divisors of k for k=1..n.

say 30.of { .isigma0_sum }    #=> OEIS: A327573

isigma_sum

isigma_sum(n, k=1)

Sum of the isigma(x,k) function, for x=1..n.

say 30.of { .isigma_sum(0) }   # OEIS: A327573
say 30.of { .isigma_sum }      # OEIS: A327566

nbdivisors

nbdivisors(n)

Returns an array with the non-bi-unitary divisors of n.

say 720.nbdivisors   #=> [3, 4, 6, 12, 15, 20, 24, 30, 36, 48, 60, 120, 180, 240]

Equivalent to the set difference between the divisors and the bi-unitary-divisors:

divisors(n) - bdivisors(n) == nbdivisors(n)

These are divisors d of n that satisfy:

gcud(n/d, d) != 1

Aliases: non_bi_unitary_divisors, nbiudivisors

nbsigma

nbsigma(n, k=1)

Returns the sum of the non-bi-unitary divisors of n, each divisor raised to the k-th power. (OEIS A319072)

say 30.of { .nbsigma }      #=> OEIS: A319072

nbsigma0

nbsigma0(n)

Returns the count of the non-bi-unitary divisors of n.

nbsigma0_sum

nbsigma0_sum(n)

Sum of number of non-bi-unitary divisors of k for k=1..n.

say 30.of { .nbsigma0_sum }    # Partial sums of OEIS A390957

nbsigma_sum

nbsigma_sum(n,k=1)

Sum of the nbsigma(x,k) function, for x = 1..n.

say 30.of { .nbsigma_sum(0) }   # Partial sums of OEIS A390957
say 30.of { .nbsigma_sum    }   # Partial sums of OEIS A319072

nedivisors

nedivisors(n)

Returns an array with the nonexponential divisors of n.

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

Equivalent to the set difference between the divisors and the exponential divisors:

divisors(n) - edivisors(n) == nedivisors(n)

Aliases: non_exponential_divisors

nesigma

nesigma(n, k=1)

Returns the sum of the nonexponential divisors of n, each divisor raised to the k-th power. (OEIS A160135)

say 30.of { .nesigma }      #=> OEIS: A160135

nesigma0

nesigma0(n)

Returns the count of the nonexponential divisors of n. (OEIS A160097)

say 30.of { .nesigma0 }     #=> OEIS: A160097

next_tau_inverse

next_tau_inverse(n, k)

Returns the next number, greater than n, that has exactly k divisors.

nidivisors

nidivisors(n)

Returns an array with the noninfinitary divisors of n.

say 240.nidivisors   #=> [2, 4, 6, 8, 10, 12, 20, 24, 30, 40, 60, 120]

Equivalent to the set difference between the divisors and the infinitary divisors:

divisors(n) - idivisors(n) == nidivisors(n)

Aliases: non_infinitary_divisors

nisigma

nisigma(n, k=1)

Returns the sum of the noninfinitary divisors of n, each divisor raised to the k-th power. (OEIS A348271)

say 30.of { .nisigma }      #=> OEIS: A348271

nisigma0

nisigma0(n)

Returns the count of the noninfinitary divisors of n. (OEIS A348341)

say 30.of { .nisigma0 }     #=> OEIS: A348341

nisigma0_sum

nisigma0_sum(n)

Sum of number of noninfinitary divisors of k for k=1..n.

say 30.of { .nisigma0_sum }    # Partial sums of OEIS A348341

nisigma_sum

nisigma_sum(n,k=1)

Sum of the nisigma(x,k) function, for x = 1..n.

say 30.of { .nisigma_sum(0) }   # Partial sums of OEIS A348341
say 30.of { .nisigma_sum    }   # Partial sums of OEIS A348271

nudivisors

nudivisors(n)

Returns an array with the non-unitary divisors of n.

say 120.nudivisors   #=> [2, 4, 6, 10, 12, 20, 30, 60]

Equivalent to the set difference between the divisors and the unitary-divisors:

divisors(n) - udivisors(n) == nudivisors(n)

These are divisors d of n that satisfy:

gcd(n/d, d) != 1

Aliases: non_unitary_divisors

nusigma

nusigma(n, k=1)

Returns the sum of the nonunitary divisors of n, each divisor raised to the k-th power. (OEIS A048146)

say 30.of { .nusigma }     #=> OEIS: A048146

nusigma0

nusigma0(n)

Returns the count of the nonunitary divisors of n. (OEIS A048105)

say 30.of { .nusigma0 }    #=> OEIS: A048105

nusigma0_sum

nusigma0_sum(n)

Sum of number of nonunitary divisors of k for k=1..n.

say 30.of { .nusigma0_sum }    # Partial sums of OEIS A048105

nusigma_sum

nusigma_sum(n,k=1)

Sum of the nusigma(x,k) function, for x = 1..n.

say 30.of { .nusigma_sum(0) }   # Partial sums of OEIS A048105
say 30.of { .nusigma_sum    }   # Partial sums of OEIS A048146

power_divisors

k.power_divisors(n)

Return an array with the k-th power divisors of n.

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

Equivalent with:

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

power_sigma

k.power_sigma(n, j=1)

Returns the sum of the k-th power divisors of n, each divisor raised to the power j.

say 30.of { 2.power_sigma(_) }        # OEIS: A035316
say 30.of { 3.power_sigma(_) }        # OEIS: A113061

Equivalent with (but much faster):

n.divisors.grep { .is_power(k) }.sum {|d| d**j }

power_sigma0

k.power_sigma0(n)

Returns the number of the k-th power divisors of n.

say 30.of { 2.power_sigma0(_) }       # OEIS: A046951

power_udivisors

k.power_udivisors(n)

Returns an array with the k-th power unitary divisors of n.

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

Equivalent with:

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

Aliases: power_unitary_divisors, unitary_power_divisors

power_usigma

k.power_usigma(n, j=1)

Returns the sum of the k-th power unitary divisors of n, each divisor raised to the power j.

Equivalent with (but much faster):

n.udivisors.grep { .is_power(k) }.sum {|d| d**j }

power_usigma0

k.power_usigma0(n)

Returns the number of the k-th power unitary divisors of n.

say 30.of { 2.power_usigma0(_) }       # OEIS: A056624

powerfree_divisors

k.powerfree_divisors(n)

Returns an array with the k-powerfree divisors of n.

say 2.powerfree_divisors(5040)      # squarefree divisors of 5040
say 3.powerfree_divisors(5040)      # cubefree divisors of 5040

powerfree_sigma

k.powerfree_sigma(n, j=1)

Returns the sum of the k-powerfree divisors of n, each divisor raised to the power j.

say 30.of { 2.powerfree_sigma(_) }        # OEIS: A048250
say 30.of { 3.powerfree_sigma(_) }        # OEIS: A073185

Equivalent with (but much faster):

n.divisors.grep { .is_powerfree(k) }.sum {|d| d**j }

powerfree_sigma0

k.powerfree_sigma0(n)

Returns the number of the k-powerfree divisors of n.

say 30.of { 3.powerfree_sigma0(_) }       # OEIS: A073184

powerfree_udivisors

k.powerfree_udivisors(n)

Returns an array with the k-powerfree unitary divisors of n.

say 2.powerfree_udivisors(5040)      # squarefree unitary divisors of 5040
say 3.powerfree_udivisors(5040)      # cubefree unitary divisors of 5040

Equivalent with (but much faster):

n.udivisors.grep { .is_powerfree(k) }
k.powerfree_divisors(n).grep {|d| gcd(n/d, d) == 1 }

powerfree_usigma

k.powerfree_usigma(n, j=1)

Returns the sum of the k-powerfree unitary divisors of n, each divisor raised to the power j.

say 20.of { 2.powerfree_usigma(_) }        # OEIS: A092261
say 20.of { 2.powerfree_usigma(_, 2) }     # OEIS: A091306

Equivalent with (but much faster):

n.udivisors.grep { .is_powerfree(k) }.sum {|d| d**j }

powerfree_usigma0

k.powerfree_usigma0(n)

Returns the number of the k-powerfree unitary divisors of n.

say 20.of { 2.powerfree_usigma0(_) }       # OEIS: A056671

powerful_divisors

k.powerful_divisors(n)

Returns an array with the k-powerful divisors of n.

powerful_sigma

k.powerful_sigma(n, j=1)

Sum of the k-powerful divisors of n, each divisor raised to the power j.

say 2.powerful_sigma(n)     # OEIS: A183097
say 3.powerful_sigma(n)     # OEIS: A385005

powerful_sigma0

k.powerful_sigma0(n)

Count the number of k-powerful divisors of n.

powerful_udivisors

k.powerful_udivisors(n)

Returns an array with the unitary k-powerful divisors of n.

Aliases: unitary_powerful_divisors, powerful_unitary_divisors

powerful_usigma

k.powerful_usigma(n, j=1)

Sum of the unitary k-powerful divisors of n, each divisor raised to the power j.

say 30.of{|n| 2.powerful_usigma(n) }    # OEIS: A360720
say 30.of{|n| 3.powerful_usigma(n) }    # OEIS: A385048

powerful_usigma0

k.powerful_usigma0(n)

Count the number of unitary k-powerful divisors of n.

pp_divisors

pp_divisors(n)

Returns an array with the perfect power divisors of n.

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

Equivalent with:

n.divisors.grep { .is_perfect_power }

Aliases: perfect_power_divisors

pp_udivisors

pp_udivisors(n)

Returns an array with the perfect power unitary divisors of n, such that each divisor d is a perfect power and gcd(n/d, d) = 1.

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

Equivalent with:

n.udivisors.grep { .is_perfect_power }

Aliases: perfect_power_udivisors, perfect_power_unitary_divisors

prev_tau_inverse

prev_tau_inverse(n, k)

Returns the previous number, smaller than n, that has exactly k divisors.

prime_power_divisors

n.prime_power_divisors

Returns the prime power divisors of n.

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

Equivalent with:

n.divisors.grep{.is_prime_power}

prime_power_sigma

n.prime_power_sigma(k=1)

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

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

Equivalent with:

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

prime_power_sigma0

prime_power_sigma0(n)

Returns the count of prime power divisors of n.

Equivalent with:

bigomega(n)

prime_power_udivisors

prime_power_udivisors(n)

Returns the unitary prime power divisors of n.

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: prime_power_unitary_divisors, unitary_prime_power_divisors

prime_power_usigma

Number.prime_power_usigma(k=1)

Returns the sum of the unitary prime power divisors of n, each divisor raised to power k.

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

prime_power_usigma0

prime_power_usigma0(n)

Returns the count of unitary prime power divisors of n.

Equivalent with:

omega(n)

prime_sigma

n.prime_sigma(k=1)

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

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

prime_sigma0

prime_sigma0(n)

Returns the count of prime divisors of n.

Equivalent with:

omega(n)

prime_udivisors

n.prime_udivisors

Returns the unique unitary prime factors of n.

Aliases: prime_unitary_divisors, unitary_prime_divisors

prime_usigma

n.prime_usigma(k=1)

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

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

prime_usigma0

prime_usigma0(n)

The number of unique unitary prime divisors of n.

proper_divisors

n.proper_divisors

Return an array with the divisors of n less than n.

say proper_divisors(6)  #=> [1, 2, 3]

proper_sigma0

n.proper_sigma0

Returns the number of proper divisors of n.

say proper_sigma0(6)    #=> 3

Aliases: proper_divisor_count

sigma

sigma(n,k=1)

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

The default k = 1 gives the usual sum-of-divisors function. When k = 0, it is equivalent to sigma0(n). Negative values of k return the reciprocal divisor-power sum as a rational value.

say sigma(12)                #=> 28
say sigma(12, 0)             #=> 6
say sigma(12, 2)             #=> 210
say sigma(12, -1).as_frac    #=> 7/3

Aliases: σ

sigma0

sigma0(n)

Returns the count of positive divisors of n.

Aliases: d, tau

sigma0_sum

sigma0_sum(n)

Sum of the sigma_0(k) function, for k = 1..n.

say 30.of { .sigma0_sum }   #=> OEIS: A006218

Aliases: tau_sum

sigma_inverse

n.sigma_inverse(k=1)

Returns an array with all the solutions x to the sigma sum of divisors function: sigma(x,k) = n.

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

Also works with arbitrary large integers:

say sigma_inverse(9325257382230393314439814176)

Aliases: inverse_sigma

sigma_inverse_len

n.sigma_inverse_len(k=1)

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

Equivalent to n.sigma_inverse(k).len, but much faster.

Aliases: inverse_sigma_len

sigma_inverse_max

n.sigma_inverse_max(k=1)

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

Equivalent to n.sigma_inverse(k).max, but much faster.

Returns nil if there are no solutions.

Aliases: inverse_sigma_max

sigma_inverse_min

n.sigma_inverse_min(k=1)

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

Equivalent to n.sigma_inverse(k).min, but much faster.

Returns nil if there are no solutions.

Aliases: inverse_sigma_min

sigma_sum

sigma_sum(n, k=1, j=0)

Returns the sum of m^j * sigma_k(m), for m=1..n, computed in O(sqrt(n)) steps.

say sigma_sum(100, 0)      #=> 482
say sigma_sum(100, 1)      #=> 8299
say sigma_sum(100, 2)      #=> 407819
say sigma_sum(100, 2, 3)   #=> 207528104123

square_divisors

square_divisors(n)

Returns the square divisors of n.

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

Equivalent with:

n.divisors.grep { .is_square }

square_sigma

n.square_sigma(k=1)

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

say 5040.square_sigma       #=> 210
say 5040.square_sigma(2)    #=> 22386

Equivalent with:

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

square_sigma0

square_sigma0(n)

Returns the count of square divisors of n.

square_udivisors

square_udivisors(n)

Returns the unitary square divisors of n, such that each divisor d is a square and gcd(n/d, d) = 1.

say 5040.square_udivisors   #=> [1, 9, 16, 144]

square_usigma

n.square_usigma(k=1)

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

say square_usigma(5040)     #=> 170
say square_usigma(5040, 2)  #=> 21074

Equivalent with:

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

square_usigma0

square_usigma0(n)

Returns the count of unitary square divisors of n.

squarefree_divisors

squarefree_divisors(n)

Returns an array with the squarefree divisors of n.

say squarefree_divisors(120)  #=> [1, 2, 3, 5, 6, 10, 15, 30]

squarefree_sigma

squarefree_sigma(n,k=1)

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

say squarefree_sigma(5040)      #=> 576
say squarefree_sigma(5040, 2)   #=> 65000

squarefree_udivisors

squarefree_udivisors(n)

Returns the unitary squarefree divisors of n.

say squarefree_udivisors(5040)   #=> [1, 5, 7, 35]

squarefree_usigma

n.squarefree_usigma(k=1)

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

say 5040.squarefree_usigma     #=> 48
say 5040.squarefree_usigma(2)  #=> 1300

squarefree_usigma0

squarefree_usigma0(n)

Returns the number of unitary squarefree divisors of n.

say squarefree_usigma0(5040)    #=> 4

squarefull_divisors

squarefull_divisors(n)

Return an array with the squarefull (2-powerful) divisors of n.

squarefull_sigma

squarefull_sigma(n, j=1)

Sum of the squarefull (2-powerful) divisors of n, each divisor raised to the power j.

squarefull_sigma0

squarefull_sigma0(n)

Count of the squarefull (2-powerful) divisors of n.

squarefull_udivisors

squarefull_udivisors(n)

Return an array with the unitary squarefull (2-powerful) divisors of n.

squarefull_usigma

squarefull_usigma(n, j=1)

Sum of the unitary squarefull (2-powerful) divisors of n, each divisor raised to the power j.

squarefull_usigma0

squarefull_usigma0(n)

Count of the unitary squarefull (2-powerful) divisors of n.

tau

Num.tau
tau(n)

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

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

Aliases: τ

tau_inverse

See tau_inverse_len below for the shared documentation of both inverse methods.

tau_inverse_len

tau_inverse(A, B, n)
tau_inverse_len(A, B, n)

Generate or count the numbers k in the range [A, B] such that tau(k) = n, where tau(k) is the number of positive divisors of k (i.e.: sigma_0(k)).

say tau_inverse(1, 100, 6)          # numbers with exactly 6 divisors in range [1, 100]
say tau_inverse(1, 2**40, 5040)     # numbers with exactly 5040 divisors in range [1, 2^40]

say tau_inverse_len(1, 100, 6)          # count the numbers with exactly 6 divisors in range [1, 100]
say tau_inverse_len(1, 2**40, 5040)     # count the numbers with exactly 5040 divisors in range [1, 2^40]

udivisors

udivisors(n)

Returns an array with the unitary divisors of n.

say 120.udivisors   #=> [1, 3, 5, 8, 15, 24, 40, 120]

These are divisors d of n that satisfy:

gcd(n/d, d) = 1

Aliases: unitary_divisors

usigma

usigma(n, k=1)

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

usigma0

usigma0(n)

Returns the count of unitary divisor of n.

Equivalently, the count of squarefree divisors of n.

say usigma0(5040)   #=> 16

Aliases: squarefree_sigma0

usigma0_sum

usigma0_sum(n)

Sum of number of unitary divisors of k for k=1..n.

say 30.of { .usigma0_sum }    #=> OEIS: A064608

usigma_inverse

n.usigma_inverse

Returns an array with all the solutions x to usigma(x) = n.

say usigma_inverse(120)      #=> [60, 87, 92, 95, 99]

Aliases: inverse_usigma

usigma_sum

usigma_sum(n, k=1)

Sum of the usigma(x,k) function, for x = 1..n.

say 30.of { .usigma_sum }       #=> OEIS: A064609
say 30.of { .usigma_sum(0) }    #=> OEIS: A064608

NUMBER THEORY: TOTIENTS & MULTIPLICATIVE FUNCTIONS

bigomega

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: Ω

bigomega_inverse

See bigomega_inverse_sum below for the shared documentation of all three inverse methods.

bigomega_inverse_len

See bigomega_inverse_sum below for the shared documentation of all three inverse methods.

bigomega_inverse_sum

bigomega_inverse(A, B, n)
bigomega_inverse_len(A, B, n)
bigomega_inverse_sum(A, B, n)

Generate, count, or sum the numbers k in the range [A, B] such that Ω(k) = n, where Ω(k) is the number of prime factors of k, counted with multiplicity.

say bigomega_inverse(1, 100, 3)
# means: (1..100).grep { .bigomega == 3 }

say bigomega_inverse_len(1, 100, 3)
# means: (1..100).count { .bigomega == 3 }

say bigomega_inverse_sum(1, 100, 3)
# means: (1..100).grep { .bigomega == 3 }.sum

bigomega_sum

bigomega_sum(n, k=0)

Returns the sum of the generalized Omega function Ω_k(j), for j=1..n.

say bigomega_sum(100, 0)      #=> 239
say bigomega_sum(100, 1)      #=> 2815
say bigomega_sum(100, 2)      #=> 66799

bphi

bphi(n)

The bi-unitary analog of Euler's totient function of n. (OEIS A116550)

cototient

cototient(n, k=1)

Generalized cototient function:

n^k - J_k(n)

where J_k(n) is the Jordan totient function.

For k=1, it equals: n - phi(n).

cototient_sum

cototient_sum(n, k=1)

Partial sums of the generalized cototient function, equivalent to:

faulhaber_sum(n, k) - totient_sum(n, k)

dconv

n.dconv(f,g)

The Dirichlet convolution of functions f and g, defined as:

Sum_{d|n} f(d) * g(n/d)

Example:

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

Aliases: dirichlet_convolution

derivative

derivative(n)

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

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

Aliases: arithmetic_derivative

dinv

n.dinv(f)

The Dirichlet inverse of a function f.

say 30.of { .dinv { .sigma } }      # Dirichlet inverse of sigma function (OEIS: A046692)

Aliases: dirichlet_inverse

dirichlet_sum

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

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

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

where F and G are the partial sums of f and g, respectively.

Computing Sum_{k=1..n} Sum_{d|k, d is squarefree} (n/d)^2:

var f = { .is_squarefree ? 1 : 0 }
var g = { .square }
var F = { .squarefree_count }
var G = { .faulhaber_sum(2) }

say 30.of {|n| dirichlet_sum(n, f, g, F, G) }   # OEIS: A321973

Computing Sum_{k=1..n} Sum_{d|k} n * mu(k/d):

var f = {|n| n }
var g = {|n| mu(n) }
var F = {|n| n*(n+1)/2 }
var G = {|n| mertens(n) }

say 30.of {|n| dirichlet_sum(n, f, g, F, G) }   # OEIS: A002088

Aliases: dirichlet_hyperbola, dsum

exp_bigomega

exp_bigomega(n,k)

Returns the result of k**bigomega(n).

exp_bigomega_sum

exp_bigomega_sum(n,k)

Partial sums of the exp_bigomega() function.

exp_bigomega_sum(n,  2)    # Sum_{k=1..n} 2^bigomega(k)
exp_bigomega_sum(n, -1)    # Sum_{k=1..n} (-1)^bigomega(k) = liouville_sum(n)

exp_mangoldt

exp_mangoldt(n)

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

say 30.of { .exp_mangoldt }     #=> OEIS: A014963

exp_mangoldt_sum

exp_mangoldt_sum(n)

Partial sums of the exponential of the Mangoldt function.

say 30.of { .exp_mangoldt_sum }     #=> OEIS: A072107
say exp_mangoldt_sum(1e10)          #=> 2220822442581396891
say exp_mangoldt_sum(1e11)          #=> 201467077843736717398

exp_omega

exp_omega(n,k)

Returns the result of k**omega(n).

exp_omega_sum

exp_omega_sum(n,k)

Partial sums of the exp_omega() function.

exp_omega_sum(n,  2)    # Sum_{k=1..n} 2^omega(k)
exp_omega_sum(n, -1)    # Sum_{k=1..n} (-1)^omega(k)

exp_squarefree_omega

exp_squarefree_omega(n,k)

Returns the result of k**omega(n) * mu(n)**2.

exp_squarefree_omega_sum

exp_squarefree_omega_sum(n,k)

Partial sums of the exp_squarefree_omega() function.

exp_squarefree_omega_sum(n,  3)     # Sum_{k=1..n} 3^omega(k) * mu(k)^2
exp_squarefree_omega_sum(n, -1)     # Sum_{k=1..n} (-1)^omega(k) * mu(k)^2 = mertens(n)

iphi

iphi(n, k=1)

Infinitary analog of Euler's phi function. (OEIS A091732)

iphi_sum

iphi_sum(n,k=1)

Partial sums of the infinitary totient function iphi(x,k) for x=1..n.

say 30.of { .iphi_sum }     #=> OEIS: A327572

is_abundant

n.is_abundant

Returns true if n is an abundant number, else false. (OEIS A005101)

say 20.by { .is_abundant }  #=> OEIS: A005101

An abundant number is smaller than the sum of its positive proper divisors.

is_amicable

is_amicable(n,m)

Returns true if the numbers n and m are "amicable", else false.

Amicable numbers are two different numbers so related that the sum of the proper divisors of each is equal to that of the other.

say is_amicable(220, 284)   #=> true

is_deficient

n.is_deficient

Returns true if n is a deficient number, else false.

A deficient number is greater than the sum of its positive proper divisors.

is_perfect

n.is_perfect

Returns true if n is a perfect number (i.e.: sigma(n) == 2*n).

is_pseudoperfect

var bool = n.is_pseudoperfect

Returns true if n is a pseudoperfect (semiperfect) number, and false otherwise.

A positive integer is pseudoperfect if it can be written as the sum of some (or all) of its proper divisors. Every perfect number is pseudoperfect, since the sum of all proper divisors equals the number itself.

say 30.by { .is_pseudoperfect }     # OEIS: A005835

Aliases: is_semiperfect

is_weird

is_weird(n)

Returns true if n is a weird number, and false otherwise.

A weird number is a positive integer that is abundant, but not pseudoperfect.

Equivalently, a weird number is an abundant number that cannot be expressed as the sum of a subset of its proper divisors.

say 10.by { .is_weird }         #=> OEIS: A006037

is_zumkeller

is_zumkeller(n)

Returns true if n is a Zumkeller number, and false otherwise.

A Zumkeller number is a positive integer whose positive divisors can be partitioned into two disjoint subsets having equal sums.

say 30.by { .is_zumkeller }     #=> OEIS: A083207

is_practical

n.is_practical

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

say 20.by { .is_practical }

is_primitive_abundant

n.is_primitive_abundant

Returns true if n is a primitive abundant number, else false. (OEIS A091191)

say 20.by { .is_primitive_abundant }    #=> OEIS: A091191

A primitive abundant number is an abundant number having no abundant proper divisor.

is_totient

n.is_totient

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

jordan_totient

jordan_totient(n,k)

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

lambda

lambda(n)

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

powmod(a, m, n) == 1

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

Alias: carmichael_lambda.

legendre_phi

n.legendre_phi(k)

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

Equivalent with:

prime(k+1).rough_count(n)

liouville

liouville(n)

The Liouville function.

Equivalent to:

(-1)**omega(n)

liouville_sum

liouville_sum(n)
liouville_sum(a,b)

Computes partial sums of the Liouville lambda function:

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

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

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

Example:

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

logarithmic_derivative

logarithmic_derivative(n)

Return the logarithmic derivative of n, defined as:

derivative(n)/n

make_coprime

n.make_coprime(k)

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

mangoldt

mangoldt(n)

The Mangoldt function. For the exponential values, see exp_mangoldt.

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)     # equivalent with: (1..100000 -> sum { .moebius })
say mertens(21, 123)    # equivalent with: (21..123 -> sum { .moebius })

mobius_range

mobius_range(a, b)

Returns an array with the Möbius values for the range a..b.

say moebius_range(7, 17) => [-1, 0, 0, 1, -1, 0, -1, 1, 1, 0, -1]

Aliases: moebius_range

moebius

mu(n)
moebius(n)

The Möbius function: μ(n). (OEIS A008683)

It is 1 when n = 1, 0 when n is divisible by the square of a prime, and (-1)^k when n is a product of k distinct primes.

say mu(1)       #=> 1
say mu(12)      #=> 0
say mu(30)      #=> -1
say mu(210)     #=> 1

Aliases: μ, mu, mobius, möbius

moebius_inverse

n.moebius_inverse(f)

The Möbius transform (inversion) of a function f, defined as:

Sum_{d|n} mu(d) * f(n/d)

Example:

say 30.of { .moebius_inverse { .psi } }   # Möbius transform of Dedekind psi function (OEIS: A063659)

Aliases: moebius_transform

nuphi

nuphi(n, k=1)

The generalized nonunitary totient function (also called nphi). (OEIS A254503)

Multiplicative with:

a(p, k) = p^k
a(p^e, k) = J_k(p^e) = (p^k - 1) * p^(k*(e-1)), for e > 1.

Example:

say 30.of { .nuphi }        # OEIS: A254503
say 30.of { .nuphi(2) }     # OEIS: A254520

nuphi_sum

nuphi_sum(n,k=1)

Partial sums of the nonunitary totient function nuphi(x,k) for x=1..n.

omega

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: ω

omega_inverse

See omega_inverse_sum below for the shared documentation of all three inverse methods.

omega_inverse_len

See omega_inverse_sum below for the shared documentation of all three inverse methods.

omega_inverse_sum

omega_inverse(A, B, n)
omega_inverse_len(A, B, n)
omega_inverse_sum(A, B, n)

Generate, count, or sum the numbers k in the range [A, B] such that ω(k) = n, where ω(k) is the number of distinct prime factors of k.

say omega_inverse(1, 100, 3)
# means: (1..100).grep { .omega == 3 }

say omega_inverse_len(1, 100, 3)
# means: (1..100).count { .omega == 3 }

say omega_inverse_sum(1, 100, 3)
# means: (1..100).grep { .omega == 3 }.sum

omega_sum

omega_sum(n, k=0)

Returns the sum of the generalized omega function ω_k(j), for j=1..n.

say omega_sum(100, 0)      #=> 171
say omega_sum(100, 1)      #=> 2298
say omega_sum(100, 2)      #=> 59962

phi

Num.phi
phi(n,k=1)

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

say Num.phi      #=> 1.618033988749...
say phi(12)      #=> 4

When called with an integer argument, phi(n) is the count of positive integers in the range 1..n that are coprime to n. When k > 1, it returns Jordan's totient function J_k(n).

say phi(12, 2)   #=> 96

Aliases: φ

phi_inverse

n.phi_inverse

Returns an array with all the solutions x to the Euler totient function: phi(x) = n.

say phi_inverse(24)     #=> [35, 39, 45, 52, 56, 70, 72, 78, 84, 90]

Aliases: inverse_phi, inverse_totient

phi_inverse_len

n.phi_inverse_len

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

Equivalent to n.phi_inverse.len, but much faster.

Aliases: inverse_phi_len, inverse_totient_len

phi_inverse_max

n.phi_inverse_max

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

Equivalent to n.phi_inverse.max, but much faster.

Returns nil if there are no solutions.

Aliases: inverse_phi_max, inverse_euler_phi_max

phi_inverse_min

n.phi_inverse_min

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

Equivalent to n.phi_inverse.min, but much faster.

Returns nil if there are no solutions.

Aliases: inverse_phi_min, inverse_euler_phi_min

pillai

pillai(n, k=1)

Pillai's arithmetical function: Sum_{k=1..n} gcd(k, n). (OEIS A018804)

Generalized as:

pillai(n,k) = Sum_{1 <= x_1, x_2, ..., x_k <= n} gcd(x_1, x_2, ..., x_k, n)^k

Which has the identity:

pillai(n,k) = Sum_{d|n} mu(n/d) * d^k * tau(d)

pillai_sum

pillai_sum(n,k=1)

Partial sums of the generalized Pillai arithmetical function pillai(n,k).

psi

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

Dedekind psi function. (OEIS A001615)

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: dedekind_psi

psi_inverse

n.psi_inverse

Returns an array with all the solutions x to Dedekind's psi function: psi(x) = n.

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

Aliases: inverse_psi, inverse_dedekind_psi

psi_inverse_len

n.psi_inverse_len

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

Equivalent to n.psi_inverse.len, but much faster.

Aliases: inverse_psi_len, inverse_dedekind_psi_len

psi_inverse_max

n.psi_inverse_max

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

Equivalent to n.psi_inverse.max, but much faster.

Returns nil if there are no solutions.

Aliases: inverse_psi_max, inverse_dedekind_psi_max

psi_inverse_min

n.psi_inverse_min

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

Equivalent to n.psi_inverse.min, but much faster.

Returns nil if there are no solutions.

Aliases: inverse_psi_min, inverse_dedekind_psi_min

psi_sum

n.psi_sum(k=1)

Returns the sum of the Dedekind psi function ψ_k(j), for j=1..n.

say psi_sum(100, 1)      #=> 7664
say psi_sum(100, 2)      #=> 400576

Aliases: dedekind_psi_sum

rad

rad(n)

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

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

Equivalent to:

n.factor.uniq.prod

Aliases: squarefree_kernel

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

sopf

sopf(n)

Sum of the distinct primes dividing n. (OEIS A008472)

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

sopf_sum

sopf_sum(n)

Partial sums of the sopf(k) function: sum of all distinct prime factors of k for k = 2..n. (OEIS A024924)

say 20.of { .sopf_sum }     #=> OEIS: A024924
say sopf_sum(1e9)           #=> 41986415795828892

sopfr

sopfr(n)

Sum of the primes dividing n (with repetition). (OEIS A001414)

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

sopfr_sum

sopfr_sum(n)

Partial sums of the sopfr(k) function: sum of primes dividing k (with repetition) for k = 2..n. (OEIS A025281)

say 20.of { .sopfr_sum }     #=> OEIS: A025281
say sopfr_sum(1e9)           #=> 41986419147434739

totient

phi(n,k=1)
totient(n,k=1)

Euler's totient function.

totient(n) counts the positive integers in the range 1..n that are coprime to n. For k > 1, it returns the generalized Jordan totient J_k(n).

say totient(12)      #=> 4
say totient(36)      #=> 12
say totient(12, 2)   #=> 96

Aliases: euler_phi, eulerphi, euler_totient

totient_range

totient_range(a, b)

Returns an array with the totient values for the range a..b.

say totient_range(7, 17)  #=> [6, 4, 6, 4, 10, 4, 12, 6, 8, 8, 16]

totient_sum

totient_sum(n, k=1)

Returns the sum of the Jordan totient function J_k(j), for j=1..n.

say totient_sum(100, 1)      #=> 3044
say totient_sum(100, 2)      #=> 280608

Aliases: phi_sum, euler_phi_sum, jordan_totient_sum

uphi

uphi(n,k=1)

The unitary totient function. (OEIS A047994)

uphi_inverse

n.uphi_inverse

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

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

Aliases: inverse_uphi

uphi_sum

uphi_sum(n,k=1)

Partial sums of the unitary totient function uphi(x,k) for x=1..n.

say 30.of { .uphi_sum }     #=> OEIS: A177754

NUMBER THEORY: SPECIAL INTEGER SETS

all_composite

all_composite(...)

Returns true if all the given values are positive composite numbers.

composite

n.composite

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

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

Aliases: nth_composite

composite_count

composite_count(n)
composite_count(a,b)

Returns the count of composite numbers <= n, or in the range a..b.

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

composite_count_lower

composite_count_lower(n)

Lower bound for composite_count(n).

composite_count_upper

composite_count_upper(n)

Upper bound for composite_count(n).

composite_lower

composite_lower(n)

Lower bound for the n-th composite number.

Aliases: nth_composite_lower

composite_sum

composite_sum(n)
composite_sum(a, b, k=1)

Returns the sum of composite numbers <= n, or in the range a..b.

say composite_sum(100)          # sum of composite numbers <= 100
say composite_sum(50, 100)      # sum of composite numbers in range 50..100

When k is specified, it returns the sum of composite numbers, each number raised to the k-th power:

say composite_sum(50, 100, 2)   # 50^2 + 51^2 + ... + 100^2

Aliases: composites_sum

composite_upper

composite_upper(n)

Upper bound for the n-th composite number.

Aliases: nth_composite_upper

composites

composites(n)
composites(a,b)

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

composites_each

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

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

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

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

Aliases: each_composite

core

core(n)

Squarefree part of n. (OEIS A007913)

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

Equivalent to PARI/GP core(n) function.

Aliases: squarefree_part

cube_part

cube_part(n)

Largest cube dividing n. (OEIS A008834)

cubefree

cubefree(n)
cubefree(a,b)

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

cubefree_count

cubefree_count(n)
cubefree_count(a,b)

Returns the count of cubefree numbers <= n, or in the range a..b.

cubefree_each

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

Iterates over the cubefree numbers <= n, or in the range a..b.

Aliases: each_cubefree

cubefree_kernel

cubefree_kernel(n)

Largest cubefree number dividing n. (OEIS A007948)

cubefree_part

cubefree_part(n)

Cubefree part of n. (OEIS A050985)

cubefree_sum

cubefree_sum(n)
cubefree_sum(a,b)

Returns the sum of cubefree numbers <= n, or in the range a..b.

cubefull

cubefull(n)
cubefull(a,b)

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

cubefull_count

cubefull_count(n)
cubefull_count(a,b)

Returns the count of cubefull numbers <= n, or in the range a..b.

cubefull_each

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

Iterates over the cubefull numbers <= n, or in the range a..b.

Aliases: each_cubefull

cubefull_part

cubefull_part(n)

Largest divisor of n that is a cubefull (3-powerful) number. (OEIS A360540)

cubefull_sum

cubefull_sum(n)
cubefull_sum(a,b)

Returns the sum of cubefull numbers <= n, or in the range a..b.

difference_of_squares

difference_of_squares(n)

Returns an array of [a,b] pairs with all the possible solutions to the equation: a^2 - b^2 = n.

say difference_of_squares(48)   #=> [[7, 1], [8, 4], [13, 11]]

Aliases: diff_of_squares

difference_of_cubes

difference_of_cubes(n)

Returns all integer solutions (x, y) satisfying:

x^3 - y^3 = n

The result is returned as an array of pairs, where each pair is itself an array [x, y]. Returns an empty array when no integer solutions exist or when n = 0.

The implementation is based on the factorization:

x^3 - y^3 = (x - y)(x^2 + xy + y^2)

and enumerates the divisors of |n|. For each divisor A = x - y, the corresponding quadratic equation is constructed, its discriminant is tested for being a perfect square, and any resulting integer solutions are recovered.

The algorithm supports both positive and negative values of n.

Examples:

say difference_of_cubes(19)      # [[-2, -3], [3, 2]]
say difference_of_cubes(1729)    # [[1, -12], [9, -10], [10, -9], [12, -1]]
say difference_of_cubes(-728)    # [[-12, -10], [-9, -1], [-8, 6], [-6, 8], [1, 9], [10, 12]]
say difference_of_cubes(30)      # []

Aliases: diff_of_cubes

fibonorial

fibonorial(n)

Returns the product of first n nonzero Fibonacci numbers F(1), ..., F(n).

is_between

n.is_between(min, max)

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

is_composite

n.is_composite

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

is_cube

n.is_cube

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

is_cubefree

n.is_cubefree

Returns true if n is not divisible by any perfect cube > 1.

is_cubefull

n.is_cubefull

Returns true if n is divisible by the cubes of all its prime factors.

is_cyclic

n.is_cyclic

Returns true when gcd(phi(n), n) = 1, where phi(n) is the Euler totient function. (OEIS A003277)

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

is_even

n.is_even

Returns true if n is an integer divisible by 2.

is_fundamental

n.is_fundamental

Returns true if n is a fundamental discriminant.

is_noncubefree

n.is_noncubefree

Returns true if n is a positive integer divisible by the cube of a prime.

is_nonpowerfree

n.is_nonpowerfree

Returns true if n is a positive integer divisible by the k-th power of a prime.

is_nonsquarefree

n.is_nonsquarefree

Returns true if n is a positive integer divisible by the square of a prime.

is_ntf

g.is_ntf(n)

Returns true if g is a nontrivial integer factor of n.

Equivalent with (assuming g and n are both integers):

g.divides(n) && g.is_between(2, n-1)

Aliases: is_nontrivial_factor

is_odd

n.is_odd

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

is_odd_composite

n.is_odd_composite

Returns true when n is an odd composite integer.

is_palindrome

n.is_palindrome(base=10)

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

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

Aliases: is_palindromic

is_pandigital

n.is_pandigital(base=10)

Returns true if n is pandigital in the specified base, and false otherwise.

A number is pandigital in base b if its representation contains every digit from 0 to b-1 at least once. Digits may appear more than once.

If no base is specified, base 10 is assumed.

say 1023456789.is_pandigital      #=> true
say 123456789.is_pandigital       #=> false
say 9876543210.is_pandigital      #=> true

Pandigitality can also be tested in other bases:

say 0b11010010.is_pandigital(2)            #=> true
say 0x123456789ABCDEF0.is_pandigital(16)   #=> true

is_smith

n.is_smith(base=10)

Returns true if n is a Smith number, and false otherwise. (OEIS A006753)

A Smith number is a composite integer whose digit sum equals the sum of the digit sums of its prime factors, counted with multiplicity.

By default, the comparison is performed in base 10. An alternative numeral base can be specified with the optional base argument.

say 30.by { .is_smith(2) }      #=> OEIS: A278909
say 30.by { .is_smith(10) }     #=> OEIS: A006753

is_pow

n.is_power
n.is_power(k)

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

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)

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

Aliases: is_power, is_perfect_power

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

is_powerfree

n.is_powerfree(k=2)

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

say 15.by { .is_powerfree(2) }      # squarefree numbers
say 15.by { .is_powerfree(3) }      # cubefree numbers

is_powerful

n.is_powerful(k=2)

Returns true when all the exponents in the prime-power factorization of n are greater than or equal to k.

is_prob_squarefree

n.is_prob_squarefree
n.is_prob_squarefree(limit)

Returns true if n is probably squarefree, checking if n is divisible by a square p^2 with p less than or equal to limit:

say is_prob_squarefree(2**512 - 1, 1e6)   # true   (probably squarefree)
say is_prob_squarefree(10**136 + 1, 1e3)  # false  (definitely not squarefree)

If n is less than limit^3 and the function returns true, then n is definitely squarefree.

If the limit parameter is omitted, multiple limits are tested internally, trying to find a square factor of n, up to limit = 10^7.

is_rough

n.is_rough(k)

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

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

is_smooth

n.is_smooth(k)

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

is_smooth_over_prod

n.is_smooth_over_prod(k)

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

is_square

n.is_square

Returns true if n is a perfect square integer.

Aliases: is_sqr, is_perfect_square

is_squarefree

n.is_squarefree

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

Aliases: is_square_free

is_squarefree_almost_prime

n.is_squarefree_almost_prime(k=2)

Returns true if n is a squarefree k-almost prime (i.e.: true iff n is the product of k distinct primes).

Equivalently, k-almost primes are numbers n that satisfy bigomega(n) == omega(n) == k.

say 20.by { .is_squarefree_almost_prime(1) }   # primes
say 20.by { .is_squarefree_almost_prime(2) }   # squarefree semiprimes
say 20.by { .is_squarefree_almost_prime(3) }   # sphenic numbers

is_squarefree_semiprime

is_squarefree_semiprime(n)

Returns true if n is a squarefree semiprime (i.e.: n = p*q, where p and q are two distinct primes).

is_squarefull

n.is_squarefull

Returns true if n is divisible by the squares of all its prime factors.

next_composite

n.next_composite

Given a nonnegative integer n, returns the next composite number greater than n.

next_composites

n.next_composites(start=4)

Returns an array with n consecutive composite numbers starting from start.

say 5.next_composites        #=> [4, 6, 8, 9, 10]
say 5.next_composites(50)    #=> [50, 51, 52, 54, 55]

next_cubefree

n.next_cubefree

Returns the next cubefree number greater than n.

next_cubefull

n.next_cubefull

Returns the next cubefull (or 2-full) number greater than n.

next_noncubefree

n.next_noncubefree

Returns the next noncubefree number greater than n.

next_nonpowerfree

n.next_nonpowerfree(k=2)

Returns the next k-nonpowerfree number greater than n.

next_nonsquarefree

n.next_nonsquarefree

Returns the next nonsquarefree number greater than n.

next_palindrome

n.next_palindrome(b=10)

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

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

next_perfect_power

next_perfect_power(n)
next_perfect_power(n,k)

Returns the next perfect power greater than n.

say next_perfect_power(1e6)     #=> 1002001

When k is given, it returns the next k-perfect-power greater than n.

say next_perfect_power(1e6, 3)  #=> 1030301

next_pow

n.next_pow(b)

Returns the next perfect power greater than n, with base b.

say 63.next_pow(2)      #=> 64
say 64.next_pow(2)      #=> 128

Equivalent with:

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

Aliases: next_power

next_powerfree

n.next_powerfree(k=2)

Returns the next k-powerfree number greater than n.

next_powerful

n.next_powerful(k=2)

Returns the next k-powerful (or k-full) number greater than n.

next_rough(n, k)

Returns the smallest k-rough integer strictly greater than n.

If n is zero, it returns 1.

next_smooth(n, k)

Returns the smallest k-smooth integer strictly greater than n.

If n is zero, it returns 1.

next_squarefree

n.next_squarefree

Returns the next squarefree number greater than n.

next_squarefull

n.next_squarefull

Returns the next squarefull (or 2-full) number greater than n.

noncubefree

noncubefree(n)
noncubefree(a,b)

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

noncubefree_count

noncubefree_count(n)
noncubefree_count(a,b)

Returns the count of noncubefree numbers <= n, or in the range a..b.

noncubefree_each

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

Iterates over the noncubefree numbers <= n, or in the range a..b.

Aliases: each_noncubefree

noncubefree_sum

noncubefree_sum(n)
noncubefree_sum(a,b)

Returns the sum of noncubefree numbers <= n, or in the range a..b.

nonpowerfree

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

Returns an array of numbers <= n, or in the range a..b, that are not k-powerfree.

say 2.nonpowerfree(100)        # numbers that are not squarefree <= 100
say 3.nonpowerfree(50, 100)    # numbers that are not cubefree in range 50..100

nonpowerfree_count

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

Returns the count of numbers <= n, or in the range a..b, that are not k-powerfree.

nonpowerfree_each

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

Iterates over the numbers <= n, or in the range a..b, that are not k-powerfree.

2.nonpowerfree_each(100, {|n| say n })         # iterate over nonsquarefree numbers <= 100
3.nonpowerfree_each(50, 100, {|n| say n })     # iterate over noncubefree numbers in range 50..100

Aliases: each_nonpowerfree

nonpowerfree_sum

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

Returns the sum of numbers <= n, or in the range a..b, that are not k-powerfree.

nonsquarefree

nonsquarefree(n)
nonsquarefree(a,b)

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

nonsquarefree_count

nonsquarefree_count(n)
nonsquarefree_count(a,b)

Returns the count of nonsquarefree numbers <= n, or in the range a..b.

nonsquarefree_each

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

Iterates over the nonsquarefree numbers <= n, or in the range a..b.

Aliases: each_nonsquarefree

nonsquarefree_sum

nonsquarefree_sum(n)
nonsquarefree_sum(a,b)

Returns the sum of nonsquarefree numbers <= n, or in the range a..b.

nth_cubefree

n.nth_cubefree

Returns the n-th cubefree number.

nth_cubefull

n.nth_cubefull

Returns the n-th cubefull (or 3-full) number.

nth_noncubefree

n.nth_noncubefree

Returns the n-th noncubefree number.

nth_nonpowerfree

n.nth_nonpowerfree(k)

Returns the n-th k-nonpowerfree number.

say nth_nonpowerfree(1e9, 2)   #=> 2550546152
say nth_nonpowerfree(1e9, 3)   #=> 5949100928
say nth_nonpowerfree(1e9, 4)   #=> 13147239114

nth_nonsquarefree

n.nth_nonsquarefree

Returns the n-th nonsquarefree number.

nth_perfect_power

nth_perfect_power(n)
nth_perfect_power(n,k)

Returns the n-th perfect power.

say nth_perfect_power(1e8)  #=> 9956760243243489
say nth_perfect_power(1e9)  #=> 997995681331086244

When k is given, it returns the n-th k-perfect-power.

nth_powerfree

nth_powerfree(n, k=2)

Returns the n-th k-powerfree number.

say nth_powerfree(1e14, 2)  #=> 164493406685659
say nth_powerfree(1e14, 3)  #=> 120205690315927
say nth_powerfree(1e14, 4)  #=> 108232323371116

nth_powerful

nth_powerful(n, k=2)

Returns the n-th k-powerful (or k-full) number.

say nth_powerful(1e4, 2)    #=> 23002083
say nth_powerful(1e4, 3)    #=> 16720797973

nth_rough

nth_rough(n, k)

Returns the n-th k-rough integer.

A number is considered k-rough (or k-resilient) if all of its prime factors are greater than or equal to k. This is the inverse of the rough_count function.

Returns NaN if n <= 0.

nth_smooth

nth_smooth(n, k)

Returns the n-th k-smooth integer (a Hamming number if k=3).

A number is considered k-smooth if all of its prime factors are less than or equal to k. This is the inverse of the smooth_count function.

Returns NaN if n <= 0 or k <= 1.

nth_squarefree

nth_squarefree(n)

Returns the n-th squarefree number. (OEIS A005117)

say nth_squarefree(1e14)    #=> 164493406685659

nth_squarefull

n.nth_squarefull

Returns the n-th squarefull (or 2-full) number.

nth_tau_inverse

nth_tau_inverse(n, k)

Returns the n-th number with exactly k divisors.

say {|n| nth_tau_inverse(n, n) }.map(1..20)    #=> OEIS: A073916

perfect_power

n.perfect_power

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

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

perfect_power_count

perfect_power_count(n)
perfect_power_count(n,k)

Returns the count of perfect powers <= n.

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

When k is provided, it returns the number of k-powers <= n.

perfect_power_sum

perfect_power_sum(n)
perfect_power_sum(n,k)

Returns the sum of perfect powers <= n.

say perfect_power_sum(10**6)      #=> 361590619
say perfect_power_sum(10**20)     #=> 333449517656248628022493884999

When k is provided, it returns the sum of k-powers <= n.

power_count

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

Returns the number of k-th power positive integers in <= n, or in the range a..b.

power_part

k.power_part(n)

Returns the largest k-th power dividing n.

say 2.power_part(n)     # largest square dividing n
say 3.power_part(n)     # largest cube dividing n

power_sum

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

Returns the sum of k-th power positive integers <= n, or in the range a..b.

powerfree

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

Returns an array with the k-powerfree <= n, or in the range a..b.

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

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 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

powerfree_each

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

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

2.powerfree_each(100, {|n| say n })          # iterate over squarefree numbers <= 100
3.powerfree_each(50, 100, {|n| say n })      # iterate over cubefree numbers in the range 50..100

Aliases: each_powerfree

powerfree_kernel

k.powerfree_kernel(n)

Returns the largest k-powerfree divisor of n.

say 2.powerfree_kernel(n)       # rad(n)
say 3.powerfree_kernel(n)       # largest cubefree number dividing n (OEIS: A007948)

powerfree_part

k.powerfree_part(n)

Returns the k-powerfree part of n.

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

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 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

powerfree_sum

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

Returns the sum of k-powerfree numbers <= n or in the range 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

powerful

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

Returns an array with the k-powerful (or k-full) numbers <= n, or in the range a..b.

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

powerful_count

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

Returns the number of k-powerful (or k-full) numbers <= n, or in the range 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

powerful_each

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

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

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 numbers in the range 50..100

Aliases: each_powerful

powerful_part

k.powerful_part(n)

The k-powerful part of n: largest k-powerful number dividing n.

say 2.powerful_part(n)      # squarefull part of n (OEIS: A057521)
say 3.powerful_part(n)      # cubefull part of n (OEIS: A360540)

powerful_sum

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

Returns the sum of k-powerful (or k-full) numbers <= n, or in the range a..b.

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

prev_composite

prev_composite(n)

Returns the previous composite number smaller than n.

prev_composites

n.prev_composites(start)

Returns an array with n consecutive decreasing composite numbers starting from start.

say 5.prev_composites(10)   #=> [10, 9, 8, 6, 4]
say 5.prev_composites(97)   #=> [96, 95, 94, 93, 92]

prev_cubefree

n.prev_cubefree

Returns the previous cubefree number smaller than n.

prev_cubefull

n.prev_cubefull

Returns the previous cubefull or (3-full) number smaller than n.

prev_noncubefree

n.prev_noncubefree

Returns the previous noncubefree number smaller than n.

prev_nonpowerfree

n.prev_nonpowerfree(k)

Returns the previous k-nonpowerfree number smaller than n.

prev_nonsquarefree

n.prev_nonsquarefree

Returns the previous nonsquarefree number smaller than n.

prev_perfect_power

prev_perfect_power(n)
prev_perfect_power(n,k)

Returns the previous perfect power smaller than n.

say prev_perfect_power(1e6)     #=> 998001

When k is given, it returns the previous k-perfect-power smaller than n.

say prev_perfect_power(1e6,3)   #=> 970299

prev_pow

n.prev_pow(b)

Returns the previous perfect power smaller than n, with base b. Returns NaN when n is <= 1.

say 65.prev_pow(2)      #=> 64
say 64.prev_pow(2)      #=> 32

Aliases: prev_power

prev_powerfree

n.prev_powerfree(k=2)

Returns the previous k-powerfree number smaller than n.

prev_powerful

n.prev_powerful(k=2)

Returns the previous k-powerful (or k-full) number smaller than n.

prev_rough(n, k)

Returns the largest k-rough integer strictly less than n.

Returns NaN if n <= 1.

prev_smooth(n, k)

Returns the largest k-smooth integer strictly less than n.

Returns NaN if n <= 1.

prev_squarefree

n.prev_squarefree

Returns the previous squarefree number smaller than n.

prev_squarefull

n.prev_squarefull

Returns the previous squarefull (or 2-full) number smaller than n.

primitive_part

n.primitive_part(f)

Returns the primitive part of 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) } }

primorial

n.primorial

Returns the product of the primes less than or equal to n.

say primorial(5)    #=>  30 (2*3*5)
say primorial(8)    #=> 210 (2*3*5*7)

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)

primorial_inflation

n.primorial_inflation

Primorial inflation of n.

Defined as:

primorial_inflation(n) = n.factor.prod {|p| primorial(p) }
primorial_inflation(n) = A108951(n)

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

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

rough_count

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

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

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

Also works with arbitrarily large n:

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

rough_divisors

k.rough_divisors(n)

Returns an array with the divisors of n that are k-rough.

Equivalent with (but more efficient):

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

rough_part

k.rough_part(n)

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

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

smooth_count

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

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

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

smooth_divisors

k.smooth_divisors(n)

Returns an array with the divisors of n that are k-smooth.

Equivalent with (but more efficient):

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

smooth_numbers

smooth_numbers(limit, [p1,p2,...])
smooth_numbers(limit, [p1,p2,...], {|n,p| ... })

Returns an array containing all the smooth numbers that factorize over the given array of primes.

say smooth_numbers(100, [2,3,5])    # 5-smooth numbers <= 100

An optional block can be provided, which is called with two arguments, n and p, where p is the current prime and n is a smooth number that is a multiple of p. If the block returns true, then the number n will be included in the returned array.

# 5-smooth numbers <= 100 that are squarefree
say smooth_numbers(100, [2,3,5], {|n,p| n.valuation(p) < 2 })

# 5-smooth numbers <= 100 that are cubefree
say smooth_numbers(100, [2,3,5], {|n,p| n.valuation(p) < 3 })

smooth_part

k.smooth_part(n)

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

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

square_part

square_part(n)

Largest square dividing n. (OEIS A008833)

say 30.of { .square_part }      #=> OEIS: A008833

squarefree

squarefree(n)
squarefree(a,b)

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

# Squarefree numbers in the interval [100, 200]
say squarefree(100, 200)

# Squarefree numbers <= 100
say squarefree(100)

squarefree_count

squarefree_count(n)
squarefree_count(a,b)

Returns the count of squarefree integers <= n, or in the range a..b.

# The number of squarefree integers between 10^5 and 10^6
say squarefree_count(10**5, 10**6) # 547132

# The number of squarefree integers <= 2^40
say squarefree_count(2**40)        # 668422917419

Aliases: square_free_count

squarefree_each

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

Iterates over the squarefree numbers in a given range.

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

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

Aliases: each_squarefree

squarefree_sum

squarefree_sum(n)
squarefree_sum(a,b)

Returns the sum of squarefree numbers <= n or in the range a..b. (OEIS A066779)

say squarefree_sum(1e12)    #=> 303963551353876732927386
say squarefree_sum(1e13)    #=> 30396355090144154315002969

squarefull

squarefull(n)
squarefull(a,b)

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

squarefull_count

squarefull_count(n)
squarefull_count(a,b)

Returns the count of squarefull numbers <= n, or in the range a..b.

squarefull_each

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

Iterates over the squarefull (or 2-full) numbers <= n, or in the range a..b.

Aliases: each_squarefull

squarefull_part

squarefull_part(n)

Largest divisor of n that is a squarefull (2-powerful) number. (OEIS A057521)

squarefull_sum

squarefull_sum(n)
squarefull_sum(a,b)

Returns the sum of squarefull numbers <= n, or in the range a..b.

squares_r

squares_r(n, k=2)

The sum of squares function r_k(n) returns the number of ways of representing n as a sum of k squares.

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: sum_of_squares_count

superprimorial

superprimorial(n)

Returns the product of first n primorials. (OEIS A006939)

say 10.of { .superprimorial }   #=> OEIS: A006939

RANDOM NUMBERS

irand

irand(n)
irand(a,b)

Returns a cryptographically-secure pseudorandom integer in the range 0..n (all inclusive), or in the range a..b (all inclusive).

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

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

The underlying CSPRNG is ISAAC-32.

miller_rabin_random

n.miller_rabin_random(k)

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

rand

rand(n)
rand(a,b)

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

random_bytes

n.random_bytes

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

random_maurer_nbit_prime

random_maurer_nbit_prime(n)

Generate a random n-bit Marurer prime.

Aliases: random_nbit_maurer_prime

random_nbit_prime

random_nbit_prime(n)

Returns a cryptographically-secure random n-bit prime.

say random_nbit_prime(128)  #=> 128-bit random prime

random_nbit_safe_prime

random_nbit_safe_prime(n)

Returns a random n-bit safe-prime.

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

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

random_nbit_strong_prime

random_nbit_strong_prime(n)

Generate a random n-bit strong prime number.

Aliases: random_strong_nbit_prime

random_ndigit_prime

random_ndigit_prime(n)

Returns a random n-digit prime.

say random_ndigit_prime(100)    # 100-digit random prime

random_prime

random_prime(n)
random_prime(a,b)

Returns a cryptographically-secure random prime <= n, or in [a,b].

say random_prime(100)       #=> Random prime in [2, 100]
say random_prime(100, 200)  #=> Random prime in [100, 200]

random_string

n.random_string

Returns a string with n random characters (bytes).

seed

seed(n)

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

urand

urand(n)
urand(a, b)

Given a positive integer n, returns a cryptographically-secure pseudorandom unsigned integer smaller than n.

say urand(100)          # pseudorandom integer in [0,99]

When two arguments are given, returns a uniform pseudorandom unsigned integer in the range [a,b] (inclusive on both ends).

say urand(100, 110)     # pseudorandom integer in [100, 110]

Both arguments must be nonnegative. If a is greater than b, then will return a pseudorandom integer in the range [b,a].

The underlying CSPRNG is ISAAC-32.

Aliases: urandomm

useed

iseed(n)
useed(n)

Takes a nonnegative integer n, which is used to seed the internal CSPRNG (currently, ISAAC-32) used for random functions, such as irand, urand and the random prime functions.

For good security, n should be at least an 128-bit random integer.

Internally, the value of n is converted into a binary string, SHA-512 hashed and expanded into a 8192-bit key.

NUMERIC PREDICATES & UTILITY

eval

eval(x,v)

Returns back x.

is_float

x.is_float

Returns true if x is stored as a floating-point number.

is_inf

x.is_inf

Returns true if x equals positive infinity (Inf).

is_int

x.is_int

Returns true if x is an integer.

is_mone

x.is_mone

Returns true if x equals -1.

is_nan

x.is_nan

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

is_neg

x.is_neg

Returns true if x is negative.

Aliases: is_negative

is_ninf

x.is_ninf

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

is_one

n.is_one

Returns true when n equals 1.

is_pos

x.is_pos

Returns true when x is a positive integer.

Aliases: is_positive

is_rat

x.is_rat

Returns true if x is a rational number.

is_zero

x.is_zero

Returns true when x equals 0.