Threaded Mode | Linear Mode

Pjwum · 07-02-2018, 12:41 PM

Hi all,

a perfect day to introduce myself. Today I'm celebrating e-Day. At least in the rest of the date format world a look at the calendar shows 2.718!

Euler found there is a limit to infinitesimal exponential growth resulting in (1+1/n)^n equals 2.718.. for large n.

There is beauty (e^iPi = -1) and magic (e = lim n->inf (n^1/n)^Prime(n).

I wondered how large n should be to meet today's event

« 1 2.718 -> n eday
« CLLCD
DO 'n' INCR
INV 1 + n ^
DUP n 1 DISP 2 DISP
eday UNTIL >
END n »
»

A simple SOLVE command should suffice but somehow I love to watch calculators at work.

Patrick

Thomas Klemm · 07-02-2018, 01:57 PM

Code:

                                                                        char

                                                            _3141592654[3141

          ],__3141[3141];_314159[31415],_3141[31415];main(){register char*

      _3_141,*_3_1415, *_3__1415; register int _314,_31415,__31415,*_31,

    _3_14159,__3_1415;*_3141592654=__31415=2,_3141592654[0][_3141592654

   -1]=1[__3141]=5;__3_1415=1;do{_3_14159=_314=0,__31415++;for( _31415

  =0;_31415<(3,14-4)*__31415;_31415++)_31415[_3141]=_314159[_31415]= -

1;_3141[*_314159=_3_14159]=_314;_3_141=_3141592654+__3_1415;_3_1415=

__3_1415    +__3141;for                 (_31415 = 3141-

           __3_1415  ;                  _31415;_31415--

           ,_3_141 ++,                  _3_1415++){_314

           +=_314<<2 ;                  _314<<=1;_314+=

          *_3_1415;_31                  =_314159+_314;

          if(!(*_31+1)                  )* _31 =_314 /

          __31415,_314                  [_3141]=_314 %

          __31415 ;* (                  _3__1415=_3_141

         )+= *_3_1415                    = *_31;while(*

         _3__1415 >=                     31415/3141 ) *

         _3__1415+= -                    10,(*--_3__1415

        )++;_314=_314                    [_3141]; if ( !

        _3_14159 && *                    _3_1415)_3_14159

        =1,__3_1415 =                    3141-_31415;}if(

        _314+(__31415                     >>1)>=__31415 )

        while ( ++ *                      _3_141==3141/314

       )*_3_141--=0                       ;}while(_3_14159

       ) ; { char *                       __3_14= "3.1415";

       write((3,1),                       (--*__3_14,__3_14

       ),(_3_14159                         ++,++_3_14159))+

      3.1415926; }                         for ( _31415 = 1;

     _31415<3141-                          1;_31415++)write(

    31415% 314-(                           3,14),_3141592654[

  _31415    ] +                           "0123456789","314"

  [ 3]+1)-_314;                           puts((*_3141592654=0

,_3141592654))                             ;_314= *"3.141592";}

One of the many ways to calculate e.

Valentin Albillo · 07-02-2018, 06:24 PM

.
Hi, Patrick:

First of all, welcome to the forum. I fully expect you'll enjoy the place and enrich us all with your own contributions. That said:

Quote:Today I'm celebrating e-Day. At least in the rest of the date format world a look at the calendar shows 2.718!

Indeed, it shows in that order where I live. To contribute to your celebration, here's a very short program I wrote many many years ago to compute e to in excess of 200 digits in an HP-15C. Here's an extract taken from my article "Long Live the HP-15C !":

This 64-step program will compute from 8 to 208 decimal places of Euler’s constant, the well-known transcendental number e = 2.71828+ . It is by no means optimized for performance but tries instead to be as short and straightforward as possible. Although you can compute more decimal places in an HP-15C, for the purposes of this article this simpler program will do nicely.

Program listing:

Very Important: steps 30 and 43 must be entered in USER mode, while all the rest must be entered out of USER mode. An u-like character must appear next to the step number only for steps 30 and 43, and no others.

01  LBL A         23  RCL A         45  FRAC
02  MATRIX 0      24  X=0?          46  CHS
03  MATRIX 1      25  ISG 2         47  RTN
04    1           26  LBL 0         48  LBL 5
05  DIM A         27  RCL A         49  FRAC
06  DIM B         28  RCL÷ I        50  RCL B
07  RESULT B      29  INT           51  INT
08  STO 2         30u STO A         52  RCL RAN#
09  STO I         31  GTO 2         53    *
10   EEX          32  RCL I         54    +
11    8           33  PSE           55   R/S
12  STO A         34  RCL MATRIX A  56  GTO 4
13   1/X          35  MATRIX 7      57  LBL 2
14  STO RAN#      36  TEST 0        58  RCL* I
15  LBL 1         37  GTO 1         59    -
16  RCL MATRIX A  38  RCL MATRIX B  60  RCL RAN#
17  RCL MATRIX B  39  RCL RAN#      61    ÷
18    +           40    *           62  STO+ A
19  RCL 2         41  FIX 8         63  RCL A
20  STO O         42  LBL 4         64  GTO 0
21  ISG I         43u RCL B
22  LBL 3         44  GTO 5

The constant e is computed using the well-known formula:

      e = 1 + 1/1! +1/2! + 1/3! + 1/4! + ...

using enough terms to achieve the desired precision. We use the powerful matrix capabilities of the HP-15C, holding the current term in a vector (matrix A), and the running sum in another vector (matrix B). Steps 01-14 dimension and initialize both matrices, as well as some indexes and ancillary constants.

Starting with 1, each term is computed by dividing the previous one by the corresponding divisor up to a maximum of 208 decimal digits, until we arrive at a term that is zero to the specified precision. This multiprecision arithmetic is done by considering each term as composed of N blocks (1-26), each holding 8 digits, the extra digits in each register being carried out to the next block. Since we may use divisors up to 125, each block is limited to 8 digits, lest the carry would make the next block larger than the 10 digits an HP-15C storage register can hold.

The addition of each multi-block term to the running total is done using matrix arithmetic in steps 16-18. Steps 19-25 increment the divisor and keep track of the first non-zero block of each term to optimize speed by avoiding unnecesary operations. Steps 26-31 and 57-64 perform the multiprecision division.

The process ends when the current term has all its blocks equal to zero. As all blocks always hold positive values, we can use an advanced matrix operation, the Row Norm, to test for finalization, because in this case the Row Norm will equal zero if and only if all block values are zero. This is tested in steps 34-37.

Once this condition is met, the result matrix which holds the running sum is scaled down for display using matrix-scalar arithmetic (steps 38-40). The computed answer is then displayed by recalling each block in turn, adding the carry from the previous block, and marking the last one negative (steps 41-56).

As you can see, though simple, this program does use some of the HP-15C’s advanced programming capabilities (as well as a trick or two), including basic matrix operations (MATRIX 1), advanced matrix functions (MATRIX 7), recall arithmetic (RCL÷ I), using registers as indexes including increment-and-test operations (ISG 2) , auto-increment-test-and-loop matrix element access (uSTO A, uRCL B), matrix arithmetic (steps 18 and 40), etc.

Usage instructions:

1) After keying in the program, make sure you’re not in complex mode (press CF 8 if in doubt), then you must commit enough storage registers to the common pool for the matrix operations. To that effect, press:

      2, f DIM (i)

