Skip to content

Commit

Permalink
proper cmpabs for acb, acf
Browse files Browse the repository at this point in the history
  • Loading branch information
@fredrik-johansson
fredrik-johansson committed Jul 20, 2024
1 parent 60d5397 commit 8645932
Show file tree
Hide file tree
Showing 2 changed files with 234 additions and 11 deletions.
103 changes: 95 additions & 8 deletions src/gr/acb.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -987,21 +987,108 @@ _gr_acb_cmp(int * res, const acb_t x, const acb_t y, const gr_ctx_t ctx)
}
}

int
_gr_arb_cmpabs(int * res, const arb_t x, const arb_t y, const gr_ctx_t ctx);

int
_gr_acb_cmpabs(int * res, const acb_t x, const acb_t y, const gr_ctx_t ctx)
{
acb_t t, u;
if (arb_is_zero(acb_imagref(x)) && arb_is_zero(acb_imagref(y)))
{
arb_srcptr a = acb_realref(x);
arb_srcptr c = acb_realref(y);

*t = *x;
*u = *y;
/* OK; ignores the context object */
return _gr_arb_cmpabs(res, a, c, ctx);
}
else
{
slong prec = ACB_CTX_PREC(ctx);
int status = GR_SUCCESS;

if (arf_sgn(arb_midref(acb_realref(t))) < 0)
ARF_NEG(arb_midref(acb_realref(t)));
arb_srcptr a = acb_realref(x);
arb_srcptr b = acb_imagref(x);
arb_srcptr c = acb_realref(y);
arb_srcptr d = acb_imagref(y);

mag_t xlo, xhi, ylo, yhi, t;

mag_init(xlo);
mag_init(xhi);
mag_init(ylo);
mag_init(yhi);
mag_init(t);

arb_get_mag_lower(xlo, a);
arb_get_mag_lower(t, b);
mag_mul_lower(xlo, xlo, xlo);
mag_mul_lower(t, t, t);
mag_add_lower(xlo, xlo, t);

arb_get_mag_lower(ylo, c);
arb_get_mag_lower(t, d);
mag_mul_lower(ylo, ylo, ylo);
mag_mul_lower(t, t, t);
mag_add_lower(ylo, ylo, t);

arb_get_mag(xhi, a);
arb_get_mag(t, b);
mag_mul(xhi, xhi, xhi);
mag_mul(t, t, t);
mag_add(xhi, xhi, t);

arb_get_mag(yhi, c);
arb_get_mag(t, d);
mag_mul(yhi, yhi, yhi);
mag_mul(t, t, t);
mag_add(yhi, yhi, t);

if (mag_cmp(xhi, ylo) < 0)
{
*res = -1;
status = GR_SUCCESS;
}
else if (mag_cmp(xlo, yhi) > 0)
{
*res = 1;
status = GR_SUCCESS;
}
else
{
arb_t t, u;

if (arf_sgn(arb_midref(acb_realref(u))) < 0)
ARF_NEG(arb_midref(acb_realref(u)));
arb_init(t);
arb_init(u);

return _gr_acb_cmp(res, t, u, ctx);
arb_mul(t, a, a, prec);
arb_addmul(t, b, b, prec);
arb_mul(u, c, c, prec);
arb_addmul(u, d, d, prec);

if ((arb_is_exact(t) && arb_is_exact(u)) || !arb_overlaps(t, u))
{
*res = arf_cmp(arb_midref(t), arb_midref(u));
status = GR_SUCCESS;
}
else
{
/* todo: worth it to do an exact computation? */
*res = 0;
status = GR_UNABLE;
}

arb_clear(t);
arb_clear(u);
}

mag_clear(xlo);
mag_clear(xhi);
mag_clear(ylo);
mag_clear(yhi);
mag_clear(t);

return status;
}
}

int
Expand Down
142 changes: 139 additions & 3 deletions src/gr/acf.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -666,13 +666,149 @@ _gr_acf_cmp(int * res, const acf_t x, const acf_t y, const gr_ctx_t ctx)
return GR_SUCCESS;
}

/* ignores ctx, so we can pass in the acf context */
int
_gr_arf_cmpabs(int * res, const arf_t x, const arf_t y, const gr_ctx_t ctx);

#include "double_extras.h"

int
_gr_acf_cmpabs(int * res, const acf_t x, const acf_t y, const gr_ctx_t ctx)
{
if (!arf_is_zero(acf_imagref(x)) || !arf_is_zero(acf_imagref(y)))
return GR_UNABLE;
arf_srcptr a = acf_realref(x);
arf_srcptr b = acf_imagref(x);
arf_srcptr c = acf_realref(y);
arf_srcptr d = acf_imagref(y);

if (arf_is_zero(b))
{
if (arf_is_zero(d))
return _gr_arf_cmpabs(res, a, c, ctx);
if (arf_is_zero(c))
return _gr_arf_cmpabs(res, a, d, ctx);
}

if (arf_is_zero(a))
{
if (arf_is_zero(d))
return _gr_arf_cmpabs(res, b, c, ctx);
if (arf_is_zero(c))
return _gr_arf_cmpabs(res, b, d, ctx);
}

if (arf_is_zero(c))
{
if (arf_is_zero(a))
return _gr_arf_cmpabs(res, b, d, ctx);
if (arf_is_zero(b))
return _gr_arf_cmpabs(res, a, d, ctx);
}

if (arf_is_zero(d))
{
if (arf_is_zero(a))
return _gr_arf_cmpabs(res, b, c, ctx);
if (arf_is_zero(b))
return _gr_arf_cmpabs(res, a, c, ctx);
}

if (ARF_IS_LAGOM(a) && ARF_IS_LAGOM(b) && ARF_IS_LAGOM(c) && ARF_IS_LAGOM(d))
{
slong aexp, bexp, cexp, dexp, xexp, yexp, exp;

aexp = arf_is_zero(a) ? WORD_MIN : ARF_EXP(a);
bexp = arf_is_zero(b) ? WORD_MIN : ARF_EXP(b);
cexp = arf_is_zero(c) ? WORD_MIN : ARF_EXP(c);
dexp = arf_is_zero(d) ? WORD_MIN : ARF_EXP(d);

/* 0.5 * 2^xexp <= |x| < sqrt(2) * 2^xexp */
xexp = FLINT_MAX(aexp, bexp);
/* 0.5 * 2^yexp <= |y| < sqrt(2) * 2^yexp */
yexp = FLINT_MAX(cexp, dexp);

if (xexp + 2 < yexp)
{
*res = -1;
return GR_SUCCESS;
}

if (xexp > yexp + 2)
{
*res = 1;
return GR_SUCCESS;
}

exp = FLINT_MAX(xexp, yexp);

double tt, xx = 0.0, yy = 0.0;
nn_srcptr xp;
slong xn;

if (aexp >= exp - 53)
{
ARF_GET_MPN_READONLY(xp, xn, a);
tt = d_mul_2exp_inrange(xp[xn - 1], aexp - exp - FLINT_BITS);
xx += tt * tt;
}

if (bexp >= exp - 53)
{
ARF_GET_MPN_READONLY(xp, xn, b);
tt = d_mul_2exp_inrange(xp[xn - 1], bexp - exp - FLINT_BITS);
xx += tt * tt;
}

if (cexp >= exp - 53)
{
ARF_GET_MPN_READONLY(xp, xn, c);
tt = d_mul_2exp_inrange(xp[xn - 1], cexp - exp - FLINT_BITS);
yy += tt * tt;
}

if (dexp >= exp - 53)
{
ARF_GET_MPN_READONLY(xp, xn, d);
tt = d_mul_2exp_inrange(xp[xn - 1], dexp - exp - FLINT_BITS);
yy += tt * tt;
}

if (xx < yy * 0.999999)
{
*res = -1;
return GR_SUCCESS;
}

if (xx * 0.999999 > yy)
{
*res = 1;
return GR_SUCCESS;
}
}

arf_struct s[5];

arf_init(s + 0);
arf_init(s + 1);
arf_init(s + 2);
arf_init(s + 3);
arf_init(s + 4);

arf_mul(s + 0, a, a, ARF_PREC_EXACT, ARF_RND_DOWN);
arf_mul(s + 1, b, b, ARF_PREC_EXACT, ARF_RND_DOWN);
arf_mul(s + 2, c, c, ARF_PREC_EXACT, ARF_RND_DOWN);
arf_mul(s + 3, d, d, ARF_PREC_EXACT, ARF_RND_DOWN);
arf_neg(s + 2, s + 2);
arf_neg(s + 3, s + 3);
arf_sum(s + 4, s, 4, 30, ARF_RND_DOWN);

*res = arf_sgn(s + 4);

arf_clear(s + 0);
arf_clear(s + 1);
arf_clear(s + 2);
arf_clear(s + 3);
arf_clear(s + 4);

*res = arf_cmpabs(acf_realref(x), acf_realref(y));
return GR_SUCCESS;
}

Expand Down

0 comments on commit 8645932

Please sign in to comment.