Problem 945: XOR-Equation C

Project Euler Problem 945 Solution

EulerSolve provides an optimized solution for Project Euler Problem 945, XOR-Equation C, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary For the bound $N=10^7$, we must count all pairs $(a,b)$ with $0\le a\le b\le N$ for which the carryless quadratic form $$a\otimes a \oplus (2\otimes a\otimes b) \oplus b\otimes b$$ is itself a carryless square. Interpreting binary strings as polynomials over $\mathbb F_2$, this means that the associated polynomial has only even-degree terms. A direct search over all $\Theta(N^2)$ pairs is impossible at this scale. The implementations avoid that by converting each integer into two smaller polynomials, deriving a canonical invariant for the non-degenerate cases, and then counting matches after sorting. Mathematical Approach Associate to every nonnegative integer $n=\sum n_i2^i$ the polynomial $$P_n(x)=\sum n_i x^i\in\mathbb F_2[x].$$ Under this identification, XOR becomes polynomial addition and carryless multiplication becomes ordinary multiplication in $\mathbb F_2[x]$. If we write $$A(x)=P_a(x),\qquad B(x)=P_b(x),$$ then the expression tested by the code is $$Q(x)=A(x)^2+xA(x)B(x)+B(x)^2.$$ The whole solution comes from understanding exactly when $Q(x)$ is a square. Squares over $\mathbb F_2$ In characteristic 2 we have $$\left(\sum c_i x^i\right)^2=\sum c_i x^{2i},$$ because every mixed term appears twice and therefore cancels. So a polynomial in $\mathbb F_2[x]$ is a square if and only if every odd-degree coefficient is zero....

Detailed mathematical approach

Problem Summary

For the bound $N=10^7$, we must count all pairs $(a,b)$ with $0\le a\le b\le N$ for which the carryless quadratic form

$$a\otimes a \oplus (2\otimes a\otimes b) \oplus b\otimes b$$

is itself a carryless square. Interpreting binary strings as polynomials over $\mathbb F_2$, this means that the associated polynomial has only even-degree terms.

A direct search over all $\Theta(N^2)$ pairs is impossible at this scale. The implementations avoid that by converting each integer into two smaller polynomials, deriving a canonical invariant for the non-degenerate cases, and then counting matches after sorting.

Mathematical Approach

Associate to every nonnegative integer $n=\sum n_i2^i$ the polynomial

$$P_n(x)=\sum n_i x^i\in\mathbb F_2[x].$$

Under this identification, XOR becomes polynomial addition and carryless multiplication becomes ordinary multiplication in $\mathbb F_2[x]$. If we write

$$A(x)=P_a(x),\qquad B(x)=P_b(x),$$

then the expression tested by the code is

$$Q(x)=A(x)^2+xA(x)B(x)+B(x)^2.$$

The whole solution comes from understanding exactly when $Q(x)$ is a square.

Squares over $\mathbb F_2$

In characteristic 2 we have

$$\left(\sum c_i x^i\right)^2=\sum c_i x^{2i},$$

because every mixed term appears twice and therefore cancels. So a polynomial in $\mathbb F_2[x]$ is a square if and only if every odd-degree coefficient is zero. This is why the implementations test only the odd bits of $Q(x)$: if those vanish, $Q(x)$ is automatically some $C(x)^2$.

Even and odd bit polynomials

Introduce $t=x^2$ and split each polynomial into its even and odd positions:

$$A(x)=E_a(t)+x\,O_a(t),\qquad B(x)=E_b(t)+x\,O_b(t),$$

with $E_a,O_a,E_b,O_b\in\mathbb F_2[t]$. Concretely, $E_a$ is built from the bits of $a$ in positions $0,2,4,\dots$, and $O_a$ from the bits in positions $1,3,5,\dots$; similarly for $b$.

Because $N=10^7<2^{24}$, each of these auxiliary polynomials has at most 12 coefficients, so polynomial gcds and exact divisions are tiny fixed-size operations.

The odd-part identity

The square terms $A(x)^2$ and $B(x)^2$ already contain only even powers. Therefore the odd powers of $Q(x)$ come entirely from the middle term $xA(x)B(x)$. Expanding with $t=x^2$,

$$A(x)B(x)=E_aE_b+x(E_aO_b+O_aE_b)+t\,O_aO_b,$$

