lattice.cpp

/*

Problem:
Given circle with center at (0,0) and integer n [1..10^8].
Find radius r [1..n] with the most number of points with integral coordinates on the circumference of the circle
Examples: n = 2, output = 1; n = 5, output = 5; n = 10000, output = 5525

References:
https://oeis.org/A046109 Number of lattice points (x,y) on the circumference of a circle of radius n with center at (0,0)
https://oeis.org/A088959 (Lowest numbers which are d-Pythagorean decomposable, essentially a solution to the problem)
https://stackoverflow.com/questions/2267146/what-is-the-fastest-integer-factorization-algorithm
https://math.stackexchange.com/questions/2069171/fast-way-to-find-pythagorean-decompositions

Solution with a precalculated lookup table (python 3, same sequence as https://oeis.org/A088959):
import bisect
n=int(input())
k=[1, 5, 25, 65, 325, 1105, 5525, 27625, 32045, 160225, 801125, 1185665, 5928325, 29641625, 48612265]
print(k[bisect.bisect(k,n)-1])

*/


#include <bits/stdc++.h>
using namespace std;

#define uint unsigned __int32
uint *factors;
// const uint MAX_F = 134217728; // 2^27

const uint MAX_F = 50000000;

void buildFactors() {
    factors = new (nothrow) uint[(MAX_F + 1) * 2]; // 4 * 2 * 2^27 = 2^30 = 1GB

    if (factors == NULL)
        return; // not able to allocate enough free memory

    int i;
    for (i = 0; i < (MAX_F + 1) * 2; i++)
        factors[i] = 0;

    // Sieve of Eratosthenese
    factors[1 * 2] = 1;
    factors[1 * 2 + 1] = 1;
    for (i = 2; i * i <= MAX_F; i++) {
        for (; factors[i * 2] && i * i <= MAX_F; i++)
            ;
        factors[i * 2] = 1;
        factors[i * 2 + 1] = i;
        for (int j = 2; i * j <= MAX_F; j++) {
            factors[i * j * 2] = i;
            factors[i * j * 2 + 1] = j;
        }
    }
    for (; i <= MAX_F; i++) {
        for (; i <= MAX_F && factors[i * 2]; i++)
            ;
        if (i > MAX_F)
            return;
        factors[i * 2] = 1;
        factors[i * 2 + 1] = i;
    }
}

uint *factor(uint x, int &factorCount) {
    if (x > MAX_F) {
        factorCount = -1;
        return NULL;
    }
    static uint tmp[70];
    uint at = x;
    int i = 0;
    while (factors[at * 2] > 1) {
        tmp[i++] = factors[at * 2];
        //cout << "at:" << at << " tmp:" << tmp[i - 1] << endl;
        at = factors[at * 2 + 1];
    }
    if (i == 0) {
        //cout << "at:" << x << " tmp:1" << endl;
        tmp[i++] = 1;
        tmp[i++] = x;
    } else {
        //cout << "at:" << at << " tmp:1" << endl;
        tmp[i++] = at;
    }
    factorCount = i;
    return tmp;
}

int a(uint n) {
    int r = 1;
    int size;
    uint *f = factor(n, size);

    for (int i=0; i<size; i++) {
        if (f[i] == 1)
            continue;

        int count = 1;
        while (i < size - 1 && f[i] == f[i + 1]) {
            count++;
            i++;
        }

        if ((f[i] % 4)==1) {
            r *= 2*count + 1;
        }
    }

    return 4*r;
}

int main() {
    time_t start = time(0);

    cout << "Loading factors lookup table" << endl;
    buildFactors();
    if (factors == NULL) {
        cerr << "Need 1GB block of free memory" << endl;
        return 0;
    }

    time_t end = time(0);

    cout << "time: " << difftime(end, start) << " sec." << endl;

    uint c = 0;

    uint maxr = 0;
    uint p0 = 0;
    uint p = 0;

    start = time(0);

    int n = 100000000;

    for (uint i=1; i<=n; i++) {
        p = a(i);
        if (p>p0) {
            p0 = p;
            maxr = i;
            cout << i << ":" << p << endl;
        }
    }

    end = time(0);
    uint diff = difftime(end, start);

    cout << "time: " << diff << " sec. " << n/diff << " rad/sec." << endl;

    delete[] factors;
}