How to resolve a function of X.

[david]

Hi All,
I'm looking for possible options to resolve a function of X given as-
f(x)= -9.42477E-06 x^2 + 0.44076 x
The total range of x is 0 to about 6000.

The function is "nearly" linear and for values of x below about 400 the x^2 term can be ignored.  I'm not looking for decimal point accuracy but the nearest integer (or two) would be nice.

In the past I would have multiplied the coefficient of the x^2 term and later divided the result to restore things but I can also do a fairly coarse lookup table and linear interpolate the intermediate values.  Is there a better approach to resolving nasty polynomials?

Cheers,
David

[JonW]

Rather than coding the full function, the only other way I have used is the reduced lookup and linear interpolation method.  Have not plotted your function but in some cases I have used a more detailed lookup in areas where gradients change the most.

[david]

Thanks for getting back.
Yes the full maths solution could get ugly and would still likely involve some rounding off. 
The lookup table and interpolation seems like an easy out for me unless someone else has a different approach.
Plot attached - very gentle curve.

Cheers,
David

[JonW]

What about breaking it down into 2 or 3 or more simple straight line equations?
One proc to compute the y=mx+c and pass the constants based on a case select of x?

[david]

2 or 3?  I was all set to throw about 64 points at it which would break it down to every 100 values and I could then linear interpolate those intermediate 100 steps.  On paper it's accurate enough and should be easy code wise.  Perhaps 64 points is an overkill...

Cheers,
David

[JonW]

The curve is "almost" linear and you did state integer or two, you can easily feed it into excel with the number of linear steps and calculate the error.  Probably be alot less than you think.
I did say more too

[david]

I'm working with a spreadsheet but haven't yet calculated error - it just looks acceptable.
You did indeed say 2 or 3 or more.  I shall check out the error.  For this particular use I'm more interested in the absolute error rather than a % error.  With enough data points it should approach an absolute error.
Thanks for your input.  It's always good to get a second opinion on these things.

Cheers,
David

[trastikata]

f(x)= -9.42477E-06 x^2 + 0.44076 x
The total range of x is 0 to about 6000.

Is there a better approach to resolving nasty polynomials?

Hello David,

I don't understand what are you trying to achieve? Do you wanna solve the equation for F(x), given x. Or to find x, given F(x)?

In first case, find F(x), given x - [f(x)= -9.42477E-06 x^2 + 0.44076 x] It is a simple math and will be faster to simply use 3 floating point operations. instead of long look-up tables.

In case you want to find x, given F(x) - this is solving a parabola formula [x = -b ± sqrt(b^2) / 2a], which takes about the same number of floating point operations.

Why do you want to complicate things with long look-up tables (100's of values) or resolving to interpolation? With this number of points solving the equation for either F(x) or x will be faster and neater.

[John Drew]

Hi David,
I'm with Trastikata. I've used both interpolation and arithmetic with floats for solving polynomials up to third order
I've mostly ended up using floating point arithmetic. Further, I've found a worthwhile accuracy and speed improvement by breaking the calculation into multiple lines.
e.g in your case of f(x)= -9.42477E-06 x^2 + 0.44076 x
ftemp1 = x * x   
ftemp2 = -9.42477E-06 * ftemp1   'fTemp2 contains the result of the first term
fTemp1 = 0.44076 * x                     'fTemp1 contains the result of the second term
Result = fTemp2 + fTemp1             'Result now contains the addition of the 2 terms

I have some very complex maths in my moon tracking program in an 18F4685 and I found breaking the line into parts resulted in more accurate results. And this was with 32 bit floats, of course if you wanted even greater accuracy you could go to a PIC24 and use Les's 64 bit arithmetic.
John

[david]

Hello Trastikata,
I'm trying to find f(x) given a value of x.   Three floating point operations sounds great but where I'm lost is the term -9.42477E-06.  How do you work with such a tiny number - do you multiply by 1,000,000, do the maths then divide again as I mentioned earlier? 

John - I was hoping you may chip in as I knew you did some complex maths for tracking and Great Circle calcs.  Taking your example, the term  ftemp2=-9.424477E-06*ftemp1  - does the compiler recognise exponents like E-06?  Have I been nodding off in class again?  That would change everything and I can understand why you and trastikata are in disbelief of my lookup table approach.
Now....the device I'm using is a modest 16F684 - will that change things? 

Many thanks for your input guys.  I think I may have been making life hard for myself.

Cheers,
David 

[trastikata]

Hello Trastikata,
I'm trying to find f(x) given a value of x.   Three floating point operations sounds great but where I'm lost is the term -9.42477E-06.  How do you work with such a tiny number - do you multiply by 1,000,000, do the maths then divide again as I mentioned earlier? 

Code Select Expand

Dim x As Word
Dim y As Float

Dim  dTemp As Dword
Dim sdTemp As SDword
Dim fTemp1 As Float
Dim fTemp2 As Float

Main:
   
'f(x)= -9.42477E-06*x^2 + 0.44076 * x   
'-9.42477E-06 = -58903 / 6250000000
'f(x) = (-58903 * x * x / 6250000000) + 0.44076 * x

dTemp = x * x                   'x^2
sdTemp = -58903 * dTemp         '-58903 * x * x; given 0 =< x =< 6000 ==> sdTemp(min) = -58903 * 6000 * 6000 = -2120508000000 ==> thus no overflow for sdTemp
fTemp1 = 0.44076 * x            '0.44076 * x
fTemp2 = sdTemp / 6250000000    '(-58903 * x * x / 6250000000)
fTemp1 = fTemp1 + fTemp2        '(-58903 * x * x / 6250000000) + 0.44076 * x


[John Drew]

G'day David,
Floats work fine in the 14 bit devices. The 18F series are better but seems the 16f684 might be a spare PIC hanging around 😁
Positron accepts the scientific notation  as written in the polynomial. The letter E is followed by the exponent. Use my suggestion and it will work.

I'm sure you know that the E-6 is just a place keeper i.e. move the decimal point to the left 6 places. It is still dealing with 9429477 (the significand) and keeping track of the exponents.
I see Trastikata has just posted. It's another approach. It will be interesting to compare the two approaches. I might do some experiments tomorrow. .
Cheers
John

[david]

Hello trastikata,
Many thanks for the explanation - yes I can follow that ok.
May I ask how you converted the 2nd line to that fraction?  (-9.42477E-06 = -58903 / 6250000000)  Why the multiply by 6250?  What about -3/318310
There used to be an excellent website that did decimal to fraction conversions with higher precision as you fed more decimal places in to it.  Sadly it vanished but just now I found this site which looks promising-
https://www.gigacalculator.com/converters/convert-decimal-to-fraction.php
I shall do some playing with the numbers but I'm sure you're right and that the maths conversion is the way to go. 

Cheers,
David

[david]

Hello John,
Well I never knew that Positron accepted scientific format - that's news to me.
So with the new decimal to fraction website and scientific notation life could be much easier than the look up table I was considering.
My thanks to both you and trastikata.

Cheers,
David

[tumbleweed]

Here's a comparison of a few different coding methods, including trastikata's.
It shows the code size and lib calls used (for the 18F... didn't change them for 16F)
Just change '$define CODE_SELECT 1' to the values 1 to 5 to see the different results.

CODE_SELECT 1 or 2 is the most straight-forward, but 3 is the smallest.

Code Select Expand

'device = 18F25K20
Device = 16F684
Declare Xtal = 8

'f(x)= -9.42477E-06*x^2 + 0.44076 * x
'The total range of x is 0 to about 6000

' set CODE_SELECT 1 thru 5
$define CODE_SELECT 1

'------------------------------------------------------
$if CODE_SELECT = 1
'18F = 870 program bytes, 30 Variable bytes
'16F = 514 program bytes, 30 Variable bytes
Dim x As Word
Dim f As Float

x = 6000
f = -9.42477E-06 * x * x + 0.44076 * x
'    rcall _uns_int16_tofl32
'    rcall __fpmult_32__
'    rcall _uns_int16_tofl32
'    rcall __fpmult_32__
'    rcall __fpsub_32__
'    rcall _uns_int16_tofl32
'    rcall __fpmult_32__
'    rcall __fpadd_32__
$endif

'------------------------------------------------------
$if CODE_SELECT = 2
'18F = 768 program bytes, 32 variable bytes
'16F = 451 program bytes, 32 variable bytes
Dim xf As Float
Dim f As Float

xf = 6000.0
f = -9.42477E-06 * xf * xf + 0.44076 * xf
'    rcall __fpmult_32__
'    rcall __fpmult_32__
'    rcall __fpsub_32__
'    rcall __fpmult_32__
'    rcall __fpadd_32__
$endif

'------------------------------------------------------
$if CODE_SELECT = 3
'18F = 720 program bytes, 28 variable bytes
'16F = 427 program bytes, 28 variable bytes
Dim xf As Float
Dim f As Float

xf = 6000.0
'f = -9.42477E-06 * xf * xf + 0.44076 * xf
f = xf * xf
f = -9.42477E-06 * f
f = f + 0.44076 * xf
'    rcall __fpmult_32__
'    rcall __fpmult_32__
'    rcall __fpsub_32__
'    rcall __fpmult_32__
'    rcall __fpadd_32__
$endif

'------------------------------------------------------
$if CODE_SELECT = 4
'18F = 960 program bytes, 40 variable bytes
'16F = 527 program bytes, 40 variable bytes
Dim x As Word
Dim xd As Dword
Dim f As Float

x = 6000
'f = -9.42477E-06 * xf * xf + 0.44076 * xf
xd = x * x
f = -9.42477E-06 * xd
f = f + 0.44076 * x
'    rcall __multiply_u3232_
'    rcall _sgn_int32_tofl32
'    rcall __fpmult_32__
'    rcall __fpsub_32__
'    rcall _uns_int16_tofl32
'    rcall __fpmult_32__
'    rcall __fpadd_32__
$endif

'------------------------------------------------------
$if CODE_SELECT = 5
'18F = 1158 program bytes, 44 variable bytes
'16F =  646 program bytes, 44 variable bytes

'-9.42477E-06 = -58903 / 6250000000
'f(x) = (-58903 * x * x / 6250000000) + 0.44076 * x

Dim x As Word

Dim  dTemp As Dword
Dim sdTemp As SDword
Dim fTemp1 As Float
Dim fTemp2 As Float

x = 6000
dTemp = x * x                   'x^2
'    rcall __multiply_u3232_
sdTemp = -58903 * dTemp         '-58903 * x * x; given 0 =< x =< 6000 ==> sdTemp(min) = -58903 * 6000 * 6000 = -2120508000000 ==> thus no overflow for sdTemp
'    rcall __multiply_u3232_
fTemp1 = 0.44076 * x            '0.44076 * x
'    rcall _uns_int16_tofl32
'    rcall __fpmult_32__
fTemp2 = sdTemp / 6250000000    '(-58903 * x * x / 6250000000)
'    rcall _sgn_int32_tofl32
'    rcall __fpdiv_32__
fTemp1 = fTemp1 + fTemp2        '(-58903 * x * x / 6250000000) + 0.44076 * x$endif
'    rcall __fpadd_32__
$endif



[John Drew]

Tumbleweed,
I recall from my experiments to squeeze the last little bit of code size and accuracy it was better not to use the construct of X=X etc but to use intermediate temporary variables.
It's 23:45 here but will get back tomorrow morning.
David has enough to go on with, but our comments may help others in future.
John

[tumbleweed]

Accuracy-wise it shouldn't make a difference, but that would depend on exactly what the equation was, the variable ranges, and how it's coded. Things like adding a small FP value to a large FP value can lose digits, so sometimes the order of things can be optimized to help that. Sometimes not.

As far as the code size goes, many times splitting it up helps the compiler not create as many temp variables to store intermediate results. That's what the difference between CODE_SELECT 2 and 3 shows.

[trastikata]

May I ask how you converted the 2nd line to that fraction?  (-9.42477E-06 = -58903 / 6250000000)  Why the multiply by 6250?  What about -3/318310
....
I shall do some playing with the numbers but I'm sure you're right and that the maths conversion is the way to go. 

I overlooked something:

- 6250000000 is off the range of 32b integers
- the OP says -9.42477E-06 and I somehow typed -9.42448E-06

Anyway, you got the idea. The decimal is converted to a an acceptable fraction. You can use this site:

https://www.calculatorsoup.com/calculators/math/decimal-to-fraction-calculator.php

If the fraction contains integers out of the range of 32b math, you can round the last digit (if some accuracy degradation is accepted) :
-9.42477E-06 = -0.00000942477 ==> round to -0.0000094248 = -94248 / 10000000000, divide the numerator and denominator by the same
number until you get fraction in any of both, thus 80 in this case and you get = -11781 / 1 250 000 000 --> 1250000000 within the range of 32b math.

[JonW]

Even if you round it to -0.0000094 * x^2 + 0.441x its <  0.1% error

[david]

What a lively discussion and my thanks to all contributors for the excellent points raised.
Key takeaways for me-
1. I'm dumber than I was 5 years ago when I would have readily used my magic decimal to fraction converter.
2. I'm really glad I just found a replacement converter after years of avoiding nasty fractions.
3. That the new converter I linked to is easier than trastikata's one.  (sorry but I think it is)
4. That Positron accepts scientific number formats.  (thanks John)
5. That JonW is right in that last post.  The function is nearly linear so the precision of the x^2 is not very important.
6. That -3/318310 is close enough to -9.42477E-6.  Actually -1/106112 is not too bad either....

Cheers to All,
David