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
mpcb_j_err(mpcb_ptr j,mpc_srcptr z,unsigned long int err_re,unsigned long int err_im)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
test_eta(void)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
main(void)145 main (void)
146 {
147 return test_eta ();
148 }
149
150