// 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; }