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