The Museum of HP Calculators

HP Forum Archive 20

[ Return to Index | Top of Index ]

wp34s program
Message #1 Posted by Gerson W. Barbosa on 14 July 2011, 10:37 a.m.

=======================================

001 LBL A 027 1/X 053 * 002 DEC X 028 3 054 RCL+ 04 003 3 029 RCL+ L 055 STOP 004 ENTER 030 1/X 056 RCL 04 005 SQRT 031 STO+ X 057 RCL* 00 006 RCL* L 032 STO+ X 058 SQRT 007 STOS 00 033 + 059 STO 04 008 4 034 DEC X 060 PSE 10 009 1/x 035 +/- 061 RCL 00 010 STO 05 036 * 062 STOP 011 * 037 ENTER 063 || 012 STO 04 038 RCL* 07 064 STO+ X 013 RCL 02 039 STO- Z 065 STO 00 014 STO 06 040 Rv 066 1 015 SIGN 041 RCL 05 067 RCL- 05 016 RCL* 05 042 STO* 07 068 SQRT 017 STO 07 043 Rv 069 DEC X 018 3 044 DSE 06 070 2 019 1/x 045 BACK 22 071 / 020 STO+ X 046 X<>Y 072 +/- 021 ENTER 047 8 073 STO 05 022 +/- 048 * 074 2 023 RCL 02 049 RCL* 01 075 STO* 01 024 INC X 050 RCL 05 076 DSE 03 025 RCL- 06 051 SQRT 077 BACK 64 026 STO+ X 052 RCL* L 078 RTN

=========================================

The curious might want to key this in and see what it does. Others please just ignore it. More details next week.

      
Re: wp34s program
Message #2 Posted by Neil Hamilton (Ottawa) on 14 July 2011, 11:06 a.m.,
in response to message #1 by Gerson W. Barbosa

... or you might want to be lazy (Virtues of a programmer) and use the assembler we have provided. :-)

Simply cut and paste the translation below into a file I'll call "barbosa.wp34s".

// George W. Barbosa's mystery program
// http://www.hpmuseum.org/cgi-sys/cgiwrap/hpmuseum/forum.cgi?read=187663#187663
//
// $ wp34s_asm.pl barbosa.wp34s -o wp34s-0.dat
//
// Emulator will now just read it in directly (provided you are in the right directory).
// Or you can load it into the calculator as well.
//
001 LBL'MYS' // Changed label so it could be more easily found via "h CAT".
002 DEC X
003 3
004 ENTER[^]
005 [sqrt]
006 RCL[times] L
007 STOS 00
008 4
009 1/x
010 STO 05
011 [times]
012 STO 04
013 RCL 02
014 STO 06
015 SIGN
016 RCL[times] 05
017 STO 07
018 3
019 1/x
020 STO+ X
021 ENTER[^]
022 +/-
023 RCL 02
024 INC X
025 RCL- 06
026 STO+ X
027 1/x
028 3
029 RCL+ L
030 1/x
031 STO+ X
032 STO+ X
033 +
034 DEC X
035 +/-
036 [times]
037 ENTER[^]
038 RCL[times] 07
039 STO- Z
040 R[v]
041 RCL 05
042 STO[times] 07
043 R[v]
044 DSE 06
045 BACK 22
046 x[<->] Y
047 8
048 [times]
049 RCL[times] 01
050 RCL 05
051 [sqrt]
052 RCL[times] L
053 [times]
054 RCL+ 04
055 STOP
056 RCL 04
057 RCL[times] 00
058 [sqrt]
059 STO 04
060 PSE 10
061 RCL 00
062 STOP
063 ||
064 STO+ X
065 STO 00
066 1
067 RCL- 05
068 [sqrt]
069 DEC X
070 2
071 /
072 +/-
073 STO 05
074 2
075 STO[times] 01
076 DSE 03
077 BACK 64
078 RTN
// CRC16: 5ACD
// Total words: 79
// Total steps: 78
            
Re: wp34s program
Message #3 Posted by Gerson W. Barbosa on 14 July 2011, 12:27 p.m.,
in response to message #2 by Neil Hamilton (Ottawa)

