hiro/arb
1
0
Fork 0
mirror of https://github.com/vale981/arb synced 2025-03-05 09:21:38 -05:00
Code Issues Projects Releases Packages Wiki Activity Actions

support partition function of bignum n (work in progress)

Browse source
This commit is contained in:
Fredrik Johansson 2014-03-03 16:32:29 +01:00
parent 2810e26428
commit 61e0c0c2ff
3 changed files with 69 additions and 46 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

6
partitions.h
View file

@ -35,9 +35,11 @@
extern "C" { extern "C" {
#endif #endif
void partitions_rademacher_bound(fmpr_t b, ulong n, ulong N); void partitions_rademacher_bound(fmpr_t b, const fmpz_t n, ulong N);
void partitions_hrr_sum_fmprb(fmprb_t x, ulong n, long N0, long N, int use_doubles); void partitions_hrr_sum_fmprb(fmprb_t x, const fmpz_t n, long N0, long N, int use_doubles);
void partitions_fmpz_fmpz(fmpz_t p, const fmpz_t n, int use_doubles);
void partitions_fmpz_ui(fmpz_t p, ulong n); void partitions_fmpz_ui(fmpz_t p, ulong n);

99
partitions/hrr_sum_fmprb.c
View file

@ -70,7 +70,7 @@ partitions_remainder_bound_log2(double n, double N)
} }
static long static long
partitions_needed_terms(ulong n) partitions_needed_terms(double n)
{ {
long N; long N;
for (N = 1; partitions_remainder_bound_log2(n, N) > 10; N++); for (N = 1; partitions_remainder_bound_log2(n, N) > 10; N++);
@ -115,7 +115,7 @@ log2_ceil(double x)
} }
static long static long
partitions_prec_bound(ulong n, long k, long N) partitions_prec_bound(double n, long k, long N)
{ {
long prec; long prec;
@ -234,24 +234,25 @@ sinh_cosh_divk_precomp(fmprb_t sh, fmprb_t ch, fmprb_t ex, long k, long prec)
void void
partitions_hrr_sum_fmprb(fmprb_t x, ulong n, long N0, long N, int use_doubles) partitions_hrr_sum_fmprb(fmprb_t x, const fmpz_t n, long N0, long N, int use_doubles)
{ {
trig_prod_t prod; trig_prod_t prod;
fmprb_t acc, C, t1, t2, t3, t4, exp1; fmprb_t acc, C, t1, t2, t3, t4, exp1;
fmpr_t bound; fmpr_t bound;
fmpz_t n24; fmpz_t n24;
long k, prec, res_prec, acc_prec, guard_bits; long k, prec, res_prec, acc_prec, guard_bits;
double Cd; double nd, Cd;
if (n <= 2) if (fmpz_cmp_ui(n, 2) <= 0)
{ {
fmprb_set_ui(x, FLINT_MAX(1, n)); abort();
return;
} }
nd = fmpz_get_d(n);
/* compute initial precision */ /* compute initial precision */
guard_bits = 2 * FLINT_BIT_COUNT(N) + 32; guard_bits = 2 * FLINT_BIT_COUNT(N) + 32;
prec = partitions_remainder_bound_log2(n, N0) + guard_bits; prec = partitions_remainder_bound_log2(nd, N0) + guard_bits;
prec = FLINT_MAX(prec, DOUBLE_PREC); prec = FLINT_MAX(prec, DOUBLE_PREC);
res_prec = acc_prec = prec; res_prec = acc_prec = prec;
@ -267,7 +268,7 @@ partitions_hrr_sum_fmprb(fmprb_t x, ulong n, long N0, long N, int use_doubles)
fmprb_zero(x); fmprb_zero(x);
/* n24 = 24n - 1 */ /* n24 = 24n - 1 */
fmpz_set_ui(n24, n); fmpz_set(n24, n);
fmpz_mul_ui(n24, n24, 24); fmpz_mul_ui(n24, n24, 24);
fmpz_sub_ui(n24, n24, 1); fmpz_sub_ui(n24, n24, 1);
@ -280,18 +281,18 @@ partitions_hrr_sum_fmprb(fmprb_t x, ulong n, long N0, long N, int use_doubles)
/* exp1 = exp(C) */ /* exp1 = exp(C) */
fmprb_exp(exp1, C, prec); fmprb_exp(exp1, C, prec);
Cd = PI * sqrt(24*n-1) / 6; Cd = PI * sqrt(24*nd-1) / 6;
for (k = N0; k <= N; k++) for (k = N0; k <= N; k++)
{ {
trig_prod_init(prod); trig_prod_init(prod);
arith_hrr_expsum_factored(prod, k, n % k); arith_hrr_expsum_factored(prod, k, fmpz_fdiv_ui(n, k));
if (prod->prefactor != 0) if (prod->prefactor != 0)
{ {
if (prec > MIN_PREC) if (prec > MIN_PREC)
prec = partitions_prec_bound(n, k, N); prec = partitions_prec_bound(nd, k, N);
prod->prefactor *= 4; prod->prefactor *= 4;
prod->sqrt_p *= 3; prod->sqrt_p *= 3;
@ -320,7 +321,7 @@ partitions_hrr_sum_fmprb(fmprb_t x, ulong n, long N0, long N, int use_doubles)
{ {
double xx, zz, xxerr; double xx, zz, xxerr;
xx = eval_trig_prod_d(prod) / (24*n - 1); xx = eval_trig_prod_d(prod) / (24*nd - 1);
zz = Cd / k; zz = Cd / k;
xx = xx * (cosh(zz) - sinh(zz) / zz); xx = xx * (cosh(zz) - sinh(zz) / zz);
@ -363,53 +364,69 @@ partitions_hrr_sum_fmprb(fmprb_t x, ulong n, long N0, long N, int use_doubles)
#define NUMBER_OF_SMALL_PARTITIONS 128 #define NUMBER_OF_SMALL_PARTITIONS 128
extern const unsigned int partitions_lookup[NUMBER_OF_SMALL_PARTITIONS]; extern const unsigned int partitions_lookup[NUMBER_OF_SMALL_PARTITIONS];
void
partitions_fmpz_fmpz(fmpz_t p, const fmpz_t n, int use_doubles)
{
if (fmpz_cmp_ui(n, NUMBER_OF_SMALL_PARTITIONS) < 0)
{
if (fmpz_sgn(n) < 0)
fmpz_zero(p);
else
fmpz_set_ui(p, partitions_lookup[fmpz_get_ui(n)]);
}
else
{
fmprb_t x;
long N;
fmprb_init(x);
N = partitions_needed_terms(fmpz_get_d(n));
partitions_hrr_sum_fmprb(x, n, 1, N, use_doubles);
if (!fmprb_get_unique_fmpz(p, x))
{
printf("not unique!\n");
fmprb_printd(x, 50);
printf("\n");
abort();
}
fmprb_clear(x);
}
}
void void
partitions_fmpz_ui(fmpz_t p, ulong n) partitions_fmpz_ui(fmpz_t p, ulong n)
{ {
fmprb_t x;
if (n < NUMBER_OF_SMALL_PARTITIONS) if (n < NUMBER_OF_SMALL_PARTITIONS)
{ {
fmpz_set_ui(p, partitions_lookup[n]); fmpz_set_ui(p, partitions_lookup[n]);
return;
} }
else
fmprb_init(x);
partitions_hrr_sum_fmprb(x, n, 1, partitions_needed_terms(n), 0);
if (!fmprb_get_unique_fmpz(p, x))
{ {
printf("not unique!\n"); fmpz_t t;
fmprb_printd(x, 50); fmpz_init(t);
printf("\n"); fmpz_set_ui(t, n);
abort(); partitions_fmpz_fmpz(p, t, 0);
fmpz_clear(t);
} }
fmprb_clear(x);
} }
void void
partitions_fmpz_ui_using_doubles(fmpz_t p, ulong n) partitions_fmpz_ui_using_doubles(fmpz_t p, ulong n)
{ {
fmprb_t x;
if (n < NUMBER_OF_SMALL_PARTITIONS) if (n < NUMBER_OF_SMALL_PARTITIONS)
{ {
fmpz_set_ui(p, partitions_lookup[n]); fmpz_set_ui(p, partitions_lookup[n]);
return;
} }
else
fmprb_init(x);
partitions_hrr_sum_fmprb(x, n, 1, partitions_needed_terms(n), 1);
if (!fmprb_get_unique_fmpz(p, x))
{ {
printf("not unique!\n"); fmpz_t t;
fmprb_printd(x, 50); fmpz_init(t);
printf("\n"); fmpz_set_ui(t, n);
abort(); partitions_fmpz_fmpz(p, t, 1);
fmpz_clear(t);
} }
fmprb_clear(x);
} }

10
partitions/rademacher_bound.c
View file

@ -38,15 +38,17 @@ _fmpr_sinh(fmpr_t y, const fmpr_t x, long prec)
/* Equation (1.8) in the paper */ /* Equation (1.8) in the paper */
void void
partitions_rademacher_bound(fmpr_t b, ulong n, ulong N) partitions_rademacher_bound(fmpr_t b, const fmpz_t n, ulong N)
{ {
fmpr_t A, B, C, t, u; fmpr_t A, B, C, t, u;
fmpz_t n1;
fmpr_init(A); fmpr_init(A);
fmpr_init(B); fmpr_init(B);
fmpr_init(C); fmpr_init(C);
fmpr_init(t); fmpr_init(t);
fmpr_init(u); fmpr_init(u);
fmpz_init(n1);
/* bound for 44*pi^2/(225*sqrt(3)) */ /* bound for 44*pi^2/(225*sqrt(3)) */
fmpr_set_si_2exp_si(A, 18695160, -24); fmpr_set_si_2exp_si(A, 18695160, -24);
@ -63,12 +65,13 @@ partitions_rademacher_bound(fmpr_t b, ulong n, ulong N)
/* B * sqrt(N/(n-1)) */ /* B * sqrt(N/(n-1)) */
fmpr_set_ui(t, N); fmpr_set_ui(t, N);
fmpr_div_ui(t, t, n - 1, FMPRB_RAD_PREC, FMPR_RND_UP); fmpz_sub_ui(n1, n, 1);
fmpr_div_fmpz(t, t, n1, FMPRB_RAD_PREC, FMPR_RND_UP);
fmpr_sqrt(t, t, FMPRB_RAD_PREC, FMPR_RND_UP); fmpr_sqrt(t, t, FMPRB_RAD_PREC, FMPR_RND_UP);
fmpr_mul(t, B, t, FMPRB_RAD_PREC, FMPR_RND_UP); fmpr_mul(t, B, t, FMPRB_RAD_PREC, FMPR_RND_UP);
/* sinh(C*sqrt(n)/N) */ /* sinh(C*sqrt(n)/N) */
fmpr_sqrt_ui(u, n, FMPRB_RAD_PREC, FMPR_RND_UP); fmpr_sqrt_fmpz(u, n, FMPRB_RAD_PREC, FMPR_RND_UP);
fmpr_div_ui(u, u, N, FMPRB_RAD_PREC, FMPR_RND_UP); fmpr_div_ui(u, u, N, FMPRB_RAD_PREC, FMPR_RND_UP);
fmpr_mul(u, C, u, FMPRB_RAD_PREC, FMPR_RND_UP); fmpr_mul(u, C, u, FMPRB_RAD_PREC, FMPR_RND_UP);
@ -83,5 +86,6 @@ partitions_rademacher_bound(fmpr_t b, ulong n, ulong N)
fmpr_clear(C); fmpr_clear(C);
fmpr_clear(t); fmpr_clear(t);
fmpr_clear(u); fmpr_clear(u);
fmpz_clear(n1);
} }