Exp(x)

[akj]

Recently I had a need for a quick and accurate exponential function and as one is not built into PureBasic and as Pow() is rather slow and inaccurate, I decided to write my own in assembler. As other people may find a use for it, I decided to publish it here:

Code: Select all

Procedure.d Exp(x.d) ; AKJ  03-Jan-09
; Returns e^x where e = 2.718281828459...
;
; Theory:  Suppose e^x can be written as 2^N * 2^t
; where N is an integer and t is real (but small).
; Then 2^N would be easy to evaluate and 2^t can be found
; using the ASM instruction F2XM1 provided -1<=t<=+1 .
;
; The problem is to determine N and t.
; To do this take logs to base 2 of the equation e^x = 2^N*2^t :
; LHS = log2(e^x) = loge(e^x)/loge(2) = x/loge(2) = x*log2(e)
; RHS = log2(2^N*2^t) = N + t
; Thus N can be chosen as the nearest integer to x*log2(e)
; whence t will be x*log2(e) - N  (therefore between -1 .. +1)
;
; This routine determines N and t by the above method and
; then computes 2^N * 2^t to give the desired result.
;
; Actually a trick (the FSCALE instruction) is used to
; provide a quick and accurate way to do this last step.
; In ASM terminology, FSCALE computes ST0 * 2^Int(ST1)
; which is ideal if ST0 = 2^t and ST1 = N .
;
Protected control.w=$037F ; FPU control word default value
EnableASM
fldcw  word [p.v_control] ; Ensure rounding is 'to the nearest'
fld qword [p.v_x] ; st0 = x
fldl2e  ; st0 = log2(e)  st1 = x
fmulp   ; st0 = x*log2(e) = z (say)
fld st0 ; st0 = z  st1 = z
frndint ; st0 = N  st1 = z
fsub st1, st0 ; st0 = N  st1 = t
fxch ; st0 = t  st1 = N
f2xm1 ; st0 = 2^t-1  st1 = N
fld1 ; st0 = 1  st1 = 2^t-1  st2 = N
faddp ; st0 = 2^t  st1 = N  (ideal)
fscale ; st0 = e^x  st1 = N
ffree st1
DisableASM
ProcedureReturn ; st0
EndProcedure

; Examples:
Debug Exp(0)     ; 1.0
Debug Exp(1)     ; 2.7182818284590451
Debug Exp(#PI)   ; 23.140692632779267
Debug Exp(43.29) ; 6318414624313423900.0
Debug Exp(-2.6)  ; 0.074273578214333877

Edit: Added additional ASM instruction ffree st1
to avoid unwanted side effects when Exp() is called many times

[Blue]

Very nice. And amazingly effficient.

And many thanks for sharing the theory that led to the algorithm.
That's always intellectually refreshng, plus almost indispensable when it comes to Assembly Language.

[Psychophanta]

Recently I had a need for a quick and accurate exponential function and as one is not built into PureBasic and as Pow() is rather slow and inaccurate, I decided to write my own in assembler. As other people may find a use for it, I decided to publish it here:

Code: Select all

Procedure.d Exp(x.d) ; AKJ  03-Jan-09
; Returns e^x where e = 2.718281828459...
;
; Theory:  Suppose e^x can be written as 2^N * 2^t
; where N is an integer and t is real (but small).
; Then 2^N would be easy to evaluate and 2^t can be found
; using the ASM instruction F2XM1 provided -1<=t<=+1 .
;
; The problem is to determine N and t.
; To do this take logs to base 2 of the equation e^x = 2^N*2^t :
; LHS = log2(e^x) = loge(e^x)/loge(2) = x/loge(2) = x*log2(e)
; RHS = log2(2^N*2^t) = N + t
; Thus N can be chosen as the nearest integer to x*log2(e)
; whence t will be x*log2(e) - N  (therefore between -1 .. +1)
;
; This routine determines N and t by the above method and
; then computes 2^N * 2^t to give the desired result.
;
; Actually a trick (the FSCALE instruction) is used to
; provide a quick and accurate way to do this last step.
; In ASM terminology, FSCALE computes ST0 * 2^Int(ST1)
; which is ideal if ST0 = 2^t and ST1 = N .
;
Protected control.w=$037F ; FPU control word default value
EnableASM
fldcw  word [p.v_control] ; Ensure rounding is 'to the nearest'
fld qword [p.v_x] ; st0 = x
fldl2e  ; st0 = log2(e)  st1 = x
fmulp   ; st0 = x*log2(e) = z (say)
fld st0 ; st0 = z  st1 = z
frndint ; st0 = N  st1 = z
fsub st1, st0 ; st0 = N  st1 = t
fxch ; st0 = t  st1 = N
f2xm1 ; st0 = 2^t-1  st1 = N
fld1 ; st0 = 1  st1 = 2^t-1  st2 = N
faddp ; st0 = 2^t  st1 = N  (ideal)
fscale ; st0 = e^x  st1 = N
DisableASM
ProcedureReturn ; st0
EndProcedure

; Examples:
Debug Exp(0)     ; 1.0
Debug Exp(1)     ; 2.7182818284590451
Debug Exp(#PI)   ; 23.140692632779267
Debug Exp(43.29) ; 6318414624313423900.0
Debug Exp(-2.6)  ; 0.074273578214333877

Why to do it like that if you can do it like this?:

Code: Select all

Macro exp(x)
  (Pow(2,x#/Log(2)))
EndMacro

Do a search before allows you to avoid 'reinventing'...
http://www.purebasic.fr/english/viewtopic.php?t=29573

[blueznl]

Pshycho, Pow() and Log() are IMHO somewhat inaccurate. There's still an old bug floating around related to this.

[Psychophanta]

Pshycho, Pow() and Log() are IMHO somewhat inaccurate. There's still an old bug floating around related to this.

Well, sorry, my way is a shorter way, and may be faster, but may be inaccurater.
Just wanted to say that this nice trick from akj might be appended to:
http://www.purebasic.fr/english/viewtopic.php?t=29573
I think it should avoid more than one thread talking about the same.

[akj]

@ Psychophanta:

I had already seen the topic defining the macro using Pow() and Log() but had rejected it, mainly because of it's inaccuracy. To see what I mean, please run this small program:

Code: Select all

Define x.d, y.d
x = 0.01 ; Arbitrarily chosen
y = Pow(2,x/Log(2))
Debug y
Debug Pow(2,x/Log(2))

On my PC the two Debug statements give different results of

Code: Select all

[Debug] 1.0100501670564135
[Debug] 1.0100501670841679

Because of this, I don't trust the Pow() method

[Psychophanta]

Yes, that's right, but you could append your tip there, as an improvement.

And anyway, take in account that for future PB versions functions like Pow() and Log(), there will probably be updated to be accurate, and then my macro will eventually be as accurate as your tip.

[akj]

I'v decided to investigate further the problem with evaluating e^x with

Code: Select all

Pow(2,x/Log(2))

I have written replacement versions of Pow() and Log() which I have called Power() and Ln() and then tried the four combinations:

Code: Select all

y = Pow(2,x/Ln(2)) : Debug y
y = Pow(2,x/Log(2)) : Debug y
y = Power(2,x/Ln(2)) : Debug y
y = Power(2,x/Log(2)) : Debug y

and as you will see (by running the program below only the combination using both Pow() and Log() gives a wrong result.

Here are the new routines Power() and Ln() and the test

Code: Select all

Procedure.d Power(x.d, y.d)
; Returns x^y
; Uses the method: y^x = 2^(y*log2(x)) = 2^z (say)
; But the instruction F2XM1 will compute 2^z only if -1<=z<=1
; Therefore let 2^z = 2^(N+t) = 2^N * 2^t
; The remaining logic is the same as for the Exp() function
;
; The returned result is always positive except for the
; special case of x<0 and y an odd integer  e.g. (-3)^5
Protected control.w=$037F ; FPU control word default value
Protected i, negate.b = #False ; For the special case mentioned above
If x=0.0: ProcedureReturn 0.0: EndIf
i = Int(y)
If x<0.0 And y-i=0.0 And i%2=1: negate = #True: EndIf
EnableASM
fldcw  word [p.v_control] ; Ensure extended double precision
fld qword [p.v_y] ; st0 = y
fld qword [p.v_x] ; st0 = x  st1 = y
fabs ; Kludge: ensure that st0 is >0
fyl2x ; st0 = y*log2(x) = z
fld st0 ; st0 = z  st1 = z
frndint ; st0 = N  st1 = z
fsub st1, st0 ; st0 = N  st1 = t
fxch    ; st0 = t  st1 = N
f2xm1   ; st0 = 2^t-1  st1 = N
fld1    ; st0 = 1  st1 = 2^t-1  st2 = N
faddp   ; st0 = 2^t  st1 = N
fscale  ; st0 = x^y  st1 = N
ffree st1 ; st0 = x^y
; Special case: If x<0 and y is an odd integer, negate the result
If negate
  fchs
EndIf
DisableASM
ProcedureReturn ; st0
EndProcedure

Procedure.d Ln(x.d)
; Returns the natural logarithm of x
; Uses the method: loge(x) = log2(x)/log2(e) = log2(x)*loge(2)
EnableASM
fldln2 ; st0 = loge(2)
fld qword [p.v_x] ; st0 = x  st1 = loge(2)
fyl2x ; st0 = log2(x)*loge(2)
DisableASM
ProcedureReturn ; st0
EndProcedure

; Examples using new Power(), Ln() and existing Pow(), Log()
Define x.d, y.d 
x = 0.01 ; Arbitrarily chosen
y = Pow(2,x/Ln(2)) : Debug y
y = Pow(2,x/Log(2)) : Debug y ; <-- Wrong result
y = Power(2,x/Ln(2)) : Debug y
y = Power(2,x/Log(2)) : Debug y

P.S. Please note that I have made an amendment to Exp() in my first post.

[Comtois]

well , i just use this one , work fine enough for me

Code: Select all

#E=2.718281828459045235
Procedure.d Exp(value.d)
  ProcedureReturn Pow(#E, value)
EndProcedure

Thank you for your ASM procedure, it could help.

[Rescator]

Do not rely on Debug behaving the same as normal code if expressions are used directly with it.

for example Debug Pow(2,x/Log(2))
you'd never use that in actual code.
Instead you'd use:
y = Pow(2,x/Log(2))
Debug y

Took me a while to learn to do this properly in my code,
especially when it comes to floats. (I kept wondering why Debug never matched the expected result)
Not entirely sure why it's like this but it's best to define the variables, put the result in a defined variable and THEN show it using Debug.
Remember, Debug is just a quick wy to test/spit out a result.
It's not like one uses: x=Debug 1+2 hehe.

[superadnim]

I've been doing something similar to what Psychophanta suggests, but compared to the original posted function, it's 4.2 times slower. Even if it's a macro, it's still doing pow and log...

Comtois' suggestion is also slower than the original function posted by akj (1.8 times slower).

But, unless this is your application's bottleneck I wouldn't go for the assembly version. Imagine if you were to port this code to another platform, or if someone else who doesn't know assembly had to do it... it would prove difficult. However, if the PB math lib were to support it, I'd be happy to use it.

What Rescator says is also true, the debug keyword is just a way of spitting out to the debug console and nothing else, it should be used carefully. Plus, with proper debugging techniques being already present in PB 4.x you'd rarely use it, unless it's for some quick test... On a real debugging scenario, you'd have CallDebugger and you'd be stepping the code while looking at the content of the variables/structures. Sometimes even the registers, specially EAX. It also helps to log down to a file, but this is useful only when you've deployed and you need to track bugs.

[Blue]

well , i just use this one , work fine enough for me [...]

Well, that's hardly the point, isn't it?
akj's aim is precision and high accuracy, and here he provides both brilliantly.
Speed is always nice too (programmers rarely complain about too much of it ), but here, it is incidental to his smart algorithm and the language he chose to implement it.

The reservations expressed by Rescator and superadnim, however, are right on and worth keeping in mind.

[Trond]

Took me a while to learn to do this properly in my code,
especially when it comes to floats. (I kept wondering why Debug never matched the expected result)
Not entirely sure why it's like this but it's best to define the variables, put the result in a defined variable and THEN show it using Debug.

Debug shows the result correctly, but it uses doubles and quads instead of float and long. "Debug expression" should be the same as "n.d = expression : debug n". If it's not, then I think there's a bug in PB.