Pathological floating point problems: Difference between revisions
m (reduce the font size for words in the Siegfried Rump's formula, added whitespace to the first formulae.) |
m (→The Chaotic Bank Society: added a comment to show the number of decimal digits in an assignment.) |
||
| Line 665: | Line 665: | ||
numeric digits digs /*have REXX use "digs" decimal digits. */ |
numeric digits digs /*have REXX use "digs" decimal digits. */ |
||
$=2.71828182845904523536028747135266249775724709369995957496696762772407663035354759457138217852516642742746639193200305992181741359662904357290033429526 |
$=2.71828182845904523536028747135266249775724709369995957496696762772407663035354759457138217852516642742746639193200305992181741359662904357290033429526 |
||
$=$ - 1 /*and subtract one 'cause that's that. */ |
$=$ - 1 /*and subtract one 'cause that's that. */ /* [↑] 150 decimal digits of e */ |
||
/* [↑] value of newly opened account. */ |
/* [↑] value of newly opened account. */ |
||
do n=1 for y /*compute the value of the account/year*/ |
do n=1 for y /*compute the value of the account/year*/ |
||
Revision as of 22:45, 9 July 2016
You are encouraged to solve this task according to the task description, using any language you may know.
Most programmers are familiar with the inexactness of floating point calculations in a binary processor.
The classic example being:
0.1 + 0.2 = 0.30000000000000004
In many situations the amount of error in such calculations is very small and can be overlooked or eliminated with rounding.
There are pathological problems however, where seemingly simple, straight-forward calculations are extremely sensitive to even tiny amounts of imprecision.
This task's purpose is to show how your language deals with such classes of problems.
A sequence that seems to converge to a wrong limit.
Consider the sequence:
- v1 = 2
- v2 = -4
- vn = 111 - 1130 / vn-1 + 3000 / (vn-1 * vn-2)
As n grows larger, the series should converge to 6 but small amounts of error will cause it to approach 100.
- Task 1
Display the values of the sequence where n = 3, 4, 5, 6, 7, 8, 20, 30, 50 & 100 to at least 16 decimal places.
n = 3 18.5
n = 4 9.378378
n = 5 7.801153
n = 6 7.154414
n = 7 6.806785
n = 8 6.5926328
n = 20 6.0435521101892689
n = 30 6.006786093031205758530554
n = 50 6.0001758466271871889456140207471954695237
n = 100 6.000000019319477929104086803403585715024350675436952458072592750856521767230266
- Task 2
The Chaotic Bank Society is offering a new investment account to their customers.
You first deposit $e - 1 where e is 2.7182818... the base of natural logarithms.
After each year, your account balance will be multiplied by the number of years that have passed, and $1 in service charges will be removed.
So ...
- after 1 year, your balance will be multiplied by 1 and $1 will be removed for service charges.
- after 2 years your balance will be doubled and $1 removed.
- after 3 years your balance will be tripled and $1 removed.
- ...
- after 10 years, multiplied by 10 and $1 removed, and so on.
What will your balance be after 25 years?
Starting balance: $e-1 Balance = (Balance * year) - 1 for 25 years Balance after 25 years: $0.0399387296732302
- Task 3, extra credit
Siegfried Rump's example. Consider the following function, designed by Siegfried Rump in 1988.
- f(a,b) = 333.75b6 + a2( 11a2b2 - b6 - 121b4 - 2 ) + 5.5b8 + a/(2b)
- compute f(a,b) where a=77617.0 and b=33096.0
- f(77617.0, 33096.0) = -0.827396059946821
Demonstrate how to solve at least one of the first two problems, or both, and the third if you're feeling particularly jaunty.
- See also;
- Floating-Point Arithmetic Section 1.3.2 Difficult problems.
360 Assembly
The system/360 hexadecimal single precision floating point format is known to its weakness
in precision. A lot of more precise formats have been added since.
A sequence that seems to converge to a wrong limit
<lang 360asm>* Pathological floating point problems 03/05/2016
PATHOFP CSECT
USING PATHOFP,R13
SAVEAR B STM-SAVEAR(R15)
DC 17F'0'
STM STM R14,R12,12(R13)
ST R13,4(R15)
ST R15,8(R13)
LR R13,R15
LE F0,=E'2'
STE F0,U u(1)=2
LE F0,=E'-4'
STE F0,U+4 u(2)=-4
LA R6,3 n=3
LA R7,U+4 @u(n-1)
LA R8,U @u(n-2)
LA R9,U+8 @u(n)
LOOPN CH R6,=H'100' do n=3 to 100
BH ELOOPN
LE F4,0(R7) u(n-1)
LE F2,=E'1130' 1130
DER F2,F4 1130/u(n-1)
LE F0,=E'111' 111
SER F0,F2 111-1130/u(n-1)
LE F2,0(R7) u(n-1)
LE F4,0(R8) u(n-2)
MER F2,F4 u(n-1)*u(n-2)
LE F6,=E'3000' 3000
DER F6,F2 3000/(u(n-1)*u(n-2))
AER F0,F6 111-1130/u(n-1)+3000/(u(n-1)*u(n-2))
STE F0,0(R9) store into u(n)
XDECO R6,PG+0 n
LE F0,0(R9) u(n)
LA R0,3 number of decimals
BAL R14,FORMATF format(u(n),'F13.3')
MVC PG+12(13),0(R1) put into buffer
XPRNT PG,80 print buffer
LA R6,1(R6) n=n+1
LA R7,4(R7) @u(n-1)
LA R8,4(R8) @u(n-2)
LA R9,4(R9) @u(n)
B LOOPN
ELOOPN L R13,4(0,R13)
LM R14,R12,12(R13)
XR R15,R15
BR R14
COPY FORMATF
LTORG
PG DC CL80' ' buffer U DS 100E
YREGS
YFPREGS
END PATHOFP</lang>
The divergence comes very soon.
- Output:
3 18.500
4 9.378
5 7.801
6 7.154
7 6.805
8 6.578
9 6.235
10 2.915
11 -111.573
12 111.905
13 100.661
14 100.040
15 100.002
16 100.000
17 100.000
18 100.000
... 100.000
Excel
A sequence that seems to converge to a wrong limit
<lang> A1: 2
A2: -4 A3: =111-1130/A2+3000/(A2*A1) A4: =111-1130/A3+3000/(A3*A2) ...</lang>
The result converges to the wrong limit!
- Output:
A 1 2 2 -4 3 18.5 4 9.378378378 5 7.801152738 6 7.154414481 7 6.806784737 8 6.592632769 9 6.449465934 10 6.348452061 11 6.274438663 12 6.218696769 13 6.175853856 14 6.142627170 15 6.120248705 16 6.166086560 17 7.235021166 18 22.06207846 19 78.57557489 20 98.34950312 21 99.89856927 22 99.99387099 23 99.99963039 24 99.99997773 25 99.99999866 26 99.99999992 27 100 ... 30 100 ... 40 100 ... 50 100 ... 100 100
Go
<lang go>package main
import (
"fmt" "math/big"
)
func main() {
sequence() bank() rump()
}
func sequence() {
// exact computations using big.Rat
var v, v1 big.Rat
v1.SetInt64(2)
v.SetInt64(-4)
n := 2
c111 := big.NewRat(111, 1)
c1130 := big.NewRat(1130, 1)
c3000 := big.NewRat(3000, 1)
var t2, t3 big.Rat
r := func() (vn big.Rat) {
vn.Add(vn.Sub(c111, t2.Quo(c1130, &v)), t3.Quo(c3000, t3.Mul(&v, &v1)))
return
}
fmt.Println(" n sequence value")
for _, x := range []int{3, 4, 5, 6, 7, 8, 20, 30, 50, 100} {
for ; n < x; n++ {
v1, v = v, r()
}
f, _ := v.Float64()
fmt.Printf("%3d %19.16f\n", n, f)
}
}
func bank() {
// balance as integer multiples of e and whole dollars using big.Int
var balance struct{ e, d big.Int }
// initial balance
balance.e.SetInt64(1)
balance.d.SetInt64(-1)
// compute balance over 25 years
var m, one big.Int
one.SetInt64(1)
for y := 1; y <= 25; y++ {
m.SetInt64(int64(y))
balance.e.Mul(&m, &balance.e)
balance.d.Mul(&m, &balance.d)
balance.d.Sub(&balance.d, &one)
}
// sum account components using big.Float
var e, ef, df, b big.Float
e.SetPrec(100).SetString("2.71828182845904523536028747135")
ef.SetInt(&balance.e)
df.SetInt(&balance.d)
b.Add(b.Mul(&e, &ef), &df)
fmt.Printf("Bank balance after 25 years: $%.2f\n", &b)
}
func rump() {
a, b := 77617., 33096.
fmt.Printf("Rump f(%g, %g): %g\n", a, b, f(a, b))
}
func f(a, b float64) float64 {
// computations done with big.Float with enough precision to give
// a correct answer.
fp := func(x float64) *big.Float { return big.NewFloat(x).SetPrec(128) }
a1 := fp(a)
b1 := fp(b)
a2 := new(big.Float).Mul(a1, a1)
b2 := new(big.Float).Mul(b1, b1)
b4 := new(big.Float).Mul(b2, b2)
b6 := new(big.Float).Mul(b2, b4)
b8 := new(big.Float).Mul(b4, b4)
two := fp(2)
t1 := fp(333.75)
t1.Mul(t1, b6)
t21 := fp(11)
t21.Mul(t21.Mul(t21, a2), b2)
t23 := fp(121)
t23.Mul(t23, b4)
t2 := new(big.Float).Sub(t21, b6)
t2.Mul(a2, t2.Sub(t2.Sub(t2, t23), two))
t3 := fp(5.5)
t3.Mul(t3, b8)
t4 := new(big.Float).Mul(two, b1)
t4.Quo(a1, t4)
s := new(big.Float).Add(t1, t2)
f64, _ := s.Add(s.Add(s, t3), t4).Float64()
return f64
}</lang>
- Output:
n sequence value 3 18.5000000000000000 4 9.3783783783783790 5 7.8011527377521617 6 7.1544144809752490 7 6.8067847369236327 8 6.5926327687044388 20 6.0435521101892693 30 6.0067860930312058 50 6.0001758466271875 100 6.0000000193194776 Bank balance after 25 years: $0.04 Rump f(77617, 33096): -0.8273960599468214
J
A sequence that seems to converge to a wrong limit.
Implementation of vn:
<lang J> vn=: 111 +(_1130 % _1&{) + (3000 % _1&{ * _2&{)</lang>
Example using IEEE-754 floating point:
<lang J> 3 21j16 ":"1] 3 4 5 6 7 8 20 30 50 100 ([,.{) (,vn)^:100(2 _4)
3 9.3783783783783861 4 7.8011527377522611 5 7.1544144809765555 6 6.8067847369419638 7 6.5926327689743687 8 6.4494659378910058 20 99.9934721906444960 30 100.0000000000000000 50 100.0000000000000000
100 100.0000000000000000</lang>
Example using exact arithmetic:
<lang J> 3 21j16 ":"1] 3 4 5 6 7 8 20 30 50 100 ([,.{) (,vn)^:100(2 _4x)
3 9.3783783783783784 4 7.8011527377521614 5 7.1544144809752494 6 6.8067847369236330 7 6.5926327687044384 8 6.4494659337902880 20 6.0360318810818568 30 6.0056486887714203 50 6.0001465345613879
100 6.0000000160995649</lang>
The Chaotic Bank Society
Let's start this example by using exact arithmetic, to make sure we have the right algorithm. We'll go a bit overboard, in representing e, so we don't have to worry too much about that.
<lang J> e=: +/%1,*/\1+i.100x
81j76":e 2.7182818284590452353602874713526624977572470936999595749669676277240766303535 21j16":+`*/,_1,.(1+i.-25),e 0.0399387296732302</lang>
(Aside: here, we are used the same mechanism for adding -1 to e that we are using to add -1 to the product of the year number and the running balance.)
Next, we will use for e, to represent the limit of what can be expressed using 64 bit IEEE 754 floating point.
<lang J> 31j16":+`*/,_1,.(1+i.-25),6157974361033%2265392166685x _2053975868590.1852178761057505</lang>
That's clearly way too low, so let's try instead using for e <lang J> 31j16":+`*/,_1,.(1+i.-25),6157974361034%2265392166685x
4793054977300.3491517765983265</lang>
So, our problem seems to be that there's no way we can express enough bits of e, using 64 bit IEEE-754 floating point arithmetic. Just to confirm:
<lang J> 1x1 2.71828
+`*/,_1,.(1+i.-25),1x1
_2.24237e9</lang>
Now let's take a closer look using our rational approximation for e:
<lang J> 21j16":+`*/,_1,.(1+i.-25),+/%1,*/\1+i.40x
0.0399387296732302 21j16":+`*/,_1,.(1+i.-25),+/%1,*/\1+i.30x 0.0399387277260840 21j16":+`*/,_1,.(1+i.-25),+/%1,*/\1+i.26x 0.0384615384615385 21j16":+`*/,_1,.(1+i.-25),+/%1,*/\1+i.25x 0.0000000000000000 21j16":+`*/,_1,.(1+i.-25),+/%1,*/\1+i.24x _1.0000000000000000</lang>
Things go haywire when our approximation for e uses the same number of terms as our bank's term. So, what does that look like, in terms of precision?
<lang J> 41j36":+/%1,*/\1+i.26x
2.718281828459045235360287471257428715 41j36":+/%1,*/\1+i.25x 2.718281828459045235360287468777832452 41j36":+/%1,*/\1+i.24x 2.718281828459045235360287404308329608</lang>
In other words, we go astray when our approximation for e is inaccurate in the 26th position after the decimal point. But IEEE-754 floating point arithmetic can only represent approximately 16 decimal digits of precision.
Siegfried Rump's example.
Again, we use exact arithmetic to see if we have the algorithm right. That said, we'll also do this in small steps, to make sure we're being exact every step of the way, and to keep from building overly long lines:
<lang J>rump=:4 :0
NB. enforce exact arithmetic add=. +&x: sub=. -&x: mul=. *&x: div=. %&x:
a=. x a2=. a mul a
b=. y b2=. b mul b b4=. b2 mul b2 b6=. b2 mul b4 b8=. b4 mul b4
c333_75=. 1335 div 4 NB. 333.75 term1=. c333_75 mul b6
t11a2b2=. 11 mul a2 mul b2 tnb6=. 0 sub b6 tn121b4=. 0 sub 121 mul b4 term2=. a2*(t11a2b2 + tnb6 + tn121b4 sub 2)
c5_5=. 11 div 2 NB. 5.5 term3=. c5_5 mul b8
term4=. a div 2 mul b
term1 add term2 add term3 add term4
)</lang>
Example use:
<lang J> 21j16": 77617 rump 33096
_0.8273960599468214</lang>
Note that replacing the definitions of add, sub, div, mul with implementations which promote to floating point gives a very different result:
<lang J> 77617 rump 33096 _1.39061e21</lang>
But given that b8 is
<lang J> 33096^8 1.43947e36</lang>
we're exceeding the limits of our representation here, if we're using 64 bit IEEE-754 floating point arithmetic.
PARI/GP
Task 1: Define recursive function V(n): <lang parigp>V(n,a=2,v=-4.)=if(n < 3,return(v));V(n--,v,111-1130/v+3000/(v*a))</lang> In order to work set precision to at least 200 digits:
\p 200: realprecision = 211 significant digits (200 digits displayed) V(50): 6.000175846627187188945614020747195469523735177... V(100): 6.0000000193194779291040868034035857150243506754369524580725927508565217672302663412282...
Task 2: Define function balance(deposit,years): <lang parigp>balance(d,y)=d--;for(n=1,y,d=d*n-1);d</lang>
Output balance(exp(1), 25):
0.039938729673230208903...
Task 3: Define function f(a,b): <lang parigp>f(a,b)=333.75*b^6+a*a*(11*a*a*b*b-b^6-121*b^4-2)+5.5*b^8+a/(2*b)</lang>
Output:
f(77617.0,33096.0): -0.827396059946821368141165...
Perl 6
The simple solution to doing calculations where floating point numbers exhibit pathological behavior is: don't do floating point calculations. :-) Perl 6 is just as susceptible to floating point error as any other C based language, however, it offers built-in rational Types; where numbers are represented as a ratio of two integers. For normal precision it uses Rats - accurate to 1/2^64, and for arbitrary precision, FatRats, which can grow as large as available memory. Rats don't require any special special setup to use. Any decimal number within its limits of precision is automatically stored as a Rat. FatRats require explicit coercion and are "sticky". Any FatRat operand in a calculation will cause all further results to be stored as FatRats. <lang perl6>say '1st: Convergent series'; my @series = 2.FatRat, -4, { 111 - 1130 / $^v + 3000 / ( $^v * $^u ) } ... *; for flat 3..8, 20, 30, 50, 100 -> $n {say "n = {$n.fmt("%3d")} @series[$n-1]"};
say "\n2nd: Chaotic banking society"; sub postfix:<!> (Int $n) { [*] 2..$n } # factorial operator my $years = 25; my $balance = sum map { 1 / FatRat.new($_!) }, 1 .. $years + 15; # Generate e-1 to sufficient precision with a Taylor series put "Starting balance, \$(e-1): \$$balance"; for 1..$years -> $i { $balance = $i * $balance - 1 } printf("After year %d, you will have \$%1.16g in your account.\n", $years, $balance);
print "\n3rd: Rump's example: f(77617.0, 33096.0) = "; sub f (\a, \b) { 333.75*b⁶ + a²*( 11*a²*b² - b⁶ - 121*b⁴ - 2 ) + 5.5*b⁸ + a/(2*b) } say f(77617.0, 33096.0).fmt("%0.16g");</lang>
- Output:
1st: Convergent series n = 3 18.5 n = 4 9.378378 n = 5 7.801153 n = 6 7.154414 n = 7 6.806785 n = 8 6.5926328 n = 20 6.0435521101892689 n = 30 6.006786093031205758530554 n = 50 6.0001758466271871889456140207471954695237 n = 100 6.000000019319477929104086803403585715024350675436952458072592750856521767230266 2nd: Chaotic banking society Starting balance, $(e-1): $1.7182818284590452353602874713526624977572470936999 After year 25, you will have $0.0399387296732302 in your account. 3rd: Rump's example: f(77617.0, 33096.0) = -0.827396059946821
Python
Task 1: Muller's sequence
Using rational numbers via standard library fractions
<lang Python>from fractions import Fraction
def muller_seq(n:int) -> float:
seq = [Fraction(0), Fraction(2), Fraction(-4)]
for i in range(3, n+1):
next_value = (111 - 1130/seq[i-1]
+ 3000/(seq[i-1]*seq[i-2]))
seq.append(next_value)
return float(seq[n])
for n in [3, 4, 5, 6, 7, 8, 20, 30, 50, 100]:
print("{:4d} -> {}".format(n, muller_seq(n)))</lang>
- Output:
3 -> 18.5 4 -> 9.378378378378379 5 -> 7.801152737752162 6 -> 7.154414480975249 7 -> 6.806784736923633 8 -> 6.592632768704439 20 -> 6.043552110189269 30 -> 6.006786093031206 50 -> 6.0001758466271875 100 -> 6.000000019319478
Task 2: The Chaotic Bank Society
Using decimal numbers with a high precision
<lang Python>from decimal import Decimal, getcontext
def bank(years:int) -> float:
"""
Warning: still will diverge and return incorrect results after 250 years
the higher the precision, the more years will cover
"""
getcontext().prec = 500
# standard math.e has not enough precision
e = Decimal('2.718281828459045235360287471352662497757247093699959574966967627724076630353547594571382178525166427427466391932003059921817413596629043572900334295260595630738132328627943490763233829880753195251019011573834187930702154089149934884167509244761460668082264800168477411853742345442437107539077744992069551702761838606261331384583000752044933826560297606737113200709328709127443747047230696977209310141692836819025515108657463772111252389784425056953696770785449969967946864454905987931636889230098793127736178215424999229576351')
decimal_balance = e - 1
for year in range(1, years+1):
decimal_balance = decimal_balance * year - 1
return(float(decimal_balance))
print("Bank balance after 25 years = ", bank(25))</lang>
- Output:
Bank balance after 25 years = 0.03993872967323021
but, still incorrectly diverging after some time, aprox. 250 years <lang Python>for year in range(200, 256, 5):
print(year, '->', bank(year))
</lang>
- Output:
200 -> 0.004999875631110097 205 -> 0.004877933277184028 210 -> 0.004761797301186607 215 -> 0.0046510626428896236 220 -> 0.004545361061789591 225 -> 0.0044443570465329246 230 -> 0.004347744257820075 235 -> 0.004255242425346535 240 -> 0.004166594632576723 245 -> 0.004081564933953891 250 -> 0.003999846590933889 255 -> -92939.78784907148
Task 3: Siegfried Rump's example
Using rational numbers via standard library fractions
<lang Python>from fractions import Fraction
def rump(generic_a, generic_b) -> float:
a = Fraction('{}'.format(generic_a))
b = Fraction('{}'.format(generic_b))
fractional_result = Fraction('333.75') * b**6 \
+ a**2 * ( 11 * a**2 * b**2 - b**6 - 121 * b**4 - 2 ) \
+ Fraction('5.5') * b**8 + a / (2 * b)
return(float(fractional_result))
print("Rump(77617, 33096) = ", rump(77617.0, 33096.0)) </lang>
- Output:
Rump(77617, 33096) = -0.8273960599468214
REXX
The REXX language uses character-based arithmetic. So effectively, it looks, feels, and tastes like decimal floating point
(implemented in software).
So, the only (minor) problem is how many decimal digits should be used to solve these pathological floating point problems.
A little extra boilerplate code was added to support the specification of how many decimal digits that should be used for the
calculations, as well how many decimal digits (past the decimal point) should be displayed.
A sequence that seems to converge to a wrong limit
<lang rexx>/*REXX pgm (pathological FP problem): a sequence that seems to converge to a wrong limit*/ parse arg digs show . /*obtain optional arguments from the CL*/ if digs== | digs=="," then digs=150 /*Not specified? Then use the default.*/ if show== | show=="," then show= 20 /* " " " " " " */ numeric digits digs /*have REXX use "digs" decimal digits. */
- = 2 4 5 6 7 8 9 20 30 50 100 /*the indices to display value of V.n */
fin=word(#, words(#) ) /*find the last (largest) index number.*/ w=length(fin) /* " " length (in dec digs) of FIN.*/ v.1= 2 /*the value of the first V element. */ v.2=-4 /* " " " " second " " */
do n=3 to fin; nm1=n-1; nm2=n-2 /*compute some values of the V elements*/
v.n=111 - 1130/v.nm1 + 3000/(v.nm1*v.nm2) /* " a value of a " element.*/
if wordpos(n, #)\==0 then say 'v.'left(n, w) "=" format(v.n, , show)
end /*n*/ /*display SHOW digs past the dec. point*/
/*stick a fork in it, we're all done. */</lang>
output when using the default inputs:
v.4 = 9.37837837837837837838 v.5 = 7.80115273775216138329 v.6 = 7.15441448097524935353 v.7 = 6.80678473692363298394 v.8 = 6.59263276870443839274 v.9 = 6.44946593379028796887 v.20 = 6.04355211018926886778 v.30 = 6.00678609303120575853 v.50 = 6.00017584662718718895 v.100 = 6.00000001931947792910
The Chaotic Bank Society
To be truly accurate, the number of decimal digits for e (the $ variable first value) should have 150 decimal
digits (or whatever is specified) as per the digs REXX variable's value, but what's currently coded will suffice
for the (default) number of years. However, it makes a difference computing the balance after sixty-five years
(when at that point, the balance becomes negative and grows increasing negative fast).
<lang rexx>/*REXX pgm (pathological FP problem): the chaotic bank society offering a new investment*/
parse arg digs use y . /*obtain optional arguments from the CL*/
if digs== | digs=="," then digs=150 /*Not specified? Then use the default.*/
if use== | use=="," then use= 20 /* " " " " " " */
if y== | y=="," then y= 25 /* " " " " " " */
numeric digits digs /*have REXX use "digs" decimal digits. */
$=2.71828182845904523536028747135266249775724709369995957496696762772407663035354759457138217852516642742746639193200305992181741359662904357290033429526
$=$ - 1 /*and subtract one 'cause that's that. */ /* [↑] 150 decimal digits of e */
/* [↑] value of newly opened account. */
do n=1 for y /*compute the value of the account/year*/
$=$*n - 1 /* " " " " " account now.*/
end /*n*/
/*display SHOW digits past the dec. pt.*/
say 'Balance after' y "years: $"format($,,use)/1 /*stick a fork in it, we're all done. */</lang> output when using the default inputs:
Balance after 25 years: $0.0399387296732302089
Siegfried Rump's example
<lang rexx>/*REXX pgm (pathological FP problem): the Siegfried Rump's example (problem dated 1988).*/ parse arg digs show . /*obtain optional arguments from the CL*/ if digs== | digs=="," then digs=150 /*Not specified? Then use the default.*/ if show== | show=="," then show= 20 /* " " " " " " */ numeric digits digs /*have REXX use "digs" decimal digits. */ a= 77617.0 /*initialize A to it's defined value.*/ b= 33096.0 /* " B " " " " */
/*display SHOW digits past the dec. pt.*/
say 'f(a,b)=' format( f(a,b), , show) /*display result from the F function.*/ exit /*stick a fork in it, we're all done. */ /*──────────────────────────────────────────────────────────────────────────────────────*/ f: procedure; parse arg a,b; a2=a**2; b2=b**2; b4=b2**2; b6=b4*b2; b8=b4**2
return 333.75*b6 + a2*(11*a2*b2 - b6 - 121*b4 - 2) + 5.5*b8 + a/(2*b)</lang>
output when using the default inputs:
f(a,b)= -0.82739605994682136814
TI-83 BASIC
A sequence that seems to converge to a wrong limit
Use the SEQ mode to enter the arithmetic progression. Note the way to set
u(1)=2 u(2)=-4
<lang ti83b> nMin=1
u(n)=111-1130/u(n-1) + 3000/(u(n-1)*u(n-2))
u(nMin)={-4;2}</lang>
The result converges to the wrong limit!
- Output:
u(20) : 100.055202 u(30) : 100 u(50) : 100 u(100) : 100
zkl
zkl doesn't have a big rational or big float library (as of this writing) but does have big ints (via GNU GMP). It does have 64 bit doubles. <lang zkl>Series:=Walker(fcn(vs){ // just keep appending new values to a list
vs.append(111.0 - 1130.0/vs[-1] + 3000.0/(vs[-1]*vs[-2])) }.fp(List(2,-4)));
series:=Series.drop(100).value;</lang> We'll use the convenient formula given in the referenced paper to create a fraction with big ints <lang zkl>var BN=Import("zklBigNum"), ten2n=BN(10).pow(64);
fcn u(n){ // use formula to create a fraction of big ints
const B=-3, Y=4; N:=BN(6).pow(n+1)*B + BN(5).pow(n+1)*Y; D:=BN(6).pow(n)*B + BN(5).pow(n)*Y; tostr(N*ten2n/D,64,32)
}
fcn tostr(bn,m,r){ // convert big int (*10^m) to float string with len r remainder, flakey
str,d:=bn.toString(), str.len()-m;
if(d<0) String(".","0"*-d,str[0,r]);
else String(str[0,d],".",str[d,r]);
}
println("1st: Convergent series"); foreach n in (T(3,4,5,6,7,8,20,30,50,100)){
"n =%3d; %3.20F %s".fmt(n,series[n-1],u(n-1)).println();
}</lang>
- Output:
Note that, at n=100, we still have diverged (at the 15th place) from the Perl6 solution and 12th place from the J solution.
1st: Convergent series n = 3; 18.50000000000000000000 18.50000000000000000000000000000000 n = 4; 9.37837837837837895449 9.37837837837837837837837837837837 n = 5; 7.80115273775216877539 7.80115273775216138328530259365994 n = 6; 7.15441448097533339023 7.15441448097524935352789065386036 n = 7; 6.80678473692481134094 6.80678473692363298394175659627200 n = 8; 6.59263276872179204702 6.59263276870443839274200277636599 n = 20; 98.34950312216535905918 6.04355211018926886777747736409754 n = 30; 99.99999999999893418590 6.00678609303120575853055404795323 n = 50; 100.00000000000000000000 6.00017584662718718894561402074719 n =100; 100.00000000000000000000 6.00000001931947792910408680340358
Chaotic banking society is just nasty so we use a five hundred digit e (the e:= text is one long line). <lang zkl>println("\n2nd: Chaotic banking society"); e:="271828182845904523536028747135266249775724709369995957496696762772407663035354759457138217852516642742746639193200305992181741359662904357290033429526059563073813232862794349076323382988075319525101901157383418793070215408914993488416750924476146066808226480016847741185374234544243710753907774499206955170276183860626133138458300075204493382656029760673711320070932870912744374704723069697720931014169283681902551510865746377211125238978442505695369677078544996996794686445490598793163688923009879312"; var en=(e.len()-1), tenEN=BN(10).pow(en); years,balance:=25, BN(e).sub(tenEN); // in place math balance=[1..years].reduce(fcn(balance,i){ balance*i - tenEN },balance); balance=tostr(balance,en,2); println("After year %d, you will have $%s in your account.".fmt(years,balance));</lang>
- Output:
2nd: Chaotic banking society After year 25, you will have $.039 in your account.
For Rump's example, multiple the formula by 10ⁿ so we can use integer math. <lang zkl>fcn rump(a,b){ b=BN(b);
b2,b4,b6,b8:=b.pow(2),b.pow(4),b.pow(6),b.pow(8); a2:=BN(a).pow(2); r:=( b6*33375 + a2*(a2*b2*11 - b6 - b4*121 - 2)*100 + b8*550 )*ten2n; r+=BN(a)*ten2n*100/(2*b); tostr(r,66,32)
} println("\n3rd: Rump's example: f(77617.0, 33096.0) = ",rump(77617,33096));</lang>
- Output:
3rd: Rump's example: f(77617.0, 33096.0) = -.82739605994682136814116509547981