Page 1 of 1
Exp(x)
Posted: Sat Jan 03, 2009 3:51 am
by 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
Posted: Sat Jan 03, 2009 8:51 am
by 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.
Re: Exp(x)
Posted: Sat Jan 03, 2009 11:43 am
by Psychophanta
akj wrote: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
Posted: Sat Jan 03, 2009 11:49 am
by blueznl
Pshycho, Pow() and Log() are IMHO somewhat inaccurate. There's still an old bug floating around related to this.
Posted: Sat Jan 03, 2009 11:51 am
by Psychophanta
blueznl wrote: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.
Posted: Sat Jan 03, 2009 2:02 pm
by 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
Posted: Sat Jan 03, 2009 2:54 pm
by 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.
Posted: Sat Jan 03, 2009 6:54 pm
by akj
I'v decided to investigate further the problem with evaluating e^x with
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.
Posted: Sat Jan 03, 2009 7:22 pm
by 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.
Posted: Sat Jan 03, 2009 11:23 pm
by 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.
Posted: Sun Jan 04, 2009 9:09 am
by 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.
Posted: Sun Jan 04, 2009 2:27 pm
by Blue
Comtois wrote: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.
Posted: Sun Jan 04, 2009 3:50 pm
by Trond
Rescator wrote:
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.