/*
et2ji : to approximate equal tempered tones by just intonation intervals
  usage:
	 et2ji #tones/oct [tolerance(cents) [show all?(Y/N)] ]
  the first argument is the number of tones per octave (may be noninteger)
  the second argument is the tolerance in cents [optional; default is 1.0]
  the third argument [optional; default is No], if present means
	 show all intermediate approximations as well as best ones
	 (warning: this may make output lines very long!)
  the output, for each tempered tone, gives tone number, cents, and
  successive JI approximations with cents error.
  NOTE: these JI intervals do NOT include all simple nearby ones.
  I got this algorithm from John Chalmers; thanks John!
*/

#include <stdio.h>
#include <math.h>

int showall=0; /* flag to show intermediate results */
static double ln1cent;  /* natural log of 1 cent */
#define cents2dec( cents ) ( exp( ln1cent*cents ) )
#define dec2cents( dec ) ( log( dec )/ln1cent )
#define rat2cents( n, d ) ( dec2cents( (double)(n)/(double)(d) ) )

void main (int argc, char *argv[])
{
double ntones=12., step, cents, maxerr, minerr=1.;
int i;
void cents2rat ();

ln1cent = log(2.)/1200.;  /* initialize to natural log of 1 cent */

/* get # of tones & tolerance from command line */
if (argc > 1) sscanf(argv[1],"%lf",&ntones);
else {
  fprintf(stderr,
	 "usage: et2ji #tones/oct [tolerance(cents) [show all?(Y/N)] ]\n");
  return;
  }
if (argc > 2) sscanf(argv[2],"%lf",&minerr);
if (argc > 3) showall = !(argv[3][0]=='N'||argv[3][0]=='n');
/* some error checking: */
ntones = fabs(ntones); if (ntones>1200.) ntones = 1200.; if(ntones<1.)return;
minerr = fabs(minerr); if (minerr<0.01) minerr = 0.01;

/* print out table, step by step */
step=1200./ntones; maxerr=step/2.;
printf("%g-tone equal temperament: just approximations to %g cent%c%s\n",
  ntones,minerr, (minerr>1.0) ? 's' : ' ',(showall) ? " (detailed)" : "" );
for ( i=1,cents=step; i<ntones; i++, cents+=step ) {
  printf("%3d (%6.1f):",i,cents);
  cents2rat(cents,maxerr,minerr);
  }
}

void cents2rat ( double cents, double maxerr, double minerr ) {
double big, small = 1., err, td;
int bn = 0, bd = 1, sn = 1, sd = 0, t;
/* bn,bd,sn,sd are big numer & denom, small numer & denom */

big = cents2dec(cents); /* decimal fraction (freq. ratio) */
do {
  while ( big >= small ) {
	 big -= small; bn += sn; bd += sd;
/* NOTE: this inner loop generates other approximations bn/bd
that are less accurate for a given complexity level */
	 if ( showall && big >= small ) {
		err = rat2cents( bn, bd ) - cents;
		if ( fabs(err) <= maxerr )
		  printf(" %d/%-d (%+.1f)", bn, bd, err);
		}
	 }
  err = rat2cents( bn, bd ) - cents;
  if ( fabs(err) <= maxerr )
	 printf(" %d/%-d (%+.1f)", bn, bd, err);
/* here, big < small, time to swap */
  td = big; big = small; small = td;
  t = bn; bn = sn; sn = t;
  t = bd; bd = sd; sd = t;
  }
while ( fabs(err) > minerr ); /* keep on til within tolerance */
putchar('\n');
}