pow(a,b)

[DiTBho]

Given two real (floating point or fixedpoint), `a` and `b`, where `b`  can be negative, how to efficiently compute the power function `pow(a, b)` ?

Mathematically, it's required to compute pow(a, b)=exp(b*ln(a))

I have to implement it in C, and then, I have to rewrite it in assembly 

edit: Typo: a^b.
a,b are real numbers
b can be negative.

[bpiphany]

Aren't Taylor series and Newton-Raphson's method the answer to most computations on irrational functions?

Also, google is your friend
https://stackoverflow.com/questions/4518011/algorithm-for-powfloat-float
https://stackoverflow.com/questions/9799041/efficient-implementation-of-natural-logarithm-ln-and-exponentiation/44232045

[Zero999]

Other than for a learning exercise, the answer normally is to use the compiler's floating/fixed point library and a powerful enough CPU, with sufficient memory, to handle it.

[newbrain]

Given two real (floating point or fixedpoint), `a` and `b`, where `n`  can be negative,

And 'n' would be? Is it 'a' or 'b'?
As it makes quite a complex difference (pun totally intended).

[T3sl4co1l]

b being negative isn't an issue, that's just fractional division rather than multiplication, and the definition of the exponent function provides continuity.  Surely, a being negative is the bigger problem -- fractional powers are imaginary but integer powers are just fine?
https://en.cppreference.com/w/c/numeric/math/pow

Tim

[SiliconWizard]

Given two real (floating point or fixedpoint), `a` and `b`, where `b`  can be negative, how to efficiently compute the power function `pow(a, b)` ?

Mathematically, it's required to compute pow(a, b)=exp(b*ln(a))

I have to implement it in C 

I suppose you're not allowed to use the C std math library, because it's all there already. And if you have access to a compliant C compiler, there's no reason you wouldn't have access to this library.
So I'm assuming it's just an exercise. Or maybe you're looking for something more efficient (in which case, have you benchmarked the standard pow() function already to know what you'd be competing against?)

You can take a look at the following implementation in openlibm:

https://github.com/JuliaMath/openlibm/blob/master/src/e_pow.c

[rstofer]

You can take a look at the following implementation in openlibm:

https://github.com/JuliaMath/openlibm/blob/master/src/e_pow.c

That code doesn't look like a homework assignment.  I can imagine a team of Applied Math PhDs sitting around discussing the boundary conditions for weeks.  Then they spend even more time coming up with the algorithm for the programmer to code.  There's some magic math in that code.

[SiliconWizard]

You can take a look at the following implementation in openlibm:

https://github.com/JuliaMath/openlibm/blob/master/src/e_pow.c

That code doesn't look like a homework assignment.  I can imagine a team of Applied Math PhDs sitting around discussing the boundary conditions for weeks.  Then they spend even more time coming up with the algorithm for the programmer to code.  There's some magic math in that code.

Definitely. But it gives the basic calculations they use, so that can serve as an inspiration for implementing the same, even if you end up with something a lot less optimized.

Now I don't know if the OP's question was actually a homework assignment. If it was, then giving them a "homework assignment" directly usable answer would not do them any favor IMHO. And I also don't know what the value of such an assignment would be. If it were a programming assignment, then I don't really see the point of using anything else than the std pow() function. But if it were an applied math assignment, then it would be interesting, but in this case I would expect the students to at least have a rough idea of what kind of approach they can use to get an approximation of this function.

[hamster_nz]

BKM might be useful. It is a bit difficult to grok though.

https://en.m.wikipedia.org/wiki/BKM_algorithm

[DiTBho]

just an exercise

Kind of personal challenge for which I will have to rewrite it in pure assembly in order to fit a very old CPU with limited resources.

Crazy, I know 

It's done for a larger hobby project, where for historical reasons I want to support all the technology I encountered during the last 25 years of my job and hobby experiences, so in the main big project there are several and different pieces of technology, including things used in the 80s, 90s, 2000s, 2010s and 2020s.

For the piece of the puzzle belonging to the 80s, I chose to use a very very old MPU to realize a height sensor, which is used in a part of the big project; in its planned scheme done on Labview as "proof of concept" between the CPU and the analogical sensor, Data linearization is done using the nonlinear curve fitting method, which requires three pow(a,b) calculations like:

h = K1*pow1(a1,b1) + K2*pow2(a2,b2) + K3*pow3(a3,b3) + q

a1,a2,a3 come from ADC sampling, so the final functionality will be height = function(ADC sampled value).
b1,b2,b3,k1,k2,k3 come from a wizard functionality offered by Labview to synthesize an interpolator.

The CPU is clocked at 2Mhz, it's very slow and limited compared to a modern MPU, there are no time-critical constraints but rather ram-rom-space constraints, since there is only 768 byte of RAM, 1Kbyte of UV-EPROM for the whole program, which should also do other stuff, like talking on the serial line, sampling data on two analog channels, properly filter data, and handling two digital lines to control a couple of motors.

The current firmware is already written in assembly. The pow(a,b) will be written in C, then translated into assembly, deep studied, and finally rewritten from scratch in assembly.

[iMo]

There is certainly a floating point library in asm (or binary) written for a similar MCU (ie from the same family but with larger rom) available.. For example a 32bit floating point lib (incl. sin, cos, tan, exp, ln, pow..) is usually ~3.5kB large with those 8bitters MCUs.. Even I doubt you would be able to cut a bigger chunk off the lib such it fits into your 1kB incl your other stuff..
In past 1kB was enough for a complete 8-9digits floating point math in vintage scientific calculators - but those MCUs were of a different arch (BCD instead of BIN)..

[T3sl4co1l]

Does LabVIEW not produce integer coefficients (a simple polynomial fit)?

If not, you will find precisely because of this, why it is cheaper to do so.  On a machine with hardware multiply, it's just a matter of pulling down parameters and chugging, on loop (Horner's method); without, it's calls/jumps to a multiplication routine (typically shift and accumulate).  If you can constrain LabVIEW to produce results of this form, you will have a much easier -- and faster -- time.

