Calculating x^y in Pbasic

bob kruggel

Can someone tell me how to perform x^y in Pbasic. I have checked the EME list of math functions and found the technique for calculating log(x), but nothing about exponentiation.
Thanks,
Bob Kruggel

TechnoRobbo

Your question is very broad so I assume you require the values of x and y to be somewhat flexible. Unless you have some very specific requirement your looking at creating an entire floating point sub-system.

For $14.95 this system has been created for us.
Check out the uM-FPU at http://www.parallax.com/detail.asp?product_id=604-00030

If you do have a specific application for which that you feel this is chip overkill, then post the specifics and maybe us users can help you code it.

▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔
Have Fun


TR

Chris Savage

Duplicate post removed.

▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔
Chris Savage
Parallax Tech Support

bob kruggel

To the member who responded. I somehow deleted your response and trying to recover it. The equation I want to calculate is:
x=(k1*e^k2)/(e^(k3*y)). k1, k2, and k3 are floating point constants. x and y are floating point variables.
BobK

Suman

The Basic stamp is a 8bit processor and doesnot support any floating point operations. Hence we cannot do any hi-end mathematical operations. Hence u need to use uMFPU which supports 32 bit floating numbers.

It connects to BS2 by SPI. And the editor fot this math co-processor is kinda easy.....
Whatever mathematical equation u write in math coprocessor's editor gets transformed into equivalent PBASIC commands.

Hope this helps,

Suman

Tracy Allen

I agree that in general that kind of problem would be handled elegantly by the Micromega µMFPU sold by Parallax. Especially if it really needs to be floating point and the constants will have to be entered at run time. On the surface of it, there are only two constants:

x=(k1*e^k2)/(e^(k3*y))
becomes
x=k4 / e^(k3*y) where k4 = k1*e^k2

That is not an especially complicated equation that involves a lot of different steps (or maybe this is just part of a longer calculation?) So there is some hope that it could be done in situ. It depends on the expected full-scale range of the variables and the value of the fixed constants.

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

bob kruggel

Tracy,
Where can I get some guidance i solving this equation? I didn't find anything about raising a number to a power at the EME site. It may be too complicated, but I'd like to give it a try.
BobK

Tracy Allen

Hi Bob,

It really does depend on the range of values and precision required. The function you will calculate with the Stamp will be wz = 2^wx, with wx between 0 and 1 and wz between 1 and 2. On the Stamp, unity will be represented as 2^15=32768, that is, with the binary radix just to the right of the MSB. So the calculation is roughly a precision of 4 1/2 decimal digits. That is the fundamental function. The rest is linear scaling. The calculation will have to keep track of the integer part of the exponent separately and apply it after the calculation of the mantissa, and that will be a shift by a power of two, or by virtual movement of the binary radix. There is also a scale factor going into the calculation, 1.4427*k3*x to scale from base e to base 2. (1.4427 = 1/ln(2)), and by your constant k3.

I think I have a PBASIC algorithm for wz = 2 ^ wx, but not on the web site. I'll try to find and extract that.

It would be possible to get more bits with double precision, but it gets complicated on the Stamp. There was recently a thread about b^x on the Propeller forum, but in that case there is a table of logs and exponentials in the ROM, and also there is a 32 bit word to play with.

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

TechnoRobbo

Here's a simple "square and multiply" algorithm for integers. Also, attached is an algorithm I used for fixed point math, it was used with rctime, upper word is the integer lower word is the mantissa. Study the integer one and the fixed point will makes sense. Simply put if the least significant bit is a 1 multiply you base times your result, if not square your exponent. Shift the exponent to the right if the result is 0 your done.

' {$STAMP BS2}
' {$PBASIC 2.5}
expo VAR Word
base VAR Word
temp VAR Word
result VAR Word
Main:
DO
· DEBUG CR,cr, "BASE:"
· DEBUGIN DEC3 Base
· DEBUG CR, "EXPONENT:"
· DEBUGIN DEC3 expo
· result=1
·
'square and multiply
· DO WHILE expo > 0
··· IF (expo & 1) THEN
····· 'base * results
····· result=base*result
····· expo = expo - 1
··· ELSE
····· 'base * base
····· base=base*base
····· expo = expo >> 1
··· ENDIF
· LOOP
· DEBUG CR,"result=",DEC5 result
LOOP
END

▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔
Have Fun


TR

Post Edited (TechnoRobbo) : 9/13/2007 12:22:31 AM GMT

Tracy Allen

Hi TR -- That's an interesting algorithm for the integer exponentiation. In the loop if x is even, then b^x becomes (b^2)^(x/2) and if x is odd, then b^x becomes b*(b^(x-1)). Fun indeed. But the one for fixed point iin the attachment is a lot more complicated and I don't get it. Could you explain a little more what is being calculated there?

---

I'm attaching a program that demos a shift/add algorithm for 2^x. Enter a number from 0 to 9999 to represent a fraction from 0.0000 to 0.9999, and it will calculate 2^x. For example, 3333 gives the cube root of 2. It uses an auxiliary table that contains the values of log2(1 + 1/(2^k)) for k=1 to 15}. The accuracy is limited by roundoff error in the table and in the calculation, but it is close to 4 digits.

Note, This is from Jorg Arndt's Algorithms for Programmers.

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

Post Edited (Tracy Allen) : 9/13/2007 6:56:31 AM GMT

TechnoRobbo

Sure, the code was from an analog input of a pressure sensor and it was measured empirically then I used excel to calculate the base and coefficient.· I simplified it and reposted it by stripping out the multiplication and I will explain the math below in Bold type:

'fixed point example

' {$STAMP BS2}
' {$PBASIC 2.5}

rct VAR Word 'rctime result
mag VAR Word(2)
temp VAR Word(4)
result VAR Word(2)
carry VAR Nib
NewVal VAR Word

Main:
DO
· DEBUG CR,CR,"Input an exponent:"
· DEBUGIN DEC3 rct· '

The code takes a base of 1.010364768 and raises it to the power of rct which is input by·the user

· mag(1)=1:mag(0)=$02A7


The mag variable holds the base in 2 words, mag(1)=1 is the integer portion and mag(0)=$02A7 is the mantissa or decimal portion which was calculated by multiplying 65536 by .010364768

· result(1)=1:result(0)=0

the result is seeded by a 1 ,· result(1)=1 is the integer and result(0)=0 is the decimal portion


· 'square and multiply

the same square and multiply math is used as was used in the integer portion except it's a 2 word by 2 word math which the stamp doesn't perform so I have to do long multiplication with carries.


· DO WHILE rct > 0

··· IF (rct & 1) THEN
····· 'base * results
····· 'r0m0
····· temp(0)=mag(0)*result(0)

····· 'r0m0+ r0m1 r1mo
····· temp(1)=(mag(0)**result(0))
····· carry=0
····· NewVal=temp(1)+ (mag(1)*result(0))
····· IF NewVal < temp(1) THEN carry=carry+1

since the stamp has no overflow flag I check for overflows and incerement the carry after evey multiplication with addition


····· temp(1)=NewVal
····· NewVal=temp(1) + (mag(0)*result(1))
····· IF NewVal < temp(1)THEN carry=carry+1
····· temp(1)=NewVal

····· 'r0m1+ r1mo+ r1m1
····· temp(2)=carry
····· carry=0
····· NewVal=temp(2)+ (mag(1)**result(0))
····· IF NewVal < temp(1)THEN carry=carry+1
····· temp(2)=NewVal
····· NewVal=temp(2)+ (mag(0)**result(1))
····· IF NewVal < temp(1)THEN carry=carry+1
····· temp(2)=NewVal
····· NewVal=temp(2)+ (mag(1)*result(1))
····· IF NewVal < temp(1)THEN carry=carry+1
····· temp(2)=NewVal

····· 'r1m1+
····· temp(3)= mag(1)**result(1)+carry

····· result(0)=temp(1):result(1)=temp(2)
····· rct = rct - 1
··· ELSE
····· 'base * base
····· 'r0m0
····· temp(0)=mag(0)*mag(0)

