diff options
Diffstat (limited to 'src/cm/p1363.c')
| -rw-r--r-- | src/cm/p1363.c | 87 |
1 files changed, 45 insertions, 42 deletions
diff --git a/src/cm/p1363.c b/src/cm/p1363.c index 119b0cb..5ead675 100644 --- a/src/cm/p1363.c +++ b/src/cm/p1363.c @@ -79,15 +79,15 @@ void p1363_free(p1363_form_t ***forms, size_t nforms) { } } -static GEN p1363_func_F(GEN z) { +static GEN p1363_func_F(GEN z, long precision) { pari_sp ltop = avma; if (gequal0(z)) { return gcopy(gen_1); } - GEN sum = cgetc(BIGDEFAULTPREC * 5); - gel(sum, 1) = real_1(BIGDEFAULTPREC * 5); - gel(sum, 2) = real_0(BIGDEFAULTPREC * 5); + GEN sum = cgetc(precision * 2); + gel(sum, 1) = real_1(precision); + gel(sum, 2) = real_0(precision); pari_printf("initial sum = %Ps\n", sum); GEN last; @@ -99,8 +99,7 @@ static GEN p1363_func_F(GEN z) { GEN quota = divis(subii(j3, j), 2); GEN quotb = divis(addii(j3, j), 2); - GEN term = gadd(gpow(z, quota, BIGDEFAULTPREC), - gpow(z, quotb, BIGDEFAULTPREC)); + GEN term = gadd(gpow(z, quota, precision), gpow(z, quotb, precision)); pari_printf("%Ps %Ps \n", sum, term); if (mod2(j) == 0) { sum = gadd(sum, term); @@ -117,37 +116,37 @@ static GEN p1363_func_F(GEN z) { return gerepilecopy(ltop, sum); } -static GEN p1363_func_fzero(GEN D, p1363_form_t *form) { +static GEN p1363_func_fzero(GEN D, p1363_form_t *form, long precision) { pari_sp ltop = avma; - GEN upper = p1363_func_F(gneg(form->theta)); - GEN lower = p1363_func_F(gsqr(form->theta)); - GEN front = gpow(form->theta, gdivgs(gen_m1, 24), BIGDEFAULTPREC * 5); + GEN upper = p1363_func_F(gneg(form->theta), precision); + GEN lower = p1363_func_F(gsqr(form->theta), precision); + GEN front = gpow(form->theta, gdivgs(gen_m1, 24), precision); GEN result = gmul(gdiv(upper, lower), front); return gerepilecopy(ltop, result); } -static GEN p1363_func_fone(GEN D, p1363_form_t *form) { +static GEN p1363_func_fone(GEN D, p1363_form_t *form, long precision) { pari_sp ltop = avma; - GEN upper = p1363_func_F(form->theta); - GEN lower = p1363_func_F(gsqr(form->theta)); - GEN front = gpow(form->theta, gdivgs(gen_m1, 24), BIGDEFAULTPREC * 5); + GEN upper = p1363_func_F(form->theta, precision); + GEN lower = p1363_func_F(gsqr(form->theta), precision); + GEN front = gpow(form->theta, gdivgs(gen_m1, 24), precision); GEN result = gmul(gdiv(upper, lower), front); return gerepilecopy(ltop, result); } -static GEN p1363_func_ftwo(GEN D, p1363_form_t *form) { +static GEN p1363_func_ftwo(GEN D, p1363_form_t *form, long precision) { pari_sp ltop = avma; - GEN upper = p1363_func_F(gpowgs(form->theta, 4)); - GEN lower = p1363_func_F(gsqr(form->theta)); - GEN front = gmul(sqrti(gen_2), - gpow(form->theta, gdivgs(gen_1, 12), BIGDEFAULTPREC * 5)); + GEN upper = p1363_func_F(gpowgs(form->theta, 4), precision); + GEN lower = p1363_func_F(gsqr(form->theta), precision); + GEN front = + gmul(sqrti(gen_2), gpow(form->theta, gdivgs(gen_1, 12), precision)); GEN result = gmul(gdiv(upper, lower), front); @@ -287,20 +286,20 @@ static void p1363_N(GEN D, p1363_form_t *form) { form->N = gerepileupto(ltop, form->N); } -static void p1363_lambda(GEN D, p1363_form_t *form) { +static void p1363_lambda(GEN D, p1363_form_t *form, long precision) { pari_sp ltop = avma; - GEN pik = mulri(mppi(BIGDEFAULTPREC), stoi(form->K)); + GEN pik = mulri(mppi(precision), stoi(form->K)); GEN quot = divrs(pik, 24); form->lambda = gerepileupto(ltop, expIr(quot)); } -static void p1363_theta(GEN D, p1363_form_t *form) { +static void p1363_theta(GEN D, p1363_form_t *form, long precision) { pari_sp ltop = avma; - GEN upper = gadd(gneg(gsqrt(D, BIGDEFAULTPREC)), gmul(form->B, gen_I())); - GEN quot = gmul(gdiv(upper, form->A), mppi(BIGDEFAULTPREC)); + GEN upper = gadd(gneg(gsqrt(D, precision)), gmul(form->B, gen_I())); + GEN quot = gmul(gdiv(upper, form->A), mppi(precision)); - form->theta = gerepileupto(ltop, gexp(quot, BIGDEFAULTPREC)); + form->theta = gerepileupto(ltop, gexp(quot, precision)); } /** @@ -310,22 +309,23 @@ static void p1363_theta(GEN D, p1363_form_t *form) { * @param D * @param forms * @param nforms - * @return + * @return The pari precision required for W_D. */ -static long p1363_bit_precision(GEN D, p1363_form_t **forms, size_t nforms) { +long p1363_bit_precision(GEN D, p1363_form_t **forms, size_t nforms) { pari_sp ltop = avma; long v0 = 32; - GEN mul = divrr(mulrr(mppi(BIGDEFAULTPREC), sqrti(D)), mplog2(BIGDEFAULTPREC)); + GEN mul = + divrr(mulrr(mppi(BIGDEFAULTPREC), sqrti(D)), mplog2(BIGDEFAULTPREC)); GEN sum = gen_0; for (size_t i = 0; i < nforms; ++i) { sum = gadd(sum, ginv(forms[i]->A)); } - long result = v0 + itos(mpceil(mulrr(mul, sum))); + long result = v0 + itos(gceil(gmul(mul, sum))); avma = ltop; - return result; + return nbits2prec(result); } -GEN p1363_invariant(GEN D, p1363_form_t *form) { +GEN p1363_invariant(GEN D, p1363_form_t *form, long precision) { pari_printf("[A,B,C] = %Pi %Pi %Pi\n", form->A, form->B, form->C); pari_sp ltop = avma; @@ -342,9 +342,9 @@ GEN p1363_invariant(GEN D, p1363_form_t *form) { pari_printf("M = %Pi\n", form->M); p1363_N(D, form); pari_printf("N = %Pi\n", form->N); - p1363_lambda(D, form); + p1363_lambda(D, form, precision); pari_printf("lambda = %Ps\n", form->lambda); - p1363_theta(D, form); + p1363_theta(D, form, precision); pari_printf("theta = %Ps\n", form->theta); GEN G = gcdii(D, stoi(3)); @@ -352,23 +352,23 @@ GEN p1363_invariant(GEN D, p1363_form_t *form) { GEN bl = mulii(form->B, form->L); - GEN lmbl = gpow(form->lambda, negi(bl), BIGDEFAULTPREC); + GEN lmbl = gpow(form->lambda, negi(bl), precision); pari_printf("lmbl = %Ps\n", lmbl); GEN mi6 = gneg(gdiv(stoi(form->I), stoi(6))); - GEN powmi6 = gpow(gen_2, mi6, BIGDEFAULTPREC); + GEN powmi6 = gpow(gen_2, mi6, precision); pari_printf("powmi6 = %Ps\n", powmi6); GEN f = gen_0; switch (form->J) { case 0: - f = p1363_func_fzero(D, form); + f = p1363_func_fzero(D, form, precision); break; case 1: - f = p1363_func_fone(D, form); + f = p1363_func_fone(D, form, precision); break; case 2: - f = p1363_func_ftwo(D, form); + f = p1363_func_ftwo(D, form, precision); break; default: pari_err_DOMAIN("p1363_invariant", "J", ">", stoi(2), @@ -377,14 +377,14 @@ GEN p1363_invariant(GEN D, p1363_form_t *form) { pari_printf("f = %Ps\n", f); - GEN fK = gpow(f, stoi(form->K), BIGDEFAULTPREC); + GEN fK = gpow(f, stoi(form->K), precision); pari_printf("fK = %Ps\n", fK); GEN in = gmul(lmbl, powmi6); pari_printf("in = %Ps\n", in); GEN inner = gmul(gmul(form->N, in), fK); pari_printf("inner = %Ps\n", inner); - GEN result = gpow(inner, G, BIGDEFAULTPREC); + GEN result = gpow(inner, G, precision); pari_printf("result = %Ps\n", result); return gerepilecopy(ltop, result); @@ -395,14 +395,17 @@ GEN p1363_poly(GEN D, p1363_form_t **forms, size_t nforms) { long v = fetch_var(); name_var(v, "t"); + long precision = p1363_bit_precision(D, forms, nforms); + GEN terms = gtovec0(gen_0, nforms); for (size_t i = 0; i < nforms; ++i) { - gel(terms, i + 1) = gsub(pol_x(v), p1363_invariant(D, forms[i])); + gel(terms, i + 1) = + gsub(pol_x(v), p1363_invariant(D, forms[i], precision)); } GEN result = gen_1; for (size_t i = 0; i < nforms; ++i) { - gmulz(result, gel(terms, i + 1), result); + result = gmul(result, gel(terms, i + 1)); } return gerepilecopy(ltop, result); } |
