Fractional Powers of 2

TomS

I need to calculate fractional powers of 2 and I can't find or derive an algorithm. I'm try to solve

2^(x/y) where both x & y are positive integers. Unfortunately, x can be 18 bits in length. y is limited to 16 bits length.

I've tried using the built-in log tables but they leave something to be desired, ie low resolution.

I'm trying to do this in asm.


Tom

rjo_

Tom...

It would be better to know exactly what you are trying to do...to give you a perfect answer.

you have to know how much accuracy you want, but for this example lets say a single a tenth in your fractional exponent is good enough.

if you resolve your exponent, x/y, into a decimal number giving the form, 2^(x.y), then

2^(x.y):=(2^x) * yconstant

you then need a lookup table corresponding to yconstant[noparse][[/noparse]1,2,... 9] array contains values 1.07177, 1.14869...1.8660

notice 2^(.1)=1.07177
2^(.2)= 1.4869
.
.
.
2^(.9)=1.8660

and that the size of your lookup table is determined by the decimal places in y

Rich

rjo_

Of course this has nothing to do with the fact that your number is stored in binary and there are tricks beyond my means for doing the same things at the bit level... which would probably serve your needs better.

deSilva

Some more suggestions:
- You could develop x/y as binary fraction, so you need not copnsider half of the terms being zero
- When looking up the anti log table, do an interpolation as you (?) used to do whith printed log tables.

rjo_

deSilva,

He was having trouble using the log table...and I've never used it either[noparse]:)[/noparse] Seems like a waste of space to me[noparse]:)[/noparse]

Rich

deSilva

It's a key to advanced algorithmes... However I have not yet seen any agorithmes on the Prop

rjo_

deSilva

We could have a good philosophical discussion about this, and it is an interesting topic... but I am too good at going "off thread" and I don't want to do that to Tom.

Let me ask... is there anyone out there that can show us how to do this at the bit level a lot better than the approach that either deSilva or I have described?


Rich

TomS

rjo_·...· I thought of doing something similar but discarded the idea.· I was going to use powers of 2 (binary) to express your y, ie 2^(1/2), 2^(1/4) ... etc and then subtract the fraction and loop . The problem is I still end up with real numbers instead of integers. I could always multiply by some constant to get rid of some of the digits to the right of the decimal point but my·result can be as high as 2^25 leaving only 6 bits magnitude in my constant.· This isn't enough to give me the final result to 12 bits.

deSilva... Thanks, but no thanks to log tables.· When I tried using them my result was off by a factor of 2.· That's with interpolation.·

My $10 calculator can do it so I know somewhere, somebody has the algorithm I need.· I'm hoping Tracy Allen·replies as he is the expert when it comes to math.



Tom

deSilva

.. and this may also enlighten the ambitious www.jjj.de/fxt/fxtbook.pdf Attention: it's 4 MB! On J

TomS

...· The ambitious also may not get much sleep in the next week.



Thanks

rjo_

deSilva,

Great link. Thanks.

Someone needs to send Jorg a Propeller[noparse]:)[/noparse]

Rich

Tracy Allen

Tom, I saw your post late last night, and my eyes were already too blurry to think about it. I am not really such an expert when it comes to math. Just persistent sometimes. The book by Jorg is quite remarkable. Especially for this problem, chapter 20, "Elementary functions with limited resources". Section 20.1.2 is "computing b^x". I wouldn't give up on the Prop ROM table (with interpolation) though, depending on how many bits are really significant.

My first reaction was to wonder about the range of the input x and y and how to scale the output. Your initial post restricted x up to 2^18 and y up to 2^16, so the result could be anywhere from 2^(2^2)=16 up to something astronomical, 2^(2^18). In a post this morning you qualified it to a result that could be as high as 2^25. If there have to be 12 bits to the right of the binary radix point, that will entail stretching the results say into one long for the integer part plus one word for the fraction. Is that necessary, or is it a question of the total number of significant bits, integer+fraction? Does it have to be fast as well as accurate?!

▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔
Tracy Allen
www.emesystems.com

TomS

Tracy,

Sorry, I meant y was 1 word. It's actual value is constant {9031 = 3 * 10K * LOG[noparse][[/noparse]10](2)}.··So what I'm trying to compute is 2^(x/9031) with x·at a high of 250_000 and a low x of 30_000. What I then tried to do is use the divide routine·in Propeller Guts to give me the integer portion and the remainder. Raising 2 to an integer exponent is easy, just shift left. Then I would need to compute 2^(remainder/9031) and multiply times the integer result (2^(i+j) = 2^i * 2^j).· The result of remainder/9031 will always be less than 1. That's where I'm at right now.· I really don't want to use spin to do this although if I must I will.· It may not be easy to do in assembly but it must be possible.



