Flawless Squares

squares.c

// Flawless Squares
// $ cc -fopenmp -O2 -o squares squares.c
// Ref: https://old.reddit.com/r/C_Programming/comments/1de40e0
// This is free and unencumbered software released into the public domain.
#include <stdio.h>

int main(void)
{
    #pragma omp parallel for schedule(dynamic, 100000)
    for (long long root = 1; root < 3037000500; root++) {
        for (int splits = 1;; splits++) {
            long long log[32];
            int       reject = 0;
            int       len    = 0;
            long long sum    = 0;
            long long accum  = 0;
            long long scale  = 1;
            long long target = root*root;
            long long digits = target;

            for (int i = 0; digits; i++, digits/=10) {
                long long digit = digits%10;
                accum += digit * scale;
                scale *= 10ull;  // overflows in some reject cases
                if ((1<<i) & splits) {
                    reject |= !digit;  // no leading zeros
                    sum += accum;
                    log[len++] = accum;
                    accum = 0;
                    scale = 1;
                }
            }
            sum     += accum;
            log[len] = accum;
            len     += !!accum;

            if (target>1 && sum==target) {
                break;  // no splits occurred
            }

            if (!reject && sum==root) {
                #pragma omp critical
                {
                    printf("%19lld = %10lld^2 = (", target, sum);
                    for (int i = len-1; i >= 0; i--) {
                        if (i < len-1) {
                            fputs(" + ", stdout);
                        }
                        printf("%lld", log[i]);
                    }
                    fputs(")^2\n", stdout);
                    fflush(stdout);
                }
                break;
            }
        }
    }
}


#if 0
// Table-based computation for [1, 10^12] derived from above.
static int num_flawless(long long min, long long max)
{
    static const int roots[] = {
             1,      9,     36,     45,     55,     82,     91,    235,    297,
           369,    379,    414,    657,    675,    703,    756,    792,    909,
           918,    945,    964,   1296,   1702,   1782,   2223,   2728,   3366,
          3646,   3682,   4132,   4906,   4950,   5050,   6832,   7191,   7272,
          7389,   7533,   7543,   7587,   7777,   7956,   8416,   8443,   8767,
          8856,   8938,   9208,   9315,   9325,   9586,   9621,   9640,   9765,
          9909,   9918,   9945,   9955,   9964,  10512,  12222,  12727,  17271,
         17344,  22222,  22231,  25354,  27414,  33669,  34327,  35442,  36558,
         36855,  41149,  41814,  42643,  45657,  46279,  48835,  49563,  49906,
         51454,  53415,  57393,  57547,  59284,  62146,  64341,  65242,  67132,
         68428,  69058,  69157,  69669,  72612,  74647,  75331,  76851,  77121,
         77778,  79516,  81118,  82045,  82656,  85941,  87571,  88768,  88812,
         88813,  89632,  91756,  93402,  94141,  94581,  97137,  99055,  99226,
         99244,  99343,  99559,  99631,  99703,  99765, 105112, 105382, 117343,
        122221, 131959, 142857, 148149, 163567, 167832, 172827, 181819, 187110,
        203274, 208495, 215829, 249832, 283924, 304930, 310339, 318682, 319176,
        329967, 333666, 336699, 338760, 339453, 344143, 349785, 351352, 352693,
        354393, 356329, 356643, 357796, 364420, 364932, 369261, 369792, 370017,
        370134, 383581, 389665, 390313, 419860, 428419, 434331, 440874, 441217,
        452133, 452241, 461539, 463284, 466830, 471537, 473967, 482608, 499321,
        499500, 500500, 512811, 514926, 524638, 532917, 533170, 534150, 536284,
        538461, 549459, 551025, 561781, 571464, 584308, 584452, 596619, 600633,
        609687, 627372, 639739, 642078, 643357, 643672, 647029, 648648, 649782,
        668332, 670033, 673651, 675037, 681318, 684099, 689778, 695071, 697834,
        719280, 736687, 739863, 758475, 758863, 784297, 790372, 791505, 791757,
        792684, 795484, 810819, 811305, 812890, 818181, 825886, 826074, 831672,
        831718, 834166, 836172, 847629, 850546, 851851, 852165, 852517, 853542,
        853867, 857143, 857151, 865242, 872146, 877501, 880749, 887824, 888363,
        889642, 890046, 906427, 906796, 907740, 908235, 910719, 920494, 921528,
        921969, 935785, 945649, 948411, 952335, 961039, 968256, 969669, 977769,
        991144, 991233, 991584, 992467, 992611, 993168, 994482, 996634, 998218,
        998298, 998704, 999036, 999055, 999082, 999091, 999208, 999244, 999297,
        999325, 999343, 999586, 999621, 999631
    };
    int nroots = sizeof(roots) / sizeof(*roots);
    int mini = 0;
    int maxi = nroots-1;

    // NOTE: these could be binary searches
    for (; mini<nroots-1; mini++) {
        long long r = roots[mini];
        if (r*r >= min) break;
    }
    for (; maxi>0; maxi--) {
        long long r = roots[maxi];
        if (r*r <= max) break;
    }

    return maxi - mini + 1;
}
#endif