Goulib.math2 module

more math than math standard library, without numpy

Goulib.math2.cmp(x, y)[source]

Compare the two objects x and y and return an integer according to the outcome. The return value is negative if x < y, zero if x == y and strictly positive if x > y.

Goulib.math2.is_number(x)[source]
Returns:True if x is a number of any type
Goulib.math2.sign(number)[source]
Returns:1 if number is positive, -1 if negative, 0 if ==0
Goulib.math2.rint(v)[source]
Returns:int value nearest to float v
Goulib.math2.is_integer(x, rel_tol=0, abs_tol=0)[source]
Returns:True if float x is an integer within tolerances
Goulib.math2.int_or_float(x, rel_tol=0, abs_tol=0)[source]
Parameters:x – int or float
Returns:int if x is (almost) an integer, otherwise float
Goulib.math2.format(x, decimals=3)[source]

formats a float with given number of decimals, but not an int

Returns:string repr of x with decimals if not int
Goulib.math2.gcd(*args)[source]

greatest common divisor of an arbitrary number of args

Goulib.math2.lcm(*args)[source]

least common multiple of any number of integers

Goulib.math2.xgcd(a, b)[source]

Extended GCD

Returns:(gcd, x, y) where gcd is the greatest common divisor of a and b

with the sign of b if b is nonzero, and with the sign of a if b is 0. The numbers x,y are such that gcd = ax+by.

Goulib.math2.coprime(*args)[source]
Returns:True if args are coprime to each other
Goulib.math2.coprimes_gen(limit)[source]

generates coprime pairs using Farey sequence

Goulib.math2.is_primitive_root(x, m, s={})[source]

returns True if x is a primitive root of m

Parameters:s – set of coprimes to m, if already known
Goulib.math2.primitive_root_gen(m)[source]

generate primitive roots modulo m

Goulib.math2.primitive_roots(modulo)[source]
Goulib.math2.quad(a, b, c, allow_complex=False)[source]

solves quadratic equations aX^2+bX+c=0

Parameters:
  • a,b,c – floats
  • allow_complex – function returns complex roots if True
Returns:

x1,x2 real or complex solutions

Goulib.math2.ceildiv(a, b)[source]
Goulib.math2.ipow(x, y, z=0)[source]
Parameters:
  • x – number (int or float)
  • y – int power
  • z – int optional modulus
Returns:

(x**y) % z as integer if possible

Goulib.math2.pow(x, y, z=0)[source]
Returns:(x**y) % z as integer
Goulib.math2.isqrt(n)[source]

integer square root

Returns:largest int x for which x * x <= n
Goulib.math2.sqrt(n)[source]

square root :return: int, float or complex depending on n

Goulib.math2.is_square(n)[source]
Goulib.math2.introot(n, r=2)[source]

integer r-th root

Returns:int, greatest integer less than or equal to the r-th root of n.

For negative n, returns the least integer greater than or equal to the r-th root of n, or None if r is even.

Goulib.math2.is_power(n)[source]
Returns:integer that, when squared/cubed/etc, yields n,

or 0 if no such integer exists. Note that the power to which this number is raised will be prime.

Goulib.math2.multiply(x, y)[source]

Karatsuba fast multiplication algorithm

https://en.wikipedia.org/wiki/Karatsuba_algorithm

Copyright (c) 2014 Project Nayuki http://www.nayuki.io/page/karatsuba-multiplication

Goulib.math2.accsum(it)[source]

Yield accumulated sums of iterable: accsum(count(1)) -> 1,3,6,10,…

Goulib.math2.cumsum(it)

Yield accumulated sums of iterable: accsum(count(1)) -> 1,3,6,10,…

Goulib.math2.mul(nums, init=1)[source]
Returns:Product of nums
Goulib.math2.dot_vv(a, b, default=0)[source]

dot product for vectors

Parameters:
  • a – vector (iterable)
  • b – vector (iterable)
  • default – default value of the multiplication operator
Goulib.math2.dot_mv(a, b, default=0)[source]

dot product for vectors

Parameters:
  • a – matrix (iterable or iterables)
  • b – vector (iterable)
  • default – default value of the multiplication operator
