//  Author: Jose Aranda, 2024.
//  Public Domain source code.
//
//  OEIS Sequence A006991
//  Make:
//  g++ -march=native -Wextra -Ofast -o p8.exe p8.cpp
//      -static

// Code modules:
// 01 INIT
// 02 Fill Vectors
// 03 Evaluate n

#pragma pack(16)

#include <cstdlib>
#include <cstdint>
#include <cstring>
#include <iostream>
#include <iomanip>
#include <array>
#include <vector>
#include <cmath>
using namespace std;

// I32 integers only 9 full digits.
typedef std::int_fast64_t I64;
typedef std::int_fast32_t I32;

// n odd:  a A1, b A2
// n even: a A3, b A4
static vector<I32> A1,A2,A3,A4;
static vector<I32> VX2;

// function is SquareFree(p)
bool isSF(const I32 p) __attribute__((const));
bool isSF(const I32 p) {
if(p>49) {
    if(p%4  ==0)return false;
    if(p%9  ==0)return false;
    if(p%25 ==0)return false;
    if(p%49 ==0)return false;
    if(p%121==0)return false;
    }

const double a= (double)p;
const double b= sqrt(a);
const I32 c= floor(b)+1;
for(I32 s= c;s>=2;s--)   {
    if(p%s>0)continue;
    const I32 ps = p/s;
    if((ps%s)==0) return false;
	}
return true;  // squarefree
}

// print one term
void toprinter(I32 qqqq,I64 n) {
cout<<setw(10)<<qqqq<<setw(12)<<n<<endl;
}

int main(int argc,char* argv[]) {

// **************
// 01 INIT module
// **************

// get line parameters
if(argc!=2) {
    cout<<"OEIS Sequence A006991"<<endl;
    cout<<"  Print All terms <= 10^Last"<<endl;
    cout<<"  Syntax: "<<argv[0]<<" {Last}"<<endl;
    return 1;
    }
I32 lasterm= atoi(argv[1]);
if(lasterm<1 || lasterm>8) {
    cout<<"  Last must be between 1,8"<<endl;
    return 1;
    }

const I32 lastprinted= (I32) pow(10,lasterm);
const I32 ASIZE= (I32) pow(10,lasterm)+2;
const I32 CLIM = (I32) ceil(sqrt(ASIZE))+1;

// A1 = (2,1, 8)    A072068 odd b
// A2 = (2,1,32)    A072069 odd a
// A3 = (4,1, 8)    A072070 even a
// A4 = (4,1,32)    A072071 even b
cout<<"Init Vector's Space"<<endl;
A1.reserve(ASIZE);A2.reserve(ASIZE);
A3.reserve(ASIZE);A4.reserve(ASIZE);
for(I32 i=0;i<ASIZE;i++) {
    A1.push_back(0);A2.push_back(0);
    A3.push_back(0);A4.push_back(0);
    };

// fill Squares
VX2.reserve(CLIM);
for(I32 i=0;i<CLIM;i++) {
    VX2.push_back(i*i);
    }

// ***************
// 02 Fill vectors
// ***************

cout<<"OEIS Sequence A006991"<<endl;
cout<<"  Computing A1, A2, A3, A4."<<endl;
cout<<"  No progress indicators will be displayed";
cout<<endl;

// explore (x,y,z) from [=0 to <CLIM].
// 4 escalar product's  (coef)*(x,y,z)
// pval variable:
// negative values are counted by variable "pval".
// pex:
// (2,3,4) 8 points, (0,3,4) 4 points
// (0,0,4) 2 points, (0,0,0) 1 point

// 3 nested for's. Added some speedup code.
// order: z,x,y. (ordered by coef)

for(I32 z= 0;z< CLIM;z++)  {
// compute the (x,y,z) space.
// variables with name t___ are temporary.
    I64 z2= (I64)VX2[z];
    I64 tz8=  (z2<<3);
    I64 tz16= (z2<<4);
    if(tz8>=ASIZE)break;
    for(I32 x= 0;x< CLIM;x++)  {
        I64 x2= (I64)VX2[x];
        I64 tx2= (x2<<1);
        I64 txz28= tx2 + tz8;
        if(txz28 >= ASIZE)break;
        for(I32 y= 0;y< CLIM;y++)  {
            I64 y2 = (I64)VX2[y];
            // 
            I64 tyz18= y2+tz8;
            if(tyz18 >= ASIZE)break;

            I64 txyz218= tyz18 + tx2;
            if(txyz218 >= ASIZE)break;
 
            I32 pval=1;    // weight value
            if(x>0)pval= pval<<1;
            if(y>0)pval= pval<<1;
            if(z>0)pval= pval<<1;
// P=(x,y,z)
// A1 n1 = (2,1, 8)*P   A072068 odd n, a== 2*b
// A2 n2 = (2,1,32)*P   A072069 odd n, b
// A3 n3 = (4,1, 8)*P   A072070 even n, a==2*b
// A4 n4 = (4,1,32)*P   A072071 even n, b

            // n1 (2,1,8)
            I64 n1= txyz218;
            if(n1<ASIZE) A1[n1]+= pval;
            
            I64 tz24= tz16 + tz8;
            
            // n2 (2,1,32)
            I64 n2= txyz218+tz24;
            if(n2<ASIZE) A2[n2]+= pval;

            // n3 (4,1,8)
            I64 n3= txyz218+tx2;
            if(n3<ASIZE) A3[n3]+= pval;
            
            // n4 (4,1,32)
            I64 n4= txyz218+tx2+tz24;
            if(n4<ASIZE) A4[n4]+= pval;
            }   // 3 nested for's closed
        }
    }

// ***********
// 03 Evaluate 
// ***********

// init interval vars.
I32 qqqq=0;     // global counter
I32 qm123=0;
I32 qm567=0;

/////////////
// Evaluate n
/////////////

for(I32 br= 1;br< ASIZE;br++)   // any integer
    {
    if(br > lastprinted)break;

// evaluate this number
    const bool bisSF= isSF(br);if(!bisSF)continue; 

// br==n and is squarefree
    I32 n= br;

// case n mod 8
    const int  nm8= n%8;

// Never
    const bool ism04= (nm8==0 || nm8==4);
    if(ism04)continue;

// Always
    const bool ism567= (nm8>=5 && nm8<=7);
    if(ism567) {
        qqqq++;
        toprinter(qqqq,br);
        qm567++;continue;
        }

// Unknown. mod8 1,2,3.
    const bool even= (n%2==0);
    I32 a,b;
// even a A3 b A4
// odd  a A1 b A2
    if(even) {
        // look i=(n/2).
        auto i= (n>>1); a =A3[i];b =A4[i];
        }
    else {
        auto i= n;      a =A1[i];b =A2[i];
        }

    // n by the Tunnell Theorem criteria.
    // and YES (0,0) match.
    if(a==(b<<1)) {
        qqqq++;
        toprinter(qqqq,br);
        qm123++;
        }
    }     // close br for

return 0;
}