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

functions for mapping to the fundamental domain

Browse source
This commit is contained in:
Fredrik Johansson 2014-10-01 19:12:18 +02:00
parent 16abd31903
commit 99cafaa7cf
8 changed files with 624 additions and 0 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

11
acb_modular.h
View file

@ -118,6 +118,17 @@ void psl2z_randtest(psl2z_t g, flint_rand_t state, long bits);
void acb_modular_transform(acb_t w, const psl2z_t g, const acb_t z, long prec);
void acb_modular_fundamental_domain_approx_d(psl2z_t g,
double x, double y, double one_minus_eps);
void acb_modular_fundamental_domain_approx_arf(psl2z_t g,
const arf_t xx, const arf_t yy, const arf_t one_minus_eps, long prec);
void acb_modular_fundamental_domain_approx(acb_t w, psl2z_t g, const acb_t z,
const arf_t one_minus_eps, long prec);
int acb_modular_is_in_fundamental_domain(const acb_t z, const arf_t tol, long prec);
#ifdef __cplusplus
}
#endif

113
acb_modular/fundamental_domain_approx.c Normal file
View file

@ -0,0 +1,113 @@
/*=============================================================================
This file is part of ARB.
ARB is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.
ARB is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with ARB; if not, write to the Free Software
Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
=============================================================================*/
/******************************************************************************
Copyright (C) 2014 Fredrik Johansson
******************************************************************************/
#include "acb_modular.h"
static int
good_enough(const acb_t z, const arf_t one_minus_eps, long prec)
{
arf_t m;
int res;
if (arf_cmpabs_2exp_si(arb_midref(acb_realref(z)), -1) > 0)
return 0;
if (arf_cmpabs_2exp_si(arb_midref(acb_imagref(z)), 0) >= 0)
return 1;
arf_init(m);
arf_mul(m, arb_midref(acb_realref(z)), arb_midref(acb_realref(z)), prec, ARF_RND_DOWN);
arf_addmul(m, arb_midref(acb_imagref(z)), arb_midref(acb_imagref(z)), prec, ARF_RND_DOWN);
res = (arf_cmp(m, one_minus_eps) >= 0);
arf_clear(m);
return res;
}
void
acb_modular_fundamental_domain_approx(acb_t w, psl2z_t g, const acb_t z,
const arf_t one_minus_eps, long prec)
{
acb_t t;
psl2z_one(g);
/* we must be in the upper half-plane */
if (!acb_is_finite(z) || arf_sgn(arb_midref(acb_imagref(z))) <= 0)
{
acb_set(w, z);
return;
}
/* too large real-value shift */
if (arf_cmpabs_2exp_si(arb_midref(acb_realref(z)), prec) > 0)
{
acb_set(w, z);
return;
}
/* y >= 1: just shift x */
if (arf_cmpabs_2exp_si(arb_midref(acb_imagref(z)), 0) >= 0)
{
arf_get_fmpz(&g->b, arb_midref(acb_realref(z)), ARF_RND_NEAR);
acb_sub_fmpz(w, z, &g->b, prec);
fmpz_neg(&g->b, &g->b);
return;
}
acb_init(t);
/* try using doubles */
if (arf_cmpabs_2exp_si(arb_midref(acb_imagref(z)), -40) > 0 &&
arf_cmpabs_2exp_si(arb_midref(acb_realref(z)), 40) < 0)
{
double zred, zimd;
zred = arf_get_d(arb_midref(acb_realref(z)), ARF_RND_DOWN);
zimd = arf_get_d(arb_midref(acb_imagref(z)), ARF_RND_DOWN);
acb_modular_fundamental_domain_approx_d(g, zred, zimd,
arf_get_d(one_minus_eps, ARF_RND_UP));
acb_modular_transform(t, g, z, prec);
if (good_enough(t, one_minus_eps, prec))
{
acb_swap(w, t);
acb_clear(t);
return;
}
}
/* try with full precision */
acb_modular_fundamental_domain_approx_arf(g,
arb_midref(acb_realref(z)), arb_midref(acb_imagref(z)),
one_minus_eps, prec);
acb_modular_transform(t, g, z, prec);
acb_swap(w, t);
acb_clear(t);
}

111
acb_modular/fundamental_domain_approx_arf.c Normal file
View file

