Summation on HP 42S

[lrdheat]

Wow! I love seeing the math employed in contrast to just doing a simple, mindless, addition loop to get to the answer.

The Free42 handles 1000! Without out of range error. The WP 34S handles it (high is displayed in the power of 10 display, but asking for a full display with left pointing yellow carrot symbol shows the exponent of 10 along with the first 34 digits of 1000! blue right pointing carrot gives final portion of 34 digits).

[Albert Chan]

Hi, Frido Bohn

For reference, 1000! ~ 4.023872600770937735437024339230039857194E+2567

To avoid overflow, sum log10 of the integers.
A straight sum will loses a few significant digits.

>>> from math import *
>>> n = 1001
>>> frac = lambda x: x - int(x)
>>> sum(log10(i) for i in xrange(2, n))
2567.6046442221304
>>> 10 ** frac(_)
4.0238726007486756

If error terms is collected, result is 1000X more accurate

>>> sum = error = 0
>>> for i in xrange(2, n):
...         v = log10(i)
...         old = sum
...         sum += v
...         error += v - (sum - old) ### collect what is not added to sum
...
>>> sum, error
(2567.6046442221304, 2.4005242238445135e-12)
>>> 10 ** (frac(sum) + error)
4.0238726007709165

Update: we could collect powers of 10, without log10()

lua> p, e = 1, 0
lua> for x = 1, 500, 2 do
        p = x * (x+1) * (1000-x) * (1001-x) * p
        if p > 1e200 then e=e+300; p=p*1e-300 end
        end
lua>
lua> p, '* 10^' .. e
4.0238726007709384e+167 * 10^2400

[lrdheat]

Interestingly, the Free42's last 5 digits of the 34 digits (30024) differs from the WP 34S's last five (26767).

[Albert Chan]

Challenge

What is formula for sum of 1 + 1/3 + 1/6 + 1/10 + ... + 1/(1 + 2 + 3 + ... + n) ?

[Valentin Albillo]

(09-25-2018 07:17 PM)Albert Chan Wrote:  Challenge

What is formula for sum of 1 + 1/3 + 1/6 + 1/10 + ... + 1/(1 + 2 + 3 + ... + n) ?

Too easy, Albert ...

V.
.

[Albert Chan]

(09-25-2018 08:01 PM)Valentin Albillo Wrote:  

(09-25-2018 07:17 PM)Albert Chan Wrote:  What is formula for sum of 1 + 1/3 + 1/6 + 1/10 + ... + 1/(1 + 2 + 3 + ... + n) ?

Too easy, Albert ...

It was actually the Huygen's test (I got that from Mathematical Universe, page 144).
Leibniz, father of calculus, need to pass this to be accepted as Huygen's student.

It were not meant to be hard ...

Edit: here is a harder one, but do not attempt until above can be solve
Find formula for sum of 1/2 + 1/5 + 1/9 + 1/14 + ... + 1/((1 + 2 + 3 + ... + n) + n)

[Frido Bohn]

(09-25-2018 04:47 PM)Dieter Wrote:  Simple. Add the (base 10) logs of 1...1000. The integer part of the sum is the tens exponent, 10^(fractional part) is the mantissa.
[...]
Dieter

Your program is really elegant as it uses the Y-register for the calculation of the logarithm as well as for decrementing n and counting in one pass.
Mine is a bit more clumsy as it uses the X-register for this purpose. The summation of the logarithm is then achieved through stack shifting and the use of LASTX.

Code:

01 LBL "FC2"
02 0
03 X<>Y
04 LBL 01
05 LOG
06 LASTX
07 Rv
08 +
09 R^
10 DSE ST X
11 GTO 01
12 X<>Y
13 IP
14 LASTX
15 FP
16 10^X
17 END

I wouldn't have made the effort to list the code if it was not for the timing of the programs. Interestingly, my code takes for 1000! about 90 seconds while the more elegant code takes about 8 seconds longer in a real HP42S.