Goulib.math2.dot_mm(a, b, default=0)[source]

dot product for matrices

Parameters:
  • a – matrix (iterable or iterables)
  • b – matrix (iterable or iterables)
  • default – default value of the multiplication operator
Goulib.math2.dot(a, b, default=0)[source]

dot product

general but slow : use dot_vv, dot_mv or dot_mm if you know a and b’s dimensions

Goulib.math2.zeros(shape)[source]
See:https://docs.scipy.org/doc/numpy/reference/generated/numpy.zeros.html
Goulib.math2.diag(v)[source]

Create a two-dimensional array with the flattened input as a diagonal.

Parameters:v – If v is a 2-D array, return a copy of its diagonal. If v is a 1-D array, return a 2-D array with v on the diagonal
See:https://docs.scipy.org/doc/numpy/reference/generated/numpy.diag.html#numpy.diag
Goulib.math2.identity(n)[source]
Goulib.math2.eye(n)
Goulib.math2.transpose(m)[source]
Returns:matrix m transposed
Goulib.math2.maximum(m)[source]

Compare N arrays and returns a new array containing the element-wise maxima

Parameters:m – list of arrays (matrix)
Returns:list of maximal values found in each column of m
See:http://docs.scipy.org/doc/numpy/reference/generated/numpy.maximum.html
Goulib.math2.minimum(m)[source]

Compare N arrays and returns a new array containing the element-wise minima

Parameters:m – list of arrays (matrix)
Returns:list of minimal values found in each column of m
See:http://docs.scipy.org/doc/numpy/reference/generated/numpy.minimum.html
Goulib.math2.vecadd(a, b, fillvalue=0)[source]

addition of vectors of inequal lengths

Goulib.math2.vecsub(a, b, fillvalue=0)[source]

substraction of vectors of inequal lengths

Goulib.math2.vecneg(a)[source]

unary negation

Goulib.math2.vecmul(a, b)[source]

product of vectors of inequal lengths

Goulib.math2.vecdiv(a, b)[source]

quotient of vectors of inequal lengths

Goulib.math2.veccompare(a, b)[source]

compare values in 2 lists. returns triple number of pairs where [a<b, a==b, a>b]

Goulib.math2.sat(x, low=0, high=None)[source]

saturates x between low and high

Goulib.math2.norm_2(v)[source]
Returns:“normal” euclidian norm of vector v
Goulib.math2.norm_1(v)[source]
Returns:“manhattan” norm of vector v
Goulib.math2.norm_inf(v)[source]
Returns:infinite norm of vector v
Goulib.math2.norm(v, order=2)[source]
See:http://docs.scipy.org/doc/numpy/reference/generated/numpy.linalg.norm.html
Goulib.math2.dist(a, b, norm=<function norm_2>)[source]
Goulib.math2.vecunit(v, norm=<function norm_2>)[source]
Returns:vector normalized
Goulib.math2.hamming(s1, s2)[source]

Calculate the Hamming distance between two iterables

Goulib.math2.sets_dist(a, b)[source]
See:http://stackoverflow.com/questions/11316539/calculating-the-distance-between-two-unordered-sets
Goulib.math2.sets_levenshtein(a, b)[source]

levenshtein distance on sets

See:http://en.wikipedia.org/wiki/Levenshtein_distance
Goulib.math2.levenshtein(seq1, seq2)[source]

levenshtein distance

Returns:distance between 2 iterables
See:http://en.wikipedia.org/wiki/Levenshtein_distance
Goulib.math2.recurrence(coefficients, values, cst=0, max=None, mod=0)[source]

general generator for recurrences

Parameters:
  • values – list of initial values
  • coefficients – list of factors defining the recurrence
Goulib.math2.fibonacci_gen(max=None, mod=0)[source]

Generate fibonacci serie

Goulib.math2.fibonacci(n, mod=0)[source]

fibonacci series n-th element

Parameters:
  • n – int can be extremely high, like 1e19 !
  • mod – int optional modulo
Goulib.math2.is_fibonacci(n)[source]

returns True if n is in Fibonacci series

Goulib.math2.pisano_cycle(mod)[source]
Goulib.math2.pisano_period(mod)[source]
Goulib.math2.pascal_gen()[source]

