Tweet
Question: If A## should equal B##, when is A##/B## >> one googol?
Answer: A=10^4932 and B=Val("1E4932")
Although B is less than the extended-precision maximum of approximately 1.2E4932, it is erroneously evaluated (in PBCC Version 5) as 2^127. Support has been notified and issue sent to R&D.
There are four distinct ways in PowerBASIC to compute powers of ten. The methodologies were found by trial and error and are listed in order of increasing accuracy:
(a) For non-integer powers two methods are used. Both PB implementations return zero rather than a denormal for x<=-4933. They are:
(a1) 10^x (only if x is not an integer) is evaluated as
Exp(x*Log(10))
(a2) Exp10(x) (whether or not x is an integer) is evaluated as
Exp2(x*Log2(10)).
(b) For integer powers two methods are also used, each involving factors of the form 10^(2^i) (0<=i<=12). If n is negative, the reciprical is taken at the end. Note that
10^27 < 2^63*2^27 < 10^28
So for 0<=n<=27, the exact powers can be generated. The two methods are:
(b1) 10^n, with n being an integer, is calculated by generating a table using:
10^(2^(0))=10, 10^(2^(i+1))=(10^(2^(i)))*(10^(2^(i)))
Then multiplication begins with the most significant bit of |n|. Multiplying in this order adds rounding error as more multiplications cannot be performed exactly.
(b2) 1En (as a literal) always equals (as long as no compile-time error occurs) Val("1E"+Format$(n)). It is calculated using a precomputed table of the accurate approximations of
10^(2^i), 0<=i<=12
Then the multiplication starting with the least significant bit of |n|. An erroneous result occurs for n=4932 or greater where 2^127 is returned. Also the compiler will generate Error 473 if a literal evaluates to less than 1E-4933 or to equal or greater than 1E+4932.
A little bit of algorithm analysis for those interested. Given:
y = 10^x
By elementary calculus:
(dy/y) = (x*ln(10))*(dx/x)
Thus the relative error in x is increased by a factor of approximately 2.3*x which partly explains the errors for non-integer values of x in the above table. If x is a rational number or can be given to higher accuracy by other means, it is best to convert it to the sum of an integer and a floating point number. Note that
x = Fix(x)+Float(x)
for all x of interest here and |Float(x)|<1. Alternatively, use CLng(x) to get a floating point term of absolute value under 0.5. In the allowed ranges, the Floating Point Unit (FPU) command !f2xm1 was designed to produce results accurate to 1.5 bits so, if x can be well-defined, splitting the agument should give excellent accuracy.
For the error analysis above, a table of the best values for powers of ten was created. As this involves almost ten thousand lines, a shorter function that gives the same results is given below:
Below is a function Val1E() that provides a run-time evaluation of integer powers of ten using the same methodology as
Val("1E"+Format$(n)) = 1En
but is considerably faster. A table of powers of five rather than ten was used so that extreme negative powers remained within extended-precision limits with just the delay of a !fscale command at the end. Please note the following (timings based on Pentium 4):
Avoids compile-time error #473 (Ivalid numeric format).
Methodology is equivalent to PB when PB gives valid results.
Takes up about 300 bytes.
Significand accuracy +/- 4 over entire range.
Up to 4000+ times more accurate and 75% faster than Exp(n*Log(10)).
Up to 2000+ times more accurate and 25% faster than Exp10(n).
Up to 36 times more accurate and four times as fast as 10^n.
Same accuracy and eleven times faster than Val("1E"+Format$(n)).
Returns zero if n<=-4951; infinity if n=>+4933.
Returns denormals for -4950 <= n <= -4932.
Equals 1En and Val("1E"+Format$(n)) for -4932 <= n <= 4931
Corrects PB error where Val("1E4932")=2^127.
Best values for -29 <= n <= +32 and exact for 0 <= n <= +27.
Use in combination with Exp10_() (see below) as:
Val1E(Fix(x))*Exp10_(Float(x)) or (perhaps better)
Val1E(CLng(x))*Exp10_(x-CLng(x)).
Lastly an assembly version for Exp10() that (perhaps because it uses !fldl2t to load Log2(10)) runs twice as fast as the PB version is given. The FPU command for raising a number to a power is !f2xm1 that returns 2^x-1 (|x|<=0.5 or 1.0 depending on the coprocessor in use). In order to use this, the following transformation is used where i is an integer and d=x*Log2(10)-i (|d|<=0.5):
10^x = 2^(x*Log2(10)) = 2^i*2^d = ((2^d-1)+1) * 2^i
Conclusions for raising ten to the power of x:
1. If just double-precision accuracy is desired, any method will work, and Exp10_(x) is the fastest.
2. If x ia an integer use Val1E(x) for its speed and accuracy. If maximum accuracy is desired and size is not a concern create a constant array within the code (using !dw pseudo-ops). The latter method, with a complete table, requires almost 97 kilobytes. A BASIC array (including zero and infinity extremes) can be dimensioned thus:
Dim MyArray(-4951 To 4933) As Static Ext At CodePtr(MyLabel)
In assembly, the starting address of the table can be accessed with:
!lea ebx,MyLabel ' Load register with equivalent address.
3. If x is not an integer but is of small maqnitude (less than 5, but preferably less than 2) use Exp10_(x).
4. If x can be split into an integer (integer) and a fractional part (delta) where the latter is more accurate than using x alone (say, x is a rational fraction or a number with known additional digits), use the product of Val1E(integer) and Exp10_(delta).
5. If x may or may not be an integer, use 10^x. Recognize however this is the least accurate method if x is a non-integer but a fairly accurate method otherwise. The compiled program makes the decision.
6. Always remember that as |x| gets larger, the relative error increases proportionately unless the programmer can use some means to circumvent it.
Note added in 1/2/09 edit:
The use of Exp2() (and probably all other power functions with a non-integer exponent) has an issue in and of itself. It appears that although PBCC is written for 32-bit machines, the Exp2() function used by PB uses a form of the FPU command f2xm1 that was written for 16-bit machines. This form limits the ST(0) value to 0<=ST(0)<=0.5. For further details see my post below and Intel, "IA-32 Intel Architecture Software Developer's Manual", Volume 3, Section 18.14.7.10.
_______________
Mark Longley-Cook
Answer: A=10^4932 and B=Val("1E4932")
Although B is less than the extended-precision maximum of approximately 1.2E4932, it is erroneously evaluated (in PBCC Version 5) as 2^127. Support has been notified and issue sent to R&D.
There are four distinct ways in PowerBASIC to compute powers of ten. The methodologies were found by trial and error and are listed in order of increasing accuracy:
(a) For non-integer powers two methods are used. Both PB implementations return zero rather than a denormal for x<=-4933. They are:
(a1) 10^x (only if x is not an integer) is evaluated as
Exp(x*Log(10))
(a2) Exp10(x) (whether or not x is an integer) is evaluated as
Exp2(x*Log2(10)).
(b) For integer powers two methods are also used, each involving factors of the form 10^(2^i) (0<=i<=12). If n is negative, the reciprical is taken at the end. Note that
10^27 < 2^63*2^27 < 10^28
So for 0<=n<=27, the exact powers can be generated. The two methods are:
(b1) 10^n, with n being an integer, is calculated by generating a table using:
10^(2^(0))=10, 10^(2^(i+1))=(10^(2^(i)))*(10^(2^(i)))
Then multiplication begins with the most significant bit of |n|. Multiplying in this order adds rounding error as more multiplications cannot be performed exactly.
(b2) 1En (as a literal) always equals (as long as no compile-time error occurs) Val("1E"+Format$(n)). It is calculated using a precomputed table of the accurate approximations of
10^(2^i), 0<=i<=12
Then the multiplication starting with the least significant bit of |n|. An erroneous result occurs for n=4932 or greater where 2^127 is returned. Also the compiler will generate Error 473 if a literal evaluates to less than 1E-4933 or to equal or greater than 1E+4932.
Code:
Maximum absolute errors in the significand (treating it as quad integer):
Maximum absolute integer power of ten: 30 100 300 1000 Limit---v
Method
(a1) Exp(n*Log(10)) = lim(x->n)(10^x) 45 164 678 3823 16637 4932
(a2) Exp10(n) = Exp2(n*Log2(10)) 56 236 459 2060 8989 4932
(b1) 10^n (n=integer) 1 2 8 35 143 4932
(b2) 1En = Val("1E"+Format$(n)) 1 1 2 3 4 4931
y = 10^x
By elementary calculus:
(dy/y) = (x*ln(10))*(dx/x)
Thus the relative error in x is increased by a factor of approximately 2.3*x which partly explains the errors for non-integer values of x in the above table. If x is a rational number or can be given to higher accuracy by other means, it is best to convert it to the sum of an integer and a floating point number. Note that
x = Fix(x)+Float(x)
for all x of interest here and |Float(x)|<1. Alternatively, use CLng(x) to get a floating point term of absolute value under 0.5. In the allowed ranges, the Floating Point Unit (FPU) command !f2xm1 was designed to produce results accurate to 1.5 bits so, if x can be well-defined, splitting the agument should give excellent accuracy.
For the error analysis above, a table of the best values for powers of ten was created. As this involves almost ten thousand lines, a shorter function that gives the same results is given below:
Code:
Function bestPowerOfTen(ByVal n As Long) As Ext
' Returns most accurate 10^n approximations.
Static LSNibble As String ' Least significant nibble (4 bits)
Local x As Ext ' Result
Dim d(0) As Dword At VarPtr(x) ' Least significant DWord
Local t As Dword ' Copy of d(0)
'
' 1E-4951 = Zero; 1E+4933 = Infinity
' 1E0 thru 1E+27 are exact.
'
' Least significant nibbles of extended-precision significands:
If LSNibble="" Then LSNibble= _
"2B279D6CD3252CAB0401D23B5EA864A9FCA1938D86386CD9BF"+ _' -4950 to -4901
"2F1A8FF782BF6CB64EA9F3860C38E15A8711D4519856FAC38E"+ _' -4900 to -4851
"125ADAD449049B7564055E215F97D9B2B15E793089268BA456"+ _' -4850 to -4801
"7EEAADC727838ED41FEA64D45AA9FBADE572FF2FB2B230C68A"+ _' -4800 to -4751
"55EAA832EAC79DDCF5783CBFB9E57A11E602238ED08A649F40"+ _' -4750 to -4701
"D8EDCF7DC7C30C7934D4948271EE2D8EFB96B2BEE3F76498ED"+ _' -4700 to -4651
"6825AD879F71598ED040087D852F53B300022F16B52BF6C079"+ _' -4650 to -4601
"A41B6FDD4600819874982AA5271644A8EDCCA5AC40157D0671"+ _' -4600 to -4551
"7012799414AC7D91E174E1DC38E963C7DCA089B2DD06CED4D2"+ _' -4550 to -4501
"3006B2D195749FB70C049FFB7525E68152578BE1FBDC3C2EA6"+ _' -4500 to -4451
"890C3ED8BB94C7516B1D478E345BE13B60048AD49F78A5F3F2"+ _' -4450 to -4401
"A0E26CFBF64A1D64D89B25BD2BAC7D000ADCBE66FB15ED45A6"+ _' -4400 to -4351
"C74D127D2741D89BACFBFB98A04BD94415EE3F764D6085AC8E"+ _' -4350 to -4301
"94D495BD4A06BA271015DC73715E57D4DC71EA0EE1303C798D"+ _' -4300 to -4251
"D8A59F3679EDC32B9A4DEE512E3443FB711111704DC3022363"+ _' -4250 to -4201
"4CC22636868DD49B11A6FBE111DC7AD4D8608D645D41156411"+ _' -4200 to -4151
"AEDD08AE2E30C38230CCF2DDCBEE745700068E9F29CF1ACF7D"+ _' -4150 to -4101
"2FAFE66F70EACE93A568AC5B9638D0D41D6CB9701E266C682F"+ _' -4100 to -4051
"BACA825B5BE66791ACB7932AD645F376CF98256C03023B641D"+ _' -4050 to -4001
"00CE56C304D8F379F6CFB5266CE2F38D4D2375AC6BA2B1BA42"+ _' -4000 to -3951
"304979A49401230E5349BBD9B2FD9322FF38AA1DCFFD8F3081"+ _' -3950 to -3901
"5A46F716538D00C5FED0C008D829308A1BD0CBE193E2A0498F"+ _' -3900 to -3851
"74FADCBE3FFF3C4C0CF6C386949E15F7971AA93ED0A52B6C16"+ _' -3850 to -3801
"7519BBA0275BE59B29C7AD89FFDD0E5BD2BE56FF302BDA1DEE"+ _' -3800 to -3751
"1D8A08EF3BA6C3A8A4A4CC267D4D879FD5A641C4CC71DE6B3C"+ _' -3750 to -3701
"B7C4A19411608D86E5BEC8DEA9BE23C38A4FA9FBA4D5B6F3CB"+ _' -3700 to -3651
"A2E23864566755A92E29BA08E3CB5BEA6C6452571042FB9CFF"+ _' -3650 to -3601
"FBE4A044ACCE6ED85FE7D4A56223CF3C23ED56A8A2B66B2979"+ _' -3600 to -3551
"8ACED8FBAA4256C76C7C79C3368138EEEAE12D832ACC5B5349"+ _' -3550 to -3501
"CB1119004271C34B5F3CFBE2A41A04557CE159C3079AC7522B"+ _' -3500 to -3451
"E2A64519C53BB2A0D46CBD045229C68E63C3C12E1978FACB6E"+ _' -3450 to -3401
"9856C5FA8E2F5E23CF98A04DEE5F3879F71D2B672B6634A1DE"+ _' -3400 to -3351
"67230A4A6FF894AD0CFF3C37B19419AC742EBD95308ACC6089"+ _' -3350 to -3301
"F386F8A08A8578D05FBD608F3BF3715A9570682F30008FA5D4"+ _' -3300 to -3251
"6BA15A49675797AD8300CB2308F3B1D9D1579387D8111D08B6"+ _' -3250 to -3201
"7A1171538AC348A5A05A2B178B1974D9FFB93CC68BE129F7C3"+ _' -3200 to -3151
"C452349E9F982B26C757D0A8B19F30815E2D45300008930860"+ _' -3150 to -3101
"E16FF308E671A270C3815B15ED008A49345756CA1F232E25A1"+ _' -3100 to -3051
"D8A45A3BA23BE30C38279F8ED60715A8E234DEE19A0493C986"+ _' -3050 to -3001
"5B29713F742230C75F3CB2D0C82ACBE7564153040CB7D8F793"+ _' -3000 to -2951
"048711197D0082E7D0CB6386A0DC30F0C28DD22322782BFFAC"+ _' -2950 to -2901
"CEBA08E223453C7C406006FB153C7027C452ED4137156FCFAE"+ _' -2900 to -2851
"2E1D9D5EFE6E1D4986C37B5AFEA2B9675FBE7451634D0C3C79"+ _' -2850 to -2801
"34045900AC386B19852B9F7A56BDC7749A82D5686F1E1B940D"+ _' -2800 to -2751
"4E26CFFF7164C38E2194C3871D8A5DCFF6FB78E9306CE9F304"+ _' -2750 to -2701
"96BA607700EACF487682FB9E5B3798A00C7645FF7C08238E97"+ _' -2700 to -2651
"4D8A45B970A1534DC302EE51A81DD4D8F34000234A192A53CE"+ _' -2650 to -2601
"6E19A876C2CA15E5D9BF68B2A2F74D9326C75663433CC3B3CF"+ _' -2600 to -2551
"19C300CFB29FB5AD045B220B2B2B115386E9898A4F2742AEA0"+ _' -2550 to -2501
"41D68EFEEBD9222908AA11682DD0A4DA52570E6B52FFB2A97C"+ _' -2500 to -2451
"E97FEAE6F90400164DA5E3BEB5F3A08568979CA5D8B9BD083A"+ _' -2450 to -2401
"C341DC4D460CD5A00430F22B9B13CA5E9F49B3FF1EDA9F1A0E"+ _' -2400 to -2351
"636FF853A8ED0CB198B123825A1938D45BE98E55B9364D041B"+ _' -2350 to -2301
"9B7DCE52345FBE1DC07179B53CF0BE111F373759BE1D00C715"+ _' -2300 to -2251
"2ABD456B5A5BA47C8563CB2D11A9BA238F27A9B308D4DEE553"+ _' -2250 to -2201
"F5793C79B16ED83413757194D8DC3AC79BE64D4119863AC7F6"+ _' -2200 to -2151
"8B6727D23F34D641B989826E9CB2367156B756E93864A93CF7"+ _' -2150 to -2101
"5BE6D1119F949941ADC2E2F787EE9671F2FD5EB526FBB90ED5"+ _' -2100 to -2051
"1D4FA153C30866F3CFE308FFEB2A0C3A9B9CBE4E9E9B308BE5"+ _' -2050 to -2001
"11A456A093C341D975F1D55EE7CFA4515FDCC1A89BA59F78A9"+ _' -2000 to -1951
"A0C9089BA237F23ED5D49ED00416890CFBB9819B5A5115B23A"+ _' -1950 to -1901
"416F75574F64A9F4C8FEA27145BD49F94DA42386CE27D81E9C"+ _' -1900 to -1851
"2F5FA379EDC7B5F15FBA02BE4E9C74B16BE230829CEFB6BA1D"+ _' -1850 to -1801
"4DE9404DE270457194152711D089F7FA53C382ED492ACA0085"+ _' -1800 to -1751
"A51DD0CF75BBD579F3C7868D000C71B5FD118FE9BA641EE130"+ _' -1750 to -1701
"0C74B9CB3B29C7CBEFEA2FB704D8A415982A6087D04D875E33"+ _' -1700 to -1651
"C2E667D82E523E5AD0866BED8AC71D8AC34FA916FE11279704"+ _' -1650 to -1601
"53004D16008A0649A8A497DA151ACD5216B385574FEEB22B63"+ _' -1600 to -1551
"0CB68E133C49BF2F38A0D482FFF7C3CFE5E9BA00C2AD8A8B2F"+ _' -1550 to -1501
"DC7197608BE9770497515AFE60F30CB38E19B9BED0C7156ED4"+ _' -1500 to -1451
"D49CEAE6BD08756AC07D417498A0049CF641D4DDCB967D8E6F"+ _' -1450 to -1401
"63CB646CF71134D6B6000E1A08EE87CC336B64DC6711970C3A"+ _' -1400 to -1351
"93863E1A4B5B9B238A233E5BF27093E1D498AC7D23B11A2AD0"+ _' -1350 to -1301
"823C32FA87C6B6CE6082E75EB9B3C700515E341A11338A68ED"+ _' -1300 to -1251
"4DA5E334B94E6BF74F70A11DE22D5ED8E56BBED8A5DC7D86FF"+ _' -1250 to -1201
"7B5EE3FB11A41119FA0C9860CE157C30CB60671DC72793CBD0"+ _' -1200 to -1151
"8533CB6A86CBDDA83E9B97D608B674160FBD8A119BD8A042DC"+ _' -1150 to -1101
"300045E52F300ED1D08E463CFA56F790A82D932A9ACB97515F"+ _' -1100 to -1051
"AA19004608B26E9BBAC0CF7DC8ED8979CF2BA90BE3459FF53B"+ _' -1050 to -1001
"711D82D8E6CBD68EB9CF389F7AD04532E5B230C781A89B62E2"+ _' -1000 to -951
"7DC07D000EEDCF79089C3CD5AA9B197EA18A5930A115E52E5A"+ _' -950 to -901
"CFF8D9B2F7159497989B2FF74D08D4DE632E2DD856C533ED9E"+ _' -900 to -851
"D815F66BA01D8A85F68A8715382F79A041DCFD5E70086C5741"+ _' -850 to -801
"1D23300C411D082D1D8A0CFF1AC826C3B3C386FD115300EA4B"+ _' -800 to -751
"D4D49A002F6227E63CF30C7C63063894816B5A5FFEB6345230"+ _' -750 to -701
"86A45B232AD6F3ADCF389FB381D0CB52F3C7C8DA5E30482BB1"+ _' -700 to -651
"E56719F5FADE229B275F5EACBED8E5A12FF7D8E75A045DD48E"+ _' -650 to -601
"A0825389BA0047864A1DEACC78B6BBA85F6260041D416004EE"+ _' -600 to -551
"E119FB600C382D5E1D8BA0EA0CA8A045D5ED04C4C0CBD08D0C"+ _' -550 to -501
"0F719F7C7D8A930CE52C642B641D2EAA1D4F6830CD52564FA1"+ _' -500 to -451
"EA5B672FFD46B9FFFBE522F34FEAADCF7D6042FEDCF3CBEB94"+ _' -450 to -401
"A978E626441A86838D45FB11DA934D0493A4E4ED23B5602ADE"+ _' -400 to -351
"A085A9CEA8644A8304649C68D8A2F7123CD97ADC3CFC7C6FBF"+ _' -350 to -301
"ADE2AADCF1E575A7CF98AEE9CBD271823264E2F253BBA853BF"+ _' -300 to -251
"75389ED9D9760C7D932A1874B1E745F3CDD8197DC7086640C4"+ _' -250 to -201
"DC575F3768AEE1FF6ACBBAC00CBFA9A825E66BAC383042F367"+ _' -200 to -151
"1563049BD8686F0C60860414E5BE9CFE793811DEA85B5FBE74"+ _' -150 to -101
"D645B225F20BAD644DC36F3CF3A82F75B67B5F9B626CFDD0C7"+ _' -100 to -51
"DE90E163C368653852A52EE8B1371FB579B9F3CFF7DC63CBAD"+ _' -50 to -1
"000000000000000000000000008A45F6E57C3041DEE90CF57C"+ _' 0 to 49
"04812608A856C75B7C7974BD5D82B6CDDCB2786CFBA10936F7"+ _' 50 to 99
"8EAC71D8A2F6E67EAC3EE5F7530008E375D4DA9323FF715EDD"+ _' 100 to 149
"8E1345307FA5DC3EA478275FDC38B6C90F45A8E215B3BA2B27"+ _' 150 to 199
"D8D45E116CEEF6C81A08A4D90BA86CE1941DC6F36FF904C75D"+ _' 200 to 249
"87A191E1F7CC33879B93022367D2B1307B9400C86F5711279E"+ _' 250 to 299
"D8BEDEE5F302649B19306382949345D1981A6089B644A6F719"+ _' 300 to 349
"712791632743759389C6E1ED49F755F6897D2B156B97193C98"+ _' 350 to 399
"6FBA68ED8221ED4D975FD8B93FE9FACBD8622F5B94E26E5349"+ _' 400 to 449
"BD4D45A0DC34DAD0000CD15F34004F2BD49415B2E3867D970C"+ _' 450 to 499
"770452E7DC9C7C3CA8EF60A9F386FEDC2EA27C86B111330C38"+ _' 500 to 549
"B2E152701168E15641597534DADCD97E681560671B16CA1B90"+ _' 550 to 599
"82E56CDC0E155742198197C07BD5975D822E175ED4DC9CBD08"+ _' 600 to 649
"9F78EE9B607119793FA0C34DDC7D675BA95B94113792687560"+ _' 650 to 699
"F306CF711F27CF7379E1DCB608EBD971638638EDEEE01DE1A2"+ _' 700 to 749
"B98ACC74193C1E5527C0B5E25BD4D971EA9F007FA936BE9BE1"+ _' 750 to 799
"D004982AA198225702B1EE68371BACFBE934CF77CF0BED04DE"+ _' 800 to 849
"A830871552B5E65F341142A4D51D86E5BB1AC340566FF1A4E2"+ _' 850 to 899
"38E90FBED049FACCA5264115A0860EDD8ED2381197115FBE52"+ _' 900 to 949
"AA93E5782326830B238ABD49B29F768AAD4D86CE2BE6B711C8"+ _' 950 to 999
"534DEE15E6CEED496BEB90E97D2F3ACBD4D6CBB16C2B1D4338"+ _' 1000 to 1049
"382EF38B5300893C90F863AC3CFB3C2CA562AD4113CFF3C981"+ _' 1050 to 1099
"BEDCB21E90304DC382E7151D990CB2BD8E33C75B17DC756E2E"+ _' 1100 to 1149
"D41D82D49E16B23C7149CB330CB64DC3CB56830CB60A41BD52"+ _' 1150 to 1199
"BA12FF30C5FB741564AD4C48FBEF238B944E55232A5B26A8AF"+ _' 1200 to 1249
"ADE162E63E5367DA8345A689CFBFA5BED2FFA416867CF3AD89"+ _' 1250 to 1299
"346C74DDCF3E15949F38A8EEF27A1D4DD6FFACBFB1B1225F6E"+ _' 1300 to 1349
"D1F2BB53086DD87D0E5349F986F1A433CF34BD10C8D4D45A2F"+ _' 1350 to 1399
"38AC79785681D5BA89B6ED1EE530456B79B2BAC8AC4055A12B"+ _' 1400 to 1449
"60BE7113BA4D9BF23ED55EA60893C1E193711D0419A4638E38"+ _' 1450 to 1499
"AADCF3468A4121532E6634D49708326B34174DA5279BDC08ED"+ _' 1500 to 1549
"08A1334A9FC2F341BE1D4D8A81E5D41E3F7463634334749822"+ _' 1550 to 1599
"FFBB22F7DC57D000268FF29C3ED4334D0CFB27415D938A0A5A"+ _' 1600 to 1649
"864E2AF2BD8A45E9793C3CFBE38DEA97528DD2AD082B686CFE"+ _' 1650 to 1699
"7D87970455B9638D045671566FFF97967D6040C3CB62EA08A2"+ _' 1700 to 1749
"B648FEBA0C75FB9EDD822B23AC7FE6089838579C57C8E94C81"+ _' 1750 to 1799
"19AC3045BA000CB194E2E30C30852BD4D23FE97BD52D4594DA"+ _' 1800 to 1849
"0812E56F82E1D0497B902A93A451A0A5AC757930C7EE9A821A"+ _' 1850 to 1899
"4512E19BB2ACB600C9CBD08BE5194DC30867D0C3CBD4D682D4"+ _' 1900 to 1949
"D0112BE5E5A7823C386C5F25FE9BED000C300829F3E1293CCB"+ _' 1950 to 1999
"1303A9FD9F9FBBAC2EAC442F7CC3345A9FEF7474DC4C4D8E5E"+ _' 2000 to 2049
"51E5F712306F7AD4B5F79B15B9F749FEFBDE2FB6327CCBD412"+ _' 2050 to 2099
"5A90F75B26689A82A82756086C37E27A5A04983F79B94D8D41"+ _' 2100 to 2149
"F3709F79B3868A00C3482FDD85640C71159F32F604153CF7DC"+ _' 2150 to 2199
"2A912E15A78A279D1986B75E9ED02375677C3CB29B2941CCA9"+ _' 2200 to 2249
"FE66BA64DC38BD971279B9BAEE9A004F68B2E5E275ED864151"+ _' 2250 to 2299
"9756F67D4DCFF71567970453A4D274FEED081907522308715D"+ _' 2300 to 2349
"56419986E1DC794D00004A972B23CFDD0826C389BAC3CF7DC0"+ _' 2350 to 2399
"862A049C9402BEEEE78EFF298A87D8930CF3A0838649BBD5B2"+ _' 2400 to 2449
"E1D4F6451985E6D97A9770004DEAC41D04D8190C3819723819"+ _' 2450 to 2499
"3EA0A193CF26F6B67522349E1A419A86CF76CF9343FF30CF75"+ _' 2500 to 2549
"E1159821D16042344E944E9CEA64D45E52341D63CA08300411"+ _' 2550 to 2599
"DC7D86BBACE1982E79B783CEE9716CBE5E19719B9AC7FE2797"+ _' 2600 to 2649
"E2BD8AE6B152B797E262604157D8D0DC7DE981DCB1D99CBBD5"+ _' 2650 to 2699
"F7064D26C7D1C8534DEA4598EFA5F70A1D8A9F3CE52F79044E"+ _' 2700 to 2749
"A5A2BA6829B6E9CFBE5A0C34EEA83C30CBD04A5257C497CF70"+ _' 2750 to 2799
"566710DC70C7E2B7CF8ED6449F7457A9B7C03A9FDD06BA64D6"+ _' 2800 to 2849
"44701AD00862A98DD4B9FCB949300875E70C8223083C70A5A4"+ _' 2850 to 2899
"987467D9386A01A5AA9FEDCB930C3879386FD1D2F6EDD8A8FB"+ _' 2900 to 2949
"25A97371D86C2F3CF85F1A0E22F7994113759BE115B2BD86AC"+ _' 2950 to 2999
"3497338FF64DC77C3A5AA97C3838E371DC782741E9C37B5F75"+ _' 3000 to 3049
"A8EA4D95FA6008ADC344DCB7DC345FB2934FEA2BDA5A46B9B2"+ _' 3050 to 3099
"A00164D45AD491674A4E2786CB60468E1151D9DDC1E98E986B"+ _' 3100 to 3149
"BE17CBB6F30455AD64556445BD0049367575EAEA02F2A4D86C"+ _' 3150 to 3199
"D93223ED89F382B5DD8D045A0274D0404A4DC82649F16FF7DC"+ _' 3200 to 3249
"63C1E1D8E33CA45F34EE9E53F898FBDC71111704D8E70DC702"+ _' 3250 to 3299
"FB3711FB9268BE9A008A4D8EB57457604D223AC3C344166B62"+ _' 3300 to 3349
"27EA45266343F302A9638D0423C12AEEDE2B115F30FA455A9A"+ _' 3350 to 3399
"4D4908E5DC89F7A52F38E193000227EA08646490B66FF6F3CB"+ _' 3400 to 3449
"1AE908E579FEDCEDC049BE6A82D522B230457523866F786CBA"+ _' 3450 to 3499
"D41594DE90C3C3A8ED86E9067DEAC081FFA83F386F483E1D82"+ _' 3500 to 3549
"F98E970A9708240972FA74D0000CF6F7498B908268E1DCF386"+ _' 3550 to 3599
"FBA945607111DAD00825B534DC826E5B970600A043B66341D4"+ _' 3600 to 3649
"F2F90349BACFB5E1B2EF230C3E97B160FF363045EB2608AAD8"+ _' 3650 to 3699
"797C34B5BBD11238300A8ED45B2267D493A87A19716B1D43F3"+ _' 3700 to 3749
"A52D04F378373796FBDC3CD1938E111F78574193049BA1DCB7"+ _' 3750 to 3799
"CFB345971B1AA0DA82AC78A9CEE523B194E63EE9E9B1538DD2"+ _' 3800 to 3849
"B9ED415689FB56F6BAE2E7D83B2FF76457979FA515BB5FF306"+ _' 3850 to 3899
"B2BA1A97974B53C745BE6E97933C71FE2756EDD2A156FD9799"+ _' 3900 to 3949
"08BE59B6C7D0453CB5225B653CE3C7A15FBE6000C3C795B1BE"+ _' 3950 to 3999
"1DC7EA49E121DD8A8B2E5E2D4EE2B3F36BEEBD99416C3A8A1D"+ _' 4000 to 4049
"CD97CFF630CB53A49E9C71D82BD453CF5F2FFFD9764525B11D"+ _' 4050 to 4099
"0A5A860E16A413F34155BD556CD5A086863A46304193E678A5"+ _' 4100 to 4149
"D0462E9CFBB530B6A8F7D4B151A8752B115FF36FF715DD02E6"+ _' 4150 to 4199
"2EE37597D5527E268A049000C75E7DCB2FBA93827CB253F2E1"+ _' 4200 to 4249
"386DD8116E1538E52EEFA13F363493C3004115AD8A478A37DC"+ _' 4250 to 4299
"483675523E9B5E6E5B970049A0CCBE25606F345A0C8FAD8A91"+ _' 4300 to 4349
"2A42A0AD8564A5EF7C6F76825AD2702373FF38E9223497E2A8"+ _' 4350 to 4399
"FE119A8649B300086349F267B5326CFFFBA55A556CFF385E17"+ _' 4400 to 4449
"11271C41EA41DC63C5743B27ADCA4926F67D0040118BD8602E"+ _' 4450 to 4499
"642F26B6487CB6041C34934CFB79BF60CE192EAC3CB2FBA9FF"+ _' 4500 to 4549
"3008D49C745341FEE5EE64D8A142663004986F163CDD897575"+ _' 4550 to 4599
"68EDC3C40D45AEA89C2B797CFF16B7D000468A0C04195F3FAC"+ _' 4600 to 4649
"30CBBD97902E23410110199459FB23413FFF387523C78332EE"+ _' 4650 to 4699
"74E11DCD9B9FF7DC86B15B15AA5A274526EA49894C8DAD4782"+ _' 4700 to 4749
"DD415AEE1BD10D493C7704704BDD4D0CB2981386982D860E6F"+ _' 4750 to 4799
"57C6381D149CF79456E5E119837D2FFD4E4E1B902AD0C7645F"+ _' 4800 to 4849
"B636753C386F67D4DCACBBA0E5B97D8E94641D04E6FE52DD87"+ _' 4850 to 4899
"1174D2B1BE5290C3413B2B2FF3CF71D4D" ' 4900 to 4932
Select Case n
Case <-4931
' Except for n=-4933 squaring 10^(n/2) is OK.
x=10^(n/2) ' To create denormals
If n=-4933 Then x=x+6E-2483 ' Only adjustment needed
x=x*x
If n<-4950 Then Function=x: Exit Function ' Zero
Case -4931 To +4931
x=Val("1E"+Format$(n)) ' Maximum significand error = 4
Case >+4931
x=1E+4096*10^(n-4096)
If n>+4932 Then Function=x: Exit Function ' Infinity
End Select
t=d(0) ' Least significant 32 bits of significand
d(0)=Val("&H"+Left$(Hex$(t,8),7)+Mid$(LSNibble,n+4951,1))
If Abs(d(0)-t)=>8 Then d(0)=d(0)-&H10*Sgn(d(0)-t) ' If carry/borrow
Function=x
End Function
Val("1E"+Format$(n)) = 1En
but is considerably faster. A table of powers of five rather than ten was used so that extreme negative powers remained within extended-precision limits with just the delay of a !fscale command at the end. Please note the following (timings based on Pentium 4):
Avoids compile-time error #473 (Ivalid numeric format).
Methodology is equivalent to PB when PB gives valid results.
Takes up about 300 bytes.
Significand accuracy +/- 4 over entire range.
Up to 4000+ times more accurate and 75% faster than Exp(n*Log(10)).
Up to 2000+ times more accurate and 25% faster than Exp10(n).
Up to 36 times more accurate and four times as fast as 10^n.
Same accuracy and eleven times faster than Val("1E"+Format$(n)).
Returns zero if n<=-4951; infinity if n=>+4933.
Returns denormals for -4950 <= n <= -4932.
Equals 1En and Val("1E"+Format$(n)) for -4932 <= n <= 4931
Corrects PB error where Val("1E4932")=2^127.
Best values for -29 <= n <= +32 and exact for 0 <= n <= +27.
Use in combination with Exp10_() (see below) as:
Val1E(Fix(x))*Exp10_(Float(x)) or (perhaps better)
Val1E(CLng(x))*Exp10_(x-CLng(x)).
Code:
Function Val1E(ByVal n As Long) As Ext
' Calculates 1En = Val("1E"+Format$(n)).
!mov eax,n ' eax=n
!mov edx,eax ' If n<0: If n=>0:
!sar edx,31 ' edx=-1=(n<0) edx=0=(n<0)
!xor eax,edx ' eax=Not(n) eax=n
!sub eax,edx ' eax=Not(n)-(-1)=|n| eax=n-0=|n|
!jz Zero ' Jump if n=0
!test eax,&HFFFFE000 ' Check if |n| greater than &H1FFF=8191
!db &H3E ' Branch taken hint
!jz PrepareFPU ' Jump if |n| within bounds
!mov eax,&H00001FFF ' eax=Min&(|n|,8191)
PrepareFPU:
!fild n ' ST(0)=n ' For !fscale at end
!fld1 ' ST(0)=1 ST(1)=n
!lea ebx,SelectedPowers_Table ' ebx=pointer to next factor
!jmp LoopEntry ' Skip over initial !add
#Align 16 ' Frequent jumps to here; Needs current version.
LoopHead:
!add ebx,10 ' Pointer to next potential factor
LoopEntry:
!shr eax,1 ' CF=next bit
!ja LoopHead ' Jump if CF=0 and ZF=0
' Here if:
' CF=1 and ZF=0 so multiply and loop
' CF=1 and ZF=1 so multiply and exit loop
' CF=0 and ZF=1 This cannot happen as n<>0.
!fld TByte [ebx] ' ST(0)=Factor ST(1)=Product ST(2)=n
!fmulp ST(1),ST(0) ' ST(0)=New_Product ST(1)=n
!jnz LoopHead
' ST(0)=Product If n<0: If n>0:
!or edx,edx ' edx=-1 ZF=0 edx=0 ZF=1
!jz AllDone ' Don't jump Jump
Zero: ' Here if n=0 with ZF=1 and empty FPU
!fld1 ' ST(0)=1.0 ST(1)=Product ST(2)=n
!jz Last ' Jump if zero
!fdivrp ST(1),ST(0) ' After pop, ST(0)=1.0/Product ST(1)=n
AllDone: ' Now ST(0)=5^n
!fscale ' ST(0)=10^n ST(1)=n
!fstp ST(1) ' After pop ST(0)=10^n
Last:
!fstp Function ' Needs current PB version.
Exit Function
#Align 2
SelectedPowers_Table:
!dw &H0000,&H0000,&H0000,&HA000, &H4001 ' 1E1/2^1=5^(2^0) Exact
!dw &H0000,&H0000,&H0000,&HC800, &H4003 ' 1E2/2^2=5^(2^1) Exact
!dw &H0000,&H0000,&H0000,&H9C40, &H4008 ' 1E4/2^4=5^(2^2) Exact
!dw &H0000,&H0000,&H2000,&HBEBC, &H4011 ' 1E8/2^8=5^(2^3) Exact
!dw &H0000,&H0400,&HC9BF,&H8E1B, &H4024 ' 1E16/2^16=5^(2^4) Exact
' The rest of the table was calculated to high precision first.
!dw &HB59E,&H2B70,&HADA8,&H9DC5, &H4049 ' 1E32/2^32=5^(2^5)
!dw &HA6D5,&HFFCF,&H1F49,&HC278, &H4093 ' 1E64/2^64=5^(2^6)
!dw &H8CE0,&H80E9,&H47C9,&H93BA, &H4128 ' 1E128/2^128=5^(2^7)
!dw &HDE8E,&H9DF9,&HEBFB,&HAA7E, &H4251 ' 1E256/2^256=5^(2^8)
!dw &H91C7,&HA60E,&HA0AE,&HE319, &H44A3 ' 1E512/2^512=5^(2^9)
!dw &H0C17,&H8175,&H7586,&HC976, &H4948 ' 1E1024/2^1024=5^(2^10)
!dw &H5DE5,&HC53D,&H3B5D,&H9E8B, &H5292 ' 1E2048/2^2048=5^(2^11)
!dw &H979B,&H8A20,&H5202,&HC460, &H6525 ' 1E4096/2^4096=5^(2^12)
End Function
10^x = 2^(x*Log2(10)) = 2^i*2^d = ((2^d-1)+1) * 2^i
Code:
Function Exp10_(ByVal x As Ext) As Ext
' Runs about twice as fast as Exp10(x).
!fldl2t ' ST(0)=Log2(10)
!fld x ' ST(0)=x ST(1)=Log2(10)
!fmulp ST(1),ST(0) ' ST(0)=x*Log2(10)
!fld ST(0) ' ST(0)=x*Log2(10) ST(1)=x*Log2(10)
!frndint ' ST(0)=i (integer) ST(1)=x*Log2(10)
!fsub ST(1),ST(0) ' ST(0)=i ST(1)=d where |d|<=0.5
!fxch ST(1) ' ST(0)=d ST(1)=i
!f2xm1 ' ST(0)=2^d-1 ST(1)=i
!fld1 ' ST(0)=1.0 ST(1)=2^d-1 ST(2)=i
!faddp ST(1),ST(0) ' ST(0)=2^d ST(1)=i
!fscale ' ST(0)=2^i*2^d=2^x ST(1)=i
!fstp ST(1) ' ST(0)=2^x
!fstp Function
End Function
1. If just double-precision accuracy is desired, any method will work, and Exp10_(x) is the fastest.
2. If x ia an integer use Val1E(x) for its speed and accuracy. If maximum accuracy is desired and size is not a concern create a constant array within the code (using !dw pseudo-ops). The latter method, with a complete table, requires almost 97 kilobytes. A BASIC array (including zero and infinity extremes) can be dimensioned thus:
Dim MyArray(-4951 To 4933) As Static Ext At CodePtr(MyLabel)
In assembly, the starting address of the table can be accessed with:
!lea ebx,MyLabel ' Load register with equivalent address.
3. If x is not an integer but is of small maqnitude (less than 5, but preferably less than 2) use Exp10_(x).
4. If x can be split into an integer (integer) and a fractional part (delta) where the latter is more accurate than using x alone (say, x is a rational fraction or a number with known additional digits), use the product of Val1E(integer) and Exp10_(delta).
5. If x may or may not be an integer, use 10^x. Recognize however this is the least accurate method if x is a non-integer but a fairly accurate method otherwise. The compiled program makes the decision.
6. Always remember that as |x| gets larger, the relative error increases proportionately unless the programmer can use some means to circumvent it.
Note added in 1/2/09 edit:
The use of Exp2() (and probably all other power functions with a non-integer exponent) has an issue in and of itself. It appears that although PBCC is written for 32-bit machines, the Exp2() function used by PB uses a form of the FPU command f2xm1 that was written for 16-bit machines. This form limits the ST(0) value to 0<=ST(0)<=0.5. For further details see my post below and Intel, "IA-32 Intel Architecture Software Developer's Manual", Volume 3, Section 18.14.7.10.
_______________
Mark Longley-Cook
Comment