Deterministic ln(x) approximation

[wolfboyft]

Hi. For my game I am looking for an approximation of the natural logarithm function that returns the exact same result on all IEEE 754 (still just all) platforms (with same rounding mode which seems garuanteed under "stock settings" LÖVE).

So that means it's got to rely on only +-*/ and sqrt. Oh, and abs, negate etc.
Here's what I have for exp(x) to show you what I mean:
https://www.desmos.com/calculator/w4blo6z6hg

Code: Select all

local function exp(x)
	local xint, xfract = modf(x)
	local exint = e ^ xint
	local exfract = 1 + xfract + (xfract ^ 2 / 2) + (xfract ^ 3 / 6) + (xfract ^ 4 / 24) -- for n = 0, 4 sum xfract^n/n!
	return exint * exfract -- e ^ (xint + xfract)
end

(Also, is e to the power of an integer (still represented through floats of course) also deterministic? I can't tell. It sounds like something that could well be implementation-specific. If it isn't, I don't mind using exponent-by-squaring or whatever.)

I have sine, cosine, etc etc already done. I'm just stuck at ln.

Thanks in advance.

[raidho36]

If you want exactly accurate math, you should use integers only. If you can confirm that addition and multiplication produce the same results on all target platforms, you can use lookup tables with linear interpolation - it will not be quite as accurate but it will produce the same results everywhere. You can be pretty sure just by doing some research online. But the only way to really confirm is to check every single combination of values you intend to use.

Generally, there's nothing deterministic about floating point math. Most if not all implementations do not honor IEEE in one way or another, particularly the hardware implementations. Most of the time basic operations will produce the same results across all platforms, but not even that is guaranteed. IIRC some nVidia GPUs will produce different results between individual operation executions, nevermind different OSes. Transcendental functions are particularly difficult to compute so virtually all of them do their own thing - that is the one thing you absolutely cannot count on.

[pgimeno]

I have sine, cosine, etc etc already done. I'm just stuck at ln.

This might help: https://en.wikipedia.org/wiki/Logarithm ... _algorithm

Edit: Quick and dirty implementation:

Code: Select all

local lntable = {
 math.log(1+2^-1),
 math.log(1+2^-2),
 math.log(1+2^-3),
 math.log(1+2^-4),
 math.log(1+2^-5),
 math.log(1+2^-6),
 math.log(1+2^-7),
 math.log(1+2^-8),
 math.log(1+2^-9),
 math.log(1+2^-10),
 math.log(1+2^-11),
 math.log(1+2^-12),
 math.log(1+2^-13),
 math.log(1+2^-14),
 math.log(1+2^-15),
}

function ln(x)
  local prod = 1
  local sum = 0
  for i = 1, 15 do
    local a = prod * (1 + 2^-i)
    if a < x then
      prod = a
      sum = sum + lntable[i]
    end
  end
  return sum
end

print(ln(1.3), math.log(1.3))

Of course, in the final implementation the table should use actual numbers, not be calculated using math.log; also the values (1 + 2^-i) themselves for i from 1 up to the desired precision (15 in the example) could be in another table. I've excluded the part where the argument is reduced to the range 1<x<2 as required by the algorithm (it actually works for up to ~2.38). The idea for the reduction is to use: ln(a*e^b) = ln(a) + b (for any b, including negatives), so you need to find the power of e that leaves a in the range 1 to 2. That's a logarithm in itself, but you only need an integer approximation; you can use e.g. binary search in an exponential table to find it.

Generally, there's nothing deterministic about floating point math. Most if not all implementations do not honor IEEE in one way or another, particularly the hardware implementations. Most of the time basic operations will produce the same results across all platforms, but not even that is guaranteed.

That sounds like FUD. All modern CPUs (including all platforms where LÖVE is supported) will use IEC559 semantics for the basic operations, excluding denormals, infinities and NaNs. Flush to Zero (FTZ) and Denormals are Zero (DAZ) is pretty common for speed reasons.

There's one major obstacle, though, and it's that the LÖVE framework doesn't guarantee a rounding mode. That's why I opened this thread: https://love2d.org/forums/viewtopic.php?f=3&t=87907 in which you spread the same FUD.

IIRC some nVidia GPUs will produce different results between individual operation executions, nevermind different OSes.

GPUs are not generally IEC559 compliant; they are happy to break IEC559 rules whenever there's room for running faster. Just don't run the calculations in the GPU if you want reproducible results.

[raidho36]

I've read plenty of blogs on the subject to know that marginally different results of floating point computations crop up virtually everywhere every now and then; one developer shared a story of spending a lot of time trying to hunt down a very subtle floating point math change between two versions of Apple software. If you need exact results every time, then it's hardly a FUD, rather it's a warning of trouble. That's why game engines use integers to this day (not having it in Lua and using it anyway in a game-specific framework is kind of an oversight; I guess you could try using integers in LuaJIT).

All modern CPUs...
Just don't run the calculations in the GPU if you want reproducible results.

All of them? Always producing exactly identical results to the least significant bit? What about older processors, are those to be listed as "not supported" by such program? What about mobile CPUs? You can't know this without either brute force checking which you can't do because of time requirements, or formally verifying the silicon (not the blueprint design) which you can't do because it's not available. You can only assume that it's true when the manufacturer specifies so and that there's no subtle bugs and whatnot. And that's not always the case, which will be critical for a program that relies on this. Nevermind software implementation of certain functions, those are subject to change without warning or even a footnote in the changelog.

People use GPUs for scientific calculations all the time, precision notwithstanding. It's a matter of programming it not to rely on exactly deterministic calculations, allowing some error bars.

[pgimeno]

What about older processors, are those to be listed as "not supported" by such program?

If LÖVE runs on a Pentium 3, you may have a point. LuaJIT works with JIT off in those cases, which compound with the slowness of the processor, means that probably there is no point for "such program" to support it.

[wolfboyft]

The idea for the reduction is to use: ln(a*e^b) = ln(a) + b (for any b, including negatives), so you need to find the power of e that leaves a in the range 1 to 2. That's a logarithm in itself, but you only need an integer approximation; you can use e.g. binary search in an exponential table to find it.

What range will b be in? What's the resolution of the table? Surely with the addition of a binary search this'll get really really slow? Couldn't I just find an implementation of what one of the CPUs I'm trying to "unify" uses and try to port it to Lua?

[wolfboyft]

Of course!

Binary searches don't have to be searching a table for a value.

[pgimeno]

What range will b be in? What's the resolution of the table?

Actually, forget that. I forgot we have math.frexp. This returns a mantissa between 0.5 and 1, and an exponent for a base of 2, i.e. it returns two numbers that are the decomposition of the input to the form m*2^x (m being the mantissa and x the exponent) with 0.5 <= m < 1. Actually if the input is negative or zero, so is m, so be careful to give an error in that case.

Multiplying the mantissa by 2 gives you a number in the 1..2 range, as required by the algorithm, so we use this instead:

original number = m*2^x = (m*2)*2^(x-1) (taking one of the twos from the power and "attaching" it to m)

Now, taking logarithms at both sides:

ln(original number) = ln(m*2^x) = ln((m*2)*2^(x-1))

... = ln(m*2) + ln(2^(x-1)) (because the logarithm of the product is the sum of the logarithms)

... = ln(m*2) + (x-1)*ln(2) (because the logarithm of a power is the exponent times the logarithm of the base)

and we can stop there. Have the value of ln(2) precalculated (edit: that's 0.6931471805599453), multiply that by the exponent minus 1, and add that to the result of the algorithm as applied to twice the mantissa. There's only a catch: if m is 0.5, m*2 is 1, and the ln algorithm is not prepared to work for that input, but we know that ln(1) = 0 so you can make that a special case and just return (x-1)*ln(2) if that happens.

[wolfboyft]

Sorry I'm late replying, but... thank you!!!
I got it done. Thank you so much :-)

Code: Select all

-- x raised to an integer is not deterministic
local function intPow(x, n) -- Exponentiation by squaring
	if n == 0 then
		return 1
	elseif n < 0 then
		x = 1 / x
		n = -n
	end
	local y = 1
	while n > 1 do
		if n % 2 == 0 then -- even
			n = n / 2
		else -- odd
			y = x * y
			n = (n - 1) / 2
		end
		x = x * x
	end
	return x * y
end

local function exp(x)
	local xint, xfract = modf(x)
	local exint = intPow(e, xint)
	local exfract = 1 + xfract + (xfract*xfract / 2) + (xfract*xfract*xfract / 6) + (xfract*xfract*xfract*xfract / 24) -- for n = 0, 4 sum xfract^n/n!
	return exint * exfract -- e ^ (xint + xfract)
end

local log
do
	local powerTable = { -- 1+2^-i
		1.5, 1.25, 1.125, 1.0625, 1.03125, 1.015625, 1.0078125, 1.00390625, 1.001953125, 1.0009765625, 1.00048828125, 1.000244140625, 1.0001220703125, 1.00006103515625, 1.000030517578125
	}
	local logTable = { -- log(1+2^-i)
		0.40546510810816438486, 0.22314355131420976486, 0.11778303565638345574, 0.06062462181643483994, 0.03077165866675368733, 0.01550418653596525448, 0.00778214044205494896, 0.00389864041565732289, 0.00195122013126174934, 0.00097608597305545892, 0.00048816207950135119, 0.00024411082752736271, 0.00012206286252567737, 0.00006103329368063853, 0.00003051711247318638
	}
	local ln2 = 0.69314718055994530942 -- log(2)
	function log(x)
		local xmant, xexp = frexp(x)
		if xmant == 0.5 then
			return ln2 * (xexp-1)
		end
		local arg = xmant * 2
		local prod = 1
		local sum = 0
		for i = 1, 15 do
			local prod2 = prod * powerTable[i]
			if prod2 < arg then
				prod = prod2
				sum = sum + logTable[i]
			end
		end
		return sum + ln2 * (xexp - 1)
	end
end

-- NOTE: Pretty big error magnification... :-/
local function pow(x, y)
	return exp(log(x)*y)
end

I wasn't sure that x^int was deterministic because the IEEE 754 standard doesn't mention that, only that x^numberInGeneral was recommended.

[wolfboyft]

Huh. Replacing that pow with

Code: Select all

local function pow(x, y)
	local yint, yfract = modf(y)
	local xyint = intPow(x, yint)
	local xyfract = exp(log(x)*yfract)
	return xyint * xyfract -- x ^ (yint + yfract)
end

leaves me with some serious error magnification too. I plotted abs(systempow(x, 1.5) - mypow(x, 1.5))

Now: https://love2d.org/imgmirrur/eMvQv86.png
Original: https://love2d.org/imgmirrur/N2VFBIV.png

It's ok tho, I probably won't even use that function much.
I'll do some research some other time.