> > Newton-Raphson seems of no use since I can't use the inverse function. > > It may be of more use than you think. For y = ln(x), the recursion > > y = y - (1 - x*exp(-y)) > > will work quickly since you can (probably) work with normalized numbers > for which the series approx for exp(-y) won't need many terms to get the > precision you want. With decent estimates this will take only a few > iterations.
That is the approach I'm leaning toward. Broadly:
1. Get a good exponentiation algorithm going. Probably split integer and fraction parts, then apply a squared Taylor expansion on the fractional part. For 24 bits of accuracy, probably end up with about 7 multiplies. (I read some stuff about binary splitting, binary reduction and other esoteric algorithms, but I fail to see how they would reduce the number of multiplies -- anyone care to explain?)
2. Given the log estimate, apply either Newton Raphson or a Halley iteration using the exponentiation above:
y = y - 2 * (exp(y)-x)/(exp(y)+x)
The Halley iteration is supposed to be cubically convergent i.e. 3x the number of bits of accuracy of the original, which according to Altivec was 5 bits. (Using NR would get 10 bits, then 20, then 40 i.e. 3 iteration; using H would get 15 bits, then 45 bits i.e. 2 iterations.) I'd estimate something like 18 multiplies.
If 2. is too slow, might have to consider some other plan...