Refined greedy computation of multiplicative sequences

Back to the main experimental results page.

Here is the source of a C++ program for investigating refined greedy algorithms in the case of completely-multiplicative sequences.

The freedom lies in choosing the values at prime indices. Many methods can be devised but only a few are implemented for the moment, feel free to add many more (please indicate what you've added or modified). The code is commented and hopefully fairly straightforward (it probably could be optimized here and there, but has been carefully checked and seems correct).

The input file must be in the format of a sequence separated by spaces like + - + (with a space after the last sign, no carriage returns). These indicate the values at the first primes: [math]\displaystyle{ x_2 }[/math], [math]\displaystyle{ x_3 }[/math], [math]\displaystyle{ x_5 }[/math]... In fact it is also possible not to specify any one value before the end by putting a 0.

There are three output files: the first one will contain the sequence found, the second may contain information on the choices made for each undetermined prime (I've commented it out for now), and the third contains in columns: n, D(n), partial sum of sequence up to n, partial sum of 2-HAP up to n, ... 11-HAP up to n (you can add more if so you wish).

Finally the constant M is the length of the sequence to be computed.

/*----------------------------------------------------------
author: Thomas Sauvaget
licence: public domain.
files: just this one.
------------------------------------------------------------
modification record: 
13-jan-2010: written from scratch. Really basic code...
14-jan-2010: implemented various methods to choose values 
             at prime indices.
----------------------------------------------------------*/


//---------------------------------------------------------
//-- libraries & preprocessor
//---------------------------------------------------------
#include <iostream> //for imput-output with terminal
#include <cmath>    //small math library
#include <algorithm> //useful for permutations and the likes
#include <vector>
#include <fstream> //for imput-output with files
#include <iomanip> //for precision control when outputing
#include <stdlib.h>
#include <string>
#include <sstream>
//#include <assert> //for end of file command


//---------------------------------------------------------
//-- namespace (standard)
//---------------------------------------------------------
using namespace std;


//---------------------------------------------------------
//-- global variables
//--------------------------------------------------------- 


int const M=3000;


int tab[M], primes[M], prtab[M], glopartsum[M];
char filenameIN[50], filenameOUT[50], filenameOUTT[50], filenameOUTTT[50];
ofstream myfileOUT;//sequence
ofstream myfileOUTT;//other data
ofstream myfileOUTTT;//yet other data
ifstream myfileIN;


//---------------------------------------------------------
//-- useful handmade functions 
//---------------------------------------------------------

//---------
//getintval 
//converts string containing a single integer to int
//---------
int getintval(string strconv){
  int intret;

  intret=atoi(strconv.c_str());

  return(intret);
}

//---------
//mypow
//integer power: a^b
//---------
int mypow(int a, int b){
  int i, ans=1;
  for(i=1;i<=b;i++){
    ans*=a;
  }
  if(b==0){
    return 1;
  }
  else{
    return ans;
  }
}


//---------
//primfact
//compute prime factors and store
//returns 2 iff n is prime
//---------
int primfact(int n){

  int a, b, c, d, ncurr;
  for(a=0;a<M;a++){
    prtab[a]=0;
  }
  ncurr=n;
  for(b=0;b<n;b++){//surely the n-th prime is greater than n already
    while(primes[b]!=0 && ncurr>1 && (ncurr%primes[b])==0){
      prtab[b]+=1;
      ncurr/=primes[b];
    }

    if(primes[b]==0){
      //exit
      b+=M;
    }
  }
  d=0;
  for(c=0;c<M;c++){
    if(prtab[c]!=0){
      d+=prtab[c];
    }
  }
  if(d==1){
    //cout<<"   "<<n<<"  is prime"<<endl;
    return 2;
  }
  if(d!=1){
    //cout<<"   "<<n<<"  is composite"<<endl;
    return 0;
  }
}



//---------
//disc
//computes discrepancy of sequence of +/-1 contained in tabb 
//over all HAPs up to len
//store signed partial sums up to len: 
//  glopartsum[0]=0, glopartsum[1]=sum of seq, glopartsum[2]=sum of 2-HAP...
//---------
int disc(int tabb[M], int len){

  int n, d, i;
  int signsum;
  int absum;
  int r;
  int beflocdisc;
  int aflocdisc;

  beflocdisc=0;
  aflocdisc=0;

  for(n=2;n<=len;n++){
    for(d=1;d<n;d++){
      signsum=0;
      absum=0;
      glopartsum[d]=0;
      r=int(floor(n/d));        
      for(i=1;i<=r;i++){
	signsum+=tabb[d*i-1];
      }
      glopartsum[d]=signsum;
      absum=abs(signsum);
      aflocdisc=max(beflocdisc,absum);
      beflocdisc=aflocdisc;
    }
  }
  return aflocdisc;
}


//---------
//lplain
//just adds the entries of tabb
//---------
int lplain(int tabb[M], int len){
  int norm, a;
  norm=0;
  for(a=0;a<len;a++){
    norm+=tabb[a];
  }
  return norm;
}


//---------
//l1norm
//---------
int l1norm(int tabb[M], int len){
  int a, norm;
  norm=0;
  for(a=0;a<len;a++){
    norm+=abs(tabb[a]);
  }
  return norm;
}


//---------
//l2norm
//---------
double l2norm(int tabb[M], int len){
  int a, snorm;
  double norm;
  snorm=0; norm=0.0;
  for(a=0;a<len;a++){
    snorm+=(tabb[a]*tabb[a]);
  }
  norm=sqrt(snorm);
  return norm;
}


//---------
//lharmw
//harmonic weighting
//---------
double lharmw(int tabb[M], int len){
  int a;
  double norm;
  norm=0;
  for(a=0;a<len;a++){
    norm+=tabb[M]/(a+1);
  }
  return norm;
}


//----------------------------------------------------------
//-- main:
//
//we compute a multiplicative sequence and its discrepancy by 
//specifying its values at first p primes
//----------------------------------------------------------
int main(int argc, char *argv[]){

  int i, k, d, u, v, w, j, c, t, imax=0, maxlength, currprime, iprec;
  int datab[M];
  int mycontinue, g, loclen, h;
  double nmplus, nmminus;
  string line, buff;
  string myplus ("+");
  string myminus ("-");
  string myzero ("0");

  //-- textlike user interface
  cout<<endl;
  cout<<" give input SEQUENCE filename:"<<endl;
  cin>>filenameIN;
  cout<<endl;


  cout<<endl;
  cout<<" give output LONG MULT SEQUENCE filename:"<<endl;
  cin>>filenameOUT;
  cout<<endl;


  cout<<endl;
  cout<<" give output CHOICE filename:"<<endl;
  cin>>filenameOUTT;
  cout<<endl;


  cout<<endl;
  cout<<" give output DISC & PARTSUM filename: "<<endl;
  cin>>filenameOUTTT;
  cout<<endl;

  //-- file creation 
  myfileOUT.open(filenameOUT,ios::out);
  myfileOUT.close();

  myfileOUTT.open(filenameOUTT,ios::out);
  myfileOUTT.close();  
  
  myfileOUTTT.open(filenameOUTTT,ios::out);
  myfileOUTTT.close();


  //fill tables with zeros
  for(i=0;i<M;i++){
    tab[i];
    datab[i]=0;
    primes[i]=0;
    prtab[i]=0;
  }

  
  //construct and store primes up to M
  //primes[a] = (a+1)-th prime
  iprec=2;
  for(i=2;i<=M;i++){
    k=0;
    for(j=1;j<=i;j++){
      if( (i%j)==0){
	k++;
      }
    }
    if(k==2){
      primes[iprec-2]=i;
      iprec+=1;
    }
  }


  //-- read the choices of values at first few primes
  myfileIN.open(filenameIN,ios::in);

  if(myfileIN.is_open()){
    
    i=0;
    getline(myfileIN,line);
    stringstream stsm(line);
    while(stsm>>buff){
      currprime=primes[i];
      if(buff.compare(myplus)==0){
	datab[currprime-1]=1;
      }
      if(buff.compare(myminus)==0){
	datab[currprime-1]=-1;
      }
      if(buff.compare(myzero)==0){
	datab[currprime-1]=0;
      }
      if( datab[currprime-1] !=1 && datab[currprime-1] !=-1 && datab[currprime-1]!=0){
      cout<<"problem: element number "<<i+1<<" is neither + nor - nor 0."<<endl;
	abort();
      }
      i+=1;
    }
    imax=i;
    cout<<" this initial sequence contains "<<imax<<" elements."<<endl;
    cout<<endl;

  }
  else{
    cout<<"there is no file with this very name"<<endl;
    cout<<" (don't forget extension like .txt or .dat in particular)"<<endl;
    abort();
  }
  myfileIN.close();



  //-- main processing:  
  //fill the rest of the sequence multiplicatively
  //and according to chosen method

  //convention: datab[i]=x_{i+1}
  //and first index is always associated to +
  datab[0]=1;

  mycontinue=0;

  while(mycontinue==0){

 
  
    for(u=2;u<M;u++){
            
      if(datab[u-1]==0){	
	
	k=0;
	k=primfact(u);
	
	
	//if prime: there's freedom, choose one method to fill it
	if(k==2){

	  //choice 1:  silly uniform -1
	  //datab[u-1]=-1;


	  //choice 2: take into account of partial sums
	  loclen=u;

	  cout<<"  adressing "<<loclen<<endl;
	  
	  //values are known below loclen, not at prime loclen, 
	  //and from loclen+1 to loclen+h. we then choose loclen 
	  //(lots of flexibility)


	  //compute h 
	  h=1;
	  for(v=loclen+1;v<M;v++){
	    if(datab[v]==0){
	      h=v-loclen;
	      v+=M;
	    }
	  }
	  cout<<" values after are known up to "<<loclen+h<<endl;
	  

	  
	  //now compute the effect of imposing datab[u-1]=+1 
	  //would have on partial sums up to loclen+h-1
          //you can change lplain to l2norm or whatever other idea
	  datab[u-1]=+1;
	  c=disc(datab,loclen+h-1);	    
	  nmplus=lplain(datab,loclen+h-1);

	  //compute the same with -1 instead
	  datab[u-1]=-1;
	  c=disc(datab,loclen+h-1);
	  nmminus=lplain(datab,loclen+h-1);
	  

	  //now choose
	  if(abs(nmplus)>abs(nmminus)){
	    datab[u-1]=-1;
	  }
	  else{
	    datab[u-1]=+1;
	  }
	  cout<<"  nmplus="<<nmplus<<", nmminus="<<nmminus<<endl;

	  /*------temporarily removed, could be useful to monitor	    
	  //record those numbers for further analysis
	  myfileOUTT.open(filenameOUTT,ios::app);
	  myfileOUTT<<u<<"\t"<<nmplus<<"\t"<<nmminus<<endl;
	  myfileOUTT.close();
	  ------------------------*/

	  //choice 3:
	  //more ideas to be added...



	  //in any case:  give info
	  cout<<"    imposing: x_"<<u<<"="<<datab[u-1]<<endl;


          //------------ in theory there's nothing to change beyond this point, except
          //             possibily the number of HAP's partial sums printed in the third file


	}
	else{//if composite then try to compute it multiplicatively
	  
	   
	  //loop on prime factors of u to see what is known so far
	  
	  w=1;
	  for(j=0;j<M;j++){
	    
	    if(prtab[j]!=0 && datab[ primes[j]-1 ]!=0){
		//then this x_j is known
	      w*=1;
	    }
	    
	    if(prtab[j]!=0 && datab[ primes[j]-1 ]==0){
	      //then this x_j not known, enough to discard
	      w*=0;
	    }
	  }
	  
	  
	  //if all x_j are known then compute x_u
	  if(w==1){
	    datab[u-1]=1;
	    for(j=0;j<M;j++){
	      if(prtab[j]!=0){
		datab[u-1]*=mypow(datab[ primes[j]-1 ],prtab[j]);
	      }
	    }
	    cout<<"      just computed x_"<<u<<"="<<datab[u-1]<<endl;    
	  }
	  else{
	    //cout<<" cannot yet compute x_"<<u<<endl;
	    //cout<<endl;
	  }
	  
	}//end of mult comp attempt
      }//end of this u
      else{
	//cout<<" already known"<<endl;
      }
    }//end of u loop
    
    //now see whether we should go on
    mycontinue=1;
    for(c=0;c<M-1;c++){
      if(datab[c]!=0){
	mycontinue*=1;
      }
      else{
	mycontinue*=0;
      }
    }
    
  }//end of mycontinue


  cout<<" ---now out, computing partial sums"<<endl;

  myfileOUT.open(filenameOUT,ios::app);
  myfileOUTTT.open(filenameOUTTT,ios::app);
  for(u=0;u<M;u++){

    if(u%50==0){
      cout<<"  computed "<<u<<" data points"<<endl;
    }

    //first output discrepancy and partial sums up to length u, and 11-th HAP
    k=disc(datab,u);
    if(u>0){
      myfileOUTTT<<u<<"\t"<<k<<"\t";
      for(d=1;d<u;d++){
	if(d<=11){
	  myfileOUTTT<<glopartsum[d]<<"\t";
	}
      }
      myfileOUTTT<<"\n";
    }

    //now output sequence value
    if(datab[u]==1){
      myfileOUT<<"+ ";
    }
    if(datab[u]==-1){
      myfileOUT<<"- ";
    }
  }
  myfileOUT.close();
  myfileOUTTT.close();

  cout<<" discrepancy="<<disc(datab,M)<<endl;

  //exit normally
  return 0;	        
}