gbadev.org forum archive

Does anybody know of a way to raise a fixed point number to a fractional exponent (without resorting to floating point)? If anyone is feeling crazy, I'm looking for a function with a prototype that looks something like this:

Code:

FIXED exp_fixed(FIXED p_base, FIXED p_power, FIXED p_frac);

p_base is the base in fixed point format.
p_power is the power to raise p_base to. May have a fractional component.
p_frac is the number of bits used for the fractional component of the base.
The fractional component of the power is defined by a macro FIXED_POWER_FRAC.

-Brendan

That's a toughie. You could use the relation:

Code:

y = a^b
ln y = b * ln a
exp(ln y) = exp(b * ln a)
y = exp(b * ln a)

Now if it's practical for you to make a LUT for ln that covers the entire domain of "a", and can also make a LUT for e^x that covers the entire domain of (b * ln a), then you'll have a really quick good approximation of your power function.

If that's not practical given your domains, then you'll probably have to calculate a series to a point of sufficient precision. Time to crack open yer old math books... :P

Dan.

That about covers it - no built in functions raise items to fractional exponents - out of curiousity, why do you need this? The only time I run into this are when I'm either doing calculus or some really crazy, in-efficient geometric transformation that should be using vector math.

I'd like a function that creates a custom gamma correction table- GBA, SP, DS, GBP/TV, VGA, and Custom. It looks like I'll have to precompute gamma correction tables and offer global brightness adjustment instead of custom gamma. (Or just brightness adjustment if I can't figure out the gamma for all the screens.)

Other than that, a fixed point fractional exponent function would be really cool. =) Run time experience to next functions based on fractional exponential curves. Sweet. =)

-Brendan

Here is the on I used in my emulators (Specificaly PCEAdvance):

Code:

;----------------------------------------------------------------------------
paletteinit; r0-r3 modified.
;called by ui.c: void map_palette(char gammavalue)
;----------------------------------------------------------------------------

stmfd sp!,{r4-r7,lr}
mov r7,#0xE0
ldr r6,=MAPPED_RGB
ldr r1,gammavalue ;gamma value = 0 -> 4
mov r4,#1024    ;pce rgb, r1=R, r2=G, r3=B
sub r4,r4,#2
nomap             ;map 0000000gggrrrbbb -> 0bbbbbgggggrrrrr
and r0,r7,r4,lsl#1 ;Red ready
bl gprefix
mov r5,r0

and r0,r7,r4,lsr#2 ;Green ready
bl gprefix
orr r5,r5,r0,lsl#5

and r0,r7,r4,lsl#4 ;Blue ready
bl gprefix
orr r5,r5,r0,lsl#10

strh r5,[r6,r4]
subs r4,r4,#2
bpl nomap

ldmfd sp!,{r4-r7,lr}
bx lr

;----------------------------------------------------------------------------
gprefix
orr r0,r0,r0,lsr#3
orr r0,r0,r0,lsr#6
;----------------------------------------------------------------------------
gammaconvert; takes value in r0(0-0xFF), gamma in r1(0-4),returns new value in r0=0x1F
;----------------------------------------------------------------------------
rsb r2,r0,#0x100
mul r3,r2,r2
rsbs r2,r3,#0x10000
; subne r2,r2,#0x2a8       ;Tweak for Gamma #4...
rsb r3,r1,#4
orr r0,r0,r0,lsl#8
mul r2,r1,r2
mla r0,r3,r0,r2
mov r0,r0,lsr#13

bx lr

C-ish version:

Code:

int GammaCorrection(int Luma, int Gamma){
int i = 0x100 - Luma;
i = i * i;
i = 0x10000 - i;
Luma = Luma | (Luma << 8);
i = i * Gamma;
Luma = (Luma * (4-Gamma)) + i;
Luma = Luma >> 13;
return(Luma);
}

Using brightness only lowers the range of the palette and you lose contrast.
_________________
I probably suck, my not is a programmer.

sgeos wrote:

Does anybody know of a way to raise a fixed point number to a fractional exponent (without resorting to floating point)?

