flint-sys 0.9.0

Bindings to the FLINT C library
Documentation 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103

/*
    Copyright (C) 2012, 2022 Fredrik Johansson

    This file is part of FLINT.

    FLINT is free software: you can redistribute it and/or modify it under
    the terms of the GNU Lesser General Public License (LGPL) as published
    by the Free Software Foundation; either version 3 of the License, or
    (at your option) any later version.  See <https://www.gnu.org/licenses/>.
*/

#include "arb.h"

/* decide |xm - ym| <= xr + yr */
int
arb_overlaps(const arb_t x, const arb_t y)
{
    arf_t t;
    mag_t a, b;
    int inexact, result;

    if (!arb_is_finite(x) || !arb_is_finite(y))
    {
        /* special cases: positive and negative infinity */
        if (arf_is_inf(arb_midref(x)) && mag_is_finite(arb_radref(x)))
        {
            if (mag_is_finite(arb_radref(y)) && (arf_is_finite(arb_midref(y)) || (arf_is_inf(arb_midref(y)) && !arf_equal(arb_midref(x), arb_midref(y)))))
                return 0;
        }

        if (arf_is_inf(arb_midref(y)) && mag_is_finite(arb_radref(y)))
        {
            if (mag_is_finite(arb_radref(x)) && (arf_is_finite(arb_midref(x)) || (arf_is_inf(arb_midref(x)) && !arf_equal(arb_midref(x), arb_midref(y)))))
                return 0;
        }

        return 1;
    }

    if (arf_equal(arb_midref(x), arb_midref(y)))
        return 1;

    arf_init(t);
    mag_init(a);
    mag_init(b);

    /* t = lower bound for |xm - ym| */
    inexact = arf_sub(t, arb_midref(x), arb_midref(y), 2 * MAG_BITS, ARF_RND_DOWN);

    /* u = lower bound for |xm - ym|; v = upper bound for xr + yr */
    arf_get_mag_lower(a, t);
    mag_add(b, arb_radref(x), arb_radref(y));

    if (mag_cmp(a, b) > 0)
    {
        result = 0;
    }
    else
    {
        /* a = upper bound for |xm - ym|; b = lower bound for xr + yr */
        arf_get_mag(a, t);

        if (inexact)
        {
            MAG_MAN(a)++;
            MAG_ADJUST_ONE_TOO_LARGE(a);
        }

        mag_add_lower(b, arb_radref(x), arb_radref(y));

        if (mag_cmp(a, b) < 0)
        {
            result = 1;
        }
        else
        {
            arf_struct u[4];

            if (arf_cmp(arb_midref(x), arb_midref(y)) >= 0)
            {
                arf_init_set_shallow(u + 0, arb_midref(x));
                arf_init_neg_shallow(u + 1, arb_midref(y));
            }
            else
            {
                arf_init_neg_shallow(u + 0, arb_midref(x));
                arf_init_set_shallow(u + 1, arb_midref(y));
            }

            arf_init_neg_mag_shallow(u + 2, arb_radref(x));
            arf_init_neg_mag_shallow(u + 3, arb_radref(y));

            arf_sum(t, u, 4, MAG_BITS, ARF_RND_DOWN);
            result = arf_sgn(t) <= 0;
        }
    }

    arf_clear(t);
    mag_clear(a);
    mag_clear(b);

    return result;
}