[Wurst] Logarithm

[Menag]

[WurstScript] Logarithm 1.2

Logarithm is a library for WurstScript which provides an easy way to calculate different logarithms.

There are no requirements except for the obligatory WurstScript StandardLib.

Code:

/**
 * A package for calculating the logarithm to different bases.
 */
package Logarithm

/**
 * Eulers number e.
 */
public constant real E = 2.718282

constant real INV_LOG_2_E = 0.6931472
constant real INV_LOG_2_10 = 0.3010300

constant real inf = Pow(2.0, 128.0)

/**
 * Calculates the logarithm to base 2 of the absolute value of the argument.
 * This is the binary logarithm or logarithm dualis.
 *
 * real a - the argument of the logarithm
 *
 * returns the result of log2|a|
 *
 * error - a = 0.0
 */
public function log2(real a) returns real
	real x = a
	real sign = 1
	real res = 0.0
	real p = 0.5
	
	if x == 0.0
		error("Cannot calculate the logarithm of 0.0!")
		return -inf
	
	if x < 0.0
		x = -x

	if x < 1
		x = 1.0 / x
		sign = -1.0
		
	while x >= 2
		res += 1
		x *= 0.5
	
	for int i = 1 to 23
		x *= x
		if x >= 2.
			x *= 0.5
			res += p
		p *= 0.5
	
	return sign * res

/**
 * This is just a wrapper for log2(a).
 *
 * real a - the argument of the logarithm
 *
 * returns the result of lb|a|
 *
 * error - a = 0
 */
public function lb(real a) returns real
	return log2(a)

/**
 * This is just a wrapper for log2(a).
 *
 * real a - the argument of the logarithm
 *
 * returns the result of ld|a|
 *
 * error - a = 0
 */
public function ld(real a) returns real
	return log2(a)

/**
 * Calculates the logarithm to base 10.
 *
 * real a - the argument of the logarithm
 *
 * returns the result of lg|a|
 *
 * error - a = 0
 */
public function lg(real a) returns real
	return log2(a) * INV_LOG_2_10

/**
 * Calculates the logarithm to base e.
 * This is the natural logarithm.
 *
 * real a - the argument of the logarithm
 *
 * returns the result of ln|a|
 *
 * error - a = 0
 */
public function ln(real a) returns real
	return log2(a) * INV_LOG_2_E

/**
 * Calculates the logarithm to an arbitrary base.
 *
 * real a - the argument of the logarithm
 * real b - the base of the logarithm
 *
 * returns the result of log_{|b|}(|a|)
 *
 * error - b = 0
 *		 - b = 1
 *		 - b = -1
 * 		 - a = 0
 */
public function log(real a, real b) returns real
	if (b == 1.0) or (b == -1.0) or (b == 0.0)
		error(b.toString() + " is not a viable base!")
		if (b == 0)
			return 0.0
		else if (a.abs() > 1)
			return inf
		else if (a.abs() < 1)
			return -inf
		return 1.0
	return log2(a) / log2(b)

Tests:

package LogarithmTest
import NoWurst
import Logarithm
import Wurstunit

function generalTest()
	real delta = 0.00001
	real r
		
	r = log(-5, 2)
	r.assertEquals(2.321928094887, delta)
		
	r = log(1, 2)
	r.assertEquals(0.0, delta)
	
@test function log2Test()
	real delta = 0.00001
	real r
	
	r = log2(2.0)
	r.assertEquals(1.0, delta)
	
	r = log2(64.0)
	r.assertEquals(6.0, delta)
	
	r = log2(1.2345)
	r.assertEquals(0.303926836480, delta)
	
	r = log2(10000.0)
	r.assertEquals(13.28771237955, delta)
	
@test function lgTest()
	real delta = 0.00001
	real r
	
	r = lg(10.0)
	r.assertEquals(1.0)
	
	r = lg(10000.0)
	r.assertEquals(4.0, delta)
	
	r = lg(1.2345)
	r.assertEquals(0.091491094267, delta)
	
@test function lnTest()
	real delta = 0.00001
	real r
		
	r = ln(E)
	r.assertEquals(1.0, delta)
	
	r = ln(148.413159)
	r.assertEquals(5.0, delta)
	
