A High-Precision Floating-Point Computational System in Win32Forth

A High-Precision Floating-Point Computational System in Win32Forth

(Adjustable few hundred decimal digits)

Ching-Tang Tseng

2013/2/12 Hamilton, New Zealand

ilikeforth@gmail.com

http://forthfortnight.blogspot.com

(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 ********************