#! perl -w
# generate c++ code for sum of integer multiples of square roots of squarefree integers
#
# Copyright 2016 Ken Takusagawa
#
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or (at
# your option) any later version.
#
# This program is distributed in the hope that it will be useful, but
# WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
# General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program.  If not, see <http://www.gnu.org/licenses/>.

#this all-integer version is 35sec versus 100sec for double
#computed with pari/gp, the fractional parts of square roots of squarefree numbers *2^64.
@squarefree=qw[7640891576956012809 13503953896175478587 4354685564936845356 8291622248878821280 11912009170470909682 2993494466006504001 5840696475078001362 11170449401992604704 13681164004195116043 16103700368098801199 2270897969802886507 6620516959819538810 10746624748413865164 12735922825918526280 14680500436340154073 1826587625549305955 7105036623409894664 8803258048406862881 10473403895298186519 13734756587530879779 15328356941775268973 16898689309937973692 1526699215303891258 3032903034902830958 4519415375026354322 7436329637833083697 8868100629346186592 10282925794625328401 14431440979901333565 15784041429090275239 2608894026076872997 5167115440072839076 7677506975318801792 10142655110352140040 11359008891344774515 12564921287827708315 14946468407256161813 16122599570302859359 1148452749468771720 2288104705630789152 3419155146083754954 5656211926141927229 6762581272094930283 7861075802853663214 10035097865011968311 11110940050264098506 14295569720461580778 15343279831676731510 16384295103918246975 1021675468707975030 2037139971618630137 4049880416483348735 5047370360242696839 6039077607422743620 8005538875866918671 9950026440786680028 11873260860885435600 12827122698685698277 13775924154384562400 15658662624704411955];

#24 2m40
#31 55 hours (mkc)
#32 144 hours
# (speed difference easily explained by 31 was running concurrently with 32, so sharing processor).
$nterms=32;
$maxx=1;
$cutoff='18446744074l';
$doshell=1;
if($doshell){
    print STDERR "DO SHELL is ON\n";
}
print << 'EOF';
#include <iostream>
#include <cstdlib>
EOF
for($i=0;$i<$nterms;++$i){
#everything done with unsigned long because signed long overflow is undefined
    print "const unsigned long s$i=$squarefree[$i]ul;\n";
}
print "const int maxx=$maxx;\n";
print << 'EOF';
int main(int argc, char**argv){
  std::cout << std::scientific;
  if(argc<2){
      std::cerr << "need first coefficient for parallelization!" << std::endl;
      return -1;
  }
  int i0=atoi(argv[1]);
  std::cerr << "got " << i0 << std::endl;
  unsigned long c0=i0*s0;
EOF
    for($i=1;$i<$nterms;++$i){

        if($i==$nterms-1 and $doshell){
            # The zero in the highest order term has been done in the previous "shell" not including the highest order term.
            # Avoid negative numbers because symmetry.
            $start=1;
            print "unsigned long c$i=c",$i-1,"+s$i;\n";
        } else {
            $start='-maxx';
            print "unsigned long c$i=c",$i-1,"-maxx*s$i;\n";
        }
        print "for(int i$i=$start;i$i<=maxx;++i$i,c$i+=s$i){\n";
}
  $i=$nterms-1;
  print << "EOF";
signed long r=labs(static_cast<signed long>(c$i));
if(r<$cutoff){
std::cout << r << " ="
EOF
for($i=0;$i<$nterms;++$i){
    print qq( << " " << i$i);
}
print << 'EOF';
 << std::endl;
}
EOF
    for($i=1;$i<$nterms;++$i){
        print "}";
}
  print "}\n";