Interesting - thanks for that. I wonder if there is a benchmark for precision?

Re speed, I have the following results, all on standard clocks:

• PBD-SF6W 78μs (hardware FP - 64 bit)
• PBD-SF2W 109μs (hardware FP - 32 bit)
• ESP32-S3 166μs (hardware FP - 32 bit) (I think)
• RP2040 394μs (software FP - 32 bit)
• Unix build on my laptop - for a laugh - 7μs (!)

There's indeed a significant problem.

void test(void) {
    float value = 1;
    int n = 1;
    while (n <= 27) {
        printf("%3d %33.30f\n", n, value);
        value += 1.0 / (1<<n);
        n++;
    }
}
int main(void) {
    test();
    return 0;
}

produces:

  1  1.000000000000000000000000000000
  2  1.500000000000000000000000000000
  3  1.750000000000000000000000000000
  4  1.875000000000000000000000000000
  5  1.937500000000000000000000000000
  6  1.968750000000000000000000000000
  7  1.984375000000000000000000000000
  8  1.992187500000000000000000000000
  9  1.996093750000000000000000000000
 10  1.998046875000000000000000000000
 11  1.999023437500000000000000000000
 12  1.999511718750000000000000000000
 13  1.999755859375000000000000000000
 14  1.999877929687500000000000000000
 15  1.999938964843750000000000000000
 16  1.999969482421875000000000000000
 17  1.999984741210937500000000000000
 18  1.999992370605468750000000000000
 19  1.999996185302734375000000000000
 20  1.999998092651367187500000000000
 21  1.999999046325683593750000000000
 22  1.999999523162841796875000000000
 23  1.999999761581420898437500000000
 24  1.999999880790710449218750000000
 25  2.000000000000000000000000000000
 26  2.000000000000000000000000000000
 27  2.000000000000000000000000000000

whereas

#!/usr/bin/env python
# -*- coding: utf-8 -*-

from array import array

def test():
    value = array('f', (1,))
    n = 1
    while n <= 27:
        print(f'{n:3d} {value[0]:33.30f}')
        value[0] += 1.0 / (1<<n)
        n += 1

test()

produces:

  1   1.00000000000000000000000000000
  2   1.50000000000000000000000000000
  3   1.75000000000000000000000000000
  4   1.87500000000000000000000000000
  5   1.93750000000000000000000000000
  6   1.96875000000000000000000000000
  7   1.98437500000000000000000000000
  8   1.99218750000000000000000000000
  9   1.99609375000000000000000000000
 10   1.99804687500000000000000000000
 11   1.99902343750000000000000000000
 12   1.99951171875000000000000000000
 13   1.99975585937500000000000000000
 14   1.99987792968750000000000000000
 15   1.99993896484375000000000000000
 16   1.99996948242187500000000000000
 17   1.99998474121093750000000000000
 18   1.99999237060546875000000000000
 19   1.99999618530273437500000000000
 20   1.99999809265136718750000000000
 21   1.99999904632568359375000000000
 22   1.99999952316284179687500000000
 23   1.99999980926513671875000000000
 24   1.99999990463256835937500000000
 25   2.00000000000000000000000000000
 26   2.00000000000000000000000000000
 27   2.00000000000000000000000000000

I would worry a little about this accuracy comparison, because the problem could be in the output formatter, too. You probably need to compare computed values to manually-constructed , exact values. These can be generated by assembling a mantissa and exponent with ldexp.

Note that in normal python, the output formatting is almost certainly being done with IEEE754 doubles, upconverted from the 32 bit floats in your array.

Note that the values you expect here should be readable by directly printing the floats as hex (get them as mem32 objects). With the suppressed 1 on IEEE754 formats, you should just see a bit walking across the hex representation, with all other bits zero.

@peterhinch I forget to mention, that for any maths calculations where accuracy is needed, but platform calculation has some limitation I'm always looking for information in any Jet Propulsion Laboratory and articles/studies about maths calculation and implementation in early microprocessors (mostly code is in Fortran, but all steps and algorithms are explained and at usually goal is eliminating calculation errors.

For example: Implementation of the Sun Position Calculation,

Wow, I looked at the paper you cited. It is based on computing using an AM9511 FPU, which is the first FPU I did a lot of work on (circa 1978). I had it interfaced to my Apple II computer, and worked in a language (I helped write) called First, which was a compiled Forth-like language. The Am9511 had some odd foibles, including not using the suppressed '1' high bit in PDP-11 and modern IEEE754 floating point formats.

@peterhinch FYI raspberrypi/pico-sdk#1255

@peterhinch yes, I knew your problem with the actual calculation wasn't the output formatter; it was just the accuracy test that was not meaningful because of that. I strongly suspect your problem is in the range reduction of the trig functions, which the pico sdk warns is not done to the usual accuracy. You may have to find an almanac calculation based on a different epoch, so that there aren't too many cycles of trig to wrap around anywhere, or at least so that any range reduction explicitly built in to the calculation doesn't have too much loss of accuracy. Single precision really doesn't work too well for calendar stuff.

Single precision really doesn't work too well for calendar stuff.

As I said above, the difference between SP and DP is about 20s. Every SP platform I tested returned identical results (within my resolution of 1s) - albeit out by 20s relative to DP. This whether FP is done in hardware or software. The one exception was RP2.

Thanks for that - I will adapt the code. In the vast majority of cases there is only one root and the numbers differ substantially. The tricky case is where the Sun or Moon is grazing the horizon (in arctic latitudes at noon or midnight). In the two hour interval two roots can occur and the numbers can be quite similar.

Thanks for that, I'll study it in detail tomorrow. That addresses one issue, improving the code.

However the second issue is the observation that RP2 seems to have lower precision than other 32-bit platforms.

Following the same way it turned out that the lmst() function introduced the deviations:

    def lmst(self, mjd):
        d = mjd - 51544.5
        t = d / 36525.0
        d_prec = 8731.875
        t_prec = 0.23906570841889116
        d = d_prec
        t = t_prec
        print('d, t:', d, t, 'deviations (%):', 100*(d-d_prec)/d_prec, 100*(t-t_prec)/t_prec)
#        lst = 280.46061837 + 360.98564736629 * d + 0.000387933 * t * t - t * t * t / 38710000
        lst = 280.46061837 + 360.98564736629 * d + 0.000387933 * t * t - 2.5833118057349522e-08* t * t * t
        lst_prec = 3152362.0102370647
#        lst = lst_prec
        res = (lst % 360) / 15.0 + self.long / 15
        res_pico = 13.40981
        lst_prec = 3152362.0102370647
        res_prec = 13.327161695342511
        print('lst, res:', lst, res, 'deviations (%):', 100*(lst-lst_prec)/lst_prec, 100*(res-res_prec)/res_prec)
#        return res_pico
        return res

With output on Pico:
lst, res: 3152363.0 13.40981 deviations (%): 2.379168e-05 0.6201418

0.6 % error, whereas it is only 0.12 % on the stm Blackpill.

Setting starting values d and t to common values we find that the following two lines introduce the deviations between platforms:

d = 8731.875
t = 0.23906571
lst = 280.46061837 + 360.98564736629 * d + 0.000387933 * t * t - 2.5833118057349522e-08* t * t * t
res = (lst % 360) / 15.0 + -2.102811634540558 / 15
res_prec = 13.327161695342511
print('res:', res, 'deviation (%):', 100*(res-res_prec)/res_prec)

Blackpill (stm): res: 13.34315 deviation (%): 0.1199252
Pico: res: 13.40981 deviation (%): 0.6201418
Further processing increases the error about 10-fold. But this is not platform dependent any more. I I feed the wrong res value into the further Python (64 bit double) computation it yields the same wrong value.

Same finding here :-)
Just noted that we have a (3152363.0 % 360) computation following which dramatically decreases precision.

I was suspicious that this could be a contribution, too. 4 LSB is a pretty big error.

I think this warrants an issue.

I'll raise an issue. The origins of the deviation of esp8266 and RP2 to stm boards are obscure to me.
I didn't find an str2float function in the RP2 SDK.

@peterhinch It may well have worked much better when it was written, probably close to mjd 51544.5, since the 't' parameter would have been very close to zero at that point. That would eliminate the roundoff error in the modular reduction. Those equations should probably be updated to a more recent epoch.

I changed in get_mjd():
mjepoch = -10957.5 # 40587 - 51544.5 # Modified Julian date of C epoch (1 Jan 70)

and in sin_alt() accordingly:
mjd1 = self.mjd + hour / 24.0

but that does not improve any further. Still I think it's a reasonable change.

I'll implement those.

I've incorporated your improvements to date in the actual code. The RP2040 and other 32-bit platforms produce rise and set times within 30s of the 64-bit results. This is (in my opinion) good enough for any likely microcontroller application. There is always the option of using a soft or hard 64-bit FPU.

On further thought I'm unsure about the change to MJD. The modified Julian date is a known concept with wider significance: changing its value could cause confusion if anyone were to extend the module.

Thank you Steward. That was the library I remembered but couldn't find now. Still I see the str2float mentioned there (as an optinal external add-on) not included in the SDK.

The issue seems to be the pico-sdk's implementation of powf... it doesn't handle 10**-11 correctly.

raspberrypi/pico-sdk#1566

@dpgeorge Should the string-to-float conversion be modified to use an ipow routine (which is only a couple lines long, if it needs to be written out) instead of pow/powf? I suspect it would be both faster and more accurate.

No, I don't think that's necessary. The point is that powf is broken and needs to be fixed (and is now fixed in MicroPython's rp2 port), regardless of whether it's used for string-to-float conversion.

And because we try to optimise for size instead of speed (for non-critical things), it's best to reuse powf instead of adding a new ipow.

@dpgeorge OK, no argument on the appropriateness of fixing powf. ipow can be really tiny, but FP conversion speed probably isn't critical enough to be worth much space.

You've lost me there for two reasons. Firstly the testing that revealed the error was done with lto==0. Secondly the round() function with default second arg yields an integer (number of decimal places is zero):

>>> round(1.5*3600)
5400
>>> round(1.75)
2
>>> round(1.33)
1

It could be that the simplification df = frac(0.5+h/24) in lstt() no longer is valid. I did it to introduce the un-reduced h/24 with full precision (role 2) into df. The 0.5 is the value that occurred always in the tests for frac(self.mjd - 51544.5). Could it be that it with is_up() has different values? Can these be re-introduced without making the error worse again?

@peterhinch No, round() with no second argument returns a float, whose fractional part is zero! That is not an int. And then, when that float (even zero) is added to the time.time(), the whole thing gets upconverted to float. Try: type(round(1.5))

@robert-hh @peterhinch OK, Robert, you're right. I am embarrassed. I thought round(1.5,0) would return the same as round(1.5), and just did a shortcut typing it. CPython does the same thing (real int if no sec ond argument), which I had never noticed.

Nice you appreciate it.