/*
 *  This program generates OEIS sequence A191837
 *  Author: Jonathan Stauduhar (jstdhr@gmail.com)
 *  Date: July, 2011
 */
#include <stdio.h>
#include <stdlib.h>
#include <limits.h>

int nCr(unsigned long* set, unsigned long nSet, unsigned long rSubset);
int is_prime(unsigned long n);

enum {FALSE,TRUE};
unsigned long* pset;
unsigned long m;
int summands=4;
int INCR=6; /* As m increases in size, INCR can be increased to 30,210,2310,... */

int main(int argc, char* argv[])
{
    unsigned long i, p, q, halfm, *indices, pset_count, total, Ysum;
    int summand_ct;
    //all integers in the sequence are 0(mod 6), so incr. m by muliples of 6
    for(m=6;m<ULONG_MAX;m+=INCR)
    {
        halfm = m/2;
        pset = (unsigned long*)malloc(sizeof(unsigned long)*halfm);
        /* Load all primes p < m/2, for which m-p is prime, into an array. 
         * Since all m's are congruent to 0 mod 6, skip 2 and 3.
         */
        for(p = 5, pset_count=0, Ysum=0, summand_ct=0; p < halfm; p+=2)
        {
            if(is_prime(p) == TRUE)
            {
                if(is_prime(m - p) == TRUE)
                {
                    /* The following two if-statements optimize the search.  
                     * Consider the following:  We need X summands.  If we 
                     * currently have fewer than X summands and their sum already 
                     * exceeds m then there can be no solution for this m, so move 
                     * to the next m.
                     * Furthermore, if we have X-1 summands that sum to Ysum <= m, 
                     * we need only continue adding new primes to the array 
                     * until Ysum+p > m -- all subsequent p's need not be considered.
                     */
                    if(Ysum + p > m) break;
                    if(summand_ct++ < summands-1) Ysum += p;
                    pset[pset_count] = p;
                    pset_count++;
                }
            }
        }
        /* Did we find enough primes to allow for a solution? */
        if(pset_count >= summands)
        {
           /* "indices" is a helper array:  Since nCr() modifies the values 
            * of the array passed to it and we don't want to modify the values 
            * of the primes, create an array of indices of the current array of 
            * primes and pass that to nCr(). The combinations of array index 
            * positions are mapped back to the array of primes when testing for 
            * solutions.
            */
            indices = (unsigned long*) calloc(pset_count, sizeof(unsigned long));
            for(i=0; i < pset_count; ++i) indices[i] = i;
            while( nCr(indices,pset_count,summands) );
            free(indices);
        }
        free(pset);
    }
}


/* 
 * Modified version of Miroslaw M. Soja's R-Combination Generator code.
 * Combinations are produced in lexicographic order.
 */
int nCr(unsigned long* set, unsigned long nSet, unsigned long rSubset)
{
    unsigned long ct, total=0;
    unsigned long indx = (nSet - 1);
    unsigned long i    = (rSubset - 1);

    if(set[0]==(nSet - rSubset))
    {
        /* If there is just one possible solution, it will be caught here. */
        for(ct = 0;ct < rSubset;ct++)
        {
            total += pset[set[ct]];
        }
        if(total == m)
        {
            printf("%d ", rSubset/2);
            printf("%lu",m);
            printf("\n");
            summands +=2;
            m = m-INCR; /* There can be multiple solutions for m, so recycle it. */
            return(FALSE);
        }
        else return (FALSE);
    }

    while( set[i]==indx && i>0 )
    {
        --i;
        --indx;
    }
    set[i]+=1;
    i++;
    
    while(i<rSubset)
    {
        set[i]=set[i-1]+1;
        ++i;
    }
    /* Does this combination provide a solution? */
    for(ct = 0;ct < rSubset;ct++)
    {
        total += pset[set[ct]];
        if(total > m) break;
    }
    if(total == m)
    {
        printf("%lu ", rSubset/2);
        printf("%lu",m);
        printf("\n");
        summands +=2;
        m = m-INCR; /* There can be multiple solutions for m, so recycle it. */
        return (FALSE);
    }
    return (TRUE) ;
}

int is_prime(unsigned long n)
{
    int is_prime = TRUE;
    unsigned long i;
    /* Check the input. */
    if (n <= 0)
    {
        printf ("Illegal input to is_prime()!\n");
        return (FALSE);
    }
    if(n == 1) is_prime = FALSE;
    else if(n == 2)
    {
        return is_prime;
    }
    else if(n%2 == 0)
    {
        is_prime = FALSE;
        return is_prime;
    }
    else
    {
        i=3;
        while(TRUE)
        {
           if( (i*i > n) || (n%i == 0) )
           {
              is_prime = FALSE;
              break;
           }
           i+=2;
        }
        if(i*i > n)
        {
            is_prime = TRUE;
            return is_prime;
        }
        else
        {
           is_prime = FALSE;
           return is_prime;
        }
    }
    return (is_prime);
}