/* -*- Mode: C; indent-tabs-mode: t; c-basic-offset: 4; tab-width: 4 -*-  */
/*
 * a098550.c
 * Copyright (C) 2015 Dev Gualtieri <gualtieri@ieee.org>
 * 
 * a098550 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.
 * 
 * a098550 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/>.
 */
 
 /*
 Produces integer sequence A098550, defined as a(n) = n if n <= 3; 
 otherwise the smallest number not occurring earlier having at least
 one common factor with a(n-2), but none with a(n-1).
 The first few elements are 1, 2, 3, 4, 9, 8, 15, 14, 5, 6, 25, 12, 35,
 16, 7, 10, 21, 20, 27, 22, 39, 11, 13, 33, 26, 45, 28, 51, 32, 17, 18,
 85, 24, 55, 34, 65, 36, 91, 30, 49, 38, 63, 19, 42, 95, 44, 57, 40, 69
 */

#include <stdio.h>
#include <stdlib.h>
#include <limits.h>

/* Prototypes */
char *strcpy(char *dest, const char *src);
void exit(int status);
/* end of prototypes */

#define limit 10000
//note that the limit must be less than one-seventh of integer_max
//see comment on b[] array, below
int i,j,m,n;
int *a; //the sequence array
int *b; //array of trials
int a_size = 0;
char fn1[64];
FILE *outdata;

//greatest common divisor function
int gcd ( int a, int b )
{
int c;
while ( a != 0 )
{
c = a; a = b%a;  b = c;
}
return b;
}

int main(int argc, char *argv[])
{

if (argc<2)
{
strcpy(fn1,"output.txt");
}
else
{
strcpy(fn1,argv[1]);
}

printf("\nOutput file selected = %s\n",fn1);

if ((outdata = fopen(fn1,"w"))==NULL)
	{printf ("\nOutput file cannot be opened.\n");
	exit (1);}
	
printf("\nInteger max = %d\n",INT_MAX);

if(limit>INT_MAX)
{
printf("Calculation exceeds integer maximum... Exiting.\n");
exit(1);
}

a = malloc(sizeof(int)*limit);

//The points lie on roughly straight lines, whose slopes are approximately
//0.467, 0.957, 1.15, 1.43, 2.40, 3.38, 5.25 and 6.20, so our flag array
//needs to be at least 6.2 larger than the sequence array

b = malloc(sizeof(int)*7*limit);

//initialize first few sequence elements
a[1] = 1;
a[2] = 2;
a[3] = 3;
a[4] = 4;
a_size = 4;

//initialize array of trial numbers

for(i=0;i<(7*limit);i++)
{
b[i]=i;
}
//zero indicates number already included in sequence
b[1] = 0;
b[2] = 0;
b[3] = 0;
b[4] = 0;

while(a_size<limit)
{

for(i=0;i<(7*limit);i++)
{

if((gcd(b[i],a[a_size])==1)&&(gcd(b[i],a[a_size-1])>1))
{
a_size++;
a[a_size] = i;
printf("\t%d",i);
b[i]=0; // flag number as included in sequence
break;
}
}
}

//print sequence
for(j=1;j<=a_size;j++)
{
//printf("\t%d",a[j]);
fprintf(outdata,"%d\n",a[j]);
}
printf("\n");

free(a);
free(b);

fclose(outdata);

return (0);
}