Here's one more example of the perils of floating-point arithmetic, due to William Kahan, designed to show that roundoff error can still bite us even when the usual suspects (cancellation, tiny divisor) do not come into play.
Let t=:3 : 'if. y=0 do. 1 else. ((^y)-1)%y end.' q=:3 : ' (|y-%:1+*:y)-%(y+%:1+*:y)' g=:3 : 't *:q y' Mathematically, t(y) is the sum of the series (y^n-1)%!n for n>:1, and in particular, t(0)=1, while q(y)=0 for all y. So if y>0, g(y)=1. However, +/ g >:i.1e4 0 Every value has been miscalculated as 0, not 1. You can get some idea of the problem by looking at require 'plot' plot q >:i.1e4 plot (];(t"0)) steps 0 1e_15 1000 Although q is almost never exactly zero, it is close. The notch in t up to 1e_16 is the problem: for values of y close to zero, ^y rounds to 1, and t(y) is 0. You probably do not expect the graph of ^ to look like this: plot (];^) steps 0 1e_15 1000 Best wishes, John ---------------------------------------------------------------------- For information about J forums see http://www.jsoftware.com/forums.htm