Yes. I've actually coded up such a thing, and am using it in the game I'm working on. The code is not simple, and ensuring that it remains numerically stable is a challenge.

I can't give you the code directly; I'm sorry; there'd be a little problem with copyrights ^^;
However... I can tell you the basic outline of the code, and where you have to look for the rest of the answers you'll need.

My power function works like this (in pseudocode):

Code:

IntPow(x, y) {
if x < 0 then return 0
z = IntLog2(x) * y
return IntPow2(z)
}

IntLog2 computes the base-2 logarithm of an integer x, returning an 8.24 fixed-point number. IntPow2 computes 2 raised to the zth power, where z is an 8.24 fixed-point number.

My implementation raises an integer x to a 16.16 fixed-point value y, and returns a result rounded to the nearest integer; if you need fractional bits in x or in the return value, that's not too hard to implement, but you'd better be up to the task of working out the math ;) There are some nasty issues in IntPow relating to retaining precision; in my version, I implemented it as assembly to ensure that I could keep all 64 bits from the multiplication.

Now, as for the uglier details of implementing IntLog2 and IntPow2...

IntLog2 is divided into two halves. First, you compute the integer portion of the logarithm using repeated shifting (i.e., find the position of the highest 1 bit, and that's the integer portion). Then, use the algorithm described in Knuth (TAoCP, of course) 1.2.2.25 to take the other bits and compute the fractional portion of the logarithm (the Knuth/Briggs algorithm). You'll need to pre-compute the additional logarithm table described in the book as well, but that can be easily done in a few lines of code (I just hacked together a Perl script to generate the table).

Code:

IntLog2(x) {
if x <= 0 then return 0
result = [position of highest 1 bit in x]
x <<= 31 - result
y = [Knuth/Briggs algorithm 1.2.2.25 performed on x]
return (result << 24) + (y >> 7)
}

IntPow2 is divided into three parts. First, it breaks apart the fixed-point number into its integer and fractional parts. Second, it uses algorithm E (Knuth 1.2.2.28, originally discovered by Feynman; note that the book has a typo in it --- it should be x = 1 - x, not x = 1 - e, in step E1) to compute 2^frac(x). Third, it corrects for rounding (not shown below), and at the same time uses a clever shift to move the decimal point in the result from algorithm E so that the fractional part happens to become the integer part, and returns it.

Code:

IntPow2(x) {
if x <= 0 then return 0
frac = (x >> 24), x &= 0xFFFFFF
if x > 0 {
x <<= 7
y = [Knuth/Feynman algorithm 1.2.2.28 performed on x]
frac++
}
return y >> (31 - frac)
}

The math behind both algorithms is really far too complex to go into full detail here; STUDY the Art of Computer Programming, Volume I if you are really planning to do this.

In reality, I don't use IntPow() directly for computing gamma, because its precision is in the wrong place. I use a variant, IntFracPow(), that raises a 1.16 fixed-point base to a 16.16 fixed-point power. This can be performed with a little extra arithmetic, and amounts to a total of four additional instructions compared to IntPow() in my version.

The code is all semi-ugly, as most deeply mathematical code is, but it runs blazingly fast (especially since I have it in IWRAM). I can perform several thousand IntPow() calls per frame with plenty of time left over for music and some minimal input-handling, which is plenty for computing gamma tables on the fly --- assuming you cache the results for the actual game after they get computed, that is ;)
_________________
Do you suppose if I put a signature here, anyone would read it? No? I didn't think so either.

gbadev.org forum archive

C/C++ > Seeking Fixed Point Fractional Power Algorithm

#33136 - sgeos - Fri Dec 31, 2004 3:02 am

#33159 - poslundc - Fri Dec 31, 2004 6:45 am

#33459 - Miked0801 - Mon Jan 03, 2005 8:44 pm

#33478 - sgeos - Tue Jan 04, 2005 3:23 am

#33484 - FluBBa - Tue Jan 04, 2005 10:47 am

#33936 - phantom-inker - Tue Jan 11, 2005 9:05 am