#include <stdio.h> #define _likely(x) (x) //Define to hints if possible in your compiler #define _unlikely(x) (x) //and worth doing #define iflike(x) if(_likely(x)) #define ifunlike(x) if(_unlikely(x)) #define uchar unsigned char #define uint unsigned int /*Types: The number of elements in a sequence is <=NMAX. Therefore, uchar The elements themselves and the scaled averages. Num The number of subsets is 2^nMAX, and so the masks for it too. Subset */ typedef uint Num; typedef int sNum; //signed num typedef uint Subset; //short is enough, but int is more efficient //The program stops when either a different average set with NMAX+1 elements is found //or the highest element from the sets beign checked exceeds HIGHESTNUM. #define NMAX 9 #define HIGHESTNUM 170 //At any moment during execution N is the number of elements in the sets beign checked. //For each such set found to be a different average set, it is checked whether an extra element //can be inserted between elements 1 and the second-lowest element so that the resulting //set is an N+1-element different average set. If so N will be set to N+1 and the process continues. static uchar N; //Table with the values of the lowest bit set in the binary representation of 1,2,3,4... 2^(NMAX-1). #define M(k,S) S,k,S //M from Mirror static const uchar TABLElb[1+(1<<(NMAX-1))]={0, M(8,M(7,M(6,M(5,M(4,M(3,M(2,1))))))), 9}; #undef M //Multipliers for the scaled averages (1+a2+... +ak)*(COMMON_MULT/k) #define COMMON_MULT (8*9*5*7) //2520 #define C(n) COMMON_MULT/n static const uint Factors[NMAX+2]={ 1, C(1),C(2),C(3),C(4),C(5),C(6),C(7),C(8),C(9),C(10) }; #undef C /*Insertar[n] is the last n-sequence found. This is the first n-sequence which admits the insertion of an extra element at the low end as explained above, the arising sequence being a different average one, namely, the first (n+1)-sequence found. So, example gratia, Insertar[2] will be {4,1}; Insertar[3] {7,5,1} and so forth. The last two sequences will be asigned but will never be used.*/ static Num Insertar[NMAX+1][NMAX]; /* Sequences are stored in variables named 'nums', from highest element to lowest. Basic Sequence: A sequence ending in a 1 with no three consecutive numbers. Only basic sequences are considered. Checked Sequence: A sequence that passed the test. */ static inline uchar next_basic(Num nums[], uchar n); //Returns (in nums) the next Basic Sequence void nextset(Num nums[NMAX]); //Returns the next Checked Sequence uchar checkseq(const Num nums[NMAX]); //Checks the sequence and returns zero or the index of the first num that has to be changed Num try_insertnum(const Num nums[NMAX], uchar n); //Sees whether a number may be inserted at the low end of the sequence void print_seq(const Num nums[]); int main(void){ Num nums[NMAX]; //The numbers of the set, from highest to lowest. Therefore, its last Num oldnmax; //element will always be 1. //getc(stdin); freopen("medias_an.txt","w",stdout); N=2; nums[0]=4; nums[1]=1; oldnmax=nums[0]; goto first_entry; while(nums[0]<=HIGHESTNUM){ ifunlike(nums[0]!=oldnmax){ oldnmax=nums[0]; fprintf(stderr,"%u ",nums[0]); } first_entry:; Num n1=try_insertnum(nums,N); ifunlike(n1!=0){ fprintf(stderr,"\n\n%u ",nums[0]); Num *pins, *pn; for(pins=Insertar[N], pn=nums; *pn!=1; *pins++=*pn++); *pins=1; *pn++=n1; N++; *pn=1; print_seq(nums); fflush(stdout); //So that the program can be interrupted at any moment if(N>NMAX) break; } nextset(nums); } return 0; } //Assumes nums has at least three elements //The number n must be 2 units less than the elements of nums. //Returns the number of elements in the sequence static inline uchar next_basic(Num nums[], uchar n){ Num *pn; pn=nums+n; nums+=2; if(*(pn-1)==*pn+1) pn--, (*pn)++; else (*pn)++; while(pn>=nums && *(pn-2)==*pn+2) pn-=2, (*pn)++; //if(pn>=nums-1 && *(pn-1)+1==*pn<<1) (*pn)++; //avoid last three elements like 2k-1, k, 1. Not worth doing pn[1]=1; return 4+(uchar)(pn-nums); } /*Returns the next Checked Sequence. Assumes nums has at least three elements. It first calls next_basic(nums) to get the next candidate. If the test checkseq for that candidate fails (which happens most often) the function checkseq will return the index k of the first number from nums whose change is unavoidable. Upon that (failure) this function clears the numbers below the k'th one and calls next_basic, whereby the k'th element, and possibly others, is increased. A call to next_basic may return a sequence with less elements than the input one. Therefore, after each call to that function, even if not after the deletion of elemens that follow a failure of checkseq(), the returned sequence has to be filled in with fresh elements in order to complete an N-element Basic Sequence. To facilitate that task next_basic returns the number of elements of the new nums. Let that number be n. Therefore, N-n=dif elements have to be inserted. The inserted elements will be the ones from Insertar[dif+1] (the lowest of which will be the closing 1 and will therfore not get "inserted" because it is already present in nums) because any (dif+1)- sequence below that one will not admit the insertion of an extra element at its lower end and is therefore not worth considering. */ void nextset(Num nums[NMAX]){ uchar n; {Num *pn; for(pn=nums+2;*pn!=1;pn++); n=(uchar)(pn-nums)-1;} redo: n=next_basic(nums,n); if(n!=N){ uchar dif=N-n; //if(n>2 && nums[n-2]==Insertar[dif+1][0]+1 && nums[n-3]==Insertar[dif+1][0]+2) goto redo; //This can never happen if(dif==1) nums[n-1]=4; else{ if(nums[0]>HIGHESTNUM) return; //Otherwise the array mbits from checkseq will not have space enough Num *pn=nums+n-1, *pins=Insertar[dif+1]; do *pn++=*pins++; while(*pins!=1); } } nums[N-1]=1; ifunlike((n=checkseq(nums))==0) return; if(n<N-2) nums[n+1]=1; goto redo; } /* The function checkseq returns 0 if the series passes the test, otherwise the index of the first num that has to be changed. The last element of the series in nums must be 1 and nums must be different from the last time the functions was called. Since successive calls almost always share the same higher values for nums, this function remembers its state from call to call, so as to test only the numbers that changed. The array mbits is declared outside the function because it is also used by try_insertnum. The program expends most of its time inside this function. */ typedef uint Bits; static Bits mbits[(COMMON_MULT*HIGHESTNUM)>>5]; //The different scaled medias computed, stored as bits. uchar checkseq(const Num _nums[NMAX]){ static Num prev_nums[NMAX]; //nums from the last time this function was called. //Its values at program startup are zero, which is right because //thus they test to false against the nums of the first call. static uint scaled[1U<<NMAX]; //The scaled averages. A value of zero marks the end. static const Subset scaled_ppios[NMAX]={0,2,6,14,30,62,126,254,510}; //Indexes into scaled. Where do //the computed scaled for each num start, so that those for highest //nums (first elements of nums) are retained between calls. uchar d; //index of first different num between prev_nums and nums sNum nums[1+NMAX]; //local non-const, signed copy of nums, with a dummy element at the beginning. //Finding of d {const Num *pnums=_nums; Num *pprev=prev_nums; while(*pnums==*pprev) pnums++,pprev++; d=(uchar)(pprev-prev_nums); pnums--; do *pprev++=*++pnums; while(*pnums!=1);} //copy of _nums. _nums will no longer be used nums[0]=0; //For the case d=0 at media=nums[d] below {sNum *pb; for(pb=nums+1; *_nums!=1; *pb++=*_nums++); *pb=1;} //Clear the scaled from the previous call above scaled_ppios[d] for(uint *p=scaled+scaled_ppios[d]; *p!=0; p++) mbits[(*p)>>5] &= ~ (1<<((*p)&31)); uint *pscaled; //Pointer into scaled; keeps increasing by one. const uchar *ptable; //Only one element is added or removed from the subset at a time uint media; const uint *esc; pscaled=scaled+scaled_ppios[d]; esc=Factors; media=nums[d]; nums[d]=-nums[d]; //If d=0 the set is set empty here so that it will be initialized to {nums[1]} esc+=(d!=0); //in the first iteration; If d!=0 the set is set to 0...010 {nums[d]} and it will ptable=TABLElb+(1<<d); //be set to 0...011 ({nums[d],nums[d+1]}) in the first iteration. for(;;){ uint m; //scaled media Bits *byte; //Retrive the number nums[*ptable] that will be added to or removed from the subset beign tested {sNum *pnum=nums+ *ptable++; //*ptable>=1 sNum k=*pnum; ifunlike(k==1) goto out0; media+=k; //If k is positive the number is beign added, otherwise it is beign removed uint test=(k>=0); esc+=test; esc-=!test; *pnum=-k;} //Hence, if the number at *pnum was added, the next time it is //encountered it will be removed, and vice-versa. m=media**esc; *pscaled=m; //pscaled++ not here, in case the test is not passed byte=mbits+(m>>5); m=1<<(m&31); if(*byte & m) break; *byte |= m; pscaled++; //test passed. Validate content of pscaled (by increasing the pointer) //The same with number one added to the subset; m=(media+1)*esc[1]; *pscaled=m; byte=mbits+(m>>5); m=1<<(m&31); if(*byte & m) break; *byte |= m; pscaled++; } *pscaled=0; Subset nmedias=(Subset)(pscaled-scaled); uchar ret=2; while(scaled_ppios[ret]<=nmedias) ret++; return ret-1; out0: *pscaled=0; return 0; } Num try_insertnum(const Num nums[NMAX], uchar n){ uint scaledthis[1<<NMAX]; Num n1, n2; const Num * const p=nums+n; //marks the end n2=1+(nums[0]-nums[1])+1; //The first 1 is nums[n-1]. The second one is because the test will be <n2 (not <=n2). n1=2; if(nums[n-2]<n2) n2=nums[n-2]; if(n>2 && nums[n-3]==nums[n-2]+1){ n1=3; if(n2==nums[n-2]) n2--; } for(;n1<n2; n1++){ uint *pscaled; Subset word, order; uint media; const uint *esc; if(mbits[n1>>5] & (1<<(n1&31))) continue; pscaled=scaledthis; word=order=0; media=n1; esc=Factors+1; do{ order++; uint m; Bits *byte; Subset testw; const Num *pnum; for(testw=1, pnum=nums; !(order & testw); testw<<=1, pnum++); ifunlike(pnum==p){order=0; break;} if(word & testw) media-=*pnum, esc--; else media+=*pnum, esc++; word^=testw; m=media**esc; *pscaled=m; byte=mbits+(m>>5); m=1<<(m&31); if(*byte & m) break; *byte |= m; pscaled++; }while(1); *pscaled=0; for(pscaled=scaledthis; *pscaled!=0; pscaled++) mbits[(*pscaled)>>5] &= ~ (1<<((*pscaled)&31)); ifunlike(order==0) return n1; } return 0; } void print_seq(const Num nums[]){ const Num *pn; for(pn=nums;*pn!=1;pn++); printf("%u\t1",1+(uint)(pn-nums)); do{pn--; printf(" %u",*pn);}while(pn>nums); fputc('\n',stdout); } /*Previous versions of this program kept also a divider in addition to a multiplier in order to avoid multiplying by, say, 7. COMMON_MULT would then be 360. In that setting, 7-element sets were checked for being multiples of 7. If they were not, they were discarded. This can be done because the average of a subset of 7 elements the sum of which is not a multiple of 7 cannot equal the average of a set of not-7 elements; as for other 7-element sets that may have the same sum, both sets would share some elements, as long as there are less than 14 elements, and those elements removed yield two other sets whose averages will be identified as equal. Not only prime numbers can be spared but also the highest prime entering a prime power. 8-element subsets would be checked for beign multiples of 2, 9-element ones for beign multiples of 3, etc. The space saved would be huge if NMAX increased a few units. But NMAX does not increase even a few units (the time it takes to find a(n) is like 2^(n^2)), and removing the dividers makes the code simpler. */