module Math

static VALUE find_fracs(VALUE mod, VALUE rv, VALUE dv)
{
  N_TYPE m[2][2], ai, maxden = R2N(rb_Integer(dv));
  double startx, x = RFLOAT_VALUE(rb_Float(rv));
  int sign = 1;

  if (maxden <= 0)
    rb_raise(rb_eArgError, "maximum denominator should be > 0");

  if (x < 0) {
    sign = -1;
    x = -x;
  }

  startx = x;

  /* initialize matrix */
  m[0][0] = m[1][1] = 1;
  m[0][1] = m[1][0] = 0;

  /* loop finding terms until denom gets too big */
  while (m[1][0] *  ( ai = (N_TYPE)x ) + m[1][1] <= maxden) {
    N_TYPE t;
    t = m[0][0] * ai + m[0][1];
    m[0][1] = m[0][0];
    m[0][0] = t;
    t = m[1][0] * ai + m[1][1];
    m[1][1] = m[1][0];
    m[1][0] = t;
    if(x==(double)ai) break;     // AF: division by zero
    x = 1/(x - (double) ai);
    if(x>(double)0x7FFFFFFF) break;  // AF: representation failure
  }

  {
    /* now remaining x is between 0 and 1/ai */
    /* approx as either 0 or 1/m where m is max that will fit in maxden */
    /* first try zero */
    VALUE num1, den1, err1, num2, den2, err2;

    num1 = N2R(sign*m[0][0]);
    den1 = N2R(m[1][0]);
    err1 = rb_float_new(startx - ((double) m[0][0] / (double) m[1][0]));

    /* now try other possibility */
    ai = (maxden - m[1][1]) / m[1][0];
    m[0][0] = m[0][0] * ai + m[0][1];
    m[1][0] = m[1][0] * ai + m[1][1];

    num2 = N2R(sign*m[0][0]);
    den2 = N2R(m[1][0]);
    err2 = rb_float_new(startx - ((double) m[0][0] / (double) m[1][0]));

    return rb_ary_new3(6, num1, den1, err1, num2, den2, err2);
  }
}

# File lib/frac.rb, line 9
def frac(float, maxden)
  arr = find_fracs(float, maxden)
  arr[2].abs > arr[5].abs ? Rational(arr[3], arr[4]) : Rational(arr[0], arr[1])
end