There are also rational approximation methods ( https://en.wikipedia.org/wiki/Pad%C3%A9_approximant ), which on an add-shift machine, might be competitive with polynomial methods -- division is a shift-subtract algorithm, taking similar time to multiplication but fewer terms may be necessary.

I suppose this isn't an excuse to wholly avoid implementing pow() -- but it can push off the need for it quite far, and indeed the same method, loaded with different coefficients, might serve to implement that very function, at least over a certain range, and the range (and selection of coefficients) can then be managed by an outside function.

How many things are you looking to get into that MPU at once, anyway?  And how much storage does it have accessible?  If practically unlimited (e.g. parallel or SPI Flash chip, SD card*), how will you access and manage it?  With so little ROM to work with, it seems likely that some kind of very basic, very general, interpreted language will be necessary.  Perhaps even just to implement access to such things.  Which means any real "thinking" (like evaluating pow()) will take very long indeed -- fast enough for a handheld calculator perhaps, but doubtful for something like a control loop.  Well, maybe like temperature control, there's a lot that can be done between seconds of samples.  But besides that...

*Probably not, at least not in the most convenient form -- FAT16 takes a few kB to implement I think, and is only good to a few hundred MB capacity (or a couple GB with some inconvenience); FAT32 takes a few kB more, and is good for pretty much whatever size is available these days.  You might be able to access the SD card as a flat memory store, with less code, but not with FAT, and so you'll have to copy data to the SD card using some raw method (disk image, hex editor, etc.).

Hmm, I wonder if you can format an SD card for FAT32 and place exactly one file on it, and have that file be at a consistent sector offset, no fragmentation, etc...

Tim

[SiliconWizard]

Given the very limited resources, doing any FP might be extremely challenging - probably impossible. I'm assuming your old CPU doesn't have any FPU. Doing FP in software with this amount of memory? I doubt it.

Anyway, it's probably what T3sl4co1l more or less suggested, but I would switch to piecewise linear interpolation instead of trying to compute those power functions. Just find a decent piecewise linear approximation of each power term. Then the only computation needed will be segment lookup, linear interpolation and additions. All this using integer/fixed-point arithmetic.

[brucehoult]

What precision is your fixed or floating point? To what precision do you need the result? Do you have an FPU? How much program code size is OK to use? How much ROM is ok to use for tables? How much working RAM can you use?

For example, if you are using any 8 bit fixed point or floating point format, and you can afford 64 KB of ROM, then you can get the most accurate possible result with one shift-and-OR and a memory load. BOOM!

[DiTBho]

Does LabVIEW not produce integer coefficients (a simple polynomial fit)?

Unfortunately no. I get irrational coefficients.

[brucehoult]

there is only 768 byte of RAM, 1Kbyte of UV-EPROM for the whole program, which should also do other stuff

Ugh.

You're clearly not going to be aiming for a fast implementation. Just getting something to fit into that size at *all* will be challenging.

It would have been useful to have this information at the outset.

We still don't know what precision and range you need to work with, either for the inputs or the result.

[DiTBho]

To what precision do you need the result?

I am trying to amplify the signal 4X in order to achieve better precision. The on board ADC is an 8bit unit working in the range of 5V (VRefL= 0V VRef_H= 5V), which means ~20mV

5V/(255) is ~20mV

But! Since the measure range is 600..1300mm, I can isolate the corresponding sub voltage range, which is around 0V..1.2V, and can be 4X amplified by an OA.

Still ~20mV, but if the analogical sensor is stable enough, this will produce better numerical results.

Do you have an FPU?

No 

How much program code size is OK to use?

At the moment I am using 768 byte of ram in total, but I am ready to prepare a voltage pump in order to program the onboard UV-Eprom (1Kbyte).

For example, if you are using any 8 bit fixed point or floating point format, and you can afford 64 KB of ROM, then you can get the most accurate possible result with one shift-and-OR and a memory load. BOOM!

64Kbyte is impossible for my system. The best I can do is 32Kbyte for both RAM and ROM.

[brucehoult]

To what precision do you need the result?

I am trying to amplify the signal 4X in order to achieve better precision. The on board ADC is an 8bit unit working in the range of 5V (VRefL= 0V VRef_H= 5V), which means ~20mV

5V/(255) is ~20mV

But! Since the measure range is 600..1300mm, I can isolate the corresponding sub voltage range, which is around 0V..1.2V, and can be 4X amplified by an OA.

Still ~20mV, but if the analogical sensor is stable enough, this will produce better numerical results.

But what numeric range do you consider these 8 bit values to be in? Is it the same for A and B? And the result?

Is B an arbitrary number, or is there only a limited set of possibilities that depend on the algorithm, not on the ADC data?

Depending on your answers, you can very quickly find that most of the results are either 0 or infinity.

[DiTBho]

But what numeric range do you consider these 8 bit values to be in?

the built-in ADC is 8bit, with an equivalent bit noise of 1 bit after filtering, that means what you read and filter is truly 20mV of resolution. Filtering costs something like 200 cycles due to filtering window.

The plan is converting the raw data from the ADC into floating-point/fixed-point, applies some interpolating function

Is it the same for A and B? And the result?

Values "a" and "b" in the exp(a,b) function come from Labview, they are both real rumbers, and they require to be floatingpoint or fixedpoint, but I cannot yet qualify how many fractional bits they need because I am still playing with Labview in order to optimize those values.

Is B an arbitrary number, or is there only a limited set of possibilities that depend on the algorithm, not on the ADC data?

There will be three "b" values, "b1", "b2", "b3", to be used in three "exp" function calls, but they are guaranteed to be constant (Labview says that), while "a" will be a kind of constant that multiply the value from the adc, kind of "a = k * adc_value"

and you have to calculate exp(a,b), so exp ( k_const * adc_value, b_const)

Depending on your answers, you can very quickly find that most of the results are either 0 or infinity.

Yup, including the LUT strategy, which seems more feasible although I haven't yet understood how to do it correctly 

[DiTBho]

This stuff is so impossible that I really wonder how Nasa did program the MPUs for the Apollo space missions. Their guidance computer has 2K words of erasable magnetic-core memory and 36 kilowords of read-only core rope memory

And this stuff has to do more tasks than just controlling the position of a vertically movable platform ... always impressed by those Big Guys!

[brucehoult]

This stuff is so impossible that I really wonder how Nasa did program the MPUs for the Apollo space missions. Their guidance computer has 2K words of erasable magnetic-core memory and 36 kilowords of read-only core rope memory

And this stuff has to do more tasks than just controlling the position of a vertically movable platform ... always impressed by those Big Guys!

I'm sure you must be overthinking all of this. It must be able to be linearised or something over the values in the possible range. And especially only talking about 8 bit precision anyway.

[DiTBho]

I'm sure you must be overthinking all of this. It must be able to be linearised or something over the values in the possible range. And especially only talking about 8 bit precision anyway.

How can you be sure if you haven't seen anything?

I wrote the above showed piece of code to acquire data and send them over the serial, so I sampled data within the range {15..1600}mm and put these samples inside a Labview's table, and when I click on the button it suggested a linearization model made with three pow(a,b), as I described above.

It means, if you want to linerize these

...
2.75   015 cm
2.55   020 cm
2.00   030 cm
1.55   040 cm
1.25   050 cm
1.15   060 cm
0.90   070 cm
0.80   080 cm
0.75   090 cm
0.65   100 cm
0.60   110 cm
0.55   120 cm
0.50   130 cm
0.44   140 cm
0.45   150 cm
...

You need to implement this: h = K1*pow1(a1,b1) + K2*pow2(a2,b2) + K3*pow3(a3,b3) + q

You have to optimize a1,b1,a2,b2,a3,b3,q, but you will always have to deal with fpow, fmul, fadd, fand fcmp

So I tried to consider only the subrange of my interest, just to see if it somehow look "locally linear".

1.15   060 cm
0.90   070 cm
0.80   080 cm
0.75   090 cm
0.65   100 cm
0.60   110 cm
0.55   120 cm
0.50   130 cm

Unfortunately, Labview says that it does not, and still suggest to synthesize two exp(a,b) functions.

So we have just moved from this

h = K1*pow1(a1,b1) + K2*pow2(a2,b2) + K3*pow3(a3,b3) + q

to this

h = K1*pow1(a1,b1) + K2*pow2(a2,b2) + q

[paf]

What is the precision/accuracy of the sensor? 

How did you get the data?

This part of the data suggests there the sensor measurements are not "very accurate":
 

Code: [Select]

...
0.50   130 cm
0.45   150 cm
0.44   140 cm
...


My suggestion is transform x into 1/x:

Code: [Select]

0,36364 15
0,39216 20
0,5 30
0,64516 40
0,8 50
0,86957 60
1,11111 70
1,25 80
1,33333 90
1,53846 100
1,66667 110
1,81818 120
2              130
2,27273 140


This way you can get a more linear relation.
 

[ledtester]

It means, if you want to linerize these

...
2.75   015 cm
2.55   020 cm
2.00   030 cm
1.55   040 cm
1.25   050 cm
1.15   060 cm
0.90   070 cm
0.80   080 cm
0.75   090 cm
0.65   100 cm
0.60   110 cm
0.55   120 cm
0.50   130 cm
0.44   140 cm
0.45   150 cm
...

You need to implement this: h = K1*pow1(a1,b1) + K2*pow2(a2,b2) + K3*pow3(a3,b3) + q

Maybe I'm misunderstanding something, but have you tried piecewise curve fitting?

For the above numbers you get a pretty good fit with quadratic curve fitting if you break up the input range into three intervals: [2.75 .. 1.55], [ 1.55 .. 1.15 ] and [ 1.15 .. 0.45 ].

[DiTBho]

Unfortunately the sensor doesn't have any good documentation.
I don't know anything about its precision/accuracy.

Code: [Select]

              _______
             /0..20Hz|
sensor _____/ OA G=1 |_____  ADC 8bit _____ raw data over serial line
            \follower|       CPU builtin
             \_______|

This is what I am doing at the moment, and how I collected data to prepare the above tables.

but have you tried piecewise curve fitting?

Not yet. Good point 

I have only tried the "NI myDAQ" curve fitting software offered by LabVIEW.
The CPU reads the ADC, and sends raw data to the computer without any manipulation.
The PC collects data and passes them to the myDAQ software.

I am also preparing a similar working scenario for Octave + GNUplot.

For the above numbers you get a pretty good fit with quadratic curve fitting if you break up the input range into three intervals: [2.75 .. 1.55], [ 1.55 .. 1.15 ] and [ 1.15 .. 0.45 ].

Nice! Thanks! Maybe breaking up the input range into several intervals is the best way to sort it out!
Thanks, again, I hadn't thought about it yet