Nice utility! I've downloaded wp34s_asm.exe and wp34s_asm.pl to the wp34s emulator folder. What should I do next? (I don't know Perl). Thanks!

Edited: 14 July 2011, 12:32 p.m.

                  
Re: wp34s program
Message #4 Posted by Neil Hamilton (Ottawa) on 14 July 2011, 1:32 p.m.,
in response to message #3 by Gerson W. Barbosa

Hi Gerson,

There is a manual in the SVN's "./doc" directory: WP34s_Assembler_Disassembler.pdf.

Basically, you don't need to know any Perl -- even if you use the Perl version. The EXE is the same thing (just SLOWER -- due to its requirement to dynamically decompress [at least the first time], and BIGGER -- same reason). Again, see the document.

To assemble an ASCII source file into a binary image (I'll show the EXE version but the Perl script is the same, just a different extension):

$ wp34s_asm.exe YOUR_ASCII_SOURCE.wp34s -o BINARY_IMAGE.dat

To disassemble a binary image back into an ASCII source file:

$ wp34s_asm.exe -dis BINARY_IMAGE.dat > RECOVERED_SOURCE.wp34s

You don't need to use the suggested file extensions but I would encourage these ones to aid in standardization. The ".dat" are eaten directly by the emulator if called "wp34s-0.dat", "wp34s-1.dat", or "wp34s-2.dat". (As of a few days ago, -3 disappeared due to flash space limitations. Can't recall the exact SVN version of this change.)

There are a bunch of options that you can use as documented in the manual.

Paul had a few cool ideas that we are trying to capture this week (or maybe the next) so it should gain some interesting new functionality soon. Most are dealing with symbolic branch targets -- you won't have to count backwards or forwards N-steps, and if you later insert a few new steps in between, you won't have to "fix up" the BACK/SKIP instructions -- these will be automatically handled. A simple example:

Loop::       RCL+ 00
             STO/ 01
             x[<->] Y
             BACK Loop
             STO 02

If you can afford the disk space, having a proper Perl package installed will result in higher performance (no decompression wait) and the doc explains about where you can get these.

I hope you find it useful! Let us know how you make out.