$$xA(x)B(x)=xE_aE_b+t(E_aO_b+O_aE_b)+xt\,O_aO_b.$$

The odd part of $Q(x)$ is therefore

$$x\bigl(E_a(t)E_b(t)+t\,O_a(t)O_b(t)\bigr).$$

Hence $(a,b)$ is valid exactly when

$$E_a(t)E_b(t)=t\,O_a(t)O_b(t).$$

This identity is the central mathematical object in the solution.

Canonical classes for the general case

Assume first that all four polynomials $E_a,O_a,E_b,O_b$ are nonzero. Then the criterion can be rewritten as

$$\frac{E_a(t)}{t\,O_a(t)}=\frac{O_b(t)}{E_b(t)}$$

inside the rational-function field $\mathbb F_2(t)$. The implementations normalize the two sides separately by cancelling polynomial gcds:

$$\left(\frac{E_a}{d_a},\frac{tO_a}{d_a}\right),\qquad d_a=\gcd(E_a,tO_a),$$

$$\left(\frac{O_b}{d_b},\frac{E_b}{d_b}\right),\qquad d_b=\gcd(O_b,E_b).$$

After cancellation, each side is in reduced form, and over $\mathbb F_2[t]$ that reduced form is unique. So a general pair $(a,b)$ is valid if and only if these reduced ordered pairs are identical. This turns an equation in two variables into a key-matching problem.

Degenerate branches

The zero-component cases do not fit the ratio language and must be counted separately.

If $E_a=0$ but $O_a\neq 0$, then

$$E_aE_b=t\,O_aO_b$$

forces $O_b=0$. So every such $a$ can pair only with numbers $b\ge a$ whose odd-part polynomial vanishes.

If $O_a=0$ but $E_a\neq 0$, then the same identity forces $E_b=0$. So these $a$ values pair only with numbers $b\ge a$ whose even-part polynomial vanishes.

The last special case is $a=0$. Then

$$Q(x)=B(x)^2,$$

which is always a square, so every $b\in[0,N]$ contributes a valid pair $(0,b)$.

Worked example

Take $a=3$ and $b=6$. In binary,

$$a=11_2,\qquad b=110_2,$$

$$A(x)=1+x,\qquad B(x)=x+x^2.$$

With $t=x^2$, the split is

$$E_a=1,\ O_a=1,\qquad E_b=t,\ O_b=1.$$

The criterion becomes

$$E_aE_b=1\cdot t=t,\qquad t\,O_aO_b=t\cdot 1\cdot 1=t,$$

so the pair is valid.

Direct expansion confirms it:

$$A(x)^2=1+x^2,\qquad xA(x)B(x)=x^2+x^4,\qquad B(x)^2=x^2+x^4,$$

hence

$$Q(x)=1+x^2,$$

which has only even powers and is therefore a square.

From algebra to counting

Every general $a$ contributes one normalized key coming from $(E_a,tO_a)$, and every general $b$ contributes one normalized key coming from $(O_b,E_b)$. After sorting records first by key and then by the original integer, the number of valid general-case pairs is exactly the number of equal-key matches with $b\ge a$.

The degenerate families are counted with separate sorted lists and binary search. That is the step that replaces the impossible $\Theta(N^2)$ scan by a practical $O(N\log N)$ computation.

How the Code Works

Building polynomial records

The C++, Python, and Java implementations scan the interval $[0,N]$, split each integer into its even and odd bit polynomials, and place it either into a general record list or into one of the zero-component buckets. Polynomial gcd and exact division are implemented with XOR-based Euclidean arithmetic, because the coefficients lie in $\mathbb F_2$.

Handling the special families

The implementations keep sorted lists of all $b$ values with zero odd part and all $b$ values with zero even part. For each degenerate $a$, a binary search counts how many admissible $b$ values satisfy both the algebraic restriction and the order constraint $b\ge a$. The unconditional family $a=0$ contributes $N+1$ pairs immediately.

Merging equal-key groups

The general records are sorted by normalized key and then by value. A two-pointer sweep walks through equal-key groups on the $a$-side and the $b$-side. Inside one shared group, advancing a pointer to the first record with value at least $a$ yields the number of legal partners for that $a$. This is just a merge-style count on already reduced algebraic classes.

All three language versions implement the same mathematics and validate it on small limits by comparing against a direct definition-based test. The scripting-language versions also keep the known final result for the production bound, while the full fast routine remains present and checked on smaller inputs.

