Queremos resolver ~\frac 1 a + \frac 1 c = \frac 2 b~. Consideremos ~x = \frac b 2~. Entonces queremos resolver $$\displaystyle \frac 1 a + \frac 1 c = \frac 1 x$$ Es bien conocido cómo factorizar esta ecuación, pero igual lo ponemos aquí:

Multiplicando por ~acx~ a ambos lados de la ecuación, obtenemos ~cx + ax = ac~, que es lo mismo que ~cx + ax - ac = 0~. Restamos ~x^2~ a ambos lados de la ecuación, obteniendo: $$\displaystyle cx + ax - ac - x^2 = - x^2$$

Ahora sacamos factor común ~a~ con los términos ~ax~ y ~-ac~, y factor común ~-x~ con los términos ~cx~ y ~-x^2~, obteniendo ~a(x - c) - x(x - c) = -x^2~, luego sacamos factor común una vez más -esta vez sería ~(x-c)~-, obteniendo ~(x-c)(a-x)=-x^2~. Multiplicando por ~-1~ a ambos lados obtenemos $$\displaystyle (x-c)(x-a)=x^2$$

Queríamos llegar esa factorización. Ahora sustituimos ~x~, quedando ~\left(\frac b 2 - c\right)\left(\frac b 2 -a\right) = \frac{b^2} 4~, que es lo mismo que $$\displaystyle (b - 2c)(b - 2a) = b^2$$

Notemos que ~b^2~ es entero, y ~(b - 2c)~ y ~(b - 2a)~ también lo son. Por tanto, ~(b - 2c)~ y ~(b - 2a)~ son divisores de ~b^2~. Entonces, solo debemos considerar las parejas de números ~d_1, d_2~, tales que ~d_1 \cdot d_2 = b^2~, (es decir, los divisores de ~b^2~), y para cada pareja, ~a~ y ~c~ toman valores únicos. Si estamos analizando la pareja ~d_1, d_2~, entonces ~d_1 = b - 2c~, y ~d_2 = b - 2a~, por lo que ~c = \frac{b - d_1}{2}~, y ~a = \frac{b - d_2}{2}~. Solo contamos los que sean enteros positivos.

Quedan varias cosas por analizar:

• La primera es acotar la cantidad de divisores de ~b^2~. Es bien conocido que la cantidad de divisores de ~n~, para ~n \le 10^{18}~, es de orden ~\mathcal{O}(\sqrt[3]{n})~. Como queremos los divisores de ~b^2~, y ~b \le 10^9~, tenemos cerca de ~10^6~ divisores en el peor caso, lo cual es suficientemente pequeño.

• Bien, ya sabemos cuántos divisores son como mucho. Ahora, cómo iterar por los divisores de ~b^2~, porque incluso cuando son pocos, iterar por todos los números chequeando si son o no divisores resulta costoso.

Veamos este subproblema de manera más general: Dado ~N~, computar los divisores de ~N~. Para hacerlo, consideremos la factorización en primos de ~N~; digamos que ~N = p_1^{a_1} \cdot p_2^{a_2}\cdot p_3^{a_3}\dots p_k^{a_k}~, donde cada ~p_i~ es un primo distinto y cada ~a_i~ es un entero positivo. Cada divisor de ~N~ se representa con exactamente los mismos primos de ~N~, y con exponentes no negativos menores o iguales a los exponentes de ~N~. Por tanto, podemos probar todas las combinaciones ~b_1, b_2, b_3, \dots, b_k~, (con ~0 \le b_i \le a_i~), y para cada combinación generar el número ~p_1^{b_1} \cdot p_2^{b_2}\cdot p_3^{b_3}\dots p_k^{b_k}~, el cual es un nuevo divisor de ~N~. Para probar todas las combinaciones podemos hacer backtracking.

• Ahí surge un problema, cómo factorizar ~b^2~. Es sencillo factorizar ~N~ en tiempo ~\mathcal{O}(\sqrt{N})~, pero en este caso ~N = b^2 \le 10^{18}~. Por suerte, la factorización de ~b^2~ es simplemente la misma factorización de ~b~, solo con los exponentes duplicados. Entonces, podemos factorizar ~b~, lo cual tomaría ~O(\sqrt{b})~ operaciones, que si nos podemos permitir.

Resumen:

• El problema se reduce a contar las soluciones enteras ~(a, c)~ de la ecuación ~(b - 2a)\cdot(b - 2c) = b^2~.
• Hay una correspondencia 1-1 entre las parejas de divisores ~(d_1, d_2)~ de ~b^2~, y las parejas de soluciones ~(a, c)~. Por tanto el problema se reduce a buscar los divisores de ~b^2~.
• Para hacerlo, factorizamos ~b~, duplicamos todos los exponentes, obtenemos ese arreglo de exponentes ~a_1, a_2, \dots, a_k~, y hacemos un backtracking probando todos los posibles arreglos de exponentes ~b_1, b_2, \dots, b_k~, con ~0 \le b_i \le a_i~ para cada ~i \in [1, k]~. Para cada uno de esos exponentes, sacamos el divisor resultante, y vemos si conduce a una pareja ~(a, c)~ válida.

Complejiad temporal: ~\mathcal{O}(\sqrt{b})~.

Código de ejemplo:

#include <bits/stdc++.h>

using namespace std;

typedef long long ll;

typedef pair<ll,ll> par;
typedef pair<ll,int> ii;
vector<par> ans;
vector<ii> v;

void factorize(int x){
    for(ll i = 2; i*i <= x; i++){
        int cnt = 0;
        while(x % i == 0){
            x /= i;
            cnt++;
        }
        if(cnt > 0) v.push_back({i,2*cnt});
    }
    if(x > 1) v.push_back({x,2});
}

ll b, n;

void bct(int id, ll d){
    if(d > b) return;
    if(id == v.size()){
        ll p = n/d;
        if(d % 2 != b % 2) return;
        if(p % 2 != b % 2) return;
        ll a = (d + b)/2;
        ll c = (p + b)/2;
        if(a > c) swap(a,c);
        ans.push_back({a,c});
        if(a != c) ans.push_back({a,c});
        return;
    }
    bct(id+1,d);
    ll cur = 1;
    for(int i = 0; i < v[id].second; i++){
        cur *= v[id].first;
        bct(id+1,d*cur);
    }
}

int main(){
    ios_base::sync_with_stdio(0); cin.tie(0);

    cin >> b;
    n = b*b;
    factorize(b);

    bct(0,1);

    sort(ans.begin(), ans.end());

    cout << ans.size() << '\n';
    for(par it:ans)
        cout << it.first << ' '<< it.second <<'\n';

    return 0;
}