Sumantra Kundu wrote:
>              congestion prob =    [{(a/b)* log (a/b)}/(log (b)] * (c/d).
> I know that its not possible to use floating point ops since the
> ...
> apps)...However, I do need **accurate** value of this quantity since

The chief advantage floating-point brings to the table is not the
ability to handle fractions, but the ability to handle wide dynamic
ranges automatically (say, when you have numbers like 6.626e-27 and
6.602e+23 in the same calculation and you don't want to futz with scale
adjustments manually).  As long as you don't need that--say, your a's
and b's are m-bit integers so you know your fractions can only be in the
range 2^-m to 2^m or zero, then it's not hard to use integers
throughout, with a fixed scaling. For example (((int2m_t)a)<<m)/b is a
2m-bit number a/b with m bits to the right of the radix point. If you
know 0 <= a/b < 1, you can cast the result back to intm_t to have an
m-bit fraction.

All you need now is an approximation to the log function. You don't need
anything as involved as the log in libm, because that's designed for
54-bit accuracy over the whole dynamic range of possible doubles, and
all you need is an approximation good enough for your purposes over
the interval of values you expect to enounter. (That does mean you need
to start by deciding how accurate is good enough for your purpose.)

A couple of techniques are: look up in a table and interpolate, or
evaluate a polynomial that fits within your accuracy requirements over
the interval you care about. If a sufficiently good fit over the whole
interval requires an ugly polynomial, you can fit a prettier polynomial
to a subinterval, and transform arguments that are outside the
subinterval into it, and transform the results back; that is especially
easy for the log function. An example of that is in sys/dev/midisyn.c
in midisyn_vol2cB(). Note there that instead of using a polynomial to
approximate log(x) only, and then computing the rest of the function
around that, it just uses a custom polynomial fit to the desired
function f: x -> 400log(x/16256)<<22 so that the only additional step
needed is to shift off the 22 fraction bits at the very end.

Another approach to covering a larger interval is to fit simple
polynomials to pieces of the interval; ISTR a B-spline could offer
an elegant way of evaluating that, but my B-spline fu has degraded
since I took the class and I could probably not give a coherent
explanation off the top of my head right now. Probably overkill for
what you need anyway.

HTH,
-Chap