@test function logTest()
	real delta = 0.00001
	real r
		
	r = log(4.0, 4.0)
	r.assertEquals(1.0, delta)
	
	r = log(10000.0, 2.0)
	r.assertEquals(13.28771237955, delta)

GitHub:
Logarithm

Changelog:
1.0 - Initial release
1.1 - Now returns the logarithm of the absolute value of the argument. Crashs the thread if 0 is passed as an argument (debugmode only!).
1.2 - Changed the behavior for invalid input. Still crashs the thread in debugmode. Added inf constant (thanks to looking_for_help).

Credits:
- The awesome penguinking

 

[Nestharus]

Maybe specifically do Ln too?

http://www.hiveworkshop.com/forums/jass-resources-412/snippet-natural-logarithm-108059/

 

[Menag]

Maybe specifically do Ln too?

http://www.hiveworkshop.com/forums/jass-resources-412/snippet-natural-logarithm-108059/

I know the methods how to implement Ln, but I don't think it's worth it since the current way is only one more multiplication then the log2...
So I don't think i will specifically do ln as long as there's no superior way to calculate it.

 

[Dr Super Good]

if x <= 0.0
return 0.0

Real mathematics tsk tsk.

At least invert x then so you get the result of RE(log2(a)). Yes it does have a real part to it so there is no need to return 0.

 

[Menag]

Real mathematics tsk tsk.

At least invert x then so you get the result of RE(log2(a)). Yes it does have a real part to it so there is no need to return 0.

Gimme NaN and inf and I'll give it the correct values. I won't give it some arbitrary value. Yes, the logarithm has a complex continuation, but that should be the task of a complex numbers package.

 

[Dr Super Good]

I won't give it some arbitrary value.

Except you currently are. 0 is pretty arbitrary. Worse is that 0 is actually the log of 1 so you will be giving it a completely wrong answer.

Thus I would strongly advise making it
RE(log2(a)) when a < 0 (ignore the imaginary part)
Thread crash when a == 0 (how JASS handles division by 0)

 

[Menag]

Thus I would strongly advise making it
RE(log2(a)) when a < 0 (ignore the imaginary part)
Thread crash when a == 0 (how JASS handles division by 0)

I'm thinking about this stuff again right now.
About the thread crash, I was thinking about that too. Wurst has some built in error function for that which crashs the thread, prints a stacktrace and an error message. But does this only in debugmode. Therefore I'd need a non-debugmode solution anyway. Again threadcrash?
For negative values. Then I'd also have to allow negative bases and I'd just say that it calculates log_{|b|}(|a|).

 

[looking_for_help]

Gimme NaN and inf and I'll give it the correct values.

I can only give you Inf Pow(2, 128) from Maths. Maybe return that in non-debug mode?

 

[Dr Super Good]

Since reals are apparently floats you should be able to get both NaN and Inf by abusing that type casting exploit thing that allows you to bitwise convert a real to an integer and back by setting an integer with the appropriate bit pattern and converting it into a real. That is of course if the type of float WC3 uses has both NaN and Inf (this behaviour is platform specific and a huge problem with using floats for anything cross platform).

 

[Menag]

Since reals are apparently floats you should be able to get both NaN and Inf by abusing that type casting exploit thing that allows you to bitwise convert a real to an integer and back by setting an integer with the appropriate bit pattern and converting it into a real. That is of course if the type of float WC3 uses has both NaN and Inf (this behaviour is platform specific and a huge problem with using floats for anything cross platform).

No. Read this.

 

[Dr Super Good]

Like I said, NaN and such values are platform specific. Not all float representations have them.

 

[PurgeandFire]

Useful snippet. Approved.

 

[Waffle]

a)this is enlightening:
http://en.wikipedia.org/wiki/IEEE_floating_point

b) war3 i think it would be safe to assume binary format (and probably.. 32bit).

Amongst the properties you will find that Inf and -Inf have a single unique value each.(highest possible exponent, all zero mantissa)

Yet there are thousands of different representations for NaN(maximum exponent, non-zero mantissa).