@ -0,0 +1,111 @@
/*=============================================================================
This file is part of ARB.
ARB is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.
ARB is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with ARB; if not, write to the Free Software
Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
=============================================================================*/
/******************************************************************************
Copyright (C) 2014 Fredrik Johansson
******************************************************************************/
#include "acb_modular.h"
void
acb_modular_fundamental_domain_approx_arf(psl2z_t g,
const arf_t xx, const arf_t yy, const arf_t one_minus_eps, long prec)
{
long i;
arf_t x, y, t;
fmpz_t n;
psl2z_one(g);
/* we must be in the upper half-plane */
if (!arf_is_finite(xx) || !arf_is_finite(yy) || arf_sgn(yy) <= 0)
return;
arf_init(x);
arf_init(y);
arf_init(t);
fmpz_init(n);
arf_set_round(x, xx, prec, ARF_RND_DOWN);
arf_set_round(y, yy, prec, ARF_RND_DOWN);
for (i = 0; i < 10 + prec / 4; i++)
{
/* too large shift */
if (arf_cmpabs_2exp_si(x, prec) > 0)
{
psl2z_one(g);
break;
}
if (fmpz_bits(&g->a) > prec || fmpz_bits(&g->b) > prec ||
fmpz_bits(&g->c) > prec || fmpz_bits(&g->d) > prec)
{
psl2z_one(g);
break;
}
/* shift */
if (arf_cmpabs_2exp_si(x, -1) > 0)
{
arf_get_fmpz(n, x, ARF_RND_NEAR);
arf_sub_fmpz(x, x, n, prec, ARF_RND_DOWN);
fmpz_submul(&g->a, &g->c, n);
fmpz_submul(&g->b, &g->d, n);
continue;
}
arf_mul(t, x, x, prec, ARF_RND_DOWN);
arf_addmul(t, y, y, prec, ARF_RND_DOWN);
/* inversion */
if (arf_cmp(t, one_minus_eps) < 0)
{
arf_div(x, x, t, prec, ARF_RND_DOWN);
arf_neg(x, x);
arf_div(y, y, t, prec, ARF_RND_DOWN);
fmpz_swap(&g->a, &g->c);
fmpz_swap(&g->b, &g->d);
fmpz_neg(&g->a, &g->a);
fmpz_neg(&g->b, &g->b);
continue;
}
/* we're good */
break;
}
if (fmpz_sgn(&g->c) < 0 || (fmpz_is_zero(&g->c) && fmpz_sgn(&g->d) < 0))
{
fmpz_neg(&g->a, &g->a);
fmpz_neg(&g->b, &g->b);
fmpz_neg(&g->c, &g->c);
fmpz_neg(&g->d, &g->d);
}
arf_clear(x);
arf_clear(y);
arf_clear(t);
fmpz_clear(n);
}

111
acb_modular/fundamental_domain_approx_d.c Normal file
View file

@ -0,0 +1,111 @@
/*=============================================================================
This file is part of ARB.
ARB is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.
ARB is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with ARB; if not, write to the Free Software
Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
=============================================================================*/
/******************************************************************************
Copyright (C) 2014 Fredrik Johansson
******************************************************************************/
#include "acb_modular.h"
static __inline__ int
d_is_ok(double x)
{
return (x > -1e15) && (x < 1e15);
}
void
acb_modular_fundamental_domain_approx_d(psl2z_t g,
double x, double y, double one_minus_eps)
{
double a, b, c, d, aa, bb, t;
int i;
a = d = 1;
b = c = 0;
for (i = 0; i < 20; i++)
{
if (!d_is_ok(x) || !d_is_ok(y) || !(y > 0.0))
{
psl2z_one(g);
return;
}
/* shift */
if (x < -0.5 || x > 0.5)
{
t = floor(x + 0.5);
x -= t;
a -= t * c;
b -= t * d;
/* too large to guarantee exactness */
if (!d_is_ok(a) || !d_is_ok(b))
{
psl2z_one(g);
return;
}
continue;
}
t = x*x + y*y;
/* can't divide by a too small number */
if (t < 1e-30)
{
psl2z_one(g);
break;
}
/* inversion */
if (t < one_minus_eps)
{
t = 1.0 / t;
x *= -t;
y *= t;
aa = a;
bb = b;
a = -c;
b = -d;
c = aa;
d = bb;
continue;
}
/* we're good */
break;
}
if (c < 0 || (c == 0 && d < 0))
{
a = -a;
b = -b;
c = -c;
d = -d;
}
fmpz_set_d(&g->a, a);
fmpz_set_d(&g->b, b);
fmpz_set_d(&g->c, c);
fmpz_set_d(&g->d, d);
}

69
acb_modular/is_in_fundamental_domain.c Normal file
View file

