Program Listing for File ntheory.h

Return to documentation for file (symengine/symengine/ntheory.h)

#ifndef SYMENGINE_NTHEORY_H
#define SYMENGINE_NTHEORY_H

#include <symengine/integer.h>

namespace SymEngine
{

// Prime Functions
int probab_prime_p(const Integer &a, unsigned reps = 25);
RCP<const Integer> nextprime(const Integer &a);

// Basic Number-theoretic functions
RCP<const Integer> gcd(const Integer &a, const Integer &b);
RCP<const Integer> lcm(const Integer &a, const Integer &b);
void gcd_ext(const Ptr<RCP<const Integer>> &g, const Ptr<RCP<const Integer>> &s,
             const Ptr<RCP<const Integer>> &t, const Integer &a,
             const Integer &b);
RCP<const Integer> mod(const Integer &n, const Integer &d);
RCP<const Integer> quotient(const Integer &n, const Integer &d);
void quotient_mod(const Ptr<RCP<const Integer>> &q,
                  const Ptr<RCP<const Integer>> &r, const Integer &a,
                  const Integer &b);
RCP<const Integer> mod_f(const Integer &n, const Integer &d);
RCP<const Integer> quotient_f(const Integer &n, const Integer &d);
void quotient_mod_f(const Ptr<RCP<const Integer>> &q,
                    const Ptr<RCP<const Integer>> &r, const Integer &a,
                    const Integer &b);
int mod_inverse(const Ptr<RCP<const Integer>> &b, const Integer &a,
                const Integer &m);

bool crt(const Ptr<RCP<const Integer>> &R,
         const std::vector<RCP<const Integer>> &rem,
         const std::vector<RCP<const Integer>> &mod);

RCP<const Integer> fibonacci(unsigned long n);

void fibonacci2(const Ptr<RCP<const Integer>> &g,
                const Ptr<RCP<const Integer>> &s, unsigned long n);

RCP<const Integer> lucas(unsigned long n);

void lucas2(const Ptr<RCP<const Integer>> &g, const Ptr<RCP<const Integer>> &s,
            unsigned long n);

RCP<const Integer> binomial(const Integer &n, unsigned long k);

RCP<const Integer> factorial(unsigned long n);

bool divides(const Integer &a, const Integer &b);

int factor(const Ptr<RCP<const Integer>> &f, const Integer &n, double B1 = 1.0);

int factor_trial_division(const Ptr<RCP<const Integer>> &f, const Integer &n);

int factor_lehman_method(const Ptr<RCP<const Integer>> &f, const Integer &n);

int factor_pollard_pm1_method(const Ptr<RCP<const Integer>> &f,
                              const Integer &n, unsigned B = 10,
                              unsigned retries = 5);

int factor_pollard_rho_method(const Ptr<RCP<const Integer>> &f,
                              const Integer &n, unsigned retries = 5);

void prime_factors(std::vector<RCP<const Integer>> &primes, const Integer &n);
void prime_factor_multiplicities(map_integer_uint &primes, const Integer &n);
// Sieve class stores all the primes upto a limit. When a prime or a list of
// prime
// is requested, if the prime is not there in the sieve, it is extended to hold
// that
// prime. The implementation is a very basic Eratosthenes sieve, but the code
// should
// be quite optimized. For limit=1e8, it is about 20x slower than the
// `primesieve` library (1206ms vs 55.63ms).
class Sieve
{

private:
    static std::vector<unsigned> _primes;
    static void _extend(unsigned limit);
    static unsigned _sieve_size;
    static bool _clear;

public:
    // Returns all primes up to the `limit` (including). The vector `primes`
    // should
    // be empty on input and it will be filled with the primes.
    static void generate_primes(std::vector<unsigned> &primes, unsigned limit);
    // Clear the array of primes stored
    static void clear();
    // Set the sieve size in kilobytes. Set it to L1d cache size for best
    // performance.
    // Default value is 32.
    static void set_sieve_size(unsigned size);
    // Set whether the sieve is cleared after the sieve is extended in internal
    // functions
    static void set_clear(bool clear);

    class iterator
    {

    private:
        unsigned _index;
        unsigned _limit;

    public:
        // Iterator that generates primes upto limit
        iterator(unsigned limit);
        // Iterator that generates primes with no limit.
        iterator();
        // Destructor
        ~iterator();
        // Next prime
        unsigned next_prime();
    };
};

RCP<const Number> bernoulli(unsigned long n);
RCP<const Number> harmonic(unsigned long n, long m = 1);
// Primitive root calculated is the smallest when n is prime.
bool primitive_root(const Ptr<RCP<const Integer>> &g, const Integer &n);
void primitive_root_list(std::vector<RCP<const Integer>> &roots,
                         const Integer &n);
RCP<const Integer> totient(const RCP<const Integer> &n);
RCP<const Integer> carmichael(const RCP<const Integer> &n);
bool multiplicative_order(const Ptr<RCP<const Integer>> &o,
                          const RCP<const Integer> &a,
                          const RCP<const Integer> &n);
int legendre(const Integer &a, const Integer &n);
int jacobi(const Integer &a, const Integer &n);
int kronecker(const Integer &a, const Integer &n);
void nthroot_mod_list(std::vector<RCP<const Integer>> &roots,
                      const RCP<const Integer> &a, const RCP<const Integer> &n,
                      const RCP<const Integer> &m);
bool nthroot_mod(const Ptr<RCP<const Integer>> &root,
                 const RCP<const Integer> &a, const RCP<const Integer> &n,
                 const RCP<const Integer> &m);
bool powermod(const Ptr<RCP<const Integer>> &powm, const RCP<const Integer> &a,
              const RCP<const Number> &b, const RCP<const Integer> &m);
void powermod_list(std::vector<RCP<const Integer>> &pows,
                   const RCP<const Integer> &a, const RCP<const Number> &b,
                   const RCP<const Integer> &m);

vec_integer_class quadratic_residues(const Integer &a);

bool is_quad_residue(const Integer &a, const Integer &p);
bool is_nth_residue(const Integer &a, const Integer &n, const Integer &mod);
// mu(n) = 1 if n is a square-free positive integer with an even number of prime
// factors
// mu(n) = −1 if n is a square-free positive integer with an odd number of prime
// factors
// mu(n) = 0 if n has a squared prime factor
int mobius(const Integer &a);
// Mertens Function
// mertens(n) -> Sum of mobius(i) for i from 1 to n
long mertens(const unsigned long a);
}
#endif