The result is (real HP42S) 4.02387322958 10^2567 for both codes. That is precise to the 5th decimal place. By using the natural logarithm (LN) instead of the decadic (LOG) the precision is increased by 1 decimal place.

Cheers!
Frido

[Gerson W. Barbosa]

(09-25-2018 08:59 PM)Albert Chan Wrote:  Edit: here is a harder one, but do not attempt until above can be solve
Find formula for sum of 1/2 + 1/5 + 1/9 + 1/14 + ... + 1/((1 + 2 + 3 + ... + n) + n)

As I am kinda lazy, I’ll just ask my HP 50g to do it for me:


'∑(k=1,n,2/(k*(k+3)))'
EVAL
COLLECT

[Albert Chan]

(09-26-2018 12:36 PM)Gerson W. Barbosa Wrote:  As I am kinda lazy, I’ll just ask my HP 50g to do it for me:

'∑(k=1,n,2/(k*(k+3)))'
EVAL
COLLECT

What does it say ? I bet it look complicated ...

Actually, the sum is very easy, once you recognized above formula

∑(k=1,n,2/(k*(k+3)))
= 2/3 * ∑(k=1,n,3/(k*(k+3)))
= 2/3 * ∑(k=1,n, 1/k - 1/(k+3))

If you tried out different values of k, the sum have at most 6 terms. All the rest cancelled out.

∑(k=1,n,2/(k*(k+3))) = 2/3 * (1 + 1/2 + 1/3 - 1/(n+1) - 1/(n+2) - 1/(n+3))

[Gerson W. Barbosa]

(09-26-2018 01:23 PM)Albert Chan Wrote:  

(09-26-2018 12:36 PM)Gerson W. Barbosa Wrote:  ...

'∑(k=1,n,2/(k*(k+3)))'
EVAL
COLLECT

What does it say ? I bet it look complicated ...

Indeed it is:

'(11*n^2+48*n+49)*n/(9*((n+3)*((n+2)*(n+1))))'

Well, that’s the price to pay...

(09-26-2018 01:23 PM)Albert Chan Wrote:  Actually, the sum is very easy, once you recognized above formula

∑(k=1,n,2/(k*(k+3)))
= 2/3 * ∑(k=1,n,3/(k*(k+3)))
= 2/3 * ∑(k=1,n, 1/k - 1/(k+3))

If you tried out different values of k, the sum have at most 6 terms. All the rest cancelled out.

∑(k=1,n,2/(k*(k+3))) = 2/3 * (1 + 1/2 + 1/3 - 1/(n+1) - 1/(n+2) - 1/(n+3))

Or

∑(k=1,n,2/(k*(k+3))) = 2/3 * (11/6 - 1/(n+1) - 1/(n+2) - 1/(n+3))

Very nice!

[Albert Chan]

(09-26-2018 10:10 AM)Frido Bohn Wrote:  The result is (real HP42S) 4.02387322958 10^2567 for both codes. That is precise to the 5th decimal place.
By using the natural logarithm (LN) instead of the decadic (LOG) the precision is increased by 1 decimal place.

To increase precision (speed too!), do 1000! in "shells", 4 numbers at a time
So, only 250 log10() calls instead of full 1000

I use my own fsum.lua code, to collect errors.

Lua> p = require "fsum" () -- goal is log10(1000!)
Lua> shell = function(x) return x * (x+1) * (1000-x) * (1001-x) end
Lua> log10 = math.log10
Lua> for x = 1,500,2 do p:add(log10(shell(x))) end
Lua> = p:total()
2567.6046442221327
Lua> p:add(-2567)
Lua> = 10 ^ p:total()
4.023872600770903

BTW, straight sum of log10 of shells produce mantissa of 4.023872600761316

[pier4r]

(09-25-2018 05:17 PM)lrdheat Wrote:  Wow! I love seeing the math employed in contrast to just doing a simple, mindless, addition loop to get to the answer.

Since years my motto is:

Who understands a bit of mathematics, computes a sum with closed formulas. (obviously not hard to compute sums)
Who doesn't, uses a for loop for the computation.