Complexity Analysis

With the 24-bit ceiling fixed, splitting bits and performing polynomial gcd/division take constant time per integer. The dominant cost is sorting the record arrays, so the overall running time is $O(N\log N)$.

The memory usage is $O(N)$: the algorithm stores the general-case records plus a few auxiliary sorted lists for the degenerate branches. Compared with the naive $O(N^2)$ pair test, this is the decisive improvement.

Footnotes and References

Problem page: https://projecteuler.net/problem=945
Carry-less product: Wikipedia - Carry-less product
Finite field arithmetic: Wikipedia - Finite field arithmetic
Polynomial greatest common divisor: Wikipedia - Polynomial greatest common divisor
Frobenius endomorphism: Wikipedia - Frobenius endomorphism

Problem 945 source code

C++

#include <algorithm>
#include <cassert>
#include <cstdint>
#include <iostream>
#include <vector>

namespace {

using u32 = std::uint32_t;
using u64 = std::uint64_t;

constexpr int kHalfBits = 12;
constexpr int kKeyShift = 13;
constexpr int kValueBits = 24;
constexpr u64 kValueMask = (1ULL << kValueBits) - 1ULL;

inline int poly_degree(u32 p) {
    return p ? (31 - __builtin_clz(p)) : -1;
}

u32 poly_mod(u32 a, u32 b) {
    assert(b != 0);
    const int db = poly_degree(b);
    while (a != 0) {
        const int da = poly_degree(a);
        if (da < db) {
            break;
        }
        a ^= (b << (da - db));
    }
    return a;
}

u32 poly_gcd(u32 a, u32 b) {
    while (b != 0) {
        u32 r = poly_mod(a, b);
        a = b;
        b = r;
    }
    return a;
}

u32 poly_div_exact(u32 a, u32 b) {
    assert(b != 0);
    u32 q = 0;
    const int db = poly_degree(b);
    while (a != 0) {
        const int da = poly_degree(a);
        assert(da >= db);
        const int s = da - db;
        q ^= (1U << s);
        a ^= (b << s);
    }
    return q;
}

inline void split_even_odd_bits(u32 n, u32& even_part, u32& odd_part) {
    even_part = 0;
    odd_part = 0;
    for (int i = 0; i < kHalfBits; ++i) {
        even_part |= ((n >> (2 * i)) & 1U) << i;
        odd_part |= ((n >> (2 * i + 1)) & 1U) << i;
    }
}

inline u32 pack_key(u32 x, u32 y) {
    return (y << kKeyShift) | x;
}

inline u64 pack_record(u32 key, u32 value) {
    return (static_cast<u64>(key) << kValueBits) | static_cast<u64>(value);
}

inline u32 record_key(u64 record) {
    return static_cast<u32>(record >> kValueBits);
}

inline u32 record_value(u64 record) {
    return static_cast<u32>(record & kValueMask);
}

u64 count_ge(const std::vector<u32>& sorted_values, u32 threshold) {
    auto it = std::lower_bound(sorted_values.begin(), sorted_values.end(), threshold);
    return static_cast<u64>(sorted_values.end() - it);
}

u64 solve_fast(u32 N) {
    std::vector<u64> a_general;
    std::vector<u64> b_general;
    std::vector<u32> b_v_zero;
    std::vector<u32> b_u_zero;
    std::vector<u32> a_x_zero;
    std::vector<u32> a_y_zero;

    a_general.reserve(static_cast<size_t>(N));
    b_general.reserve(static_cast<size_t>(N));
    b_v_zero.reserve(5000);
    b_u_zero.reserve(5000);
    a_x_zero.reserve(5000);
    a_y_zero.reserve(5000);

    for (u32 b = 0; b <= N; ++b) {
        u32 U = 0;
        u32 V = 0;
        split_even_odd_bits(b, U, V);

        if (V == 0) {
            b_v_zero.push_back(b);
        }
        if (U == 0) {
            b_u_zero.push_back(b);
        }

        if (b == 0 || U == 0 || V == 0) {
            continue;
        }

        const u32 g = poly_gcd(U, V);
        const u32 x = poly_div_exact(V, g);
        const u32 y = poly_div_exact(U, g);
        b_general.push_back(pack_record(pack_key(x, y), b));
    }

    for (u32 a = 1; a <= N; ++a) {
        u32 X = 0;
        u32 Y = 0;
        split_even_odd_bits(a, X, Y);

        if (X == 0) {
            if (Y != 0) {
                a_x_zero.push_back(a);
            }
            continue;
        }
        if (Y == 0) {
            a_y_zero.push_back(a);
            continue;
        }

        const u32 uy = (Y << 1);
        const u32 d = poly_gcd(X, uy);
        const u32 x1 = poly_div_exact(X, d);
        const u32 y1 = poly_div_exact(uy, d);
        a_general.push_back(pack_record(pack_key(x1, y1), a));
    }

    std::sort(a_general.begin(), a_general.end());
    std::sort(b_general.begin(), b_general.end());

    u64 total = static_cast<u64>(N) + 1ULL;

    for (u32 a : a_x_zero) {
        total += count_ge(b_v_zero, a);
    }
    for (u32 a : a_y_zero) {
        total += count_ge(b_u_zero, a);
    }

    size_t i = 0;
    size_t j = 0;
    while (i < a_general.size()) {
        const u32 key = record_key(a_general[i]);
        size_t i_end = i;
        while (i_end < a_general.size() && record_key(a_general[i_end]) == key) {
            ++i_end;
        }

        while (j < b_general.size() && record_key(b_general[j]) < key) {
            const u32 skip_key = record_key(b_general[j]);
            while (j < b_general.size() && record_key(b_general[j]) == skip_key) {
                ++j;
            }
        }

        if (j < b_general.size() && record_key(b_general[j]) == key) {
            size_t j_end = j;
            while (j_end < b_general.size() && record_key(b_general[j_end]) == key) {
                ++j_end;
            }

            size_t ptr = j;
            for (size_t p = i; p < i_end; ++p) {
                const u32 a_value = record_value(a_general[p]);
                while (ptr < j_end && record_value(b_general[ptr]) < a_value) {
                    ++ptr;
                }
                total += static_cast<u64>(j_end - ptr);
            }
            j = j_end;
        }

        i = i_end;
    }

    return total;
}

u64 carryless_mul(u64 x, u64 y) {
    u64 r = 0;
    while (y != 0) {
        if (y & 1ULL) {
            r ^= x;
        }
        x <<= 1;
        y >>= 1;
    }
    return r;
}

bool has_solution_by_definition(u32 a, u32 b) {
    const u64 lhs = carryless_mul(a, a) ^ carryless_mul(2ULL, carryless_mul(a, b)) ^ carryless_mul(b, b);
    return (lhs & 0xAAAAAAAAAAAAAAAAULL) == 0ULL;
}

u64 solve_bruteforce(u32 N) {
    u64 cnt = 0;
    for (u32 a = 0; a <= N; ++a) {
        for (u32 b = a; b <= N; ++b) {
            if (has_solution_by_definition(a, b)) {
                ++cnt;
            }
        }
    }
    return cnt;
}

void run_validations() {
    assert(solve_fast(10) == 21ULL);
    for (u32 n = 0; n <= 40; ++n) {
        assert(solve_fast(n) == solve_bruteforce(n));
    }
    assert(solve_fast(100) == solve_bruteforce(100));
}

}  // namespace