BTW: All the programs in the SVN "./library" directory have a ".wp34s" counterpart that can be directly assembled by this tool (assuming I haven't made any transcription errors :-). These can be used as examples for writing your own compatible code. We know the syntax is a bit weird and unforgiving. To help, the tool can write an opcode syntax guide using the "-syntax FILENAME" switch. This should help a bit.

Cheers...

                        
Re: wp34s program
Message #5 Posted by Gerson W. Barbosa on 14 July 2011, 5:52 p.m.,
in response to message #4 by Neil Hamilton (Ottawa)

Hello Neil,

Thanks! I've installed the latest version of Strawberry Perl:

>wp34s_asm.pl mpgm.wp34s -o wp34s-0.dat

// CRC16: 5ACD // Total words: 79 // Total steps: 78

>wp34s_asm.pl -dis wp34s-0.dat > test.wp34s

This reverses to the original listing, as expected, but I don't know how to properly load the .dat file into the emulator. I tried the LOAD command, but this is what I got:

001 LBL'MYS'
002 DEC X
003 3
004 ENTER[^]
005 ???
006 ???
007 ???
...

Most likely I've done something wrong. I'll read through the manual later. For the time being I am not interesting in transferring programs to the real calculator, as I don't have a suitable USB-to-RS232 converter yet.

Regards,

Gerson.

                              
Re: wp34s program
Message #6 Posted by Paul Dale on 14 July 2011, 8:46 p.m.,
in response to message #5 by Gerson W. Barbosa

Firmware version mismatch between the assembler and the calculator?

- Pauli

                                    
Re: wp34s program
Message #7 Posted by Gerson W. Barbosa on 14 July 2011, 9:00 p.m.,
in response to message #6 by Paul Dale

I'll check that, thanks! Emulator build date is 2011-7-14 (0).

Gerson.

                                          
Re: wp34s program
Message #8 Posted by Paul Dale on 14 July 2011, 9:23 p.m.,
in response to message #7 by Gerson W. Barbosa

The current emulator is broken too :-(

- Pauli

                              
Re: wp34s program
Message #9 Posted by Neil Hamilton (Ottawa) on 14 July 2011, 9:25 p.m.,
in response to message #5 by Gerson W. Barbosa

Hi Gerson,

No, your commands are correct. The fact that you got the same CRC16 as I did means that the assembly did work. Mail me your DAT file and I'll see if I can figure out what might have occurred. You can mail it to the account attached to this message.

As Paul pointed out, firmware versions are VERY important as the op-code numbers sometimes migrate with the ever-improving state of the development. The assembler can handle this IFF you have the correct op-code tables. (There is info in the manual about this but it can be tricky to create the alternate tables.) However, looking at your message, the first op-code that is wacky is the square root (step 005) and it hasn't changed in a while (0x020a). The version numbers we need to track this are the ones reported by "h X.FCN|3|ENTER" (in other words, the VERS command) from the calculator and the SVN number of the assembler. Tack the -svn switch onto the command line and tell me what number it reports. It is likely:

$ wp34s_asm.pl  library/gerson.wp34s -o wp34s-0.dat -svn
// Opcode SVN version: 1210
// CRC16: 5ACD
// Total words: 79
// Total steps: 78

Just for interest, if you have an MD-5 or SHA1 tool (look on Wiki), the hashes for the DAT file should be:

$ sha1sum wp34s-0.dat
7b5213cf6e9176c282cc1933f578891a89b4bbf6  wp34s-0.dat
$ md5sum wp34s-0.dat
9f875de2bc41422e3c1a639c84103cfc  wp34s-0.dat

If you have such a hash tool, let me know what value(s) you got. Either one will do -- no need for both.

If you make sure that the "wp34s-0.dat" file is present in the same directory you started the emulator from, the emulator will automatically read that specific file into emulated flash page 0 (wp34s-1.dat would be read into flash page 1, and -2 into page 2).

Though the current Assembler documentation mentions flash page 3, the latest SVN versions of the calculator have run out flash space -- so only 0-2 for now. :-(

What directory the emulator is starting in depends on how you start it. In Windows, you can create links to programs so you can start them from a shortcut. I am pretty sure you can edit the properties of that shortcut to tell it what directory to start in. (I am running Linux and don't recall the exact sequence off the top of my head.) Make sure that the DAT file is in that directory. (Going from memory here so YMMV.)

Once you start the emulator, because I switched your label from A to 'MYS', the program will be visible with the calculator's "h CAT" function. Using "h CAT", search for MYS (since you have used the flash page 0 file name, wp34s-0.dat, the CAT function should show it as being in "PG 00") and hit XEQ to run your program directly from emulated flash. There are other ways to do this but this is probably the easiest.

I see Bruce Bergman has kindly built a web front end for the tool here. You could give this a try as well.

I am still working kinks out but I have also run a fair bit of code through on both Windows and Linux with few issues (other than SVN bumps *evolving* the internal opcode numbers!! :-). For example, I had a Great Circle program that, as a last step, converted nautical miles to km but after one SVN bump of the calculator, it was converting astronomical units to km! It changed the distance from Boston to Vancouver to something out of this world -- literally!

This is why the assembler can read in opcode definitions from virtually any SVN revision after 1133 (or some such number).

                                    
Re: wp34s program
Message #10 Posted by Gerson W. Barbosa on 14 July 2011, 10:03 p.m.,
in response to message #9 by Neil Hamilton (Ottawa)

Hello Neil,

Here are what I have:

>wp34s_asm.pl mpgm.wp34s -o wp34s-0.dat -svn

// Opcode SVN version: 1199 // CRC16: 5ACD // Total words: 79 // Total steps: 78

VERS on the emulator returns only

34S 2.0
PAULI, WALTEr

no build date (like my hardware wp34s which was flashed two or three weeks ago). It is Marcus's stack lift fix as of this afternoon. I replaced only the files wp34s.exe and wp34sgui.exe inside the emulator folder.

Thanks for your concern.

Gerson.

                                          
Re: wp34s program
Message #11 Posted by Neil Hamilton (Ottawa) on 14 July 2011, 10:31 p.m.,
in response to message #10 by Gerson W. Barbosa

Hi Gerson,

It looks like your emulator might be broken as Paul suggested. Mine is the latest on the tree and VERS is showing "34s 2.0 1220" (meaning SVN 1220). The actual complete SVN file set is now at 1222.

The (lack of) SVN number your emulator is showing is mighty peculiar! I would be leery of just updating a few files since this thing is typically built as a set and there may be things in other files that need to be there as well (DLLs, etc.)

My emulator reads in the flash image of your program just fine and if I execute a "PRCL 0", to swap the flash into RAM, I can scroll through your program with no trouble. (I don't know how to *run* the program -- do I seed it with something in X? -- but I can see it. :-)

If I were you I would do a clean SVN checkout to a new directory and try it all again. (I am *not* an SVN expert so there may be cleaner ways of getting a clean checkout but I don't know how to do that. I just use "svn up".)

Let me know what happens.

Cheers...

                                                
Re: wp34s program
Message #12 Posted by Gerson W. Barbosa on 14 July 2011, 11:36 p.m.,
in response to message #11 by Neil Hamilton (Ottawa)

Hello Neil,

I am not familiarized with SVN. Previously I'd dowload the whole build when is was available. I'll try to fix things tomorrow. I realize the assembler is a great tool, so it's worth to learn how to use it. Thanks again!

Regarding the program, it takes two arguments from the stack (see Jim Horn's post just below). The equivalent C++ program is 33 lines long (I'll post it later).

Cheers,

Gerson.

Edited: 14 July 2011, 11:38 p.m.

                                                      
Re: wp34s program
Message #13 Posted by Neil Hamilton (Ottawa) on 15 July 2011, 5:31 a.m.,
in response to message #12 by Gerson W. Barbosa

Quote:
I am not familiarized with SVN. Previously I'd dowload the whole build when is was available. I'll try to fix things tomorrow. I realize the assembler is a great tool, so it's worth to learn how to use it. Thanks again!

Ahh! I think the root of the problem is becoming clear!

Many of us obsessive ones are running on the bleeding edge of the techo-trinity's every twitch and whim with this project. When I look in the latest ZIP bundle, I see that the newest file puts it at about SVN 1136. We are currently on SVN 1223. The emulator in the bundle is hopelessly out-of-date with the current assembler. Unfortunately, you will likely have to be up near the latest to make use of the assembler (recall that opcodes "evolve") -- at least for now. These targets are moving fairly rapidly. The assembler was not reliable until about SVN 1161 -- and that one was set up to use that vintage of op-code map.

[WARNING: Using the bundled version will also be a problem with the assembler that Bruce Bergman is hosting as well. His CGI script is simply running the very same assembler from SVN and it is still assuming you are at or near the latest SVN build!]

Not to worry! If you want to pursue the assembler route (and get a more modern version of the calculator!), the best thing you can do is get an SVN client. If you are on Windows, check out TortoiseSVN -- it is quite easy to use and hooks right into Explorer. Use the download address from here. TortoiseSVN needs the https address and a clean directory (and an optional build number -- this is where you want to go).

There are some very stable builds that should get you going and I am sure we can give you some "landing points" for those builds. Normally you would right click on TortoiseSVN to tell it to get the latest -- but in your case you want to tell to get revision XYZ (TBD).

Basically everything was incredibly stable just before the latest series of serial functions went in. If you are still game to try this, I will find a good SVN number and let you know. SVN is not that daunting, and has the wonderful advantage that, if you don't like the version you are running, you just tell it to get another specific version. You can easily backward or forward.

Note that you do get quite a few files with the SVN build but I am assuming that a 140MB and a slightly more involved set of paths won't deter you! :-)

                                                            
Re: wp34s stability
Message #14 Posted by Marcus von Cube, Germany on 15 July 2011, 12:49 p.m.,
in response to message #13 by Neil Hamilton (Ottawa)

Quote:
Basically everything was incredibly stable just before the latest series of serial functions went in.

We had some issues with stack lift and y^ lately which have nothing to do with the serial I/O stuff. The problem is flash space. Fixing bugs and adding features at the same time may lead to some inconsistencies. This will definitely be fixed!

There are sometimes SVN updates which simply don't work but are necessary for Pauli and me to exchange the source files. We are working on the same files sometimes and need to exchange them often. That's what SVN is all about.

                                                                  
Re: wp34s stability
Message #15 Posted by Neil Hamilton (Ottawa) on 15 July 2011, 1:09 p.m.,
in response to message #14 by Marcus von Cube, Germany

Quote:
There are sometimes SVN updates which simply don't work but are necessary for Pauli and me to exchange the source files. We are working on the same files sometimes and need to exchange them often. That's what SVN is all about.

Sorry, Marcus. You said it better than I. It was not a criticism at all. I was just pointing out that there are "islands of stability" if Gerson wanted to give it another try. The last stable island was just before the serial effort began... not that it was necessarily related.

We understand that you have to break a few eggs to make an omelette. :-)

                                          
Re: wp34s program
Message #16 Posted by Paul Dale on 15 July 2011, 12:07 a.m.,
in response to message #10 by Gerson W. Barbosa

No version results from using the web interface to download the subversion files. You need a proper subversion client.

- Pauli

      
Re: wp34s program
Message #17 Posted by Jim Horn on 14 July 2011, 6:52 p.m.,
in response to message #1 by Gerson W. Barbosa

Interesting - it appears to compute pi with the number of iterations entered in X and the follow-up iterations in Y or such. When run, it returns the value found. Pressing R/S repeatedly shows number pairs that straddle pi as found. Anything above 14 in X gives pi to the resolution of the

I have to note that RPN programming with all the WP-34S features is remarkably powerful and remarkably write-only. It's not obvious to me yet what algorithm you're using. But it is fun to dig into...

            
Re: wp34s program
Message #18 Posted by Gerson W. Barbosa on 15 July 2011, 11:34 a.m.,
in response to message #17 by Jim Horn

Quote:
it appears to compute pi with the number of iterations entered in X and the follow-up iterations in Y or such.

Quite exactly!

Quote:
I have to note that RPN programming with all the WP-34S features is remarkably powerful and remarkably write-only.

I haven't read all through the manual yet, so it is possible the program can be shrunk even more. Minor improvements can always be made. Before posting, for instance, I managed to make it one step shorter so that the three columns were the same height. The following is much easier to read and will compile in Dev-C++ :

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

int main(int argc, char *argv[]) { int i, j, mi, mj, n = 3; double a, c, p, h, g, s, x; scanf("%d %d", &mi,&mj); h = 3*sqrt(3); g = h/4; x = 1/4.; for (i = 1; i <= mi; i++) { p = x; a = 2/3.; c = -a; for (j = 1; j < mj;) { c *= 1-1/(2.*j)-4/(2.*j++ +3); a -= c*p; p *= x; } s = g+8*n*x*sqrt(x)*a; g = sqrt(g*h); printf("%12d %.16f %.16f %.16f\n",n,g,h,s); h *= 2*g/(g+h); x = (1-sqrt(1-x))/2; n += n; } printf("\n"); system("pause"); return 0; }

            
Re: wp34s program
Message #19 Posted by Paul Dale on 15 July 2011, 7:25 p.m.,
in response to message #17 by Jim Horn

I'd figured out the PI but not the followup portion of the program...

Quote:
I have to note that RPN programming with all the WP-34S features is remarkably powerful and remarkably write-only. It's not obvious to me yet what algorithm you're using. But it is fun to dig into...

I think this is more a function of RPN programming than the WP 34S feature set :-)

- Pauli

                  
Re: wp34s program
Message #20 Posted by Gerson W. Barbosa on 15 July 2011, 10:04 p.m.,
in response to message #19 by Paul Dale

Quote:
I'd figured out the PI but not the followup portion of the program...
1 ENTER 20 A  ->  3.141592653589793 

That's Newton's approximation of pi. In 1666 he obtained 3.1415926535897928 carrying out his calculations to 17 significant digits and evaluating the first 20 terms of his series.

6 ENTER 2 A     ->  3.14903810568          Same as above, with only the first two terms of the series
          R/S   ->  2.59807621135          perimeter of the inscribed triangle (diameter of the circumference = 1)
                    5.19615242271          perimeter of the circumscribed triangle                   
          R/S   ->  3.14172962149          sum of the area of he inscribed hexagon and the area of the sectors comprised between, the latter 
                                           ones approximated with the first two terms of Newton's series
          R/S   ->  3                      perimeter of the inscribed hexagon
                    3.46410161514          perimeter of circumscribed hexagon
          .
          .                                and so on, until the 96-gon (3*2^6)
          .

R/S -> 3.14159265360 area of the 96-gon and the approximated area of the inner sectors R/S -> 3.14103195089 Archimedes' lower bound R/S -> 3.14271459965 Archimedes' upper bound R/S -> 2 end of iterations

The number of correct digits yielded by the approximation of the inner sectors through the Newton's series appear to be

d = (n + 1)*lg(4^a) + 1.5

where

n = number of terms of the Newton's series
a = number of duplications of the sides of the polygons, starting with triangles
lg: base-10 logarithm

Thus, if the calculations are carried out with 1000+ digits, we should expect 1001 correct digits when n = 82 and a = 20 ( 83*lgt(4^20) + 1.5 ). The drawback is the extraction of 40 square roots which might take some time.

I've come up with a variation of Archimedes' method that will produce about 1.8 digits per iteration. Unlike Newton's method, no Calculus needed, but it will require a cube root extraction in the last iteration (van Ceulen could have used it to compute pi to 70 digits in the 15th century :-). Here are some results obtained on Free42 Decimal:

       6  3.141460911773498978633978   
      12  3.141590947963119992966256
      24  3.141592628123496238671042
      48  3.141592653196343012723174
      96  3.141592653583662847516782
     192  3.141592653589697518409744
     384  3.141592653589791743099584
     768  3.141592653589793215098622

Regards,

Gerson.

Edited: 15 July 2011, 11:12 p.m.

                        
Re: wp34s program
Message #21 Posted by Gerson W. Barbosa on 16 July 2011, 11:27 p.m.,
in response to message #20 by Gerson W. Barbosa

Quote:
I've come up with a variation of Archimedes' method that will produce about 1.8 digits per iteration.

It is essentially the Archimedes' method, only an attempt has been made to improve the approximations. (He did it this way).

Using trigonometry, the Archimedes' lower and upper bounds for pi are, respectively:

a = x*sin(pi/x)

b = x*tan(pi/x)

where x=3*2^k, k=1, 2,...

Let's consider the following limit:

        x*tan(pi/x) - pi
  lim   ---------------- = 2
x->inf  pi - x*sin(pi/x)

then the following weighted arithmetic mean of the bounds gives successively better approximations to pi at each iteration:

             2*a + b
   wam   =  ---------       ~   pi
                3
Likewise, the following weighted geometric mean appears to hold:
   wgm   =  (a^2*b)^(1/3)   ~   pi   

also

         wam(x) - pi     9
  lim   ------------- = ---
x->inf   wgm(x) - pi     4

this leads to

         9*wgm(x) - 4*wam(x)
   pi = ---------------------,  when x equals infinity
                  5
or
         27*(a^2*b)^(1/3) - 4*(2*a + b)   
   pi ~ --------------------------------
                       15

The demonstration may not be rigorously correct, nevertheless we can use the result to write the following wp34s program:

001 LBL 77
002 STO 00         
003 3        
004 SQRT        
005 STO+ X       
006 FILL            
007 4       
008 RCL/ L                   
009 /                                     
010 *                                
011 SQRT                                     
012 DSE 00
013 SKIP 01   
014 SKIP 06
015 STO 01
016 || 
017 STO+ X                     
018 FILL
019 RCL 01
020 BACK 10   
021 RCL* X
022 RCL L
023 Rv
024 *
025 CUBERT
026 2
027 7
028 *
029 RCL Z
030 RCL+ X
031 RCL+ Z
032 4
033 *
034 -
035 1
036 5
037 /
038 RTN

1 XEQ 77 -> 3.141460911773500 2 XEQ 77 -> 3.141590947963119 3 XEQ 77 -> 3.141592628123497 4 XEQ 77 -> 3.141592653196343 5 XEQ 77 -> 3.141592653583662 6 XEQ 77 -> 3.141592653589791

An equivalent C++ program is:
#include <stdio.h>
#include <stdlib.h>
#include <math.h>

int main(int argc, char *argv[]) { char i, m; int n = 6; double gm, hm, wam, wgm; // start with inscribed and circumscribed scanf("%d", &m); // hexagons (diameter of the circumference = 1) system("cls"); // hm: harmonic mean hm = 2*sqrt(3); // gm: geometric mean gm = 3*hm/4; // 2*sqrt(3): perimeter of the circumscribed printf("\ n-gon\t wgm\t wam\t\t(9*wgm-4*wam)/5 \n\n"); // hexagon for (i = 1; i <= m; i++) // gm(3/4*2*sqrt(3),2*sqrt(3)) = 3, this is the { // perimeter of the inscribed hexagon gm = sqrt(gm*hm); // wam: weighted arithmetic mean of hm and gm wam = (2*gm+hm)/3; // wgm: weighted geometric mean of hm and gm wgm = pow(gm*gm*hm,1/3.); // weights of gm and hm: 2 and 1, respectively printf("%12d %.16f %.16f %.16f\n",n,wgm,wam,(9*wgm-4*wam)/5); // the weighted arithmetic mean of wgm and wam hm *= 2*gm/(gm+hm); // (weights 9 and -4, respectively) will approximate n += n; // pi at a rate of 1.8 correct digits per iteration. } printf("\n"); system("pause"); return 0; }

Output for nine iterations:
       n-gon          wgm                 wam           (9*wgm-4*wam)/5

6 3.1473451902649443 3.1547005383792515 3.1414609117734984 12 3.1419279179993587 3.1423491305446571 3.1415909479631199 24 3.1416132628330500 3.1416390562199923 3.1415926281234960 48 3.1415939364016969 3.1415955404083897 3.1415926531963425 96 3.1415927336837215 3.1415928338087955 3.1415926535836625 192 3.1415926585943867 3.1415926648502488 3.1415926535896972 384 3.1415926539025598 3.1415926542935209 3.1415926535897909 768 3.1415926536093401 3.1415926536337748 3.1415926535897922 1536 3.1415926535910144 3.1415926535925416 3.1415926535897927

=======================================

Update:

9*a*c 7*(a - 4*c) 10*pi ~ --------- - ----------- 4*a - c 3

where

a = lower bound c = next lower bound

avoids the cube root exraction and gives out the same accuracy (1.8 decimal places per iteration).

=======================================

Update:

Found a relationship between the previous formulas and had W|A solve the equation (solve (p-(27*(a^2*b)^(1/3)-4*(2*a+b))/15 )/((9*a*c/(4*a-c)-7/3*(a-4*c))/10-p)==128/19 for p) and simplify the result:

p ~ (513*(a^2*b)^(1/3) + 24*a*(-(288*a)/(c - 4*a) - 97) - 76*b + 1792*c)/2205

where

a = lower bound b = upper bound c = the next lower bound

This gives 2.4 decimal digits per iteration.

Archimedes' 96-gon is enough for 15 correct digits:

3.14159265358979259

The following procedure can be used to check the formula without running all the iterations:

n:=50;
a:=3*2^n*sin(pi/(3*2^n));
b:=3*2^n*tan(pi/(3*2^n));
c:=3*2^(n+1)*sin(pi/(3*2^(n+1)));
p:=(513*(a^2*b)^(1/3)+24*a*(-(288*a)/(c-4*a)-97)-76*b+1792*c)/2205;

More details here.

Edited: 22 July 2011, 4:58 p.m. after one or more responses were posted

                              
Re: wp34s program
Message #22 Posted by Gerson W. Barbosa on 20 July 2011, 1:47 a.m.,
in response to message #21 by Gerson W. Barbosa

FWIW, the TurboBCD program below is equivalent to the 38-step wp34s program above, except for the computation of p inside the loop. Sqrt(x) and Cbrt(x) are necessary because Turbo Pascal BCD lacks these functions. Notice that only eight iterations suffice for an approximation of pi to 17 digits (as the number of iterations grows, the algorithm yields 1.8 digits per iteration).

Program Quick_Archimedes;

var a, b, p: real; i, m: byte; n: integer;

function Sqrt(x: Real): real;

var s, t: Real;

begin if x<>0 then begin s:=x/2; repeat t:=s; s:=(s+x/s)/2 until s=t; Sqrt:=s end else Sqrt:=0 end;

function Cbrt(x: Real): real;

var c, t: Real;

begin if x<>0 then begin c:=x/2; repeat t:=c; c:=(2*c+x/(c*c))/3 until Abs(c-t)<1e-16; Cbrt:=c end else Cbrt:=0 end;

begin ReadLn(m); n:=6; b:=2*Sqrt(3); a:=3*b/4; ClrScr; Writeln(' n-gon',' ':11,'a',' ':21,'b',' ':21,'p',#10); for i:=1 to m-1 do begin a:=Sqrt(a*b); p:=(27*Cbrt(a*a*b)-4*(2*a+b))/15; WriteLn(n:5,a:22:17,b:22:17,p:22:17); b:=2*a*b/(a+b); n:=2*n end; a:=Sqrt(a*b); p:=(27*Cbrt(a*a*b)-4*(2*a+b))/15; WriteLn(n:5,a:22:17,b:22:17,p:22:17) end.

Output for eight iterations:

n-gon a b p

6 3.00000000000000001 3.46410161513775460 3.14146091177349899 12 3.10582854123024916 3.21539030917347249 3.14159094796312001 24 3.13262861328123821 3.15965994209750050 3.14159262812349626 48 3.13935020304686722 3.14608621513143498 3.14159265319634303 96 3.14103195089050965 3.14271459964536831 3.14159265358366287 192 3.14145247228546209 3.14187304997982388 3.14159265358969753 384 3.14155760791185767 3.14166274705684855 3.14159265358979177 768 3.14158389214831843 3.14161017660468956 3.14159265358979324

=======================================

Update:

TurboBCD program using the new formula above, which avoids the cube root extraction whilst keeping the same accuracy (1.8 places per iteration).

Program Quick_Archimedes;

var a, b, c, p: real; i, m: byte; n: integer;

function Sqrt(x: Real): real;

var s, t: Real;

begin if x<>0 then begin s:=x/2; repeat t:=s; s:=(s+x/s)/2 until s=t; Sqrt:=s end else Sqrt:=0 end;

begin ReadLn(m); n:=6; b:=2*Sqrt(3); a:=3*b/4; ClrScr; Writeln(' n-gon',' ':11,'a',' ':21,'b',' ':21,'p',#10); for i:=1 to m-1 do begin a:=Sqrt(a*b); WriteLn(n:5,a:22:17,b:22:17); b:=2*a*b/(a+b); n:=2*n end; c:=Sqrt(a*b); p:=((9*a*c/(4*a-c)-7/3*(a-4*c))/10); WriteLn(n:5,c:22:17,b:22:17,p:22:17) end.

Output for nine iterations
 n-gon           a                     b                     p

6 3.00000000000000001 3.46410161513775460 12 3.10582854123024916 3.21539030917347249 24 3.13262861328123821 3.15965994209750050 48 3.13935020304686722 3.14608621513143498 96 3.14103195089050965 3.14271459964536831 192 3.14145247228546209 3.14187304997982388 384 3.14155760791185767 3.14166274705684855 768 3.14158389214831843 3.14161017660468956 1536 3.14159046322805012 3.14159703432152617 3.14159265358979327

=======================================

Edited: 21 July 2011, 12:49 p.m.

      
Re: wp34s program
Message #23 Posted by Bruce Bergman on 14 July 2011, 8:04 p.m.,
in response to message #1 by Gerson W. Barbosa

Gerson, as a present for you, I've just created an online version of the assembler. Check above in the forum postings for an announcement. :-)

Thanks,

Bruce

            
Re: wp34s program
Message #24 Posted by Gerson W. Barbosa on 14 July 2011, 8:14 p.m.,
in response to message #23 by Bruce Bergman

Hey, thank you! Just in time for my birthday next week :-)

Gerson.


[ Return to Index | Top of Index ]

Go back to the main exhibit hall