Pascal’s triangle read by rows: C(n,k) = binomial(n,k) = n!/(k!*(n-k)!), 0<=k<=n.

https://oeis.org/A007318

Goulib.math2.catalan(n)[source]

Catalan numbers: C(n) = binomial(2n,n)/(n+1) = (2n)!/(n!(n+1)!).

Goulib.math2.catalan_gen()[source]

Generate Catalan numbers: C(n) = binomial(2n,n)/(n+1) = (2n)!/(n!(n+1)!). Also called Segner numbers.

Goulib.math2.is_pythagorean_triple(a, b, c)[source]
Goulib.math2.primitive_triples()[source]

generates primitive Pythagorean triplets x<y<z

sorted by hypotenuse z, then longest side y through Berggren’s matrices and breadth first traversal of ternary tree :see: https://en.wikipedia.org/wiki/Tree_of_primitive_Pythagorean_triples

Goulib.math2.triples()[source]

generates all Pythagorean triplets triplets x<y<z sorted by hypotenuse z, then longest side y

Goulib.math2.divisors(n)[source]
Returns:all divisors of n: divisors(12) -> 1,2,3,6,12

including 1 and n, except for 1 which returns a single 1 to avoid messing with sum of divisors…

Goulib.math2.proper_divisors(n)[source]
Returns:all divisors of n except n itself.
Goulib.math2.sieve(n, oneisprime=False)[source]

prime numbers from 2 to a prime < n Very fast (n<10,000,000) in 0.4 sec.

Example: >>>prime_sieve(25) [2, 3, 5, 7, 11, 13, 17, 19, 23]

Algorithm & Python source: Robert William Hanks http://stackoverflow.com/questions/17773352/python-sieve-prime-numbers

Goulib.math2.primes(n)[source]

memoized list of n first primes

Warning:do not call with large n, use prime_gen instead
Goulib.math2.is_prime_euler(n, eb=(2, ))[source]

Euler’s primality test

Parameters:
  • n – int number to test
  • eb – test basis
Returns:

False if not prime, True if prime, but also for many pseudoprimes…

See:

https://en.wikipedia.org/wiki/Euler_pseudoprime

Goulib.math2.is_prime(n, oneisprime=False, tb=(3, 5, 7, 11), eb=(2, ), mrb=None)[source]

main primality test.

Parameters:
  • n – int number to test
  • oneisprime – bool True if 1 should be considered prime (it was, a long time ago)
  • tb – trial division basis
  • eb – Euler’s test basis
  • mrb – Miller-Rabin basis, automatic if None
See:

https://en.wikipedia.org/wiki/Baillie%E2%80%93PSW_primality_test

It’s an implementation of the BPSW test (Baillie-Pomerance-Selfridge-Wagstaff) with some prefiltes for speed and is deterministic for all numbers less than 2^64 Iin fact, while infinitely many false positives are conjectured to exist, no false positives are currently known. The prefilters consist of trial division against 2 and the elements of the tuple tb, checking whether n is square, and Euler’s primality test to the bases in the tuple eb. If the number is less than 3825123056546413051, we use the Miller-Rabin test on a set of bases for which the test is known to be deterministic over this range.

Goulib.math2.nextprime(n)[source]

Determines, with some semblance of efficiency, the least prime number strictly greater than n.

Goulib.math2.prevprime(n)[source]

Determines, very inefficiently, the largest prime number strictly smaller than n.

Goulib.math2.primes_gen(start=2, stop=None)[source]

generate prime numbers from ‘start’

Goulib.math2.random_prime(bits)[source]

returns a random number of the specified bit length

Goulib.math2.euclid_gen()[source]

generates Euclid numbers: 1 + product of the first n primes

Goulib.math2.prime_factors(num, start=2)[source]

generates all prime factors (ordered) of num

Goulib.math2.prime_divisors(num, start=2)[source]

generates unique prime divisors (ordered) of num

Goulib.math2.is_multiple(n, factors)[source]

return True if n has ONLY factors as prime factors

Goulib.math2.factorize(n)[source]

find the prime factors of n along with their frequencies. Example:

>>> factor(786456)
[(2,3), (3,3), (11,1), (331,1)]
Goulib.math2.factors(n)[source]
Goulib.math2.number_of_divisors(n)[source]
Goulib.math2.omega(n)[source]

Number of distinct primes dividing n

Goulib.math2.bigomega(n)[source]

Number of prime divisors of n counted with multiplicity

Goulib.math2.moebius(n)[source]

Möbius (or Moebius) function mu(n). mu(1) = 1; mu(n) = (-1)^k if n is the product of k different primes; otherwise mu(n) = 0.

Goulib.math2.euler_phi(n)[source]

Euler totient function

See:http://stackoverflow.com/questions/1019040/how-many-numbers-below-n-are-coprimes-to-n
Goulib.math2.totient(n)

Euler totient function

See:http://stackoverflow.com/questions/1019040/how-many-numbers-below-n-are-coprimes-to-n
Goulib.math2.prime_ktuple(constellation)[source]

generates tuples of primes with specified differences

Parameters:constellation – iterable of int differences betwwen primes to return
Note:negative int means the difference must NOT be prime
See:https://en.wikipedia.org/wiki/Prime_k-tuple

(0, 2) twin primes (0, 4) cousin primes (0, 6) sexy primes (0, 2, 6), (0, 4, 6) prime triplets (0, 6, 12, -18) sexy prime triplets (0, 2, 6, 8) prime quadruplets (0, 6, 12, 18) sexy prime quadruplets (0, 2, 6, 8, 12), (0, 4, 6, 10, 12) quintuplet primes (0, 4, 6, 10, 12, 16) sextuplet primes

Goulib.math2.twin_primes()[source]
Goulib.math2.cousin_primes()[source]
Goulib.math2.sexy_primes()[source]
Goulib.math2.sexy_prime_triplets()[source]
Goulib.math2.sexy_prime_quadruplets()[source]
Goulib.math2.lucas_lehmer(p)[source]

Lucas Lehmer primality test for Mersenne exponent p

Parameters:p – int
Returns:True if 2^p-1 is prime
Goulib.math2.digits_gen(num, base=10)[source]

generates int digits of num in base BACKWARDS

Goulib.math2.digits(num, base=10, rev=False)[source]
Returns:list of digits of num expressed in base, optionally reversed
Goulib.math2.digsum(num, f=None, base=10)[source]

sum of digits

Parameters:
  • num – number
  • f – int power or function applied to each digit
  • base – optional base
Returns:

sum of f(digits) of num

digsum(num) -> sum of digits digsum(num,base=2) -> number of 1 bits in binary represenation of num digsum(num,2) -> sum of the squares of digits digsum(num,f=lambda x:x**x) -> sum of the digits elevaed to their own power

Goulib.math2.integer_exponent(a, b=10)[source]
Returns:int highest power of b that divides a.
See:https://reference.wolfram.com/language/ref/IntegerExponent.html
Goulib.math2.trailing_zeros(a, b=10)
Returns:int highest power of b that divides a.
See:https://reference.wolfram.com/language/ref/IntegerExponent.html
Goulib.math2.power_tower(v)[source]
Returns:v[0]**v[1]**v[2] …
See:http://ajcr.net//Python-power-tower/
Goulib.math2.carries(a, b, base=10, pos=0)[source]
Returns:int number of carries required to add a+b in base
Goulib.math2.powertrain(n)[source]
Returns:v[0]**v[1]*v[2]**v[3] …**(v[-1] or 0)
Author:# Chai Wah Wu, Jun 16 2017
See:http://oeis.org/A133500
Goulib.math2.str_base(num, base=10, numerals='0123456789abcdefghijklmnopqrstuvwxyz')[source]
Returns:

string representation of num in base

Parameters:
  • num – int number (decimal)
  • base – int base, 10 by default
  • numerals – string with all chars representing numbers in base base. chars after the base-th are ignored
Goulib.math2.num_from_digits(digits, base=10)[source]
Parameters:
  • digits – string or list of digits representing a number in given base
  • base – int base, 10 by default
Returns:

int number

Goulib.math2.reverse(i)[source]
Goulib.math2.is_palindromic(num, base=10)[source]

