1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182
|
#ifndef MPU_UTILITY_H
#define MPU_UTILITY_H
#include <math.h>
#include <gmp.h>
#include "ptypes.h"
extern int get_verbose_level(void);
extern void set_verbose_level(int level);
extern gmp_randstate_t* get_randstate(void);
extern void init_randstate(unsigned long seed);
extern void clear_randstate(void);
extern void mpz_isaac_urandomb(mpz_t rop, int nbits);
extern void mpz_isaac_urandomm(mpz_t rop, mpz_t n);
extern UV irand64(int nbits);
extern NV drand64(void);
extern UV is_power(mpz_t n, UV a);
extern UV prime_power(mpz_t prime, mpz_t n);
extern int is_primitive_root(mpz_t a, mpz_t b, int nprime);
/* tdiv_r is faster, but we'd need to guarantee the input is positive */
#define mpz_mulmod(r, a, b, n, t) \
do { mpz_mul(t, a, b); mpz_mod(r, t, n); } while (0)
#undef mpz_divmod
extern int mpz_divmod(mpz_t r, mpz_t a, mpz_t b, mpz_t n, mpz_t t);
/* s = sqrt(a) mod p */
extern int sqrtmod(mpz_t s, mpz_t a, mpz_t p);
extern int sqrtmod_t(mpz_t s, mpz_t a, mpz_t p,
mpz_t t, mpz_t t2, mpz_t b, mpz_t g); /* 4 tmp variables */
extern unsigned long modinverse(unsigned long a, unsigned long p);
#if __GNU_MP_VERSION < 5
/* Older versions left out a normalization step */
extern void gcdext(mpz_t g, mpz_t s, mpz_t t, const mpz_t a, const mpz_t b);
#else
#define gcdext(g,s,t,a,b) mpz_gcdext(g,s,t,a,b)
#endif
extern int chinese(mpz_t ret, mpz_t lcm, mpz_t *a, mpz_t *m, int items);
extern UV mpz_order_ui(UV r, mpz_t n, UV limit);
extern void mpz_arctan(mpz_t r, unsigned long base, mpz_t pow, mpz_t t1, mpz_t t2);
extern void mpz_arctanh(mpz_t r, unsigned long base, mpz_t pow, mpz_t t1, mpz_t t2);
extern void mpz_product(mpz_t* A, UV a, UV b);
extern void poly_mod_mul(mpz_t* px, mpz_t* py, UV r, mpz_t mod, mpz_t t1, mpz_t t2, mpz_t t3);
extern void poly_mod_pow(mpz_t *pres, mpz_t *pn, mpz_t power, UV r, mpz_t mod);
extern void poly_mod(mpz_t *pres, mpz_t *pn, UV *dn, mpz_t mod);
extern void polyz_mod(mpz_t *pres, mpz_t *pn, long *dn, mpz_t mod);
extern void polyz_set(mpz_t* pr, long* dr, mpz_t* ps, long ds);
extern void polyz_print(const char* header, mpz_t* pn, long dn);
extern void polyz_mulmod(mpz_t* pr, mpz_t* px, mpz_t *py, long *dr, long dx, long dy, mpz_t mod);
extern void polyz_div(mpz_t *pq, mpz_t *pr, mpz_t *pn, mpz_t *pd,
long *dq, long *dr, long dn, long dd, mpz_t NMOD);
extern void polyz_pow_polymod(mpz_t* pres, mpz_t* pn, mpz_t* pmod,
long *dres, long dn, long dmod,
mpz_t power, mpz_t NMOD);
extern void polyz_gcd(mpz_t* pres, mpz_t* pa, mpz_t* pb, long* dres, long da, long db, mpz_t MODN);
extern void polyz_root_deg1(mpz_t root, mpz_t* pn, mpz_t NMOD);
extern void polyz_root_deg2(mpz_t root1, mpz_t root2, mpz_t* pn, mpz_t NMOD);
/* Find roots of a polynomial mod a prime, slightly modified. */
/* We will stop if we've found at least maxroots unique roots. */
extern void polyz_roots_modp(mpz_t** roots, long *nroots, long maxroots,
mpz_t *pP, long dP, mpz_t NMOD);
/* Solve x^2 + |D|y^2 = p */
extern int cornacchia(mpz_t x, mpz_t y, mpz_t D, mpz_t p);
/* Solve x^2 + |D|y^2 = 4p */
extern int modified_cornacchia(mpz_t x, mpz_t y, mpz_t D, mpz_t p);
/* return a class poly (Hilbert [type 1] or Weber [type 2]) */
extern UV poly_class_poly(IV D, mpz_t**T, int* type);
extern const char* poly_class_type_name(int type);
/* return a 0 terminated list of all D's sorted by degree */
extern IV* poly_class_degrees(int insert_1s);
/* List of class polynomial indices in order */
extern int* poly_class_nums(void);
/* Given a class poly index, return the degree and fill in (if not null):
* D the discriminant number
* T the polynomial coefficients
* type the poly type: 1 Hilber, 2 Weber
*/
extern UV poly_class_poly_num(int i, int *D, mpz_t**T, int* type);
#define BITS2DIGS(bits) ceil(bits/3.3219281)
#define DIGS2BITS(digs) ceil(digs*3.3219281)
extern void mpf_log(mpf_t logx, mpf_t x);
extern void mpf_exp(mpf_t expx, mpf_t x);
extern void mpf_pow(mpf_t powx, mpf_t b, mpf_t x);
extern void mpf_root(mpf_t rootx, mpf_t x, mpf_t n);
extern void mpf_agm(mpf_t r, mpf_t a, mpf_t b);
extern UV logint(mpz_t n, UV base);
#if defined(FUNC_isqrt) || defined(FUNC_is_perfect_square)
static UV isqrt(UV n) {
UV root;
#if BITS_PER_WORD == 32
if (n >= UVCONST(4294836225)) return UVCONST(65535);
#else
if (n >= UVCONST(18446744065119617025)) return UVCONST(4294967295);
#endif
root = (UV) sqrt((double)n);
while (root*root > n) root--;
while ((root+1)*(root+1) <= n) root++;
return root;
}
#endif
#if defined(FUNC_gcd_ui) || defined(FUNC_lcm_ui)
static UV gcd_ui(UV x, UV y) {
UV t;
if (y < x) { t = x; x = y; y = t; }
while (y > 0) {
t = y; y = x % y; x = t; /* y1 <- x0 % y0 ; x1 <- y0 */
}
return x;
}
#endif
#ifdef FUNC_lcm_ui
static UV lcm_ui(UV x, UV y) {
/* Can overflow if lcm(x,y) > 2^64 (e.g. two primes each > 2^32) */
return x * (y / gcd_ui(x,y));
}
#endif
#ifdef FUNC_is_perfect_square
/* Return 0 if n is not a perfect square. Set sqrtn to int(sqrt(n)) if so.
* See: http://mersenneforum.org/showpost.php?p=110896
*/
static int is_perfect_square(UV n, UV* sqrtn)
{
UV m;
m = n & 127;
if ((m*0x8bc40d7d) & (m*0xa1e2f5d1) & 0x14020a) return 0;
/* This cuts out another 80%: */
m = n % 63; if ((m*0x3d491df7) & (m*0xc824a9f9) & 0x10f14008) return 0;
/* m = n % 25; if ((m*0x1929fc1b) & (m*0x4c9ea3b2) & 0x51001005) return 0; */
m = isqrt(n);
if (n != m*m) return 0;
if (sqrtn != 0) *sqrtn = m;
return 1;
}
#endif
#ifdef FUNC_mpz_logn
static double mpz_logn(mpz_t n)
{
long exp;
double logn = mpz_get_d_2exp(&exp, n);
logn = log(logn) + (log(2) * exp);
return logn;
}
#endif
#ifdef FUNC_mpz_log2
static double mpz_log2(mpz_t n)
{
long exp;
double logn = mpz_get_d_2exp(&exp, n);
logn = exp + log(logn)/log(2);
return logn;
}
#endif
#endif
|