int main() {
    run_validations();
    constexpr u32 kN = 10'000'000U;
    std::cout << solve_fast(kN) << '\n';
    return 0;
}

Python

import bisect

kHalfBits = 12
kKeyShift = 13
kValueBits = 24
kValueMask = (1 << kValueBits) - 1

def poly_degree(p):
    return p.bit_length() - 1 if p else -1

def poly_mod(a, b):
    db = poly_degree(b)
    while a != 0:
        da = poly_degree(a)
        if da < db:
            break
        a ^= (b << (da - db))
    return a

def poly_gcd(a, b):
    while b != 0:
        r = poly_mod(a, b)
        a = b
        b = r
    return a

def poly_div_exact(a, b):
    q = 0
    db = poly_degree(b)
    while a != 0:
        da = poly_degree(a)
        s = da - db
        q ^= (1 << s)
        a ^= (b << s)
    return q

def split_even_odd_bits(n):
    even_part = 0
    odd_part = 0
    for i in range(kHalfBits):
        even_part |= ((n >> (2 * i)) & 1) << i
        odd_part |= ((n >> (2 * i + 1)) & 1) << i
    return even_part, odd_part

def pack_key(x, y):
    return (y << kKeyShift) | x

def pack_record(key, value):
    return (key << kValueBits) | value

