/*-------------- Telecommunications & Signal Processing Lab ---------------
                             McGill University

Routine:
  void RSratio (double Sratio, int Ir, double *Ns, double *Ds, long int NsMax,
                long int DsMax);

Purpose:
  Find a rational approximation to a given value

Description:
  This procedure finds a rational approximation to a given value Sratio,
    Sratio = Ns / Ds.
  Ns and Ds are doubles which take on integer values if Sratio can be
  represented exactly as a ratio of integers.  Otherwise, Ds is set to one and
  Ns is set equal to Sratio.

  In addition, if possible, Ns and Ds are modified by multiplying each by the
  same value, such that the resultant Ns is divisible by Ir.

Parameters:
   -> double Sratio
      Input value to be represented
   -> int Ir
      Input value.  If possible the output value Ns will be such that Ir is
      a factor of Ns.
  <-  double *Ns
      Output numerator value for the rational representation
  <-  double *Ds
      Output denominator value for the rational representation
      Integral part of the result
   -> long int NsMax
      Maximum value for Ns
   -> long int DsMax
      Maximum value for Ds

Author / revision:
  P. Kabal  Copyright (C) 1996
  $Revision: 1.6 $  $Date: 1996/08/14 19:04:02 $

-------------------------------------------------------------------------*/

static char rcsid[] = "$Id: RSratio.c 1.6 1996/08/14 AFsp-V2R1 $";

#include <libtsp.h>
#include "ResampAudio.h"

void
RSratio (Sratio, Ir, Ns, Ds, NsMax, DsMax)

     double Sratio;
     int Ir;
     double *Ns;
     double *Ds;
     long int NsMax;
     long int DsMax;

{
  long int Lgcd, LNs, LDs, M;
  double Sra;

  if (Sratio <= 0.0 || Ir <= 0)
    UThalt ("Error: Invalid Sratio and/or Ir");

  /* Find a rational approximation to Sratio, Sratio = Ns / Ds */
  MSratio (Sratio, &LNs, &LDs, 0.0, NsMax, DsMax);
  Sra = (double) LNs / LDs;

  if (Sratio != Sra) {

    /* Inexact */
    *Ds = 1.0;
    *Ns = Sratio;
    printf (" Sampling ratio: %g\n", Sratio);
  }

  else {

    /* Exact ratio of integers */
    /* For ordinary interpolation Ns = Ir.  In other cases, we multiply Ns and
       Ds by M such that the resultant Ns is divisible by Ir.
    */
    Lgcd = FNgcd (LNs, (long int) Ir);
    M = Ir / Lgcd;
    if (LDs <= DsMax/M && LNs <= NsMax/M) {
      *Ds = M * LDs;
      *Ns = M * LNs;
    }
    else {
      *Ds = LDs;
      *Ns = LNs;
    }
    printf (" Sampling ratio: %ld/%ld\n", LNs, LDs);
  }

  return;
}