····· 'r0m0+ r0m1 r1mo
····· temp(1)=(mag(0)**mag(0))
····· carry=0
····· NewVal=temp(1)+ (mag(1)*mag(0))
····· IF NewVal < temp(1)THEN carry=carry+1
····· temp(1)=NewVal
····· NewVal=temp(1) + (mag(0)*mag(1))
····· IF NewVal < temp(1)THEN carry=carry+1
····· temp(1)=NewVal

····· 'r0m1+ r1mo+ r1m1
····· temp(2)=carry
····· carry=0
····· NewVal=temp(2)+ (mag(1)**mag(0))
····· IF NewVal < temp(1)THEN carry=carry+1
····· temp(2)=NewVal
····· NewVal=temp(2)+ (mag(0)**mag(1))
····· IF NewVal < temp(1)THEN carry=carry+1
····· temp(2)=NewVal
····· NewVal=temp(2)+ (mag(1)*mag(1))
····· IF NewVal < temp(1)THEN carry=carry+1
····· temp(2)=NewVal

····· 'r1m1+
····· temp(3)= mag(1)**mag(1)+carry

····· mag(0)=temp(1):mag(1)=temp(2)
····· rct = rct >> 1
··· ENDIF
· LOOP

I now convert the·least significant word back into decimal by multiplying it by 1000 and displaying the leading·zeros with dec4


· Newval=1000**result(0)
· DEBUG CR,"Final result=", DEC result(1),".",DEC4 Newval,CR


LOOP
END

▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔
Have Fun


TR

bob kruggel

To all the recent posts,

The floating point processor is looking better all the time. Making Pbasic jump through hoops is tough when Pbasic doesn't have legs.

BobK

Tracy Allen

TR, Thanks for clarifying that. I'd been thrown off by the connection with RCTIME and the pressure sensor and didn't focus in on the double precision calculation. Nice work on that, by the way! What a pain with the carry, huh?

Bob, you can't go wrong with the µFPU coprocessor. It will save you a lot of headaches. The time you might want to code it directly is if it is not going to change much in a dedicated system, and it is a relatively simple equation. The Stamp has strong legs, but they don't go easily in certain directions.

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

TechnoRobbo

Bob, If you really like a good challenge, run the Iditarod every year,eat your salsa hot and spicy and never want to leave your home office and spend time with your family then the by all means write the math. On the other hand if you actually want to accomplish something in a reasonable amount of time then buy the math coprocessor.

Personally I'm a math geek and a huge Lt. Columbo fan so I'll do the first.

▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔
Have Fun


TR

Post Edited (TechnoRobbo) : 9/13/2007 10:51:05 AM GMT

Tracy Allen

TR,
The base of b = 1.010364768 seems pretty small, but when combined with the integer result from RCTIME it grows, uhhh, exponentially, so I'm guessing that the values of RCTIME you have to deal with are constrained to low values, right? It turns our that b^1075 = 65171 will fit in one word, but b ^ 1076 = 65847 will not. And b^65535 is googleplexian.

The same problem can be reworked as a power of 2.
b^x = 2^(0.01487623751 * x)

In Stampese, using (0.01487623751 * 65536 = 975)

xInteger = 975 ** x ' over the range x=1 to x=1075, xInteger can be as high as 14
xRemainder = x - (xInteger * 975)

And from that using binary division we can find the binary fractional part, xRemainder/975, and then use the algorithm I posted above for 2^x, when x is in the interval {0,1}. The final answer is then a shift by the power of 2 contained in xInteger. All that without resorting to double precision, to get the same 3 or 4 digits in the answer. Or, there might be double precision involved, but only a shift by power of 2. The point is that the function 2^x in the interval 0<=x<1 is cannonical and can be scaled to solve any problem involving b^x with any base b and over any range of x.

But not. The thing about undertaking these math problems on the Stamp especially is the sense of scale, the hard wall at 65535. Bob, it is also something to be kept in mind with the coprocessor. There are problems that are so "stiff" that even a coprocessor will flake out.

I'll also lay claim to being a math geek,. Add long bike rides to the hot salsa and bumbling along slowly to solutions like good Lt. Columbo.

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

TechnoRobbo

Yes, your right Tracy - they were small values and after all that work we ended up with a lookup table. That goes to a very important point that sometimes the math is not necessary if there are a limited number of scenarios.

▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔
Have Fun


TR