2) Now, enter the number of 8-digit blocks you want to use, from 1 (8 digits) to 26 (208 digits). For example, if you want to compute 200 decimal digits of e, you must specify 200/8 = 25 blocks. Start the program by pressing either

      f PRGM, R/S or GSB A or (in User mode) A

While running, it will briefly show each succesive divisor used (2, 3, ...), then once the computation is over, it will display each block of 8 decimal digits (with an initial “0.”, in order to preserve leading zeros at the left end of the block), starting with decimals 1st-8th. The very last block will be marked negative, to signal the end of the output.

Let’s see an example: to compute the first 24 decimal digits of e, we specify 24/8 = 3 blocks, and proceed like this (in USER mode):

     3, A -> (2.00000000) [divisor = 2]
          ... [after 2’26”]
          -> (25.00000000) [divisor = 25]
          -> 0.71828182 [decimals 1st- 8th]
      R/S -> 0.84590452 [decimals 9th-16th]
      R/S -> -0.35360274 [decimals 17th-24th]

so, after adding a “2.” at the front and writing down all 8-digit blocks (that is, minus the initial “0.” or “-0.”) in their proper order, we finally get:

      e = 2.71828182 84590452 35360274

where, due to the accumulation of rounding errors during the process, the last block of the computed answer comes out as “3536 0274” while correct is “3536 0287”, so we have an error of 13 ulps (units in the last place).

3) Should you need to display the output again, press: GSB 4

For 208 decimal digits, the final result is:

      e = 2.71828182 84590452 35360287 47135266 24977572
            47093699 95957496 69676277 24076630 35354759
            45713821 78525166 42742746 63919320 03059921
            81741359 66290435 72900334 29526059 56307381
            32328627 94349076 32338298 80753195 25101901
            15738281

The last block of the computed answer comes out as “1573 8281” while correct is “1573 8341”, so the error is 60 ulps, and thus after rounding the last two places we’ve got 206 decimals fully correct.

