A High-Precision Floating-Point Computational System in Win32Forth
(Adjustable few hundred decimal digits)
Ching-Tang Tseng
2013/2/12 Hamilton, New Zealand
(1). General description
2013/02/12
Once upon a time, I’ve got a high-precision floating point system. I typed
PRINT 1000*0.228
This was the answer:
227.9999999999999999757
Uh hoh….. What would be your reaction? To me, unsatisfied, it should be 228. This system was not good enough. If I want a system done right, I have to design it myself.
For more than 30 years, FORTH is the one and only one programming language I preferred to use. Most of the time, When I used FORTH, I was concentrated on mathematic calculation application. During the past 4 years, ABC FORTH has been built by me. This is a BASIC programming style system in FORTH, easy for using in math, and is designed by FORTH high level definitions only. I built it step by step, a lot of public domain beautiful FORTH codes, created by the FORTH sages and able men of the past, had been used in this system.
4 years ago, I would like to do a non-profit research by FORTH. It is about the topic: To forecast the big earthquake a few days early. Signal analysis and math calculation with hardware device operation are needed in the system. And may be, to collect much enough data by automatic internet connection function should be included as well. But I have no money to buy any kinds of software or hardware. Win32Forth system and my personal technique are all I have, during these days.
After 4 years development, ABC FORTH performance is good enough for math using and for my personal research application. There were many articles related to ABC FORTH had been posted on this Blog. They were written in Chinese. All of them were going to let you know what I have done in the ABC FORTH?
There are six independent math systems in the ABC FORTH already. They are:
(1) traditional integer number
(2) traditional floating point number
(3) complex number
(4) fraction number
(5) big integer number
(6) big floating point number
When I developed my ABC FORTH, there was no supporting from anywhere or anybody. I have to do all testing by myself. Then, many practical application programs had been implemented out simultaneously. For the purpose of paying all my attention on the system development, I have no time and have no interesting to introduce all achievements in public domain.
Recently, a prototype of many digits high precision floating point arithmetic function has been carried out and has been added into my ABC FORTH. Now, I am able to use BASIC style program to deal with all math problems for high precision purpose. This is just a prototype system temporarily. I am going to improve it from now on. And wish it could be a nice tool for high precision math study in near future.
This is the first time to introduce ABC FORTH in English. Here, I would like to post some testing program for big floating point number computation and their output results instead of my poor English discussion. The following testing will be posted step by step but not in once.
Wish you enjoy it.
In this system, the precision of big floating point number is adjustable, independent with hardware math co-processor. A few hundred or thousand digits precision could be adjusted out by just to set the value of one variable named “CellsSetting” only. I am running Win32Forth in W7.
In this system, I am concerned about the accuracy only. As far as how much memory needed? How fast the running speed should be? Both of them are not concerned by me. When I set this system to be one hundred and twelve digits precision, I am able to get the output on the screen almost just after touch ENTER. Speeds are acceptable and response does make sense. Except for special condition, 112 digits precision will be a fixed setting in the program.
What is an adjustable precision function? For example, you are able to get the value of sine and let his precision higher up to more than thousand digits. This is such a value of sin ( 1 ).
ok
BIGF1 BIGFSIN BF.
8.414709848078965066525023216302989996225630607983
71065672751709991910404391239668948639743543052695
85434903790792067429325911892099189888119341032772
92124094807919558267666069999077640119784087827325
66347484802870298656157017962455394893572924670127
08648628105338203056137721820386844966776167426623
90133827533979567642555654779639897648243286902756
96429120630058303651523031278255289853264851398193
45213597095596206217211481444178105760107567413664
80550089167266058041400780623930703718779562612888
04636081734524656391420252404187763420749206952007
71334780981427902145268255663208233521544160916442
09058929870224733844604489723713979912740819247250
48855487311931035068190815153260745739291118331962
82150897348688114214528382298651257016673840744551
92375614322129060592482739703681801585630905432667
84643107531263812173256701985601106836028901895019
42151616655191791451720046686595971691072197805885
40646001994013701405309580855205280525317113323054
61638363601816994797150048515079398383039567816794
81612214022089169871097439312119047662675566086294
39084
X10^ -1
ok
According to the primitive function FSIN in the same Win32Forth system, we have
20 SIGDIGITS ! ok
1E0 FSIN FE. 841.47098480789644800E-3 ok
Directly manipulation above is just let the system seems to be the same as a calculator only. Here is a very fundamental programmable example. Let’s get the square root of the numbers from 1 to 10, and all of them, the precision are more than 100 digits.
Before I post the next testing material, these output results will be kept here.
And, I am sure, if I let my system to do this calculation: 1000 * 0.228. The output should be 228 exactly. But I would like to show it as 2.28 X10^ 2.
This is a very high precision floating point system.
2 BIGREALS AA BB
: T34
11 0
DO
H{{ AA = I>BIGF ( I ) }}H
H{{ BB = SQRT ( AA ) }}H
CR ." SQRT( " I . ." )="
BB BF.
LOOP
;
ok
T34
SQRT( 0 )=
0.0
X10^ 0
SQRT( 1 )=
1.0
X10^ 0
SQRT( 2 )=
1.414213562373095048801688724209698078569671875376
94807317667973799073247846210703885038753432764157
273501384623
X10^ 0
SQRT( 3 )=
1.732050807568877293527446341505872366942805253810
38062805580697945193301690880003708114618675724857
5675626141415
X10^ 0
SQRT( 4 )=
2.0
X10^ 0
SQRT( 5 )=
2.236067977499789696409173668731276235440618359611
52572427089724541052092563780489941441440837878227
496950817615
X10^ 0
SQRT( 6 )=
2.449489742783178098197284074705891391965947480656
67012843269256725096037745731502653985943310464023
4818594601226
X10^ 0
SQRT( 7 )=
2.645751311064590590501615753639260425710259183082
45018036833445920106882323028362776039288647454361
0615064578338
X10^ 0
SQRT( 8 )=
2.828427124746190097603377448419396157139343750753
89614635335947598146495692421407770077506865528314
5470027692461
X10^ 0
SQRT( 9 )=
3.0
X10^ 0
SQRT( 10 )=
3.162277660168379331998893544432718533719555139325
21682685750485279259443863923822134424810837930029
5187347284152
X10^ 0
Go on to the next section
\ *************************************************
To be continued.
(2). I/O format
2013/02/16
At the beginning, when I started to develop this system. I knew. The most difficulties of coding are dealing with the problem of how to define a few I/O instructions? They are helped to input and output the big floating point numbers into and out of the system. All of their function must work normally, no matter the system is under interpreted or compiled mode.
For output, I made my mind, let all big floating point numbers to be printed out with only one kind of format. This format has been showed above already. 50 characters per line are fixed. The decimal point of fraction part is always to be fixed at the position after the first mantissa integer number. Follow the mantissa, another line is used for ten to the power. It has been showed as is started with X10^ and ended with another minus or positive integer.
Two fundamental words for output are:
BF. ( addr exp -- ) “Big Floating point number Dot”
addr means the beginning address of mantissa
exp means the exponent of big floating point number
BF@. ( addr -- ) “Big Floating point number Fetch Dot”
addr means the beginning address of a big floating point number
As for input, it is more complicated than output. Lots of facts you have to consider. Lots of conditions you have to face to. Lots of frustration I have suffered while coding was progressed. Of course, at last, there were a lot of instructions I have finished. All of them are helped to convert many kinds of number’s format into a big floating point number, and then let my system is able to keep going to have successful operation on it.
For the time being, I would like to introduce only one kind of the input format to you. This format is good enough to let me easy to finish the discussion of this article and is apt to remember for the future using.
No matter how simple or how many or how big the input number is, this format is always to be represented as following:
S” 1.0”X10^ 1 POWER
S” 1.23” X10^ 0 POWER
S” -9.87654321” X10^ -397 POWER
S” 9.999999999999999999999999999999999999999999999” X10^ 0 POWER
The rules are:
(1) Mantissa is always quoted by one pair of instructions S” ………”
(2) The format of mantissa is always in terms of the same format as output. The decimal point should be place at the position following the first digit. Leading negative mark is permit. Number is delimited by a double quotation mark, has been showed as above.
(3) X10^ is one instruction. It means “times ten to the power of”.
(4) Exponent is a general integer number.
(5) All of all are ended with another instruction named POWER.
(6) If the input number is too big, too long or longer than one line width in one screen, this kind of number could be cut into many sections, other two instructions have to use. They are FIRST-SECTION and NEXT-SECTION.
A typical BASIC style simple program and its output result after execution can tell you how did this program to reach these rules?
In this program, “(“and “)”, both of them are unnecessary everywhere. They are using for the purpose of easy to identify there.
3 BIGREALS AA BB CC
: T40
BASIC
10 LET H{ AA = ( S" 1.23" X10^ 1 POWER ) + ( S" 4.56" X10^ 1 POWER ) }H
:: H{ BB = ( S" 7.89" X10^ 5 POWER ) }H
:: H{ CC = AA / BB + AA }H
20 RUN CC BF.
30 LET H{ CC = NEGATE ( CC ) }H
40 RUN CC BF.
50 LET H{ CC = ABS ( CC ) }H
60 RUN CC BF.
70 END
;
ok
T40
5.790007338403041825095057034220532319391634980988
59315589353612167300380228136882129277566539923954
3726235741444
X10^ 1
-5.79000733840304182509505703422053231939163498098
85931558935361216730038022813688212927756653992395
43726235741444
X10^ 1
5.790007338403041825095057034220532319391634980988
59315589353612167300380228136882129277566539923954
3726235741444
X10^ 1
ok
As to the data structure, for the present, only BIGREALS has been used in this article.
Except BIGREALS,
BIGREAL, BIGFVARIABLE, BIGFVARIABLES, BIGFVALUE, BIGFVALUES, BIGFCONSTANT had been built as well.
No need to introduce you the others now.
Go on to the next section.
\ *******************************************************
To be continued
(3). Demonstration of Application
2013/2/18
\ Error function program:
\ MIMI: mini mini tiny register for checking of lower bound
\ XXR : x Argument register for input
\ YYR : y = f(x) Accumulator for output
\ WKR : Working register for working in one term
\ NR1 : Numerator register 1
\ NR2 : Numerator register 2
\ NR3 : Numerator register 3
\ DR1 : Denominator register 1
\ DR2 : Denominator register 2
\ DR3 : Denominator register 3
10 BIGREALS MIMI XXR YYR WKR NR1 NR2 NR3 DR1 DR2 DR3
CellsSetting 10 * NEGATE VALUE TINY \ to be fixed while loading
\ CONSTANT 10 was changed for ERFC as follows.
\ x=12 --> use 4, x=14 --> use 6, x=16 --> use 8, x=17 --> use 10
: InputData
BASIC
10 LET H{ MIMI = S" 1.0" X10^ TINY POWER }H
20 LET H{ XXR = S" 3.0" X10^ 1 POWER }H \ 17.0 ~~> huge #, OK for erfc
\ 30 PRINT " MIMI = :"
\ 40 RUN MIMI BF.
50 PRINT " X = :"
60 RUN XXR BF.
\ 70 PRINT " ERF(x) = :"
80 PRINT " ERFC(x) = :"
90 END
;
\ erf(x) = (2/sqrt(pi)) * ( integration(0-->x) of (exp(-t^2) dt ) )
\ = (2x/sqrt(pi)) * [ 1 - x^2/1!*3 + x^4/2!*5 - x^6/3!*7 + ..... ]
\ 0 <= x <= 11 OK for erf.
\ 4 <= x <= 11 OK for erfc, if needed, erfc(x)=1-erf(x)
: ERF
\ H{{ XXR }}H BIGF! \ ( addr exp -- addr exp )
BASIC
10 RUN InputData
20 LET H{ NR1 = BIGF1 }H
:: H{ DR1 = BIGF0 }H
:: H{ DR2 = BIGF1 }H
:: H{ NR2 = BIGF1 }H
:: H{ YYR = BIGF1 }H
30 LET H{ DR1 = DR1 + BIGF1 }H
:: H{ DR2 = DR2 * DR1 }H
:: H{ NR1 = NR1 * XXR * XXR }H
:: H{ NR2 = BIGF2 * DR1 + BIGF1 }H
:: H{ WKR = NR1 / ( DR2 * NR2 ) }H
:: H{ YYR = YYR - WKR }H
40 LET H{ DR1 = DR1 + BIGF1 }H
:: H{ DR2 = DR2 * DR1 }H
:: H{ NR1 = NR1 * XXR * XXR }H
:: H{ NR2 = BIGF2 * DR1 + BIGF1 }H
:: H{ WKR = NR1 / ( DR2 * NR2 ) }H
:: H{ YYR = YYR + WKR }H
50 IF H{ WKR > MIMI }H THEN -30
60 LET H{ YYR = YYR * ( BIGF2 * XXR ) / SQRT ( PI ) }H
\ 65 let h{ yyr = bigf1 - yyr }h \ for erfc(x<11)
70 RUN YYR BF.
80 END
\ YYR \ ( addr exp -- addr exp )
;
\ erfc(x) = [1/(exp(x*x)*sqrt(pi)*x)] * [ 1 - 1/2x^2 + 1*3/(2x^2)^2
\ - 1*3*5/(2x^2)^3 + ..... ]
\ x = 17 --> huge x, OK for erfc
: ERFC
BASIC
10 RUN InputData
20 LET H{ NR1 = BIGF1 }H
:: H{ NR2 = BIGF1 }H
:: H{ DR1 = BIGF2 * XXR * XXR }H
:: H{ DR2 = DR1 }H
:: H{ YYR = BIGF1 - ( NR2 / DR2 ) }H
30 LET H{ DR2 = DR2 * DR1 }H
:: H{ NR1 = NR1 + BIGF2 }H
:: H{ NR2 = NR2 * NR1 }H
:: H{ WKR = NR2 / DR2 }H
:: H{ YYR = YYR + WKR }H
40 LET H{ DR2 = DR2 * DR1 }H
:: H{ NR1 = NR1 + BIGF2 }H
:: H{ NR2 = NR2 * NR1 }H
:: H{ WKR = NR2 / DR2 }H
:: H{ YYR = YYR - WKR }H
\ 45 RUN WKR BF. \ See divergent point from MIMI
50 IF H{ WKR > MIMI }H THEN -30
60 LET H{ DR3 = EXP ( XXR * XXR ) }H
:: H{ DR3 = DR3 * SQRT ( PI ) }H
:: H{ DR3 = DR3 * XXR }H
70 LET H{ YYR = YYR / DR3 }H
80 RUN YYR BF.
90 END
;
\S
\ *********************************************************
\ Testing results:
erf
X = :
0.0
X10^ 0
ERF(x) = :
0.0
X10^ 0
erf
X = :
1.0
X10^ -1
ERF(x) = :
1.124629160182848922032750717439683832216962991597
02547534494144817599244025531359170221174136405196
9494684444739
X10^ -1
erf
X = :
2.0
X10^ -1
ERF(x) = :
2.227025892104784541401390068001438163882690384302
27605620935023888363674271912182870353723778740700
3625084786272
X10^ -1
erf
X = :
3.0
X10^ -1
ERF(x) = :
3.286267594591274276389140478667565511699180962626
75822609143139853710579890618085895322222697478706
2285426775476
X10^ -1
erf
X = :
4.0
X10^ -1
ERF(x) = :
4.283923550466684551036038453201724441218629285225
90383495086346113333368505671802422726158110538627
1012415050084
X10^ -1
erf
X = :
5.0
X10^ -1
ERF(x) = :
5.204998778130465376827466538919645287364515757579
63700058805725647193521716853570914788218734787757
0329661243901
X10^ -1
erf
X = :
1.0
X10^ 0
ERF(x) = :
8.427007929497148693412206350826092592960669979663
02908459937897834717254096010841261983325348144888
4541582615366
X10^ -1
erf
X = :
2.0
X10^ 0
ERF(x) = :
9.953222650189527341620692563672529286108917970400
60076738352326200437280719995177367629008019680680
487939328715
X10^ -1
erf
X = :
3.0
X10^ 0
ERF(x) = :
9.999779095030014145586272238704176796201522929126
00750342761045157057543316379867732183745349184683
712856553818
X10^ -1
erf
X = :
4.0
X10^ 0
ERF(x) = :
9.999999845827420997199811478403265131159514278547
46410808831657095005786958973188745927234986543819
418211441985
X10^ -1
erf
X = :
5.0
X10^ 0
ERF(x) = :
9.999999999984625402055719651498116565146166211098
81949685276620069312085944079608635441308535281856
140283421212
X10^ -1
erf
X = :
6.0
X10^ 0
ERF(x) = :
9.999999999999999784802632875010868834066496008126
15369522485938311457899472107948943662761515072149
437615991014
X10^ -1
erf
X = :
7.0
X10^ 0
ERF(x) = :
9.999999999999999999999581617439222058560138598977
61000677499703825861875398479583656469504534922860
172790264017
X10^ -1
erf
X = :
8.0
X10^ 0
ERF(x) = :
9.999999999999999999999999999887757028270170729200
32111556829720906568070835521036614058191866173481
9237285022
X10^ -1
erf
X = :
9.0
X10^ 0
ERF(x) = :
9.999999999999999999999999999999999995862968253486
18976194609653263747540428980389591497158310489836
00184378937
X10^ -1
erf
X = :
1.0
X10^ 1
ERF(x) = :
9.999999999999999999999999999999999999999999979115
12416237455242999213691611093205070508043921132553
652299156938
X10^ -1
erf
X = :
1.1
X10^ 1
ERF(x) = :
9.999999999999999999999999999999999999999999999999
99998559133974903630731655919342161305061590478873
27341811488
X10^ -1
erf
X = :
1.2
X10^ 1
ERF(x) = : \ ***** Beware! this is a wrong result! *****
1.000000000000000000000000000000000000000000000000 00028964232747623645928111862273111180446503830761
9596656136387
X10^ 0
\ **************************************************
\ When 11<x<17, their precision are not the same as others.
erfc
MIMI = :
1.0
X10^ -48
X = :
1.2
X10^ 1
ERFC(x) = :
1.356261169205904212780306156590417572666782233287
92211903650655213650444305226488837167452692414329
9360066305547
X10^ -64
erfc
MIMI = :
1.0
X10^ -72
X = :
1.4
X10^ 1
ERFC(x) = :
3.037229847750311665115172806783328447911896669343
42259515007060199595805908881107988253347572678975
2872977727062
X10^ -87
erfc
MIMI = :
1.0
X10^ -96
X = :
1.6
X10^ 1
ERFC(x) = :
2.328485751571530693364872854573442597534396948094
94802151649509557369523737790080619238332625392233
1153810817962
X10^ -113
\ ***************************************************
erfc
X = :
1.7
X10^ 1
ERFC(x) = :
1.021228015094260881145599235077652994401543730587
57117255350655550943197600942423095274220475440447
1071116583101
X10^ -127
erfc
X = :
1.8
X10^ 1
ERFC(x) = :
6.082369231816399307668466715702274949587610076867
57685159278589190532744511115561184378302587387186
5094174183213
X10^ -143
erfc
X = :
1.9
X10^ 1
ERFC(x) = :
4.917722839256475446413297625239608170930672249989
25903631521698356756256672056877874774045924679965
6966106818571
X10^ -159
erfc
X = :
2.0
X10^ 1
ERFC(x) = :
5.395865611607900928934999167905345604088272670923
60528347010378491301876315666993333543074240234371
03862266973
X10^ -176
erfc
X = :
2.1
X10^ 1
ERFC(x) = :
8.032453871022455669021356947138268888967875692513
37885278602378401460562630878875250682660497529814
4290532553308
X10^ -194
erfc
X = :
2.5
X10^ 1
ERFC(x) = :
8.300172571196522752044012769513722768714223191886
69959605444858200128083858466509582535907094157923
0337303216692
X10^ -274
erfc
X = :
3.0
X10^ 1
ERFC(x) = :
2.564656203756111600033397277501447146548889722778
61705412259958618423869477919735075745592460023192
5697543021348
X10^ -393
erfc
X = :
4.0
X10^ 1
ERFC(x) = :
1.896961059966276509268278259713415434936907563929
18618346283475290041180520511188660525669077676004
1365306022708
X10^ -697
erfc
X = :
5.0
X10^ 1
ERFC(x) = :
2.070920778841656048448447875165788792932250920995
39968376623553375989816942513882393211712997100144
289795083031
X10^ -1088
erfc
X = :
1.0
X10^ 2
ERFC(x) = :
6.405961424921732039021339148586394148214414399460
33805776710765024890255482950583122794586657987149
7097903402585
X10^ -4346
erfc
X = :
2.0
X10^ 2
ERFC(x) = :
4.689359295983601498035699776992316770422384630544
86566387750769583162476707139525208321549676277950
7240566295404
X10^ -17375
\ Testing was stopped at x = 200.
\ If x>200, every x must take a long period of time for testing, owing to the exp(x*x) evaluation when x*x > 40000.
(4). Acknowledgement
These days, I am very sad for the disappearance of FSL (Forth Scientific Library) website. This is the first time I have the opportunity to say many thanks to all contributors of FSL in English. In my Chinese posted articles, I was always to introduce Chinese all the FSL program contributors and to say thanks in Chinese for FSL to all readers. I gave the name “FORTH sages and able man of the past” to all contributors of FSL, they helped to build Forth Scientific Library. And so let me easy to develop my ABC FORTH. Even now this system has not been finished yet. I would like to say thanks first to some key FORTH sages for this article.
Thanks to:
Leonard Francis Zettel, his big number program “Forth Scientific Library Algorithm #47” had been used in my system.
There are many very important persons else I have to say thanks to them, especially, all contributors of Win32Forth system. And if let us get back to the early years, Both Forth sages Charles H. Moore and Michael Perry, had to be included also. Their code “Tiny BASIC compiler” had been used in my system as well, and it was the fundamentals of my ABC FORTH system.
I believe. You are coming here for the technical reasons not for English. And I prefer to practice FORTH coding rather than to prepare English article. I posted this article by way of Tom Zimmer’s editor WinEd first, reviewed it on site later. It is inevitable to modify something of this article from time to time.
2013/03/17 modified
2016/06/09 modified
\ ******************* THE END ********************