Tom

rjo_

Tom,

If this is a time sensitive encryption problem... you can still use a look up table that you can generate externally and hold on a memory card of some sort... and then just parallelize the heck out of it until you get your time constraint satisfied.

Rich

Graham Stabler

2^(x/y) is the same as the yth root of 2^x so perhaps you can find assembly for each of these processes.

e.g. 2^(1/2) = sqrt(2^1) = sqrt(2)

Graham

TomS

I have decided to try again the idea I earlier discarded.· When I divide x by y, I will shift left to get the integer power of two and then decompose the remainder into a binary fraction. I will make up a table of the fractional powers of 2, say to 2^(1/65536). I can drop the 1 to the left of the binary point? and express the rest as a binary fraction. Then using the binary fraction from the remainder decomposition I will multiply the result·for every bit position which has a "1' by the corresponding table entry. Tests using Excel to accomplish this show I should be able to maintain 0.01% accuracy over my intended range of x. Hopefully I will have this coded this weekend.



Thanks for all the ideas & I'll post a reply with my results

Tom


Post Edited (TomS) : 8/24/2007 4:03:27 PM GMT

Tracy Allen

Tom,

That is precisely what is contained in the table in the ROM.

The remainder comes out as a number from 0 to 9030, which implies a fraction from (y//9031) / 9031 in {0/9031 to 9030/9031}. That needs to be renormalized to a binary fraction, that is to a fraction like f11/2048 or f16/65536. The values f11 or f16 provide the closest approximation to the target, specifically, f11/2048 <= (y//9031)/9031 < (f11+1)/2048. I choose 11 bits because that is the number of words in the ROM table. You could also extend it to, say, f16, and use the additional 5 bits for interpolation in the ROM table. The values in the ROM table are 16 bit words, so there is grist there for interpolation.

The code to find the values f11 or f16 is the same code you would use to find the FRQx values for the COG Counters.

pasm

'  general proper fraction, approx y/x by f/(2^n).
' enter with y and x: x > y, x < 2^31, y <= 2^31
' exit with f: f/(2^n) =< y/x =< (f+1) / (2^n)
fractyxb              ' enter here , n contains number of bits
    mov f,#0    ' result
:loop
    shl y,#1        ' double the value of y
    cmpsub y, x   wc,wr  ' if y>=x then y:=y-x, carry:=1, else a unchanged, carry:=0
    rcl f,#1        ' shift in carry, double the estimate of f, appoximation f/(2^n) =~ y/x
    djnz n,#:loop    ' do all bits
fractyxb_ret    ret

spin:

PRI fractyxb (y, x, b) : f            ' calculate f = y/x * 2^b
'' b is number of bits
'' same as f: y/x = f / (2^b) < y/x < (f+1)/(2^b), 
'' that is, f over a power of two is a close appoximation to the original fraction.
  repeat b                      
    y <<= 1
    f <<= 1
    if y => x    '  
      y -= x
      f++

Once you have the 11 bit value, shift is left one and use it as a RDWORD to look up the value in the ROM table. Remember that the value going in is in {0,1} while the value coming out is in {1,2} with the most significant bit (the 17th bit) implied. The power of two in you integer part then determines the final position of the binary radix point.

For interpolation, the two words, Word(f11) and Word(f11+1). Then use the additional bits from f16 (the ones to the right of f11), to interpolate.

▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔
Tracy Allen
www.emesystems.com

Tracy Allen

I think that one difficulty we have with this is to realize that the exponential and the logarithm are periodic functions in much the same way as sine and cosine. It is easy to see when looking at a piece of semilog graph paper. The number base, (2, e, 10 etc) provides the reference point for the periodicity, but they are all essentially the same. Base to base convertion is a simple multiplier. Here is a graph of a base 2 exponential function in blue on the left axis and the periodicity superimposed in yellow on the right axis. The function in yellow is the value in between the integer values on the x axis: 1.x, 2.x, 3.x and so on, but divided by the integer value, so it is p.x/p for p an integer and 0<=x<1. The values over one period suffice to determine the whole function, whether by table lookup or by computation.



▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔
Tracy Allen
www.emesystems.com

TomS

Tracy,

Thanks, I'll give it a go.



Tom

rokicki

Note that by approximating 30K * log_10 2 by 9031 you are already throwing away a certain amount of
accuracy. To pick a number in the middle of your range, 140K, your approximation of 9031 would
give you a result of 215397638.96 where the correct answer is 215443469.00. It's not in the twelfth
bit but it's getting pretty close. To fix this I would just do a scale and divide by as many bits of
9030.89986991943 similarly scaled as you want to use.

TomS

I've already scaled it up so that the denominator is just under 2^16.

It's accurate enough at the values I'll use it.

Tom

TomS

Aaargh, the @*&^** manual and the Propeller Guts document both have the same typo in the Anti-Log routines.

The following line from the manual is wrong!

table_antilog long $C000 'anti-log table base

It should read:

table_antilog long $D000 'anti-log table base

I cut and pasted from the manual and it burned me. At least now I know my procedure works, although without interpolation the error is around 5%. I'll now try interpolation. More to come.

Tom

deSilva

Figure 1-5 on p. 31 in the Manual is correct.

And to my disappointment this typo is not yet in the otherwise excellent supplement/erratum www.parallax.com/dl/docs/prod/prop/PMv1.0Supplement-v1.1.pdf

Post Edited (deSilva) : 8/25/2007 1:35:22 PM GMT

TomS

Yes, it is correct. As is the heading in the appendix. The problem is I didn't double check the code.

I think a sticky thread with all know manual errors would be a good idea. Surely there are more than this one.

deSilva

Good idea! Parallax probably thought they found most -> Supplement v1.1

Tracy Allen

That's a good catch about the typo. For an additional reference, it is very worthwhile to check out Cam Thompson's implementations in float32 and float32full, from the object library. Those libraries use the tables without interpolation for the log and power functions.

I have trouble with some of the terminology in the appendix in the Propeller manual, where it seems to me that terms for "exponent" and "mantissa" are turned on their heads. For example, it says on page 422, "The anti-log table is comprised of 2,048 unsigned words which are each the lower 16-bits of a 17-bit mantissa (the 17th bit is implied and must be set separately)." I would argue that the index for the lookup in that table is actually the mantissa. When we take the logarithm of a number to a certain base, we come up with,
logb x = characteristic.mantissa
the characteristic being the integer part and the mantissa is the fractional part. The anitlog is then,
x = b^(characteristic.mantissa) = b^(characteristic) * b^(mantissa)
Using the HUB table, we are looking up the value b^(mantissa). Mantissa is the index, not the value. We have b=2, so the b^characteristic amounts to a shift operation, right or left.

Tom, I am surprised by the error of 5%. That sounds terrible. 5% of what?

▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔
Tracy Allen
www.emesystems.com

TomS

Tracy,

I'm investigating the error right now. I must have made a mistake in my code. I'm checking it against what my TI-36X calculator says using the y^x function. I bought this calculator several years ago to replace the same model that I dropped an broke. It does fractional exponents just fine. Earlier this year I bought another for work. I tried yesterday to check my programming and it doesn't work with fractional exponents. I guess they changed the code on me.

Tom

TomS

I found the error. When n = 11, f must be shifted left 5 bits, not 1 bit. When n = 16 no shift is required.

TomS

With interpolation I'm getting 15 bit accuracy.

Tom

Fred Hawkins

Tom,
Can you post a minimal code example that documents how to use this table as well as the correct manual citations and the incorrect ones?

Fred

Tracy Allen

Hi Tom,

That's good news about the success.

What n and f are you referring to? I don't understand the error.

In the routine on page 422 of the Prop manual, it starts with a 21 bit exponent, which consists of 5 integer bits followed by 16 bits to the right of the radix point. The code separates out the fractional part so that it represents a number in {0,1} and so that the answer will end up in {1,2}. Then it shifts the 16 fraction bits right by 4 (losing 4 of them--no interpolation), then of the 12 remaining bits it knocks down the least significant bit so that it becomes a word offset into the antilog table, which is achieved by adding the $D000 base address for RDWORD. The 16 bit value retrieved from the table is then left justified in the 32 bit long, except that a 1 is tacked on in bit 31. That represents the integer 1 in the mapping from {0,1} to {1,2}. Then the original exponent is effectively subtracted from 15 (xor $1FF) and that becomes the number of bits that the value is shifted right for the final answer. If the original exponent happened to be zero, then the final answer ends up with an integer part of 1.

I still object to calling the value that is read from the table a mantissa. Maybe it will become one in floating point (which in programming terms is always reduced to exponent.mantissa), but that does not make sense in the integer fixed point calculation. Maybe the text is focused on using this table in the context of floating point.

▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔
Tracy Allen
www.emesystems.com