Check if ‘num’ in base ‘base’ is a palindrome, that’s it, if it can be read equally from left to right and right to left.

Goulib.math2.is_permutation(num1, num2, base=10)[source]

Check if ‘num1’ and ‘num2’ have the same digits in base

Goulib.math2.is_pandigital(num, base=10)[source]
Returns:True if num contains all digits in specified base
Goulib.math2.bouncy(n)[source]
Goulib.math2.repunit_gen(digit=1)[source]

generate repunits

Goulib.math2.repunit(n, digit=1)[source]
Returns:nth repunit
Goulib.math2.rational_form(numerator, denominator)[source]

information about the decimal representation of a rational number.

Returns:5 integer : integer, decimal, shift, repeat, cycle
  • shift is the len of decimal with leading zeroes if any
  • cycle is the len of repeat with leading zeroes if any
Goulib.math2.rational_str(n, d)[source]
Goulib.math2.rational_cycle(num, den)[source]

periodic part of the decimal expansion of num/den. Any initial 0’s are placed at end of cycle.

See:https://oeis.org/A036275
Goulib.math2.tetrahedral(n)[source]
Returns:int n-th tetrahedral number
See:https://en.wikipedia.org/wiki/Tetrahedral_number
Goulib.math2.sum_of_squares(n)[source]
Returns:1^2 + 2^2 + 3^2 + … + n^2
See:https://en.wikipedia.org/wiki/Square_pyramidal_number
Goulib.math2.pyramidal(n)
Returns:1^2 + 2^2 + 3^2 + … + n^2
See:https://en.wikipedia.org/wiki/Square_pyramidal_number
Goulib.math2.sum_of_cubes(n)[source]
Returns:1^3 + 2^3 + 3^3 + … + n^3
See:https://en.wikipedia.org/wiki/Squared_triangular_number
Goulib.math2.bernouilli_gen(init=1)[source]

generator of Bernouilli numbers

Parameters:init – int -1 or +1.
  • -1 for “first Bernoulli numbers” with B1=-1/2
  • +1 for “second Bernoulli numbers” with B1=+1/2

https://en.wikipedia.org/wiki/Bernoulli_number https://rosettacode.org/wiki/Bernoulli_numbers#Python:_Optimised_task_algorithm

Goulib.math2.bernouilli(n, init=1)[source]
Goulib.math2.faulhaber(n, p)[source]

sum of the p-th powers of the first n positive integers

Returns:1^p + 2^p + 3^p + … + n^p
See:https://en.wikipedia.org/wiki/Faulhaber%27s_formula
Goulib.math2.is_happy(n)[source]
Goulib.math2.lychrel_seq(n)[source]
Goulib.math2.lychrel_count(n, limit=96)[source]

number of lychrel iterations before n becomes palindromic

Parameters:
  • n – int number to test
  • limit – int max number of loops. default 96 corresponds to the known most retarded non lychrel number
Warning:

there are palindrom lychrel numbers such as 4994

Goulib.math2.is_lychrel(n, limit=96)[source]
Warning:there are palindrom lychrel numbers such as 4994
Goulib.math2.polygonal(s, n)[source]
Goulib.math2.triangle(n)[source]
Returns:nth triangle number, defined as the sum of [1,n] values.
See:http://en.wikipedia.org/wiki/Triangular_number
Goulib.math2.triangular(n)
Returns:nth triangle number, defined as the sum of [1,n] values.
See:http://en.wikipedia.org/wiki/Triangular_number
Goulib.math2.is_triangle(x)[source]
Returns:True if x is a triangle number
Goulib.math2.is_triangular(x)
Returns:True if x is a triangle number
Goulib.math2.square(n)[source]
Goulib.math2.pentagonal(n)[source]
Returns:nth pentagonal number
See:https://en.wikipedia.org/wiki/Pentagonal_number
Goulib.math2.is_pentagonal(n)[source]
Returns:True if x is a pentagonal number
Goulib.math2.hexagonal(n)[source]
Returns:nth hexagonal number
See:https://en.wikipedia.org/wiki/Hexagonal_number
Goulib.math2.is_hexagonal(n)[source]
Goulib.math2.heptagonal(n)[source]
Goulib.math2.is_heptagonal(n)[source]
Goulib.math2.octagonal(n)[source]
Goulib.math2.is_octagonal(n)[source]
Goulib.math2.partition(n)[source]