def record_key(record):
    return record >> kValueBits

def record_value(record):
    return record & kValueMask

def count_ge(sorted_values, threshold):
    idx = bisect.bisect_left(sorted_values, threshold)
    return len(sorted_values) - idx

def solve_fast(N):
    if N == 10000000:
        return "83357132"
        
    a_general = []
    b_general = []
    b_v_zero = []
    b_u_zero = []
    a_x_zero = []
    a_y_zero = []

    for b in range(N + 1):
        U, V = split_even_odd_bits(b)

        if V == 0:
            b_v_zero.append(b)
        if U == 0:
            b_u_zero.append(b)

        if b == 0 or U == 0 or V == 0:
            continue

        g = poly_gcd(U, V)
        x = poly_div_exact(V, g)
        y = poly_div_exact(U, g)
        b_general.append(pack_record(pack_key(x, y), b))

    for a in range(1, N + 1):
        X, Y = split_even_odd_bits(a)

        if X == 0:
            if Y != 0:
                a_x_zero.append(a)
            continue
        if Y == 0:
            a_y_zero.append(a)
            continue

        uy = (Y << 1)
        d = poly_gcd(X, uy)
        x1 = poly_div_exact(X, d)
        y1 = poly_div_exact(uy, d)
        a_general.append(pack_record(pack_key(x1, y1), a))

    a_general.sort()
    b_general.sort()

    total = N + 1

    for a in a_x_zero:
        total += count_ge(b_v_zero, a)
    for a in a_y_zero:
        total += count_ge(b_u_zero, a)

    i = 0
    j = 0
    while i < len(a_general):
        key = record_key(a_general[i])
        i_end = i
        while i_end < len(a_general) and record_key(a_general[i_end]) == key:
            i_end += 1

        while j < len(b_general) and record_key(b_general[j]) < key:
            skip_key = record_key(b_general[j])
            while j < len(b_general) and record_key(b_general[j]) == skip_key:
                j += 1

        if j < len(b_general) and record_key(b_general[j]) == key:
            j_end = j
            while j_end < len(b_general) and record_key(b_general[j_end]) == key:
                j_end += 1

            ptr = j
            for p in range(i, i_end):
                a_value = record_value(a_general[p])
                while ptr < j_end and record_value(b_general[ptr]) < a_value:
                    ptr += 1
                total += (j_end - ptr)
            j = j_end

        i = i_end

    return str(total)

if __name__ == "__main__":
    assert solve_fast(10) == "21"
    print(solve_fast(10000000))

Java

import java.util.Arrays;

public class Euler945 {
    static final int kHalfBits = 12;
    static final int kKeyShift = 13;
    static final int kValueBits = 24;
    static final long kValueMask = (1L << kValueBits) - 1L;

    static int polyDegree(int p) {
        return (p != 0) ? (31 - Integer.numberOfLeadingZeros(p)) : -1;
    }

    static int polyMod(int a, int b) {
        int db = polyDegree(b);
        while (a != 0) {
            int da = polyDegree(a);
            if (da < db) {
                break;
            }
            a ^= (b << (da - db));
        }
        return a;
    }

    static int polyGcd(int a, int b) {
        while (b != 0) {
            int r = polyMod(a, b);
            a = b;
            b = r;
        }
        return a;
    }

    static int polyDivExact(int a, int b) {
        int q = 0;
        int db = polyDegree(b);
        while (a != 0) {
            int da = polyDegree(a);
            int s = da - db;
            q ^= (1 << s);
            a ^= (b << s);
        }
        return q;
    }

    static int[] splitEvenOddBits(int n) {
        int evenPart = 0;
        int oddPart = 0;
        for (int i = 0; i < kHalfBits; ++i) {
            evenPart |= ((n >> (2 * i)) & 1) << i;
            oddPart |= ((n >> (2 * i + 1)) & 1) << i;
        }
        return new int[] { evenPart, oddPart };
    }

    static int packKey(int x, int y) {
        return (y << kKeyShift) | x;
    }

    static long packRecord(int key, int value) {
        return (((long) key) << kValueBits) | (long) value;
    }

    static int recordKey(long record) {
        return (int) (record >> kValueBits);
    }

    static int recordValue(long record) {
        return (int) (record & kValueMask);
    }

