A (one dimensional) cellular automaton is a function1 F : Σ → Σ with the property that there is a K > 0 such that F (x)i depends only on the 2K + 1 coordinates xi−K , xi−K+1, . . . , xi−1, xi, xi+1, . . . , xi+K . A periodic point of σ is any x such that σ^p (x) = x for some p ∈ N, and a periodic point of F is any x such that F^q (x) = x for some q ∈ N. Given a cellular automaton F, a point x ∈ Σ is jointly periodic if there are p, q ∈ N such that σ^p (x) = F^q (x) = x, that is, it is a periodic point under both functions.
This project aims to explore the nature of one-dimensional Cellular Automata, in the hope of finding the structure of cellular automata through its periodic points.
License: MIT
ubuntu2004
/*1* Copyright (C) 2004 Bryant Lee2*3* This file is part of FPeriod.4*5* FPeriod is free software; you can redistribute it and/or modify6* it under the terms of the GNU General Public License as published by7* the Free Software Foundation; either version 2 of the License, or8* (at your option) any later version.9*10* FPeriod is distributed in the hope that it will be useful,11* but WITHOUT ANY WARRANTY; without even the implied warranty of12* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the13* GNU General Public License for more details.14*15* You should have received a copy of the GNU General Public License16* along with FPeriod; if not, write to the Free Software17* Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA18*/1920/* Main function21* FPeriod program finds points that are periodic under the block code22* induced by F23*24* Written by: Bryant Lee25* Date: 10/29/0426*/2728#define USEARRAY 029#define DEBUG 03031//Standard C++32#include <string.h>33#include <iostream>34#include <iomanip>35#include <fstream> //for file ops3637//C38#include <stdlib.h> //for atoi39#include <math.h> //for pow4041//STL42#include <vector>43#include <list>44#include <map>45#include <hash_map>4647//local48#include "StringOps.h"49#include "Comp.h"50#include "StorageKey.h"51#include "StorageVal.h"52#include "RecNode.h"5354//function definitions55void convertPeriodList(const string &strInput, vector<int> &periodList);56void exhaustiveExp(unsigned int shiftSize, const vector<int> &periodList,57string &prefix, Comp *comp, bool btrunc,58list<RecNode*> &recList,59list<RecNode*> &leastRecList,60vector<string> &commandbuffer);61void exhaustiveExpTwoShift(const vector<int> &periodList,62string &prefix, Comp *comp, bool btrunc,63list<RecNode*> &recList,64list<RecNode*> &leastRecList,65vector<string> &commandbuffer);66void clearStorageNodes(vector<StorageVal*> &table);67void clearRecList(list<RecNode*> &recList);68void output(ofstream & fout, ofstream &fout_preperiod, Comp * comp,69list<RecNode*> &recList, list<RecNode*> &leastRecList);70void truncoutput(ofstream & fout, Comp * comp, list<RecNode*> &leastRecList);71void findDivisors(unsigned int x, vector<int> &divisorsList);72bool leastPeriod(byte *word, unsigned int currPeriod,73const vector<int> &divisorsList);74bool leastPeriod(unsigned int word, unsigned int currPeriod,75const vector<int> &divisorsList);7677void outcontrol(string &prefix,78Comp * comp,79bool btrunc,80list<RecNode*> &recList,81list<RecNode*> &leastRecList,82vector<string> &commandbuffer);8384/*85* GLOBAL variables86*/87int UINTSIZE = 32;8889//boolean that determines whether output is given after each k90//or only at end91bool OUTPUT_AFTER_EACH = false;9293//main94int main() {95unsigned int shiftSize; //# symbols in shift96string prefix; //prefix for output files97string trunc; //"trunc" or "no trunc": determine whether to truncate output98bool btrunc; //true if truncate output is on99string strInput, strInput2, strInput3; //temp variables100vector<int> periodList; //list of periods101unsigned int i = 0, temp; //counter102list<RecNode*> recList; //list of records103list<RecNode*> leastRecList; //list of records for points of least period k104Comp *comp = NULL; //composition to apply105string outoption;106107//hold the command strings108vector<string> commandbuffer;109int commandbufferit;110111//grab all the input from cin and put in an array of strings112while(!cin.eof()) {113std::getline(cin, strInput);114commandbuffer.push_back(strInput);115}116commandbufferit = 0;117118//get "output after each" or "output at end"119outoption = commandbuffer[commandbufferit++];120trim(outoption);121transform(outoption.begin(), outoption.end(), outoption.begin(), tolower);122123if(outoption == "output after each")124OUTPUT_AFTER_EACH = true;125else if(outoption == "output at end")126OUTPUT_AFTER_EACH = false;127else {128cout << "Error on line 1: invalid value for output mode\n";129goto cleanup;130}131132//get truncate value133trunc = commandbuffer[commandbufferit++];134trim(trunc);135transform(trunc.begin(), trunc.end(), trunc.begin(), tolower);136137if(trunc == "trunc")138btrunc = true;139else if(trunc == "no trunc")140btrunc = false;141else {142cout << "Error on line 2: invalid value for truncate\n";143goto cleanup;144}145146//get filename prefix147prefix = commandbuffer[commandbufferit++];148trim(prefix);149150//get N, the number of symbols in the full shift151strInput = commandbuffer[commandbufferit++];152trim(strInput);153shiftSize = atoi(strInput.c_str());154155//get set of K, the periodicities to be tested156strInput = commandbuffer[commandbufferit++];157convertPeriodList(strInput, periodList);158159//Allocate composition object160comp = new Comp(shiftSize);161162//get "Begin Definitions"163strInput = commandbuffer[commandbufferit++];164trim(strInput);165166//make it lowercase167transform(strInput.begin(), strInput.end(), strInput.begin(), tolower);168169if(strInput.compare("begin definitions") != 0) {170cout << "ERROR: Expecting \"BEGIN DEFINITIONS\"\n";171goto cleanup;172}173174//read definitions175while(1) {176strInput = commandbuffer[commandbufferit++];177178//strInput2 is lowercased and trimmed179strInput2.resize(strInput.length());180transform(strInput.begin(), strInput.end(), strInput2.begin(), tolower);181trim(strInput2);182183//strInput3 is lowercased and NOT TRIMMED184strInput3.resize(strInput.length());185transform(strInput.begin(), strInput.end(), strInput3.begin(), tolower);186187if(strInput2 == "" || (strInput2[0] == '/' && strInput2[1] == '/')) {188//ignore blank lines and allow comment out189}190else if(strInput2 == "begin commands") {191break; //BREAK192}193else {194i = 0;195196while(i < strInput.length()) {197if(strInput[i] == '=') {198break;199}200201i++;202}203204if(strInput2[0] == 'o' && strInput2[1] == 'p' && strInput2[2] == 't') {205temp = strInput3.find("opt");206temp += 3;207208comp->addFunc(strInput.substr(i+1, strInput.length() - (i+1)),209strInput.substr(temp, i-temp), true);210}211else {212comp->addFunc(strInput.substr(i+1, strInput.length() - (i+1)),213strInput.substr(0,i), false);214}215}216}217218//read commands219while(commandbufferit < (int) commandbuffer.size()) { //read until end220strInput = commandbuffer[commandbufferit++];221trim(strInput);222223//ignore blank lines and allow comment out with "//"224if(strInput != "" && !(strInput[0] == '/' && strInput[1] == '/')) {225226//set composition227if(!comp->setComposition(strInput)) {228cout << "ERROR: Command \"" << strInput229<< "\" references undefined function\n";230goto cleanup;231}232233//check maxshiftperiod234int maxshiftperiod = 0;235for(i = 0; i < periodList.size(); i++)236maxshiftperiod = (maxshiftperiod < periodList[i] ?237periodList[i] : maxshiftperiod);238239//do experiment240if(shiftSize == 2241&& sizeof(unsigned int) == 4242&& !USEARRAY243&& maxshiftperiod <= 32) //if shift period goes > 32 use array244exhaustiveExpTwoShift(periodList, prefix, comp, btrunc,245recList, leastRecList,246commandbuffer);247else248exhaustiveExp(shiftSize, periodList,249prefix, comp, btrunc, recList, leastRecList,250commandbuffer);251252if(!OUTPUT_AFTER_EACH)253outcontrol(prefix, comp, btrunc,254recList, leastRecList, commandbuffer);255256//clear records (delete the nodes with ptrs in lists)257clearRecList(recList);258clearRecList(leastRecList);259recList.clear();260leastRecList.clear();261}262}263264//cleanup and exit265cleanup:266//delete the nodes with ptrs in lists267clearRecList(recList);268clearRecList(leastRecList);269270//delete composition (and the functions in the Comp object's storage)271delete comp;272273return 0;274}275276/*277* This function does the output. It chooses between full output278* and truncated output accordingly.279*/280281void outcontrol(string &prefix,282Comp * comp,283bool btrunc,284list<RecNode*> &recList,285list<RecNode*> &leastRecList,286vector<string> &commandbuffer) {287unsigned int commandbufferit=0;288289string test;290291ofstream fout;292string fname = prefix + comp->returnName() + ".txt";293fout.open(fname.c_str(), ios::trunc | ios::out);294295if(fout.is_open()) {296297//HEADER298fout << "FPeriod 3.67\n\n";299300//copy of the input301fout << "INPUT\n";302fout << "--------------------\n";303for(commandbufferit = 0; commandbufferit < commandbuffer.size(); commandbufferit++) {304test = commandbuffer[commandbufferit];305trim(test);306if(test.compare("") != 0) {307fout << commandbuffer[commandbufferit] << "\n";308}309}310fout << "--------------------\n";311fout << "\n";312313if(btrunc)314truncoutput(fout, comp, leastRecList);315else {316ofstream fout_preperiod; //file for preperiods317string fname_preperiod = prefix + comp->returnName() +318"_preperiods.txt";319fout_preperiod.open(fname_preperiod.c_str(), ios::trunc | ios::out);320321if(fout_preperiod.is_open()) {322323//HEADER324fout_preperiod << "FPeriod 3.67\n\n";325326//copy of the input327fout_preperiod << "INPUT\n";328fout_preperiod << "--------------------\n";329for(commandbufferit = 0; commandbufferit < commandbuffer.size(); commandbufferit++) {330test = commandbuffer[commandbufferit];331trim(test);332if(test.compare("") != 0) {333fout_preperiod << commandbuffer[commandbufferit] << "\n";334}335}336fout_preperiod << "--------------------\n";337fout_preperiod << "\n";338339output(fout, fout_preperiod, comp, recList, leastRecList);340341fout_preperiod.close();342}343else {344cout << "ERROR: Could not create file " << fname_preperiod << "\n";345}346}347348fout.close();349}350else {351cout << "ERROR: Could not create file " << fname << "\n";352}353}354355/*356* Truncated output357*358* Outputs only a table for least period k with the column P^(1/k) i.e.359*360* K P^(1/k)361* 1 1.000362* 2 1.4112363* 3 1.4432364* etc.365*/366367void truncoutput(ofstream & fout, Comp * comp, list<RecNode*> &leastRecList) {368list<RecNode*>::const_iterator it;369370fout << "Truncated output\n\n";371372//Output373fout << "COMMAND:\n";374fout << comp->returnName() << "\n";375fout << "WITH DEFINITIONS:\n";376comp->printDefinitions(fout);377fout << "\n";378379//Summary info380fout << "LEAST PERIOD K\n\n";381fout << "K P^(1/k)\n";382383//fill on right side, show in fixed notation384fout << setiosflags(ios::left | ios::fixed);385386for(it = leastRecList.begin(); it != leastRecList.end(); it++) {387fout << setw(3) << (*it)->currPeriod << " ";388389if((*it)->pending) {390fout << "PENDING" << "\n";391}392else {393fout << setw(8) << setprecision(4) << pow((*it)->numPeriodic, ((double)1)/((double)(*it)->currPeriod)) << "\n";394}395}396}397398/*399* Full output400*/401402void output(ofstream & fout, ofstream & fout_preperiod,403Comp * comp, list<RecNode*> &recList,404list<RecNode*> &leastRecList) {405list<RecNode*>::const_iterator it;406407/* Output file */408409//fill on right side, show in fixed notation410fout << setiosflags(ios::left | ios::fixed);411412fout << "Full output\n\n";413414fout << "COMMAND:\n";415fout << comp->returnName() << "\n";416fout << "WITH DEFINITIONS:\n";417comp->printDefinitions(fout);418fout << "\n";419420/* Preperiod file */421422//fill on right side, show in fixed notation423fout_preperiod << setiosflags(ios::left | ios::fixed);424425fout_preperiod << "\nPreperiod output\n\n";426427fout_preperiod << "COMMAND:\n";428fout_preperiod << comp->returnName() << "\n";429fout_preperiod << "WITH DEFINITIONS:\n";430comp->printDefinitions(fout_preperiod);431fout_preperiod << "\n";432433/* output file */434fout << "PERIOD K (NOT LEAST)\n\n";435436/* preperiod file */437fout_preperiod << "PERIOD K (NOT LEAST)\n\n";438439/* output file */440//Summary info441fout << "K FracPdc P^(1/k) L^(1/k) P L Not-P AvgPd AvgPrepd MaxPrepd\n";442443for(it = recList.begin(); it != recList.end(); it++) {444fout << setw(3) << (*it)->currPeriod << " ";445446if((*it)->pending) {447fout << "PENDING" << "\n";448}449else {450fout << setw(8) << setprecision(6) << (*it)->fractionPeriodic() << " ";451fout << setw(8) << setprecision(4) << pow((*it)->numPeriodic, ((double)1)/((double)(*it)->currPeriod)) << " ";452fout << setw(8) << setprecision(4) << pow((*it)->largestorbit, ((double)1)/((double)(*it)->currPeriod)) << " ";453fout << setw(9) << (*it)->numPeriodic << " ";454fout << setw(9) << (*it)->largestorbit << " ";455fout << setw(9) << (*it)->numNonPeriodic() << " ";456fout << setw(12) << setprecision(2) << (*it)->avgPeriod() << " ";457fout << setw(12) << setprecision(2) << (*it)->avgPreperiod() << " ";458fout << setw(9) << (*it)->maxPreperiod() <<"\n";459}460}461462//orbit and period multiplicity463fout << "\n";464fout << "K Period Mult(#orbits) Mult(#perpoints) Mult(#points)\n";465for(it = recList.begin(); it != recList.end(); it++) {466fout << setw(3) << (*it)->currPeriod << " ";467468if((*it)->pending) {469fout << "PENDING" << "\n";470}471else {472map<unsigned long long, unsigned long long>::const_iterator pitorbit, pitpoint, pitperpoint;473bool pfirst = true;474475//note that the period and orbit maps have the same periods in the same orders476// so we can iterate through both at once477for(pitorbit = (*it)->orbitMap.begin(),pitpoint = (*it)->periodMap.begin();478pitorbit != (*it)->orbitMap.end();479pitorbit++, pitpoint++) {480if(!pfirst) {481fout << " ";482}483else {484pfirst = false;485}486487fout << setw(11) << pitorbit->first << " ";488fout << setw(13) << pitorbit->second << " "; //orbits489490if((pitperpoint = (*it)->periodicPeriodMap.find(pitorbit->first))491!= (map<unsigned long long, unsigned long long>::const_iterator) (*it)->periodicPeriodMap.end()) {492fout << setw(16) << pitperpoint->second << " "; //periodic points493}494else {495fout << setw(16) << 0 << " "; //periodic points496}497498fout << setw(13) << pitpoint->second << "\n"; //points499}500}501}502503/* Preperiod file */504505//Preperiod multiplicity506fout_preperiod << "K Preperiod Mult\n";507for(it = recList.begin(); it != recList.end(); it++) {508fout_preperiod << setw(3) << (*it)->currPeriod << " ";509510if((*it)->pending) {511fout_preperiod << "PENDING" << "\n";512}513else {514515map<unsigned long long, unsigned long long>::const_iterator pit;516bool pfirst = true;517518for(pit = (*it)->preperiodMap.begin(); pit != (*it)->preperiodMap.end();519pit++) {520if(!pfirst) {521fout_preperiod << " ";522}523else {524pfirst = false;525}526527fout_preperiod << setw(11) << pit->first << " ";528fout_preperiod << setw(13) << pit->second << "\n";529}530}531}532533/* output file */534fout << "\n";535fout << "LEAST PERIOD K\n\n";536537/* preperiod file */538fout_preperiod << "\n";539fout_preperiod << "LEAST PERIOD K\n\n";540541/* output file */542//Summary info543fout << "K FracPdc P^(1/k) L^(1/k) P L Not-P AvgPd AvgPrepd MaxPrepd\n";544545for(it = leastRecList.begin(); it != leastRecList.end(); it++) {546fout << setw(3) << (*it)->currPeriod << " ";547548if((*it)->pending) {549fout << "PENDING" << "\n";550}551else {552fout << setw(8) << setprecision(6) << (*it)->fractionPeriodic() << " ";553fout << setw(8) << setprecision(4) << pow((*it)->numPeriodic, ((double)1)/((double)(*it)->currPeriod)) << " ";554fout << setw(8) << setprecision(4) << pow((*it)->largestorbit, ((double)1)/((double)(*it)->currPeriod)) << " ";555fout << setw(9) << (*it)->numPeriodic << " ";556fout << setw(9) << (*it)->largestorbit << " ";557fout << setw(9) << (*it)->numNonPeriodic() << " ";558fout << setw(12) << setprecision(2) << (*it)->avgPeriod() << " ";559fout << setw(12) << setprecision(2) << (*it)->avgPreperiod() << " ";560fout << setw(9) << (*it)->maxPreperiod() <<"\n";561}562}563564//orbit and period multiplicity565fout << "\n";566fout << "K Period Mult(#perpoints) Mult(#points)\n";567for(it = leastRecList.begin(); it != leastRecList.end(); it++) {568fout << setw(3) << (*it)->currPeriod << " ";569570if((*it)->pending) {571fout << "PENDING" << "\n";572}573else {574575map<unsigned long long, unsigned long long>::const_iterator pit, pitperpoint;576bool pfirst = true;577578for(pit = (*it)->periodMap.begin();579pit != (*it)->periodMap.end();580pit++) {581if(!pfirst) {582fout << " ";583}584else {585pfirst = false;586}587588fout << setw(11) << pit->first << " ";589590if((pitperpoint = (*it)->periodicPeriodMap.find(pit->first))591!= (map<unsigned long long, unsigned long long>::const_iterator) (*it)->periodicPeriodMap.end()) {592fout << setw(16) << pitperpoint->second << " "; //periodic points593}594else {595fout << setw(16) << 0 << " "; //periodic points596}597598fout << setw(13) << pit->second << "\n"; //points599}600}601}602603/* Preperiod file */604//Preperiod multiplicity605fout_preperiod << "K Preperiod Mult\n";606for(it = leastRecList.begin(); it != leastRecList.end(); it++) {607fout_preperiod << setw(3) << (*it)->currPeriod << " ";608609if((*it)->pending) {610fout_preperiod << "PENDING" << "\n";611}612else {613614map<unsigned long long, unsigned long long>::const_iterator pit;615bool pfirst = true;616617for(pit = (*it)->preperiodMap.begin(); pit != (*it)->preperiodMap.end();618pit++) {619if(!pfirst) {620fout_preperiod << " ";621}622else {623pfirst = false;624}625626fout_preperiod << setw(11) << pit->first << " ";627fout_preperiod << setw(13) << pit->second << "\n";628}629}630}631632//DONE633}634635void exhaustiveExp(unsigned int shiftSize, const vector<int> &periodList,636string &prefix, Comp *comp, bool btrunc,637list<RecNode*> &recList,638list<RecNode*> &leastRecList,639vector<string> &commandbuffer) {640unsigned int i = 0, j = 0; //k = 0; //counters641unsigned int currPeriod;642643list<StorageVal*> orbitList;644list<bool> leastPeriodList;645646list<RecNode*>::const_iterator recIterator, leastRecIterator;647648vector<StorageVal*> table;649650vector<int> divisorsList;651652if(DEBUG)653cout << "Array" << endl;654655//make a list of recNodes and leastRecNodes656//each list has one node per K value657for(i = 0; i < periodList.size(); i++) {658recList.push_back(new RecNode(periodList[i], pow(shiftSize,periodList[i])));659leastRecList.push_back(new RecNode(periodList[i], 0));660}661662//initialize recIterator and leastRecIterator663recIterator = recList.begin();664leastRecIterator = leastRecList.begin();665666for(i = 0; i < periodList.size(); i++) {667currPeriod = periodList[i];668byte itWord[currPeriod]; //iterating word669670//set table size671table.resize(pow(shiftSize,currPeriod));672673//record keeping674RecNode *record = *recIterator;675RecNode *leastRecord = *leastRecIterator;676677//zero the itWord678for(j = 0; j < currPeriod; j++) {679itWord[j] = 0;680}681682//Find divisors683findDivisors(currPeriod, divisorsList);684685while(true) {686byte cWord[currPeriod]; //curr word687byte output[currPeriod]; //helper word688unsigned long long period = 0, preperiod = 0;689StorageVal *cStore;690StorageVal *rt;691bool halt = false;692693while(halt != true &&694table[StorageKey(itWord,currPeriod).hash(shiftSize)] != NULL)695{696halt = iterateWord(itWord, currPeriod, shiftSize);697}698699if(halt == true)700break; //exit loop701702//Set current word703copyWord(cWord, itWord, currPeriod);704705cStore = new StorageVal(0);706707StorageKey sk1 = StorageKey(cWord, currPeriod);708709table[sk1.hash(shiftSize)] = cStore;710711orbitList.push_back(cStore);712if(leastPeriod(cWord, currPeriod, divisorsList))713leastPeriodList.push_back(true);714else715leastPeriodList.push_back(false);716717//DEBUG718//cout << "-- Start orbit --\n";719//printArray(cWord, currPeriod);720721//Perform experiment722j = 1;723while(true) {724comp->image(cWord, currPeriod, output);725copyWord(cWord, output, currPeriod);726727//DEBUG728//printArray(cWord, currPeriod);729730StorageKey sk2 = StorageKey(cWord, currPeriod);731732rt = table[sk2.hash(shiftSize)];733if(rt != NULL) {734break;735}736else {737cStore = new StorageVal(j);738739table[sk2.hash(shiftSize)] = cStore;740orbitList.push_back(cStore);741742if(leastPeriod(cWord, currPeriod, divisorsList))743leastPeriodList.push_back(true);744else745leastPeriodList.push_back(false);746}747748j++;749}750751//DEBUG752//cout << "-- Returned to seen --\n";753754if(rt->imageNum == -1) {755//we hit a previously found point756period = rt->period;757preperiod = rt->preperiod + j;758759list<StorageVal*>::iterator it;760list<bool>::iterator lpit = leastPeriodList.begin();761for(it = orbitList.begin(); it != orbitList.end(); it++) {762(*it)->imageNum = -1;763(*it)->period = period;764(*it)->preperiod = preperiod;765766if(*lpit == true) {767leastRecord->recordPeriod(period, preperiod);768leastRecord->numTrials++;769}770771lpit++;772773record->recordPeriod(period,preperiod);774775preperiod--;776777//NOTE: none of these points can be periodic778}779}780else {781//we circled back to an active point782period = j - rt->imageNum;783preperiod = rt->imageNum;784785record->recordOrbit(period);786787list<StorageVal*>::iterator it;788list<bool>::iterator lpit = leastPeriodList.begin();789for(it = orbitList.begin(); it != orbitList.end(); it++) {790(*it)->imageNum = -1;791(*it)->period = period;792(*it)->preperiod = preperiod;793794if(*lpit == true) {795leastRecord->recordPeriod(period, preperiod);796leastRecord->numTrials++;797if(preperiod == 0) {798leastRecord->recordPeriodicPeriod(period);799leastRecord->numPeriodic++;800}801}802803lpit++;804805record->recordPeriod(period, preperiod);806807if(preperiod > 0) {808preperiod--;809}810else {811record->recordPeriodicPeriod(period);812record->numPeriodic++;813}814}815}816817//cout << "Clearing list\n";818orbitList.clear();819leastPeriodList.clear();820}821822record->pending = false;823leastRecord->pending = false;824825//advance recIterator and leastRecIterator826recIterator++;827leastRecIterator++;828829//cout << "Clearing table\n";830clearStorageNodes(table);831table.clear();832833if(OUTPUT_AFTER_EACH) {834outcontrol(prefix, comp, btrunc,835recList, leastRecList, commandbuffer);836}837}838}839840void exhaustiveExpTwoShift(const vector<int> &periodList,841string &prefix, Comp *comp, bool btrunc,842list<RecNode*> &recList,843list<RecNode*> &leastRecList,844vector<string> &commandbuffer) {845unsigned int i = 0, j = 0; //k = 0; //counters846unsigned int currPeriod;847unsigned long long numTrials;848849list<StorageVal*> orbitList;850list<bool> leastPeriodList;851vector<StorageVal*> table;852853vector<int> divisorsList;854855list<RecNode*>::const_iterator recIterator, leastRecIterator;856857//DEBUG858if(DEBUG)859cout << "Special case 2-shift" << endl;860861//make a list of recNodes and leastRecNodes862for(i = 0; i < periodList.size(); i++) {863recList.push_back(new RecNode(periodList[i], pow(2, periodList[i])));864leastRecList.push_back(new RecNode(periodList[i], 0));865}866867//initialize recIterator868recIterator = recList.begin();869leastRecIterator = leastRecList.begin();870871for(i = 0; i < periodList.size(); i++) {872currPeriod = periodList[i];873unsigned int itWord = 0; //iterating word874875numTrials = (unsigned long long) pow(2, periodList[i]);876877table.resize(numTrials);878879RecNode *record = *recIterator;880RecNode *leastRecord = *leastRecIterator;881882//Find divisors883findDivisors(currPeriod, divisorsList);884885while(true) {886unsigned int cWord = 0; //curr word887unsigned int output = 0; //helper word888unsigned long long period = 0, preperiod = 0;889StorageVal *cStore;890StorageVal* rt;891892while(itWord != numTrials &&893table[itWord] != NULL) {894itWord++;895}896897if(itWord == numTrials)898break; //exit loop899900//Set current word901cWord = itWord;902903cStore = new StorageVal(0);904table[cWord] = cStore;905orbitList.push_back(cStore);906907if(leastPeriod(cWord, currPeriod, divisorsList))908leastPeriodList.push_back(true);909else910leastPeriodList.push_back(false);911912//DEBUG913//cout << "-- start orbit --\n";914//printUInt(cWord, 4);915916//Perform experiment917j = 1;918while(true) {919comp->image(&cWord, currPeriod, &output, 1);920921cWord = output;922923//DEBUG924//printUInt(cWord, 4);925926rt = table[cWord];927if(rt != NULL) {928break;929}930else {931cStore = new StorageVal(j);932table[cWord] = cStore;933orbitList.push_back(cStore);934935if(leastPeriod(cWord, currPeriod, divisorsList))936leastPeriodList.push_back(true);937else938leastPeriodList.push_back(false);939}940941j++;942}943944//DEBUG945//cout << "-- end orbit --\n";946947if(rt->imageNum == -1) {948//we hit a previously found point949period = rt->period;950preperiod = rt->preperiod + j;951952list<StorageVal*>::iterator it;953list<bool>::iterator lpit = leastPeriodList.begin();954for(it = orbitList.begin(); it != orbitList.end(); it++) {955(*it)->imageNum = -1;956(*it)->period = period;957(*it)->preperiod = preperiod;958959if(*lpit == true) {960leastRecord->recordPeriod(period, preperiod);961leastRecord->numTrials++;962}963964lpit++;965966record->recordPeriod(period,preperiod);967968preperiod--;969970//NOTE: none of these points can be periodic971}972}973else {974//we circled back to an active point975period = j - rt->imageNum;976preperiod = rt->imageNum;977978record->recordOrbit(period);979980list<StorageVal*>::iterator it;981list<bool>::iterator lpit = leastPeriodList.begin();982for(it = orbitList.begin(); it != orbitList.end(); it++) {983(*it)->imageNum = -1;984(*it)->period = period;985(*it)->preperiod = preperiod;986987if(*lpit == true) {988leastRecord->recordPeriod(period, preperiod);989leastRecord->numTrials++;990if(preperiod == 0) {991leastRecord->recordPeriodicPeriod(period);992leastRecord->numPeriodic++;993}994}995996lpit++;997998record->recordPeriod(period, preperiod);9991000if(preperiod > 0) {1001preperiod--;1002}1003else {1004record->recordPeriodicPeriod(period);1005record->numPeriodic++;1006}1007}1008}10091010//cout << "Clearing list\n";1011orbitList.clear();1012leastPeriodList.clear();1013}10141015record->pending = false;1016leastRecord->pending = false;10171018//advance recIterator and leastRecIterator1019recIterator++;1020leastRecIterator++;10211022//cout << "Clearing table\n";1023clearStorageNodes(table);1024table.clear();10251026if(OUTPUT_AFTER_EACH) {1027outcontrol(prefix, comp, btrunc,1028recList, leastRecList, commandbuffer);1029}1030}1031}10321033void convertPeriodList(const string &strInput,1034vector<int> &periodList) {1035vector<string> strList; //all terms1036vector<string> subList; //vector representation of an "[a,b]" term1037unsigned int strListLength = 0, i = 0;10381039split(strInput,';', strList);1040strListLength = strList.size();10411042for(i = 0; i < strListLength; i++) {1043string *temp = &(strList[i]);10441045trim(*temp);10461047if((*temp)[0] == '[') {1048int j, low, high;10491050temp->erase(0,1); //remove '['1051temp->erase(temp->length() - 1, 1); //remove ']'10521053split(*temp, ',', subList);10541055trim(subList[0]);1056trim(subList[1]);10571058low = atoi(subList[0].c_str());1059high = atoi(subList[1].c_str());10601061for(j = low; j <= high; j++) {1062periodList.push_back(j);1063}10641065subList.clear();1066}1067else {1068periodList.push_back(atoi(temp->c_str()));1069}1070}1071}10721073/*1074* Delete the nodes pointed to by ptrs in recList1075*/1076void clearRecList(list<RecNode*> &recList) {1077list<RecNode*>::iterator it;10781079for(it = recList.begin(); it != recList.end(); it++) {1080delete (*it);1081}1082}10831084/*1085* Delete the nodes pointed to by the values in the table1086*/1087void clearStorageNodes(vector<StorageVal*> &table) {1088vector<StorageVal*>::iterator it;10891090for(it = table.begin(); it != table.end(); it++) {1091delete (*it);1092}1093}10941095/*1096* Find divisors of X, excluding X itself1097*/1098void findDivisors(unsigned int x, vector<int> &divisorsList) {1099unsigned int i = 2;1100unsigned int halfX = x/2;11011102divisorsList.clear();11031104if(x > 1) {1105divisorsList.push_back(1);1106}11071108for(i = 2; i <= halfX; i++) {1109if(x % i == 0) {1110divisorsList.push_back(i);1111}1112}1113}11141115/*1116* returns true if this word has least period p1117*/1118bool leastPeriod(byte *word, unsigned int currPeriod,1119const vector<int> &divisorsList) {1120bool result = false;1121unsigned int i, j, k;11221123if(currPeriod == 1)1124return true;11251126//for every divisor d1127for(i = 0; i < divisorsList.size(); i++) {1128unsigned int d = divisorsList[i];11291130result = false;11311132//check in chunks of length d1133for(j = 0; j < currPeriod; j+= d) {1134//check that this chunk is same as other chunks1135for(k = 0; k < d; k++) {1136if(word[k + j] != word[k]) {1137result = true;1138break;1139}1140}11411142if(result == true)1143break;1144}11451146if(result == false)1147break;1148}11491150return result;1151}11521153/*1154* returns true if this word has least period p1155*/1156bool leastPeriod(unsigned int word, unsigned int currPeriod,1157const vector<int> &divisorsList) {1158bool result = false;1159unsigned int i, j, k;11601161if(currPeriod == 1)1162return true;11631164//for every divisor d1165for(i = 0; i < divisorsList.size(); i++) {1166unsigned int d = divisorsList[i];11671168result = false;11691170//check in chunks of length d1171for(j = 0; j < currPeriod; j+= d) {1172//check that this chunk is same as other chunks1173for(k = 0; k < d; k++) {1174if(((word >> (k + j)) & 1) != ((word >> k) & 1)) {1175result = true;1176break;1177}1178}11791180if(result == true)1181break;1182}11831184if(result == false)1185break;1186}11871188return result;1189}119011911192119311941195119611971198119912001201120212031204120512061207