@ -0,0 +1,69 @@
/*=============================================================================
This file is part of ARB.
ARB is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.
ARB is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with ARB; if not, write to the Free Software
Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
=============================================================================*/
/******************************************************************************
Copyright (C) 2014 Fredrik Johansson
******************************************************************************/
#include "acb_modular.h"
int
acb_modular_is_in_fundamental_domain(const acb_t z, const arf_t tol, long prec)
{
arb_t t;
arb_init(t);
/* require re(w) + 1/2 >= 0 */
arb_set_ui(t, 1);
arb_mul_2exp_si(t, t, -1);
arb_add(t, t, acb_realref(z), prec);
arb_add_arf(t, t, tol, prec);
if (!arb_is_nonnegative(t))
{
arb_clear(t);
return 0;
}
/* require re(w) - 1/2 <= 0 */
arb_set_ui(t, 1);
arb_mul_2exp_si(t, t, -1);
arb_sub(t, acb_realref(z), t, prec);
arb_sub_arf(t, t, tol, prec);
if (!arb_is_nonpositive(t))
{
arb_clear(t);
return 0;
}
/* require |w| >= 1 - tol, i.e. |w| - 1 + tol >= 0 */
acb_abs(t, z, prec);
arb_sub_ui(t, t, 1, prec);
arb_add_arf(t, t, tol, prec);
if (!arb_is_nonnegative(t))
{
arb_clear(t);
return 0;
}
arb_clear(t);
return 1;
}

142
acb_modular/test/t-fundamental_domain_approx.c Normal file
View file

@ -0,0 +1,142 @@
/*=============================================================================
This file is part of ARB.
ARB is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.
ARB is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with ARB; if not, write to the Free Software
Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
=============================================================================*/
/******************************************************************************
Copyright (C) 2014 Fredrik Johansson
******************************************************************************/
#include "acb_modular.h"
int main()
{
long iter;
flint_rand_t state;
printf("fundamental_domain_approx....");
fflush(stdout);
flint_randinit(state);
for (iter = 0; iter < 100000; iter++)
{
fmpq_t x, y;
psl2z_t g;
arf_t one_minus_eps, tol;
acb_t z, w, w2;
arb_t t;
long prec;
fmpq_init(x);
fmpq_init(y);
psl2z_init(g);
acb_init(z);
acb_init(w);
acb_init(w2);
arf_init(one_minus_eps);
arf_init(tol);
arb_init(t);
/* pick an exact point in the upper half plane */
fmpq_randtest(x, state, 1 + n_randint(state, 10));
do {
fmpq_randtest(y, state, 1 + n_randint(state, 10));
} while (fmpz_sgn(fmpq_numref(y)) <= 0);
/* pick a tolerance */
arf_set_ui_2exp_si(tol, 1, -(long) n_randint(state, 10));
/* now increase the precision until convergence */
for (prec = 32; ; prec *= 2)
{
if (prec > 16384)
{
printf("FAIL (no convergence)\n");
printf("x = "); fmpq_print(x); printf("\n\n");
printf("y = "); fmpq_print(y); printf("\n\n");
printf("z = "); acb_printd(z, 50); printf("\n\n");
printf("w = "); acb_printd(w, 50); printf("\n\n");
printf("w2 = "); acb_printd(w2, 50); printf("\n\n");
printf("g = "); psl2z_print(g); printf("\n\n");
abort();
}
arb_set_fmpq(acb_realref(z), x, prec);
arb_set_fmpq(acb_imagref(z), y, prec);
arf_set_ui_2exp_si(one_minus_eps, 1, -prec / 4);
arf_sub_ui(one_minus_eps, one_minus_eps, 1, prec, ARF_RND_DOWN);
arf_neg(one_minus_eps, one_minus_eps);
acb_modular_fundamental_domain_approx(w, g, z, one_minus_eps, prec);
acb_modular_transform(w2, g, z, prec);
if (!psl2z_is_correct(g) || !acb_overlaps(w, w2))
{
printf("FAIL (incorrect transformation)\n");
printf("x = "); fmpq_print(x); printf("\n\n");
printf("y = "); fmpq_print(y); printf("\n\n");
printf("z = "); acb_printd(z, 50); printf("\n\n");
printf("w = "); acb_printd(w, 50); printf("\n\n");
printf("w2 = "); acb_printd(w2, 50); printf("\n\n");
printf("g = "); psl2z_print(g); printf("\n\n");
abort();
}
/* success */
if (acb_modular_is_in_fundamental_domain(w, tol, prec))
break;
}
/* check that g^(-1) * w contains x+yi */
psl2z_inv(g, g);
acb_modular_transform(w2, g, w, 2 + n_randint(state, 1000));
if (!arb_contains_fmpq(acb_realref(w2), x) ||
!arb_contains_fmpq(acb_imagref(w2), y))
{
printf("FAIL (inverse containment)\n");
printf("x = "); fmpq_print(x); printf("\n\n");
printf("y = "); fmpq_print(y); printf("\n\n");
printf("z = "); acb_printd(z, 50); printf("\n\n");
printf("w = "); acb_printd(w, 50); printf("\n\n");
printf("w2 = "); acb_printd(w2, 50); printf("\n\n");
printf("g = "); psl2z_print(g); printf("\n\n");
abort();
}
fmpq_clear(x);
fmpq_clear(y);
psl2z_clear(g);
acb_clear(z);
acb_clear(w);
acb_clear(w2);
arf_clear(one_minus_eps);
arf_clear(tol);
arb_clear(t);
}
flint_randclear(state);
flint_cleanup();
printf("PASS\n");
return EXIT_SUCCESS;
}

