Pathological floating point problems
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
FreeBASIC
<lang freebasic>' FB 1.05.0 Win64
' As FB's native types have only 64 bit precision at most we need to use the ' C library, GMP v6.1.0, for arbitrary precision arithmetic
- Include Once "gmp.bi"
mpf_set_default_prec(640) 640 bit precision, enough for this exercise
Function v(n As UInteger, prev As __mpf_struct, prev2 As __mpf_struct) As __mpf_struct
Dim As __mpf_struct a, b, c mpf_init(@a) : mpf_init(@b) : mpf_init(@c) If n = 0 Then mpf_set_ui(@a, 0UL) If n = 1 Then mpf_set_ui(@a, 2UL) If n = 2 Then mpf_set_si(@a, -4L) If n < 3 Then Return a mpf_ui_div(@a, 1130UL, @prev) mpf_mul(@b, @prev, @prev2) mpf_ui_div(@c, 3000UL, @b) mpf_ui_sub(@b, 111UL, @a) mpf_add(@a, @b, @c) mpf_clear(@b) mpf_clear(@c) Return a
End Function
Function f(a As Double, b As Double) As __mpf_Struct
Dim As __mpf_struct temp1, temp2, temp3, temp4, temp5, temp6, temp7, temp8 mpf_init(@temp1) : mpf_init(@temp2) : mpf_init(@temp3) : mpf_init(@temp4) mpf_init(@temp5) : mpf_init(@temp6) : mpf_init(@temp7) : mpf_init(@temp8) mpf_set_d(@temp1, a) a mpf_set_d(@temp2, b) b mpf_set_d(@temp3, 333.75) 333.75 mpf_pow_ui(@temp4, @temp2, 6UL) b ^ 6 mpf_mul(@temp3, @temp3, @temp4) 333.75 * b^6 mpf_pow_ui(@temp5, @temp1, 2UL) a^2 mpf_pow_ui(@temp6, @temp2, 2UL) b^2 mpf_mul_ui(@temp7, @temp5, 11UL) 11 * a^2 mpf_mul(@temp7, @temp7, @temp6) 11 * a^2 * b^2 mpf_sub(@temp7, @temp7, @temp4) 11 * a^2 * b^2 - b^6 mpf_pow_ui(@temp4, @temp2, 4UL) b^4 mpf_mul_ui(@temp4, @temp4, 121UL) 121 * b^4 mpf_sub(@temp7, @temp7, @temp4) 11 * a^2 * b^2 - b^6 - 121 * b^4 mpf_sub_ui(@temp7, @temp7, 2UL) 11 * a^2 * b^2 - b^6 - 121 * b^4 - 2 mpf_mul(@temp7, @temp7, @temp5) (11 * a^2 * b^2 - b^6 - 121 * b^4 - 2) * a^2 mpf_add(@temp3, @temp3, @temp7) 333.75 * b^6 + (11 * a^2 * b^2 - b^6 - 121 * b^4 - 2) * a^2 mpf_set_d(@temp4, 5.5) 5.5 mpf_pow_ui(@temp5, @temp2, 8UL) b^8 mpf_mul(@temp4, @temp4, @temp5) 5.5 * b^8 mpf_add(@temp3, @temp3, @temp4) 333.75 * b^6 + (11 * a^2 * b^2 - b^6 - 121 * b^4 - 2) * a^2 + 5.5 * b^8 mpf_mul_ui(@temp4, @temp2, 2UL) 2 * b mpf_div(@temp5, @temp1, @temp4) a / (2 * b) mpf_add(@temp3, @temp3, @temp5) 333.75 * b^6 + (11 * a^2 * b^2 - b^6 - 121 * b^4 - 2) * a^2 + 5.5 * b^8 + a / (2 * b) mpf_clear(@temp1) : mpf_clear(@temp2) : mpf_clear(@temp4) : mpf_clear(@temp5) mpf_clear(@temp6) : mpf_clear(@temp7) : mpf_clear(@temp8) Return temp3
End Function
Dim As Zstring * 60 z Dim As __mpf_struct result, prev, prev2 ' We cache the two previous results to avoid recursive calls to v For i As Integer = 1 To 100
result = v(i, prev, prev2) If (i >= 3 AndAlso i <= 8) OrElse i = 20 OrElse i = 30 OrElse i = 50 OrElse i = 100 Then gmp_sprintf(@z,"%53.50Ff",@result) express result to 50 decimal places Print "n ="; i , z End If prev2 = prev prev = result
Next
mpf_clear(@prev) : mpf_clear(@prev2) note : prev = result
Dim As __mpf_struct e, balance, ii, temp mpf_init(@e) : mpf_init(@balance) : mpf_init(@ii) : mpf_init(@temp) mpf_set_str(@e, "2.71828182845904523536028747135266249775724709369995", 10) e to 50 decimal places mpf_sub_ui(@balance, @e, 1UL)
For i As ULong = 1 To 25
mpf_set_ui(@ii, i) mpf_mul(@temp, @balance, @ii) mpf_sub_ui(@balance, @temp, 1UL)
Next
Print Print "Chaotic B/S balance after 25 years : "; gmp_sprintf(@z,"%.16Ff",@balance) express balance to 16 decimal places Print z mpf_clear(@e) : mpf_clear(@balance) : mpf_clear(@ii) : mpf_clear(@temp)
Print Dim rump As __mpf_struct rump = f(77617.0, 33096.0) gmp_sprintf(@z,"%.16Ff", @rump) express rump to 16 decimal places Print "f(77617.0, 33096.0) = "; z
Print Print "Press any key to quit" Sleep</lang>
- Output:
n = 3 18.50000000000000000000000000000000000000000000000000 n = 4 9.37837837837837837837837837837837837837837837837838 n = 5 7.80115273775216138328530259365994236311239193083573 n = 6 7.15441448097524935352789065386036202438123383819727 n = 7 6.80678473692363298394175659627200908762327670780193 n = 8 6.59263276870443839274200277636599482655298231773461 n = 20 6.04355211018926886777747736409754013318771500000612 n = 30 6.00678609303120575853055404795323970583307231443837 n = 50 6.00017584662718718894561402074719546952373517709933 n = 100 6.00000001931947792910408680340358571502435067543695 Chaotic B/S balance after 25 years : 0.0399387296732302 f(77617.0, 33096.0) = -0.8273960599468214
Fortran
Problems arise because the floating-point arithmetic as performed by digital computers has only an oblique relationship to the arithmetic of Real numbers: many axia are violated, even if only by a little, and in certain situations. Most seriously, only a limited precision is available even if the floating-point variables are declared via such words as "REAL". Actions such as subtraction (of nearly-equal values) can in one step destroy many or all the digits of accuracy of the value being developed.
Fortran's only "built-in" assistance in this is provided via the ability to declare variables DOUBLE PRECISION, and on some systems, QUADRUPLE PRECISION is available. Earlier systems such as the IBM1620 supported decimal arithmetic of up to ninety-nine decimal digits (via hardware!), specified via control cards, but this sort of ability has been abandoned. Special "bignumber" arithmetic routines can be written supporting floating-point (or integer, or rational) arithmetic of hundreds or thousands or more words of storage per number, but this is not a standard arrangement.
Otherwise, a troublesome calculation might be recast into a different form that avoids a catastrophic loss of precision, probably after a lot of careful and difficult analysis and exploration and ingenuity.
Here, no such attempt is made. In the spirit of Formula Translation, this is a direct translation of the specified formulae into Fortran, with single and double precision results on display. There is no REAL*16 option, nor the REAL*10 that some systems allow to correspond to the eighty-bit floating-point format supported by the floating-point processor. The various integer constants cause no difficulty and I'm not bothering with writing them as <integer>.0 - the compiler can deal with this. The constants with fractional parts happen to be exactly represented in binary so there is no fuss over 333.75 and 333.75D0 whereas by contrast 0.1 and 0.1D0 are not equal. Similarly, there is no attempt to rearrange the formulae, for instance to have A**2 * B**2 replaced by (A*B)**2, nor worry over B**8 where 33096**8 = 1.439E36 and the largest possible single-precision number is 3.4028235E+38, in part because arithmetic within an expression can be conducted with a greater dynamic range. Most of all, no attention has been given to the subtractions...
This would be F77 style Fortran, except for certain conveniences offered by F90, especially the availability of generic functions such as EXP whose type is determined by the type of its parameter, rather than having to use EXP and DEXP for single and double precision respectively, or else... The END statement for subroutines and functions names the routine being ended, a useful matter to have checked. <lang Fortran> SUBROUTINE MULLER
REAL*4 VN,VNL1,VNL2 !The exact precision and dynamic range
REAL*8 WN,WNL1,WNL2 !Depends on the format's precise usage of bits.
INTEGER I !A stepper.
WRITE (6,1) !A heading.
1 FORMAT ("Muller's sequence should converge to six...",/
1 " N Single Double")
VNL1 = 2; VN = -4 !Initialise for N = 2.
WNL1 = 2; WN = -4 !No fractional parts yet.
DO I = 3,36 !No point going any further.
VNL2 = VNL1; VNL1 = VN !Shuffle the values along one place.
WNL2 = WNL1; WNL1 = WN !Ready for the next term's calculation.
VN = 111 - 1130/VNL1 + 3000/(VNL1*VNL2) !Calculate the next term.
WN = 111 - 1130/WNL1 + 3000/(WNL1*WNL2) !In double precision.
WRITE (6,2) I,VN,WN !Show both.
2 FORMAT (I3,F12.7,F21.16) !With too many fractional digits.
END DO !On to the next term.
END SUBROUTINE MULLER !That was easy. Too bad the results are wrong.
SUBROUTINE CBS !The Chaotic Bank Society.
INTEGER YEAR !A stepper.
REAL*4 V !The balance.
REAL*8 W !In double precision as well.
V = 1; W = 1 !Initial values, without dozy 1D0 stuff.
V = EXP(V) - 1 !Actual initial value desired is e - 1..,
W = EXP(W) - 1 !This relies on double-precision W selecting DEXP.
WRITE (6,1) !Here we go.
1 FORMAT (///"The Chaotic Bank Society in action..."/"Year")
WRITE (6,2) 0,V,W !Show the initial deposit.
2 FORMAT (I3,F16.7,F28.16)
DO YEAR = 1,25 !Step through some years.
V = V*YEAR - 1 !The specified procedure.
W = W*YEAR - 1 !The compiler handles type conversions.
WRITE (6,2) YEAR,V,W !The current balance.
END DO !On to the following year.
END SUBROUTINE CBS !Madness!
REAL*4 FUNCTION SR4(A,B) !Siegfried Rump's example function of 1988.
REAL*4 A,B
SR4 = 333.75*B**6
1 + A**2*(11*A**2*B**2 - B**6 - 121*B**4 - 2)
2 + 5.5*B**8 + A/(2*B)
END FUNCTION SR4
REAL*8 FUNCTION SR8(A,B) !Siegfried Rump's example function.
REAL*8 A,B
SR8 = 333.75*B**6 !.75 is exactly represented in binary.
1 + A**2*(11*A**2*B**2 - B**6 - 121*B**4 - 2)
2 + 5.5*B**8 + A/(2*B)!.5 is exactly represented in binary.
END FUNCTION SR8
PROGRAM POKE
REAL*4 V !Some example variables.
REAL*8 W !Whose type goes to the inquiry function.
WRITE (6,1) RADIX(V),DIGITS(V),"single",DIGITS(W),"double"
1 FORMAT ("Floating-point arithmetic is conducted in base ",I0,/
1 2(I3," digits for ",A," precision",/))
WRITE (6,*) "Single precision limit",HUGE(V)
WRITE (6,*) "Double precision limit",HUGE(W)
WRITE (6,*)
CALL MULLER
CALL CBS
WRITE (6,10)
10 FORMAT (///"Evaluation of Siegfried Rump's function of 1988",
1 " where F(77617,33096) = -0.827396059946821")
WRITE (6,*) "Single precision:",SR4(77617.0,33096.0)
WRITE (6,*) "Double precision:",SR8(77617.0D0,33096.0D0) !Must match the types.
END</lang>
Output
Floating-point numbers in single and double precision use the "implicit leading one" binary format on this system: there have been many variations across different computers across the years. One can write strange routines that will test the workings of arithmetic (and other matters) so as to determine the situation on the computer of the moment, but F90 introduced special "inquiry" routines that reveal certain details as standard. This information could be used to make choices amongst calculation paths and ploys appropriate for different results, at of course a large expenditure in thought to produce a compound scheme that will (should?) work correctly on a variety of computers. No such effort has been made here!
Fifty-three binary digits corresponds to 15·95 decimal digits: there is no simple conversion so the usual ploy is to show additional decimal digits, knowing that the lower-order digits will be fuzz due to the binary/decimal conversion. The "Muller" sequence has for its fourth term 9.3783783783783790 - note that this is a recurring sequence, and its precision is less than the displayed sixteen decimal digits (seventeen digits of "precision" are on show) - when trying for maximum accuracy, converting a binary value to decimal adds confusion.
Floating-point arithmetic is conducted in base 2 24 digits for single precision 53 digits for double precision Single precision limit 3.4028235E+38 Double precision limit 1.797693134862316E+308 Muller's sequence should converge to six... N Single Double 3 18.5000000 18.5000000000000000 4 9.3783779 9.3783783783783790 5 7.8011475 7.8011527377521688 6 7.1543465 7.1544144809753334 7 6.8058305 6.8067847369248113 8 6.5785794 6.5926327687217920 9 6.2355156 6.4494659340539329 10 2.9135900 6.3484520607466237 11-111.7097931 6.2744386627281159 12 111.8982391 6.2186967685821628 13 100.6615448 6.1758538558153901 14 100.0406036 6.1426271704810063 15 100.0024948 6.1202487045701588 16 100.0001526 6.1660865595980994 17 100.0000076 7.2350211655349312 18 100.0000000 22.0620784635257934 19 100.0000000 78.5755748878722358 20 100.0000000 98.3495031221653591 21 100.0000000 99.8985692661829034 22 100.0000000 99.9938709889027848 23 100.0000000 99.9996303872863450 24 100.0000000 99.9999777306794897 25 100.0000000 99.9999986592166863 26 100.0000000 99.9999999193218088 27 100.0000000 99.9999999951477605 28 100.0000000 99.9999999997082796 29 100.0000000 99.9999999999824638 30 100.0000000 99.9999999999989342 31 100.0000000 99.9999999999999289 32 100.0000000 99.9999999999999858 33 100.0000000 100.0000000000000000 34 100.0000000 100.0000000000000000 35 100.0000000 100.0000000000000000 36 100.0000000 100.0000000000000000 The Chaotic Bank Society in action... Year 0 1.7182819 1.7182818284590453 1 0.7182819 0.7182818284590453 2 0.4365637 0.4365636569180906 3 0.3096912 0.3096909707542719 4 0.2387648 0.2387638830170875 5 0.1938238 0.1938194150854375 6 0.1629429 0.1629164905126252 7 0.1406002 0.1404154335883767 8 0.1248016 0.1233234687070137 9 0.1232147 0.1099112183631235 10 0.2321472 0.0991121836312345 11 1.5536194 0.0902340199435798 12 17.6434326 0.0828082393229579 13 228.3646240 0.0765071111984525 14 3196.1047363 0.0710995567783357 15 47940.5703125 0.0664933516750352 16 767048.1250000 0.0638936268005637 1713039817.0000000 0.0861916556095821 18**************** 0.5514498009724775 19**************** 9.4775462184770731 20**************** 188.5509243695414625 21**************** 3958.5694117603707127 22**************** 87087.5270587281556800 23**************** 2003012.1223507476970553 24**************** 48072289.9364179447293282 25**************** 1201807247.4104485511779785 Evaluation of Siegfried Rump's function of 1988 where F(77617,33096) = -0.827396059946821 Single precision: -1.1805916E+21 Double precision: -1.1805916E+21
None of the results are remotely correct! In the absence of a Fortran compiler supporting still higher precision (such as quadruple, or REAL*16) only two options remain: either devise multi-word high-precision arithmetic routines and try again with even more more brute-force, or, analyse the calculation with a view to finding a way to avoid the loss of accuracy with calculations conducted in the available precision.
Alternatively, do not present various intermediate results such as might give rise to doubts, nor yet entertain any doubts, just declare the answer to be what appears, and move on. In a letter from F.S. Acton, "A former student of mine now hands out millions of dollars for computation ... and he dismally estimates that 70% of the "answers" are worthless because of poor analysis and poor programming."
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
Icon and Unicon
Icon and Unicon support large integers. Used for scaling the intermediates. Task 1 includes an extra step, 200 iterations, to demonstrate a closer convergence.
<lang unicon>#
- Pathological floating point problems
procedure main()
sequence() chaotic()
end
- First task, sequence convergence
link printf procedure sequence()
local l := [2, -4]
local iters := [3, 4, 5, 6, 7, 8, 20, 30, 50, 100, 200]
local i, j, k
local n := 1
write("Sequence convergence")
# Demonstrate the convergence problem with various precision values
every k := (100 | 300) do {
n := 10^k
write("\n", k, " digits of intermediate precision")
# numbers are scaled up using large integer powers of 10
every i := !iters do {
l := [2 * n, -4 * n]
printf("i: %3d", i)
every j := 3 to i do {
# build out a list of intermediate passes
# order of scaling operations matters
put(l, 111 * n - (1130 * n * n / l[j - 1]) +
(3000 * n * n * n / (l[j - 1] * l[j - 2])))
}
# down scale the result to a real
# some precision may be lost in the final display
printf(" %20.16r\n", l[i] * 1.0 / n)
}
}
end
- Task 2, chaotic bank of Euler
procedure chaotic()
local euler, e, scale, show, y, d
write("\nChaotic Banking Society of Euler")
# format the number for listing, string form, way overboard on digits
euler :=
"2718281828459045235360287471352662497757247093699959574966967627724076630353_
547594571382178525166427427466391932003059921817413596629043572900334295260_ 595630738132328627943490763233829880753195251019011573834187930702154089149_ 934884167509244761460668082264800168477411853742345442437107539077744992069_ 551702761838606261331384583000752044933826560297606737113200709328709127443_ 747047230696977209310141692836819025515108657463772111252389784425056953696_ 770785449969967946864454905987931636889230098793127736178215424999229576351_ 482208269895193668033182528869398496465105820939239829488793320362509443117_ 301238197068416140397019837679320683282376464804295311802328782509819455815_ 301756717361332069811250996181881593041690351598888519345807273866738589422_ 879228499892086805825749279610484198444363463244968487560233624827041978623_ 209002160990235304369941849146314093431738143640546253152096183690888707016_ 768396424378140592714563549061303107208510383750510115747704171898610687396_ 9655212671546889570350354"
# precise math with long integers, string form just for pretty listing e := integer(euler)
# 1000 digits after the decimal for scaling intermediates and service fee scale := 10^1000
# initial deposit - $1 d := e - scale
# show balance with 16 digits
show := 10^16
write("Starting balance: $", d * show / scale * 1.0 / show, "...")
# wait 25 years, with only a trivial $1 service fee
every y := 1 to 25 do {
d := d * y - scale
}
# show final balance with 4 digits after the decimal (truncation)
show := 10^4
write("Balance after ", y, " years: $", d * show / scale * 1.0 / show)
end</lang>
- Output:
prompt$ time unicon -s patho.icn -x Sequence convergence 100 digits of intermediate precision i: 3 18.5000000000000000 i: 4 9.3783783783783790 i: 5 7.8011527377521620 i: 6 7.1544144809752490 i: 7 6.8067847369236330 i: 8 6.5926327687044380 i: 20 6.0435521101892680 i: 30 6.0067860930312060 i: 50 6.0001758466271870 i: 100 99.9999999999998400 i: 200 100.0000000000000000 300 digits of intermediate precision i: 3 18.5000000000000000 i: 4 9.3783783783783790 i: 5 7.8011527377521610 i: 6 7.1544144809752490 i: 7 6.8067847369236320 i: 8 6.5926327687044380 i: 20 6.0435521101892680 i: 30 6.0067860930312060 i: 50 6.0001758466271870 i: 100 6.0000000193194780 i: 200 6.0000000000000000 Chaotic Banking Society of Euler Starting balance: $1.718281828459045... Balance after 25 years: $0.0399 real 0m0.075s user 0m0.044s sys 0m0.020s
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.
Kotlin
<lang scala>// version 1.0.6
import java.math.*
const val LIMIT = 100
val con480 = MathContext(480) val bigTwo = BigDecimal(2) val bigE = BigDecimal("2.71828182845904523536028747135266249775724709369995") // precise enough!
fun main(args: Array<String>) {
// v(n) sequence task
val c1 = BigDecimal(111)
val c2 = BigDecimal(1130)
val c3 = BigDecimal(3000)
var v1 = bigTwo
var v2 = BigDecimal(-4)
var v3: BigDecimal
for (i in 3 .. LIMIT) {
v3 = c1 - c2.divide(v2, con480) + c3.divide(v2 * v1, con480)
println("${"%3d".format(i)} : ${"%19.16f".format(v3)}")
v1 = v2
v2 = v3
}
// Chaotic Building Society task
var balance = bigE - BigDecimal.ONE
for (year in 1..25) balance = balance.multiply(BigDecimal(year), con480) - BigDecimal.ONE
println("\nBalance after 25 years is ${"%18.16f".format(balance)}")
// Siegfried Rump task
val a = BigDecimal(77617)
val b = BigDecimal(33096)
val c4 = BigDecimal("333.75")
val c5 = BigDecimal(11)
val c6 = BigDecimal(121)
val c7 = BigDecimal("5.5")
var f = c4 * b.pow(6, con480) + c7 * b.pow(8, con480) + a.divide(bigTwo * b, con480)
val c8 = c5 * a.pow(2, con480) * b.pow(2, con480) - b.pow(6, con480) - c6 * b.pow(4, con480) - bigTwo
f += c8 * a.pow(2, con480)
println("\nf(77617.0, 33096.0) is ${"%18.16f".format(f)}")
}</lang>
- Output:
3 : 18.5000000000000000 4 : 9.3783783783783784 5 : 7.8011527377521614 6 : 7.1544144809752494 7 : 6.8067847369236330 8 : 6.5926327687044384 9 : 6.4494659337902880 10 : 6.3484520566543571 11 : 6.2744385982163279 12 : 6.2186957398023978 13 : 6.1758373049212301 14 : 6.1423590812383559 15 : 6.1158830665510808 16 : 6.0947394393336811 17 : 6.0777223048472427 18 : 6.0639403224998088 19 : 6.0527217610161522 20 : 6.0435521101892689 21 : 6.0360318810818568 22 : 6.0298473250239019 23 : 6.0247496523668479 24 : 6.0205399840615161 25 : 6.0170582573289876 26 : 6.0141749145508190 27 : 6.0117845878713337 28 : 6.0098012392984846 29 : 6.0081543789122289 30 : 6.0067860930312058 31 : 6.0056486887714203 32 : 6.0047028131881752 33 : 6.0039159416664605 34 : 6.0032611563057406 35 : 6.0027161539543513 36 : 6.0022624374405593 37 : 6.0018846538818819 38 : 6.0015700517342190 39 : 6.0013080341649643 40 : 6.0010897908901841 41 : 6.0009079941545271 42 : 6.0007565473053508 43 : 6.0006303766028389 44 : 6.0005252586505718 45 : 6.0004376772265183 46 : 6.0003647044182955 47 : 6.0003039018761868 48 : 6.0002532387368678 49 : 6.0002110233741743 50 : 6.0001758466271872 51 : 6.0001465345613879 52 : 6.0001221091522881 53 : 6.0001017555560260 54 : 6.0000847948586303 55 : 6.0000706613835716 56 : 6.0000588837928413 57 : 6.0000490693458029 58 : 6.0000408907870884 59 : 6.0000340754236785 60 : 6.0000283960251310 61 : 6.0000236632422855 62 : 6.0000197192908008 63 : 6.0000164326883272 64 : 6.0000136938694348 65 : 6.0000114115318177 66 : 6.0000095095917616 67 : 6.0000079246472413 68 : 6.0000066038639788 69 : 6.0000055032139253 70 : 6.0000045860073981 71 : 6.0000038216699107 72 : 6.0000031847228971 73 : 6.0000026539343389 74 : 6.0000022116109709 75 : 6.0000018430084630 76 : 6.0000015358399141 77 : 6.0000012798662675 78 : 6.0000010665549954 79 : 6.0000008887956715 80 : 6.0000007406629499 81 : 6.0000006172190487 82 : 6.0000005143491543 83 : 6.0000004286242585 84 : 6.0000003571868566 85 : 6.0000002976556961 86 : 6.0000002480464011 87 : 6.0000002067053257 88 : 6.0000001722544322 89 : 6.0000001435453560 90 : 6.0000001196211272 91 : 6.0000000996842706 92 : 6.0000000830702242 93 : 6.0000000692251858 94 : 6.0000000576876542 95 : 6.0000000480730447 96 : 6.0000000400608703 97 : 6.0000000333840583 98 : 6.0000000278200485 99 : 6.0000000231833736 100 : 6.0000000193194779 Balance after 25 years is 0.0399387296732302 f(77617.0, 33096.0) is -0.8273960599468214
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 show y . /*obtain optional arguments from the CL*/
if digs== | digs=="," then digs=150 /*Not specified? Then use the default.*/
if show== | show=="," then show= 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*/
@baf= 'Balance after' /*display SHOW digits past the dec. pt.*/ say @baf y "years: $"format($, , show) / 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