1 1.1 mrg /* teta -- test file for the Dedekind eta function. 2 1.1 mrg 3 1.1 mrg Copyright (C) 2022 INRIA 4 1.1 mrg 5 1.1 mrg This file is part of GNU MPC. 6 1.1 mrg 7 1.1 mrg GNU MPC is free software; you can redistribute it and/or modify it under 8 1.1 mrg the terms of the GNU Lesser General Public License as published by the 9 1.1 mrg Free Software Foundation; either version 3 of the License, or (at your 10 1.1 mrg option) any later version. 11 1.1 mrg 12 1.1 mrg GNU MPC is distributed in the hope that it will be useful, but WITHOUT ANY 13 1.1 mrg WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS 14 1.1 mrg FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for 15 1.1 mrg more details. 16 1.1 mrg 17 1.1 mrg You should have received a copy of the GNU Lesser General Public License 18 1.1 mrg along with this program. If not, see http://www.gnu.org/licenses/ . 19 1.1 mrg */ 20 1.1 mrg 21 1.1 mrg #include <math.h> 22 1.1 mrg #include "mpc-tests.h" 23 1.1 mrg 24 1.1 mrg static void 25 1.1 mrg mpcb_j_err (mpcb_ptr j, mpc_srcptr z, unsigned long int err_re, 26 1.1 mrg unsigned long int err_im) 27 1.1 mrg /* Compute the j-function in z; err_re and err_im have the same meaning 28 1.1 mrg as in mpcb_eta_err. */ 29 1.1 mrg { 30 1.1 mrg const mpfr_prec_t p = mpc_get_prec (z); 31 1.1 mrg mpc_t z2; 32 1.1 mrg mpcb_t eta, eta2, sixteen; 33 1.1 mrg 34 1.1 mrg mpc_init2 (z2, p); 35 1.1 mrg mpcb_init (eta); 36 1.1 mrg mpcb_init (eta2); 37 1.1 mrg mpcb_init (sixteen); 38 1.1 mrg 39 1.1 mrg /* Compute f1 (z) = eta (z/2) / eta (z). */ 40 1.1 mrg mpcb_eta_err (eta, z, err_re, err_im); 41 1.1 mrg mpc_div_2ui (z2, z, 1, MPC_RNDNN); 42 1.1 mrg mpcb_eta_err (eta2, z2, err_re, err_im); 43 1.1 mrg mpcb_div (eta, eta2, eta); 44 1.1 mrg 45 1.1 mrg /* Compute gamma2 = (f1^24 + 16) / f1^8. */ 46 1.1 mrg mpcb_set_ui_ui (sixteen, 16, 0, p); 47 1.1 mrg mpcb_pow_ui (eta, eta, 8); 48 1.1 mrg mpcb_pow_ui (eta2, eta, 3); 49 1.1 mrg mpcb_add (eta2, eta2, sixteen); 50 1.1 mrg mpcb_div (eta2, eta2, eta); 51 1.1 mrg 52 1.1 mrg /* Compute j = gamma2^3. */ 53 1.1 mrg mpcb_pow_ui (j, eta2, 3); 54 1.1 mrg 55 1.1 mrg mpc_clear (z2); 56 1.1 mrg mpcb_clear (eta); 57 1.1 mrg mpcb_clear (eta2); 58 1.1 mrg mpcb_clear (sixteen); 59 1.1 mrg } 60 1.1 mrg 61 1.1 mrg 62 1.1 mrg static int 63 1.1 mrg test_eta (void) 64 1.1 mrg { 65 1.1 mrg const mpfr_prec_t p = 300; 66 1.1 mrg mpc_t z, eta; 67 1.1 mrg mpcb_t j; 68 1.1 mrg mpfr_t fuzz; 69 1.1 mrg mpz_t re_z, tmp; 70 1.1 mrg long int re, im; 71 1.1 mrg int ok; 72 1.1 mrg 73 1.1 mrg mpc_init2 (z, p); 74 1.1 mrg mpc_init2 (eta, p); 75 1.1 mrg mpcb_init (j); 76 1.1 mrg mpfr_init2 (fuzz, 2*p); 77 1.1 mrg 78 1.1 mrg mpfr_set_si (mpc_realref (z), -1, MPFR_RNDN); 79 1.1 mrg mpfr_set_ui (mpc_imagref (z), 163, MPFR_RNDD); 80 1.1 mrg mpfr_sqrt (mpc_imagref (z), mpc_imagref (z), MPFR_RNDD); 81 1.1 mrg mpc_div_2ui (z, z, 1, MPC_RNDNN); 82 1.1 mrg 83 1.1 mrg /* Check whether mpc_eta_fund avoids an infinite loop. */ 84 1.1 mrg mpc_eta_fund (eta, z, MPC_RNDNN); 85 1.1 mrg 86 1.1 mrg /* The error is bounded by 1/2 ulp in the real and 3 ulp in the 87 1.1 mrg imaginary part, see algorithms.tex. */ 88 1.1 mrg mpcb_j_err (j, z, 1, 6); 89 1.1 mrg 90 1.1 mrg /* Check whether j ((-1+sqrt(-163))/2) equals -262537412640768000. 91 1.1 mrg Rounding is impossible since the result is exact, and the imaginary 92 1.1 mrg part is 0; for a quick and dirty check, add the non-representable 93 1.1 mrg number 0.1 + 1.1 I and use the precisions that just work. */ 94 1.1 mrg mpfr_set_ui (fuzz, 1, MPFR_RNDN); 95 1.1 mrg mpfr_div_ui (fuzz, fuzz, 10, MPFR_RNDN); 96 1.1 mrg mpfr_add (mpc_realref (j->c), mpc_realref (j->c), fuzz, MPFR_RNDN); 97 1.1 mrg mpfr_add (mpc_imagref (j->c), mpc_imagref (j->c), fuzz, MPFR_RNDN); 98 1.1 mrg mpfr_add_ui (mpc_imagref (j->c), mpc_imagref (j->c), 1, MPFR_RNDN); 99 1.1 mrg ok = mpcb_can_round (j, 291, 234, MPC_RNDNN); 100 1.1 mrg mpcb_round (z, j, MPC_RNDNN); 101 1.1 mrg mpz_init (re_z); 102 1.1 mrg mpz_init_set_str (tmp, "-262537412640768000", 10); 103 1.1 mrg mpfr_get_z (re_z, mpc_realref (z), MPFR_RNDN); 104 1.1 mrg im = mpfr_get_si (mpc_imagref (z), MPFR_RNDN); 105 1.1 mrg ok &= (!mpz_cmp (re_z, tmp) && im == 1); 106 1.1 mrg if (!ok) { 107 1.1 mrg printf ("Error for -163:\n"); 108 1.1 mrg MPC_OUT (z); 109 1.1 mrg mpz_out_str (stdout, 10, re_z); 110 1.1 mrg printf ("\n"); 111 1.1 mrg } 112 1.1 mrg mpz_clear (tmp); 113 1.1 mrg mpz_clear (re_z); 114 1.1 mrg 115 1.1 mrg /* Check whether mpc_eta_fund (I) avoids an infinite loop. */ 116 1.1 mrg mpc_set_ui_ui (z, 0, 1, MPC_RNDNN); 117 1.1 mrg mpc_eta_fund (eta, z, MPC_RNDNN); 118 1.1 mrg 119 1.1 mrg /* Check whether j (I) equals 1728. */ 120 1.1 mrg mpcb_j_err (j, z, 0, 0); 121 1.1 mrg mpfr_add (mpc_realref (j->c), mpc_realref (j->c), fuzz, MPFR_RNDN); 122 1.1 mrg mpfr_add (mpc_imagref (j->c), mpc_imagref (j->c), fuzz, MPFR_RNDN); 123 1.1 mrg mpfr_add_ui (mpc_imagref (j->c), mpc_imagref (j->c), 1, MPFR_RNDN); 124 1.1 mrg ok &= mpcb_can_round (j, 292, 282, MPC_RNDNN); 125 1.1 mrg mpcb_round (z, j, MPC_RNDNN); 126 1.1 mrg re = mpfr_get_si (mpc_realref (z), MPFR_RNDN); 127 1.1 mrg im = mpfr_get_si (mpc_imagref (z), MPFR_RNDN); 128 1.1 mrg ok &= (re == 1728 && im == 1); 129 1.1 mrg if (!ok) { 130 1.1 mrg printf ("Error for -4:\n"); 131 1.1 mrg MPC_OUT (z); 132 1.1 mrg printf ("%li\n", re); 133 1.1 mrg } 134 1.1 mrg 135 1.1 mrg mpc_clear (eta); 136 1.1 mrg mpc_clear (z); 137 1.1 mrg mpcb_clear (j); 138 1.1 mrg mpfr_clear (fuzz); 139 1.1 mrg 140 1.1 mrg return !ok; 141 1.1 mrg } 142 1.1 mrg 143 1.1 mrg 144 1.1 mrg int 145 1.1 mrg main (void) 146 1.1 mrg { 147 1.1 mrg return test_eta (); 148 1.1 mrg } 149 1.1 mrg 150
Indexes created Sat Feb 28 05:31:39 UTC 2026