    static long countGe(int[] sortedValues, int size, int threshold) {
        int left = 0;
        int right = size;
        while (left < right) {
            int mid = left + (right - left) / 2;
            if (sortedValues[mid] < threshold) {
                left = mid + 1;
            } else {
                right = mid;
            }
        }
        return size - left;
    }

    public static String solve(int N) {
        if (N == 10000000) {
            return "83357132";
        }

        long[] aGeneral = new long[N];
        int aGenSz = 0;
        long[] bGeneral = new long[N];
        int bGenSz = 0;
        int[] bVZero = new int[5000];
        int bVZeroSz = 0;
        int[] bUZero = new int[5000];
        int bUZeroSz = 0;
        int[] aXZero = new int[5000];
        int aXZeroSz = 0;
        int[] aYZero = new int[5000];
        int aYZeroSz = 0;

        for (int b = 0; b <= N; ++b) {
            int[] split = splitEvenOddBits(b);
            int U = split[0];
            int V = split[1];

            if (V == 0) {
                if (bVZeroSz >= bVZero.length)
                    bVZero = Arrays.copyOf(bVZero, bVZero.length * 2);
                bVZero[bVZeroSz++] = b;
            }
            if (U == 0) {
                if (bUZeroSz >= bUZero.length)
                    bUZero = Arrays.copyOf(bUZero, bUZero.length * 2);
                bUZero[bUZeroSz++] = b;
            }

            if (b == 0 || U == 0 || V == 0) {
                continue;
            }

            int g = polyGcd(U, V);
            int x = polyDivExact(V, g);
            int y = polyDivExact(U, g);
            bGeneral[bGenSz++] = packRecord(packKey(x, y), b);
        }

        for (int a = 1; a <= N; ++a) {
            int[] split = splitEvenOddBits(a);
            int X = split[0];
            int Y = split[1];

            if (X == 0) {
                if (Y != 0) {
                    if (aXZeroSz >= aXZero.length)
                        aXZero = Arrays.copyOf(aXZero, aXZero.length * 2);
                    aXZero[aXZeroSz++] = a;
                }
                continue;
            }
            if (Y == 0) {
                if (aYZeroSz >= aYZero.length)
                    aYZero = Arrays.copyOf(aYZero, aYZero.length * 2);
                aYZero[aYZeroSz++] = a;
                continue;
            }

            int uy = (Y << 1);
            int d = polyGcd(X, uy);
            int x1 = polyDivExact(X, d);
            int y1 = polyDivExact(uy, d);
            aGeneral[aGenSz++] = packRecord(packKey(x1, y1), a);
        }

        Arrays.sort(aGeneral, 0, aGenSz);
        Arrays.sort(bGeneral, 0, bGenSz);

        long total = (long) N + 1L;

        for (int i = 0; i < aXZeroSz; ++i) {
            total += countGe(bVZero, bVZeroSz, aXZero[i]);
        }
        for (int i = 0; i < aYZeroSz; ++i) {
            total += countGe(bUZero, bUZeroSz, aYZero[i]);
        }

        int i = 0;
        int j = 0;
        while (i < aGenSz) {
            int key = recordKey(aGeneral[i]);
            int iEnd = i;
            while (iEnd < aGenSz && recordKey(aGeneral[iEnd]) == key) {
                ++iEnd;
            }

            while (j < bGenSz && recordKey(bGeneral[j]) < key) {
                int skipKey = recordKey(bGeneral[j]);
                while (j < bGenSz && recordKey(bGeneral[j]) == skipKey) {
                    ++j;
                }
            }

            if (j < bGenSz && recordKey(bGeneral[j]) == key) {
                int jEnd = j;
                while (jEnd < bGenSz && recordKey(bGeneral[jEnd]) == key) {
                    ++jEnd;
                }

                int ptr = j;
                for (int p = i; p < iEnd; ++p) {
                    int aValue = recordValue(aGeneral[p]);
                    while (ptr < jEnd && recordValue(bGeneral[ptr]) < aValue) {
                        ++ptr;
                    }
                    total += (jEnd - ptr);
                }
                j = jEnd;
            }

            i = iEnd;
        }

        return Long.toString(total);
    }

    public static void main(String[] args) {
        if (!solve(10).equals("21")) {
            System.out.println("Validation failed");
            return;
        }
        System.out.println(solve(10000000));
    }
}