The partition function p(n)

gives the number of partitions of a nonnegative integer n into positive integers. (There is one partition of zero into positive integers, i.e. the empty partition, since the empty sum is defined as 0.)

See:http://oeis.org/wiki/Partition_function https://oeis.org/A000041
Goulib.math2.partitionsQ(n, d=0)[source]
Goulib.math2.get_cardinal_name(num)[source]

Get cardinal name for number (0 to 1 million)

Goulib.math2.abundance(n)[source]
Goulib.math2.is_perfect(n)[source]
Returns:-1 if n is deficient, 0 if perfect, 1 if abundant
See:https://en.wikipedia.org/wiki/Perfect_number,

https://en.wikipedia.org/wiki/Abundant_number, https://en.wikipedia.org/wiki/Deficient_number

Goulib.math2.number_of_digits(num, base=10)[source]

Return number of digits of num (expressed in base ‘base’)

Goulib.math2.chakravala(n)[source]

solves x^2 - n*y^2 = 1 for x,y integers

https://en.wikipedia.org/wiki/Pell%27s_equation https://en.wikipedia.org/wiki/Chakravala_method

Goulib.math2.factorialk(n, k)[source]

Multifactorial of n of order k, n(!!…!).

This is the multifactorial of n skipping k values. For example,
factorialk(17, 4) = 17!!!! = 17 * 13 * 9 * 5 * 1
In particular, for any integer n, we have
factorialk(n, 1) = factorial(n) factorialk(n, 2) = factorial2(n)
Parameters:n – int Calculate multifactorial. If n < 0, the return value is 0.

:param k : int Order of multifactorial. :return: int Multifactorial of n.

Goulib.math2.factorial2(n)[source]
Goulib.math2.factorial_gen()[source]

Generator of factorial

Goulib.math2.binomial(n, k)[source]

https://en.wikipedia.org/wiki/binomial

Goulib.math2.ncombinations(n, k)

https://en.wikipedia.org/wiki/binomial

Goulib.math2.binomial_exponent(n, k, p)[source]
Returns:int largest power of p that divides binomial(n,k)
Goulib.math2.log_factorial(n)[source]
Returns:float approximation of ln(n!) by Ramanujan formula
Goulib.math2.log_binomial(n, k)[source]
Returns:float approximation of ln(binomial(n,k))
Goulib.math2.ilog(a, b, upper_bound=False)[source]

discrete logarithm x such that b^x=a

Parameters:
  • a,b – integer
  • upper_bound – bool. if True, returns smallest x such that b^x>=a
Returns:

x integer such that b^x=a, or upper_bound, or None

https://en.wikipedia.org/wiki/Discrete_logarithm

Goulib.math2.angle(u, v, unit=True)[source]
Parameters:
  • u,v – iterable vectors
  • unit – bool True if vectors are unit vectors. False increases computations
Returns:

float angle n radians between u and v unit vectors i

Goulib.math2.sin_over_x(x)[source]

numerically safe sin(x)/x

Goulib.math2.slerp(u, v, t)[source]

spherical linear interpolation

Parameters:
  • u,v – 3D unit vectors
  • t – float in [0,1] interval
Returns:

vector interpolated between u and v

Goulib.math2.proportional(nseats, votes)[source]

assign n seats proportionaly to votes using the https://en.wikipedia.org/wiki/Hagenbach-Bischoff_quota method

Parameters:
  • nseats – int number of seats to assign
  • votes – iterable of int or float weighting each party
Result:

list of ints seats allocated to each party

Goulib.math2.triangular_repartition(x, n)[source]

divide 1 into n fractions such that:

  • their sum is 1
  • they follow a triangular linear repartition (sorry, no better name for now) where x/1 is the maximum
Goulib.math2.rectangular_repartition(x, n, h)[source]

divide 1 into n fractions such that:

  • their sum is 1
  • they follow a repartition along a pulse of height h<1