While the Warcraft 3's JASS2 engine apparently does not treat these values differently.. most modern CPU-s (least ones you'd wanna play games on) do have FPU-s that comply with the IEEE 32bit binary float standard.

Thus if you would test for NaN-s manually you should still be able to detect them fairly reliably..

(a whole separate issue would be if you'd want to risk using NaN detection in any logic that effects actual gameplay.. rare setups tend to make themselves visible when you got thousands of different systems running the code.)

EDIT: http://en.wikipedia.org/wiki/NaN
if you actually want to deal with NaN's this is at minimum.. educational.

 

[Menag]

Ehm...? We know, but wc3 doesn't give us the possibility to use NaN and inf. Just test it.

 

[MapperMalte]

Why not use Taylor Polynomials to calculate ln? They converge pretty quickly

 

[Menag]

Why not use Taylor Polynomials to calculate ln? They converge pretty quickly

not. With ln(1+h) and 0 <= h <= 1 and 0 < theta < 1 we have
|R_n(h)| = |1/(n+1) * h^(n+1)/(1+theta*h)^(n+1)|
<= h^(n+1)/(n+1)
We now want this to be smaller or equal to, let's say 10^-5.
With h = 1 we have the worst case and obtain
n+1 = 10^5. Huh. Pretty quickly. To make this more obvious, let's calculate it a bit
ln(2) - (1 - 1/2 + 1/3 - 1/4 + 1/5 - 1/6) = 0.0764805
doing this till 50 gives us a difference of
0.00990002
And now. The real fast aspect: n = 10000 yields:
0.0000499975.
That's pretty fast.

 

[Dark-Zalor]

/**
 * A package for calculating the logarithm to different bases.
 */
library Logarithm

/**
 * Eulers number e.
 */
globals
  private constant real E = 2.718282

  private constant real INV_LOG_2_E = 0.6931472
  private constant real INV_LOG_2_10 = 0.3010300

  private constant real inf = Pow(2.0, 128.0)
endglobals

/**
 * Calculates the logarithm to base 2 of the absolute value of the argument.
 * This is the binary logarithm or logarithm dualis.
 *
 * real a - the argument of the logarithm
 *
 * returns the result of log2|a|
 *
 * error - a = 0.0
 */
private function log2 takes real a returns real
  local real x = a
  local real sign = 1
  local real res = 0.0
  local real p = 0.5
  local integer I = 1
 
  if x == 0.0 then
    call BJDebugMsg("Cannot calculate the logarithm of 0.0!")
    return -inf
  endif

  if x < 0.0 then
    set x = -x
  endif

  if x < 1 then
    set x = 1.0 / x
    set sign = -1.0
  endif

  loop
    exitwhen x >= 2
    set res = res + 1
    set x = x * 0.5
  endloop

  // for int i = 1 to 23
  loop
    exitwhen I >= 23
    set x = x * x
    if x >= 2. then
      set x = x * 0.5
      set res = res + p
    endif
    set p = p * 0.5
    set I = 1 + 1
  endloop

  return sign * res
endfunction

/**
 * This is just a wrapper for log2(a).
 *
 * real a - the argument of the logarithm
 *
 * returns the result of lb|a|
 *
 * error - a = 0
 */
public function lb takes real a returns real
  return log2(a)
endfunction

/**
 * This is just a wrapper for log2(a).
 *
 * real a - the argument of the logarithm
 *
 * returns the result of ld|a|
 *
 * error - a = 0
 */
public function ld takes real a returns real
  return log2(a)
endfunction

/**
 * Calculates the logarithm to base 10.
 *
 * real a - the argument of the logarithm
 *
 * returns the result of lg|a|
 *
 * error - a = 0
 */
public function lg takes real a returns real
  return log2(a) * INV_LOG_2_10
endfunction

/**
 * Calculates the logarithm to base e.
 * This is the natural logarithm.
 *
 * real a - the argument of the logarithm
 *
 * returns the result of ln|a|
 *
 * error - a = 0
 */
public function ln takes real a returns real
  return log2(a) * INV_LOG_2_E
endfunction

endlibrary

As vJass