18
acb_modular/transform.c
View file

@ -46,6 +46,24 @@ void acb_modular_transform(acb_t w, const psl2z_t g, const acb_t z, long prec)
acb_inv(w, w, prec);
acb_neg(w, w);
}
else if (0)
{
acb_t t, u;
acb_init(t);
acb_init(u);
acb_set_fmpz(t, b);
acb_addmul_fmpz(t, z, a, prec);
acb_set_fmpz(u, d);
acb_addmul_fmpz(u, z, c, prec);
acb_div(w, t, u, prec);
acb_clear(t);
acb_clear(u);
}
else
{
/* (az+b)/(cz+d) = (re+im*i)/den where

49
doc/source/acb_modular.rst
View file

@ -82,6 +82,55 @@ Modular transformations
w = g z = \frac{az+b}{cz+d}.
.. function:: void acb_modular_fundamental_domain_approx_d(psl2z_t g, double x, double y, double one_minus_eps)
.. function:: void acb_modular_fundamental_domain_approx_arf(psl2z_t g, const arf_t x, const arf_t y, const arf_t one_minus_eps, long prec)
Attempts to determine a modular transformation *g* that maps the
complex number `x+yi` to the fundamental domain or just
slightly outside the fundamental domain, where the target tolerance
(not a strict bound) is specified by *one_minus_eps*.
The inputs are assumed to be finite numbers, with *y* positive.
Uses floating-point iteration, repeatedly applying either
the transformation `z \gets z + b` or `z \gets -1/z`. The iteration is
terminated if `|x| \le 1/2` and `x^2 + y^2 \ge 1 - \varepsilon` where
`1 - \varepsilon` is passed as *one_minus_eps*. It is also terminated
if too many steps have been taken without convergence, or if the numbers
end up too large or too small for the working precision.
The algorithm can fail to produce a satisfactory transformation.
The output *g* is always set to *some* correct modular transformation,
but it is up to the user to verify a posteriori that *g* maps `x+yi`
close enough to the fundamental domain.
.. function:: void acb_modular_fundamental_domain_approx(acb_t w, psl2z_t g, const acb_t z, const arf_t one_minus_eps, long prec)
Attempts to determine a modular transformation *g* that maps the
complex number `z` to the fundamental domain or just
slightly outside the fundamental domain, where the target tolerance
(not a strict bound) is specified by *one_minus_eps*. It also computes
the transformed value `w = gz`.
This function first tries to use
:func:`acb_modular_fundamental_domain_approx_d` and checks if the
result is acceptable. If this fails, it calls
:func:`acb_modular_fundamental_domain_approx_arf` with higher precision.
Finally, `w = gz` is evaluated by a single application of *g*.
The algorithm can fail to produce a satisfactory transformation.
The output *g* is always set to *some* correct modular transformation,
but it is up to the user to verify a posteriori that `w` is close enough
to the fundamental domain.
.. function:: int acb_modular_is_in_fundamental_domain(const acb_t z, const arf_t tol, long prec)
Returns nonzero if it is certainly true that `|z| \ge 1 - \varepsilon` and
`|\operatorname{Re}(z)| \le 1/2 + \varepsilon` where `\varepsilon` is
specified by *tol*. Returns zero if this is false or cannot be determined.
The Dedekind eta function
-------------------------------------------------------------------------------