▼
Posts: 2,761
Threads: 100
Joined: Jul 2005
=======================================
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.
▼
Posts: 255
Threads: 22
Joined: May 2011
... 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/cgisys/cgiwrap/hpmuseum/forum.cgi?read=187663#187663
//
// $ wp34s_asm.pl barbosa.wp34s o wp34s0.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
▼
Posts: 2,761
Threads: 100
Joined: Jul 2005
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.
▼
Posts: 255
Threads: 22
Joined: May 2011
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 "wp34s0.dat", "wp34s1.dat", or "wp34s2.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 Nsteps, 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...
▼
Posts: 2,761
Threads: 100
Joined: Jul 2005
Hello Neil,
Thanks! I've installed the latest version of Strawberry Perl:
>wp34s_asm.pl mpgm.wp34s o wp34s0.dat
// CRC16: 5ACD
// Total words: 79
// Total steps: 78
>wp34s_asm.pl dis wp34s0.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 USBtoRS232 converter yet.
Regards,
Gerson.
▼
Posts: 3,229
Threads: 42
Joined: Jul 2006
Firmware version mismatch between the assembler and the calculator?
 Pauli
▼
Posts: 2,761
Threads: 100
Joined: Jul 2005
I'll check that, thanks! Emulator build date is 2011714 (0).
Gerson.
▼
Posts: 3,229
Threads: 42
Joined: Jul 2006
The current emulator is broken too :(
 Pauli
Posts: 255
Threads: 22
Joined: May 2011
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 opcode numbers sometimes migrate with the everimproving state of the development. The assembler can handle this IFF you have the correct opcode 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 opcode 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.FCN3ENTER" (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 wp34s0.dat svn
// Opcode SVN version: 1210
// CRC16: 5ACD
// Total words: 79
// Total steps: 78
Just for interest, if you have an MD5 or SHA1 tool (look on Wiki), the hashes for the DAT file should be:
$ sha1sum wp34s0.dat
7b5213cf6e9176c282cc1933f578891a89b4bbf6 wp34s0.dat
$ md5sum wp34s0.dat
9f875de2bc41422e3c1a639c84103cfc wp34s0.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 "wp34s0.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 (wp34s1.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 02 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, wp34s0.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).
▼
Posts: 2,761
Threads: 100
Joined: Jul 2005
Hello Neil,
Here are what I have:
>wp34s_asm.pl mpgm.wp34s o wp34s0.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.
▼
Posts: 255
Threads: 22
Joined: May 2011
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...
▼
Posts: 2,761
Threads: 100
Joined: Jul 2005
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.
▼
Posts: 255
Threads: 22
Joined: May 2011
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 techotrinity'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 outofdate 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 opcode 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! :)
▼
Posts: 3,283
Threads: 104
Joined: Jul 2005
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.
▼
Posts: 255
Threads: 22
Joined: May 2011
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. :)
Posts: 3,229
Threads: 42
Joined: Jul 2006
No version results from using the web interface to download the subversion files. You need a proper subversion client.
 Pauli
Posts: 151
Threads: 8
Joined: Oct 2009
Interesting  it appears to compute pi with the number of iterations entered in X and the followup 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 WP34S features is remarkably powerful and remarkably writeonly. It's not obvious to me yet what algorithm you're using. But it is fun to dig into...
▼
Posts: 2,761
Threads: 100
Joined: Jul 2005
Quote:
it appears to compute pi with the number of iterations entered in X and the followup iterations in Y or such.
Quite exactly!
Quote:
I have to note that RPN programming with all the WP34S features is remarkably powerful and remarkably writeonly.
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 DevC++ :
#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 *= 11/(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 = (1sqrt(1x))/2;
n += n;
}
printf("\n");
system("pause");
return 0;
}
Posts: 3,229
Threads: 42
Joined: Jul 2006
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 WP34S features is remarkably powerful and remarkably writeonly. 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
▼
Posts: 2,761
Threads: 100
Joined: Jul 2005
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 96gon (3*2^6)
.
R/S > 3.14159265360 area of the 96gon 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: base10 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.
▼
Posts: 2,761
Threads: 100
Joined: Jul 2005
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("\ ngon\t wgm\t wam\t\t(9*wgm4*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*wgm4*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:
ngon wgm wam (9*wgm4*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 WA solve the equation (solve (p(27*(a^2*b)^(1/3)4*(2*a+b))/15 )/((9*a*c/(4*ac)7/3*(a4*c))/10p)==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' 96gon 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)/(c4*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
▼
Posts: 2,761
Threads: 100
Joined: Jul 2005
FWIW, the TurboBCD program below is equivalent to the 38step 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(ct)<1e16;
Cbrt:=c
end
else
Cbrt:=0
end;
begin
ReadLn(m);
n:=6;
b:=2*Sqrt(3);
a:=3*b/4;
ClrScr;
Writeln(' ngon',' ':11,'a',' ':21,'b',' ':21,'p',#10);
for i:=1 to m1 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:
ngon 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(' ngon',' ':11,'a',' ':21,'b',' ':21,'p',#10);
for i:=1 to m1 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*ac)7/3*(a4*c))/10);
WriteLn(n:5,c:22:17,b:22:17,p:22:17)
end.
Output for nine iterations
ngon 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.
Posts: 553
Threads: 57
Joined: Sep 2006
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
▼
Posts: 2,761
Threads: 100
Joined: Jul 2005
Hey, thank you! Just in time for my birthday next week :)
Gerson.