Obviously I am in the second group. I do for loop and then I benchmark things to be faster as I am lazy.

[Thomas Klemm]

(09-26-2018 04:09 PM)pier4r Wrote:  I do for loop and then I benchmark things to be faster as I am lazy.

Or then you just have the compiler figure out the closed formula for you.

Cheers
Thomas

---

Here's the code for \(\sum_{x=1}^n x^2 - 3 x = \frac{1}{3} (n - 4) n (n + 1)\).

[Ángel Martin]

(09-27-2018 03:47 AM)Thomas Klemm Wrote:  Or then you just have the compiler figure out the closed formula for you.

Interesting talk, the presenter's English is certainly pre-compiled ;-)

[pier4r]

(09-27-2018 03:47 AM)Thomas Klemm Wrote:  Or then you just have the compiler figure out the closed formula for you.

Cheers
Thomas

---

Here's the code for \(\sum_{x=1}^n x^2 - 3 x = \frac{1}{3} (n - 4) n (n + 1)\).

That's impressive

[Thomas Okken]

(09-27-2018 10:37 PM)pier4r Wrote:  That's impressive

It's nice as long as it works...

Free42 is at version 2.0.21a on iOS, and 2.0.21 on all other platforms, because the latest iOS SDK broke part of the number display logic, causing it to display unwanted leading zeroes in certain cases in HEX and OCT modes. This is code that has worked fine on more than half a dozen platforms, compiled with many different compilers, over the course of more than 13 years... and suddenly an optimizer bug in the latest iOS compiler caused it to break.

Don't get me wrong, optimizers are great when they work, but when they don't, they can cause malfunctions in the weirdest places. Keep in mind that in the greater scheme of things, Free42 is a pretty simple piece of software, where bugs tend to be noticed quickly. I shudder to think what obscure optimizer bugs could do in really complex information systems!

The person holding that talk apparently works for a high-speed trading company, where speed is everything. Good for him, but in the places where I have worked, and work today, I would rather disable all optimizations everywhere, and only enable them in those places where they make a significant and predictable difference. Having the compiler surprise me by going Mathematica on me and translating code to a closed-form equivalent strikes me as scary rather than exciting.

[pier4r]

(09-28-2018 12:33 AM)Thomas Okken Wrote:  The person holding that talk apparently works for a high-speed trading company, where speed is everything. Good for him, but in the places where I have worked, and work today, I would rather disable all optimizations everywhere, and only enable them in those places where they make a significant and predictable difference. Having the compiler surprise me by going Mathematica on me and translating code to a closed-form equivalent strikes me as scary rather than exciting.

This is also a good point. Stability and predictability.

[ijabbott]

(09-28-2018 12:33 AM)Thomas Okken Wrote:  Don't get me wrong, optimizers are great when they work, but when they don't, they can cause malfunctions in the weirdest places. Keep in mind that in the greater scheme of things, Free42 is a pretty simple piece of software, where bugs tend to be noticed quickly. I shudder to think what obscure optimizer bugs could do in really complex information systems!

Question is, how do you know the optimizer is at fault? It could be at fault, but more often than not it is invalid assumptions of the programmer that is at fault. (No accusations intended here. I'm just stating the general case.)

[Albert Chan]

(09-28-2018 11:20 PM)ijabbott Wrote:  Question is, how do you know the optimizer is at fault? It could be at fault, but more often than not it is invalid assumptions of the programmer that is at fault. (No accusations intended here. I'm just stating the general case.)

I agreed that many times we don't know the answer.
What assumptions the optimizer were using ? Who knows ...

The fully optimized code might still work.
But, once in a while, it may hit a wrong assumption, causing undefined behavior.
If lucky, the rare cases only affect some rounding errors. Again, you don't know.

With so many unknowns, you are better off without too much optimization.
Example: I will never buy a self-driving car (optimized away the driver)

[Valentin Albillo]

Quote:With so many unknowns, you are better off without too much optimization.

This^

Quote: I will never buy a self-driving car

Never is a long, long time.

Have a nice weekend.
V.
.