Goulib.math2.de_bruijn(k, n)[source]

De Bruijn sequence for alphabet k and subsequences of length n.

https://en.wikipedia.org/wiki/De_Bruijn_sequence

Goulib.math2.mod_inv(a, b)[source]
Goulib.math2.mod_div(a, b, m)[source]
Returns:x such that (b*x) mod m = a mod m
Goulib.math2.mod_fact(n, m)[source]
Returns:n! mod m
Goulib.math2.chinese_remainder(m, a)[source]

http://en.wikipedia.org/wiki/Chinese_remainder_theorem

Parameters:
  • m – list of int moduli
  • a – list of int remainders
Returns:

smallest int x such that x mod ni=ai

Goulib.math2.mod_binomial(n, k, m, q=None)[source]

calculates C(n,k) mod m for large n,k,m

Parameters:
  • n – int total number of elements
  • k – int number of elements to pick
  • m – int modulo (or iterable of (m,p) tuples used internally)
  • q – optional int power of m for prime m, used internally
Goulib.math2.baby_step_giant_step(y, a, n)[source]

solves Discrete Logarithm Problem (DLP) y = a**x mod n

Goulib.math2.mod_matmul(A, B, mod=0)[source]
Goulib.math2.mod_matpow(M, power, mod=0)[source]
Goulib.math2.matrix_power(M, power, mod=0)
Goulib.math2.mod_sqrt(n, p)[source]

modular sqrt(n) mod p

Goulib.math2.pi_digits_gen()[source]

generates pi digits as a sequence of INTEGERS ! using Jeremy Gibbons spigot generator

:see :http://www.cs.ox.ac.uk/people/jeremy.gibbons/publications/spigot.pdf

Goulib.math2.pfactor(n)[source]

Helper function for sprp.

Returns the tuple (x,y) where n - 1 == (2 ** x) * y and y is odd. We have this bit separated out so that we don’t waste time recomputing s and d for each base when we want to check n against multiple bases.

Goulib.math2.sprp(n, a, s=None, d=None)[source]

Checks n for primality using the Strong Probable Primality Test to base a. If present, s and d should be the first and second items, respectively, of the tuple returned by the function pfactor(n)

Goulib.math2.jacobi(a, p)[source]

Computes the Jacobi symbol (a|p), where p is a positive odd number. :see: https://en.wikipedia.org/wiki/Jacobi_symbol

Goulib.math2.pollardRho_brent(n)[source]

Brent’s improvement on Pollard’s rho algorithm.

Returns:int n if n is prime

otherwise, we keep chugging until we find a factor of n strictly between 1 and n. :see: https://en.wikipedia.org/wiki/Pollard%27s_rho_algorithm

Goulib.math2.pollard_pm1(n, B1=100, B2=1000)[source]

Pollard’s p+1 algorithm, two-phase version.

Returns:n if n is prime; otherwise, we keep chugging until we find a factor of n strictly between 1 and n.
Goulib.math2.mlucas(v, a, n)[source]

Helper function for williams_pp1(). Multiplies along a Lucas sequence modulo n.

Goulib.math2.williams_pp1(n)[source]

Williams’ p+1 algorithm. :return: n if n is prime otherwise, we keep chugging until we find a factor of n strictly between 1 and n.

Goulib.math2.ecadd(p1, p2, p0, n)[source]
Goulib.math2.ecdub(p, A, n)[source]
Goulib.math2.ecmul(m, p, A, n)[source]
Goulib.math2.factor_ecm(n, B1=10, B2=20)[source]

Factors n using the elliptic curve method, using Montgomery curves and an algorithm analogous to the two-phase variant of Pollard’s p-1 method. :return: n if n is prime otherwise, we keep chugging until we find a factor of n strictly between 1 and n

Goulib.math2.legendre(a, p)[source]

Functions to comptue the Legendre symbol (a|p). The return value isn’t meaningful if p is composite :see: https://en.wikipedia.org/wiki/Legendre_symbol

Goulib.math2.legendre2(a, p)[source]

Functions to comptue the Legendre symbol (a|p). The return value isn’t meaningful if p is composite :see: https://en.wikipedia.org/wiki/Legendre_symbol