Regards.
V.
.

Pjwum · 07-03-2018, 11:53 AM

Thank you for showing tremendous ideas to approximate e.

The real fun in programming seems to come from how you deal with given technical restrictions. I remember a nice review about the HP-25 published in Creative Computing outlining how these limitations strengthen strategic thinking. Article
I should give up my Prime and go back to HP-25 or TI-57...

Regarding my little program to find n for 2.718 there are interesting forensic effects. In fact, the numerical solver delivers another n. Graphing (1+1/n)^n and zooming around 4822.55 a pronounced sawtooth alternates above and below eday.

Dieter · 07-03-2018, 09:37 PM

(07-03-2018 11:53 AM)Pjwum Wrote: Regarding my little program to find n for 2.718 there are interesting forensic effects. In fact, the numerical solver delivers another n. Graphing (1+1/n)^n and zooming around 4822.55 a pronounced sawtooth alternates above and below eday.

Of course the is no such sawtooth, mathematically. What you see is the result of roundoff errors, especially when calculating 1+1/n. After the addition of 1 only 8 out of 12 digits of the original 1/n value are left. Maybe you can try rewriting the equation as exp[n · ln1+x(1/n)].

BTW the exact result is n = 4821,66555538...

Dieter

Valentin Albillo · 07-03-2018, 11:25 PM

.
Hi again, Patrick:

(07-03-2018 11:53 AM)Pjwum Wrote: Thank you for showing tremendous ideas to approximate e.

You're welcome. I wrote the HP-15C program as part of my series of articles "Long Live the ..." intended as a celebration of classic HP calcs and sometimes particular Modules or ROMs (such as the Advantage module for the HP-41C or the Math ROM for the HP-71B).

The C program posted by Mr. Klemm was obviously taken from the book "Obfuscated C and other mysteries", page 260. The program was written by Roemer B. Lievaart and it won the 1989 Obfuscated C Code Contest prize for Best Layout. It's a very interesting book, I can heartily recommend it.

Quote:The real fun in programming seems to come from how you deal with given technical restrictions. I remember a nice review about the HP-25 published in Creative Computing outlining how these limitations strengthen strategic thinking. Article

Completely agree, and thanks a lot for the link to the HP-25 article, I also agree with its author and really enjoyed my HP-25 (my very first HP calc) back then at the end of the 70's.

If you liked that article, I really think that you'd also like to read my "Long Live the HP-25" article where, among other interesting things, you'll find two HP-25 programs implementing both 3rd-order and 4th-order Runge-Kutta differential equations solvers in less than 49 steps and leaving enough steps to define the differential equation to solve. They make 3 and 4 calls to the equation's definition, respectively, in a machine without subroutines.

That's overcoming technical restrictions for you !

Regards.
V.
.

Pjwum · 07-04-2018, 06:46 PM

Quote:Of course the is no such sawtooth, mathematically. What you see is the result of roundoff errors, especially when calculating 1+1/n. After the addition of 1 only 8 out of 12 digits of the original 1/n value are left. Maybe you can try rewriting the equation as exp[n · ln1+x(1/n)].

BTW the exact result is n = 4821,66555538...

Yes, correct interpretation. The ramp is the exponent, the step the new rounded 1+1/n.
And your result is dead on; today I found n = 4821.66555538178034.. but can't go beyond that.

Quote:I wrote the HP-15C program as part of my series of articles "Long Live the ..." intended as a celebration of classic HP calcs

Cool! Thank you..

Thomas Klemm · 07-04-2018, 07:06 PM

(07-03-2018 11:25 PM)Valentin Albillo Wrote: The C program posted by Mr. Klemm was obviously taken from the book "Obfuscated C and other mysteries", page 260.

Nah, I found it in Pi Unleashed by Arndt, Jörg, Haenel, Christoph on page 36.

Quote:It's a very interesting book, I can heartily recommend it.

Pi Unleashed is another book that we full-heartedly recommend.

Later on page 85 they show this small program that calculates e to 15,000 places:

Code:

/* note: N=15000, LEN=87700 >= 1.4*N*log10(N), 84700=LEN-N/5 */

a[87700],b,c=87700,d,e=1e4,f=1e5,h;

main(){for(;b=c--,b>84700;h=printf("%05d",e+d/f),e=d%=f)

for(;--b;d+=f*(h?a[b]:e),a[b]=d%b,d/=b);}

Both programs use the spigot algorithm to calculate e.

Cheers
Thomas

PS: Here's a link to Roemer B. Lievaart wining contribution of 1989 for the best layout.

Thomas Klemm · 07-06-2018, 09:27 PM

Related programs in the General Software Library: