An Exclusive Integral

It is easy to show that if $0 \le x \le 1$ and $0 \le y \le 1$ , then $x \wedge y + x \wedge (1 - y) = 1.$

Let $I = \iint_C x \wedge y \: dx \: dy.$ Using the transformation $y \mapsto 1 - y$ , we get $I = \iint_C x \wedge (1 - y) \: dx \: dy.$ Then $\begin{aligned} 2I &= \iint_C x \wedge y \: dx \: dy + \iint_C x \wedge (1 - y) \: dx \: dy \\ &= \iint_C (x \wedge y + x \wedge (1 - y)) \: dx \: dy \\ &= \iint_C 1 \: dx \: dy \\ &= \frac{\pi}{4}, \end{aligned}$ so $I = \pi/8$ .

The closest integer to $2^{20} \cdot \pi/8$ is 411775.

The function $f(x, y) = x \wedge y$ is neither differentiable nor continuous, and has no closed form anti-derivative. Despite all of this, the function is, in fact, integrable.

Before we can make any progress, we need to get some idea of what this function looks like, at least in $[0, 1] \times [0, 1]$ . It's clear that at all such points, $x \wedge y$ is bounded in $[0, 1$ ), but we need to have some idea of what the distribution is like if we're going to integrate it. Here's a rendering of what the function looks like for each bit, starting with the $1/2$ bit, then the $1/4$ bit, and so on. The bluer a pixel is, the closer it is in value, for the acknowledged bits, to $1$ , while the whiter pixels are closer to $0$ in value.

Rendering

The rendering should make obvious what we can also derive analytically - that $x \wedge y$ has fractal behavior on the coordinate plane. Each binary bit below $2^0$ is set for exactly half of the area of $[0, 1] \times [0, 1]$ , but over increasingly finer-grained checkerboards for each bit. This enables us to answer a useful question:

$\int\limits_0^1 \int\limits_0^1 x \wedge y\:dx\:dy = \frac{1}{2}(\frac{1}{2^1} + \frac{1}{2^2} + \frac{1}{2^3} + ...) = \frac{1}{2}$

More directly, we can actually compute the integral of any square in any fractal layer, by computing the minimum value of $x \wedge y$ based on the square's coordinates and dimensions, and then using the geometric series we discovered to compute the total value across all remaining bits. However, the domain of integration - a circle - cannot be fully divided into exactly these fractal squares. This is okay, though, because we don't need to compute $I$ exactly - we just need to get it accurate to $20$ bits.

We can borrow a technique used by cartographers here to get $I$ arbitrarily accurate. We can pick an $n \geq 0$ , divide the unit square into $2^n \times 2^n$ sub-squares. For each one, we know how to compute the true integral, $V$ , over that square. We add a value to the total integral volume based on a simple condition for each square

• If the square is totally within the circular domain, add $V$

• If the square is partially within, and partially without, add $\frac{1}{2} V$

• Otherwise, add $0$

To get $20$ bits of accuracy, we need to ensure that the total volume of the boundary squares in our calculation does not exceed $2^{-20}$ . For a given $n$ , the approximate total volume of the boundary squares is:

$\frac{1}{2} \cdot \pi \cdot \left( \left( \frac{1}{2} \right)^2 - \left( \frac{1}{2} - \frac{1}{2^n} \right)^2 \right)$

To get this expression below $2^{-20}$ , $n$ needs to be $21$ or greater. However, we can't evaluate $2^{21} \times 2^{21}$ squares directly - it would take days, if not weeks, to do $2^{42}$ calculations like that. Fortunately, we don't have to - we can, instead, recursively evaluate the unit square $[0, 1] \times [0, 1]$ , splitting it into four equal quadrants for further evaluation when the square is not strictly within or without the circle. In this way, we will make a number of calculations that is linear in the number of boundary squares, which is linear in $2^n$ , rendering the calculation tractable. Here's a rendering of the recursive breakdown, where red overlays are outside, green are inside, and yellow are border:

Rendering

Running the code below up through $n = 21$ (R = 22 in the code) yields the answer $411775$ , though the error bounds cross the fractional $0.5$ mark, so we can't be certain the rounded answer is correct. Running an additional step, with $n = 22$ , gives us tight enough error bounds to be completely confident in the answer, $\boxed{411775}$ .

 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

import java.util.ArrayList;
import java.util.List;

public class Exclusive {

    private static final int R = 24;
    private static final int P = (1 << 20);

    private static double border = 0.0;
    private static double volume(long numx, long numy, int power, int r) {
        if (r == 0) {
            double fv = fullVolume(numx, numy, power);
            border += fv / 2.0;
            return fv / 2.0;
        }

        long x2 = numx * 2;
        long y2 = numy * 2;
        int p2 = power + 1;

        List<Boolean> bools = new ArrayList<Boolean>();
        bools.add(inCircle(x2, y2, p2)); // lower-left
        bools.add(inCircle(x2, y2 + 2, p2)); // upper-left
        bools.add(inCircle(x2 + 2, y2, p2)); // lower-right
        bools.add(inCircle(x2 + 2, y2 + 2, p2)); // upper-right
        bools.add(inCircle(x2 + 1, y2 + 1, p2)); // center

        double ret;
        if (bools.stream().allMatch(x -> x)) {
            ret = fullVolume(numx, numy, power);
        } else if (bools.stream().allMatch(x -> !x)) {
            ret = 0.0;
        } else {
            double v1 = volume(x2, y2, p2, r-1);
            double v2 = volume(x2 + 1, y2, p2, r-1);
            double v3 = volume(x2, y2 + 1, p2, r-1);
            double v4 = volume(x2 + 1, y2 + 1, p2, r-1);
            ret = v1 + v2 + v3 + v4;
        }
        return ret;
    }

    private static double fullVolume(long numx, long numy, int power) {
        long num = 2 * (numx ^ numy) + 1;
        double den = Math.pow(2.0, 3*power + 1);
        return num / den;
    }

    private static boolean inCircle(long numx, long numy, int power) {
        // (x-1/2)^2 + (y-1/2)^2 < (1/2)^2
        // (nx/d - 1/2)^2 + (yx/d - 1/2)^2 < 1/4
        // (2nx - d)^2 + (2yx - d)^2 < d^2
        long x1 = 2*numx - (1L << power);
        long y1 = 2*numy - (1L << power);
        boolean ret = (x1*x1 + y1*y1) < (1L << 2*power);
        return ret;
    }

    public static void main(String[] args) {
        for (int r = 1; r <= R; r++) {
            border = 0.0;
            double v = volume(0, 0, 0, r);
            System.out.println(v * P + " +- " + border * P);
        }
    }

}