Interval arithmetic by Robert Elton Maas:

Demonstration of my newly-written software for interval
arithmetic with rational endpoints, including basic
arithmetic functions (plus difference times quotient)
and Newton's method for computing square root, in
equivalent form where arithmetic mean of guess and
num/guess is taken as good approximation for geometric
mean which is sqrt(num), adapted for interval
arithmetic. Here I set up the initial values of each
value from its dependencies, then manually cycle
through all the iterators until I have nearly 80
significant digits. The calculation is sqrt(3 +
2*sqrt(2)) - sqrt(2), where the sqrt(2) calculation is
shared between the two uses for it. Mathematically the
final result is an exact representable number, which
you'll see if you skip down a ways:

*       (setq guess (interval-to-sqrt-firstguess+bounds 2))

5/3
*       (setq sqrt2 (interval+guess-to-sqrt-bounds 2 guess))

(6/5 5/3)
*       (show-interval sqrt2)
1.2d0
1.6666666666666667d0
1.[2..7]
NIL
*       (setq sqrt8 nil)

NIL
* (setq sqrt8 (update-intervals-times sqrt8 2 sqrt2))
2.4d0
3.3333333333333335d0
2.[4..E]
(12/5 10/3)
*       (setq surdsum nil)

NIL
* (setq surdsum (update-intervals-plus surdsum 3 sqrt8))
5.4d0
6.333333333333333d0
5.[4..E]
(27/5 19/3)
*       (setq guess (interval-to-sqrt-firstguess+bounds surdsum))

12/5
*       (setq sqrtsurd (interval+guess-to-sqrt-bounds surdsum guess))

(9/4 95/36)
*       (show-interval sqrtsurd)
2.25d0
2.638888888888889d0
2.[2..7]
NIL
* sqrtsurd sqrt2

(9/4 95/36)
*
(6/5 5/3)
* (show-interval sqrtsurd)
2.25d0
2.638888888888889d0
2.[2..7]
NIL
* (show-interval sqrt2)
1.2d0
1.6666666666666667d0
1.[2..7]
NIL
*       (setq result nil)

NIL
* (setq result (update-intervals-diff result sqrtsurd sqrt2))
0.5833333333333334d0
1.4388888888888889d0
0.[5..F]
(7/12 259/180)
* (setq guessres (ratio-interval-to-simplest-ratio-within result))

1
* (setq sqrt2 (interval-newton-sqrt 2 sqrt2))
1.2d0
1.6666666666666667d0
1.[2..7]
[-----------------------------!!!!!!-----------------------------------]
1.3953488372093024d0
1.4333333333333333d0
1.3[9..E]
[------------------------------!!!!!-----------------------------------]
1.4d0
1.4285714285714286d0
1.4[0..3]
(7/5 10/7)
* (setq sqrt8 (update-intervals-times sqrt8 2 sqrt2))
2.4d0
3.3333333333333335d0
2.[4..E]
[------------------------------!!!!!-----------------------------------]
2.8d0
2.857142857142857d0
2.8[0..6]
(14/5 20/7)
* (setq surdsum (update-intervals-plus surdsum 3 sqrt8))
5.4d0
6.333333333333333d0
5.[4..E]
[------------------------------!!!!!-----------------------------------]
5.8d0
5.857142857142857d0
5.8[0..6]
(29/5 41/7)
* (setq sqrtsurd (interval-newton-sqrt surdsum sqrtsurd))
2.25d0
2.638888888888889d0
2.[2..7]
[----------------------!!!!!!!!!!!!!-----------------------------------]
2.3727272727272726d0
2.4444444444444446d0
2.3[7..F]
[---------------------------!!!!!!!!-----------------------------------]
2.4d0
2.4404761904761907d0
2.4[0..5]
(12/5 205/84)
* (setq result (update-intervals-diff result sqrtsurd sqrt2))
0.5833333333333334d0
1.4388888888888889d0
0.[5..F]
[-------------------------------!!!!!!!--------------------------------]
0.9714285714285714d0
1.0404761904761906d0
0.9[7..F]
(34/35 437/420)
* (setq guessres (ratio-interval-to-simplest-ratio-within result))

1
* (setq sqrt2 (interval-newton-sqrt 2 sqrt2))
1.4d0
1.4285714285714286d0
1.4[0..3]
[----------------------------------!-----------------------------------]
1.4141414141414141d0
1.4142857142857144d0
1.4141[4..J]
[----------------------------------!-----------------------------------]
1.4141414141414141d0
1.4142857142857144d0
1.4141[4..J]
(140/99 99/70)
* (setq sqrt8 (update-intervals-times sqrt8 2 sqrt2))
2.8d0
2.857142857142857d0
2.8[0..6]
[----------------------------------!-----------------------------------]
2.8282828282828283d0
2.8285714285714287d0
2.828[2..6]
(280/99 99/35)
* (setq surdsum (update-intervals-plus surdsum 3 sqrt8))
5.8d0
5.857142857142857d0
5.8[0..6]
[----------------------------------!-----------------------------------]
5.828282828282828d0
5.828571428571428d0
5.828[2..6]
(577/99 204/35)
* (setq sqrtsurd (interval-newton-sqrt surdsum sqrtsurd))
2.4d0
2.4404761904761907d0
2.4[0..5]
[--------------!!!!!!!!!!!!!!!!!!!!!-----------------------------------]
2.408144405193102d0
2.420238095238095d0
2.4[0..3]
[--------------------!!!!!!!!!-----------------------------------------]
2.411703239289446d0
2.4166666666666665d0
2.41[1..7]
(2308/957 29/12)
* (setq result (update-intervals-diff result sqrtsurd sqrt2))
0.9714285714285714d0
1.0404761904761906d0
0.9[7..F]
[--------------------------!!!!!!--------------------------------------]
0.9974175250037319d0
1.0025252525252526d0
0.99[7..D]
(66817/66990 397/396)
* (setq guessres (ratio-interval-to-simplest-ratio-within result))

1
* (setq sqrt2 (interval-newton-sqrt 2 sqrt2))
1.4141414141414141d0
1.4142857142857144d0
1.4141[4..J]
[----------------------------------!-----------------------------------]
1.4142135605326258d0
1.4142135642135643d0
1.41421356[0..5]
[----------------------------------!-----------------------------------]
1.4142135605326258d0
1.4142135642135643d0
1.41421356[0..5]
(27720/19601 19601/13860)
* (setq sqrt8 (update-intervals-times sqrt8 2 sqrt2))
2.8282828282828283d0
2.8285714285714287d0
2.828[2..6]
[----------------------------------!-----------------------------------]
2.8284271210652516d0
2.8284271284271285d0
2.82842712[1..9]
(55440/19601 19601/6930)
* (setq surdsum (update-intervals-plus surdsum 3 sqrt8))
5.828282828282828d0
5.828571428571428d0
5.828[2..6]
[----------------------------------!-----------------------------------]
5.828427121065252d0
5.828427128427128d0
5.82842712[1..9]
(114243/19601 40391/6930)
* (setq sqrtsurd (interval-newton-sqrt surdsum sqrtsurd))
2.411703239289446d0
2.4166666666666665d0
2.41[1..7]
[-----------------------------------!----------------------------------]
2.4141849529780566d0
2.4142421736318833d0
2.4141[8..F]
[-----------------------------------!----------------------------------]
2.4142011834319526d0
2.4142259429024135d0
2.4142[0..3]
(408/169 6826079/2827440)
* (setq result (update-intervals-diff result sqrtsurd sqrt2))
0.9974175250037319d0
1.0025252525252526d0
0.99[7..D]
[-----------------------------------!----------------------------------]
0.9999876192183884d0
1.0000123823697875d0
0.9999[8..C]
(2342311/2342340 3260078687/3260038320)
* (setq guessres (ratio-interval-to-simplest-ratio-within result))

1
* (setq sqrt2 (interval-newton-sqrt 2 sqrt2))
1.4142135605326258d0
1.4142135642135643d0
1.41421356[0..5]
[----------------------------------!-----------------------------------]
1.4142135623730951d0
1.4142135623730951d0
1.41421356237309504[7..A]
[----------------------------------!-----------------------------------]
1.4142135623730951d0
1.4142135623730951d0
1.41421356237309504[7..A]
(1086679440/768398401 768398401/543339720)
* (setq sqrt8 (update-intervals-times sqrt8 2 sqrt2))
2.8284271210652516d0
2.8284271284271285d0
2.82842712[1..9]
[----------------------------------!-----------------------------------]
2.8284271247461903d0
2.8284271247461903d0
2.82842712474619009[5..A]
(2173358880/768398401 768398401/271669860)
* (setq surdsum (update-intervals-plus surdsum 3 sqrt8))
5.828427121065252d0
5.828427128427128d0
5.82842712[1..9]
[----------------------------------!-----------------------------------]
5.82842712474619d0
5.82842712474619d0
5.82842712474619009[5..A]
(4478554083/768398401 1583407981/271669860)
* (setq sqrtsurd (interval-newton-sqrt surdsum sqrtsurd))
2.4142011834319526d0
2.4142259429024135d0
2.4142[0..3]
[----------------------------------!-----------------------------------]
2.414213561579007d0
2.4142135631671833d0
2.41421356[1..4]
[----------------------------------!-----------------------------------]
2.4142135620573204d0
2.4142135626888694d0
2.414213562[0..7]
(80782/33461 52982414452241/21946034630520)
* (setq result (update-intervals-diff result sqrtsurd sqrt2))
0.9999876192183884d0
1.0000123823697875d0
0.9999[8..C]
[----------------------------------!-----------------------------------]
0.9999999996842254d0
1.0000000003157745d0
0.999999999[6..E]
(18180690365179/18180690370920 16863297923707194757841/16863297918382193798520)
* (setq guessres (ratio-interval-to-simplest-ratio-within result))

1
* (setq sqrt2 (interval-newton-sqrt 2 sqrt2))
1.4142135623730951d0
1.4142135623730951d0
1.41421356237309504[7..A]
[----------------------------------!-----------------------------------]
1.4142135623730951d0
1.4142135623730951d0
1.414213562373095048801688724209698078[0..B]
[----------------------------------!-----------------------------------]
1.4142135623730951d0
1.4142135623730951d0
1.414213562373095048801688724209698078[0..B]
(1670005488191150880/1180872205318713601 1180872205318713601/835002744095575440)
* (setq sqrt8 (update-intervals-times sqrt8 2 sqrt2))
2.8284271247461903d0
2.8284271247461903d0
2.82842712474619009[5..A]
[----------------------------------!-----------------------------------]
2.8284271247461903d0
2.8284271247461903d0
2.82842712474619009760337744841939615[6..9]
(3340010976382301760/1180872205318713601 1180872205318713601/417501372047787720)
* (setq surdsum (update-intervals-plus surdsum 3 sqrt8))
5.82842712474619d0
5.82842712474619d0
5.82842712474619009[5..A]
[----------------------------------!-----------------------------------]
5.82842712474619d0
5.82842712474619d0
5.82842712474619009760337744841939615[6..9]
(6882627592338442563/1180872205318713601 2433376321462076761/417501372047787720)
* (setq sqrtsurd (interval-newton-sqrt surdsum sqrtsurd))
2.4142135620573204d0
2.4142135626888694d0
2.414213562[0..7]
[----------------------------------!-----------------------------------]
2.414213562373095d0
2.414213562373095d0
2.414213562373095048[2..E]
[----------------------------------!-----------------------------------]
2.414213562373095d0
2.414213562373095d0
2.414213562373095048[5..B]
(3166815962/1311738121
 3191952483600556543231906081/1322150009157834778483586640)
* (setq result (update-intervals-diff result sqrtsurd sqrt2))
0.9999999996842254d0
1.0000000003157745d0
0.999999999[6..E]
[----------------------------------!-----------------------------------]
1.0d0
1.0d0
0.999999999999999999[7..D]
(1095304930569773971854289559/1095304930569773972079348240
 1561290197076369738695283398195708598377064481/15612901970763697383744760114988
82417429890640)
* (setq guessres (ratio-interval-to-simplest-ratio-within result))

1
* (setq sqrt2 (interval-newton-sqrt 2 sqrt2))
1.4142135623730951d0
1.4142135623730951d0
1.414213562373095048801688724209698078[0..B]
[----------------------------------!-----------------------------------]
1.4142135623730951d0
1.4142135623730951d0
1.414213562373095048801688724209698078569671875376948073176679737990732478[3..6]
[----------------------------------!-----------------------------------]
1.4142135623730951d0
1.4142135623730951d0
1.414213562373095048801688724209698078569671875376948073176679737990732478[3..6]
(3944126127469278527968910146598237760/2788918330588564181308597538924774401
 2788918330588564181308597538924774401/1972063063734639263984455073299118880)
* (setq sqrt8 (update-intervals-times sqrt8 2 sqrt2))
2.8284271247461903d0
2.8284271247461903d0
2.82842712474619009760337744841939615[6..9]
[----------------------------------!-----------------------------------]
2.8284271247461903d0
2.8284271247461903d0
2.828427124746190097603377448419396157139343750753896146353359475981464956[7..C]
(7888252254938557055937820293196475520/2788918330588564181308597538924774401
 2788918330588564181308597538924774401/986031531867319631992227536649559440)
* (setq surdsum (update-intervals-plus surdsum 3 sqrt8))
5.82842712474619d0
5.82842712474619d0
5.82842712474619009760337744841939615[6..9]
[----------------------------------!-----------------------------------]
5.82842712474619d0
5.82842712474619d0
5.828427124746190097603377448419396157139343750753896146353359475981464956[7..C]
(16255007246704249599863612909970798723/2788918330588564181308597538924774401
 5747012926190523077285280148873452721/986031531867319631992227536649559440)
* (setq sqrtsurd (interval-newton-sqrt surdsum sqrtsurd))
2.414213562373095d0
2.414213562373095d0
2.414213562373095048[5..B]
[----------------------------------!-----------------------------------]
2.414213562373095d0
2.414213562373095d0
2.414213562373095048801688724209698078[3..8]
[----------------------------------!-----------------------------------]
2.414213562373095d0
2.414213562373095d0
2.414213562373095048801688724209698078[4..7]
(4866752642924153522/2015874949414289041
 11585259391867585946511862767537871507692906187141930561/4798771563721829525434
879762059756187376114916224347680)
* (setq result (update-intervals-diff result sqrtsurd sqrt2))
1.0d0
1.0d0
0.999999999999999999[7..D]
[----------------------------------!-----------------------------------]
1.0d0
1.0d0
0.9999999999999999999999999999999999999[1..J]
(3975432528847853791284206486836718264843744325717055919/39754325288478537912842
06486836718265189613786940194080
 1338338197837095844130945439064043907024031450098297617793624047027247167548672
8860987972161/133833819783709584413094543906404390690759372714330250164448105020
77150463062548147387739680)
* (setq guessres (ratio-interval-to-simplest-ratio-within result))

1
* (setq sqrtsurd (interval-newton-sqrt surdsum sqrtsurd))
2.414213562373095d0
2.414213562373095d0
2.414213562373095048801688724209698078[4..7]
[----------------------------------!-----------------------------------]
2.414213562373095d0
2.414213562373095d0
2.414213562373095048801688724209698078569671875376948073176679737990732478[3..6]
[----------------------------------!!----------------------------------]
2.414213562373095d0
2.414213562373095d0
2.414213562373095048801688724209698078569671875376948073176679737990732478[3..7]
(6733044458057842709277507685523012161/2788918330588564181308597538924774401
 16027949695982172840402177156513573335102057071633379140051638962864595121/6638
994141109541574736124039057831663823073029360974515611089182412349840)
* (setq result (update-intervals-diff result sqrtsurd sqrt2))
1.0d0
1.0d0
0.9999999999999999999999999999999999999[1..J]
[----------------------------------!!----------------------------------]
1.0d0
1.0d0
0.999999999999999999999999999999999999999999999999999999999999999999999999[8..D]
(5499922827526179381859858156795820014911822025822860217658921196079790879/54999
22827526179381859858156795820014911822025822860217658921196079790880
 1851561245681048118705335426678444074181234761664440285203531957696284269057890
8189829678694334707090512339121/185156124568104811870533542667844407418123476166
44402852035319576962842685817926795506475249041654478288445840)
* (setq guessres (ratio-interval-to-simplest-ratio-within result))

1

Yes, the correct mathematical answer is exactly 1. How
quickly did you figure that out?
Proof:
(1 + sqrt(2))**2 = 1 + 2*sqrt(2) + sqrt(2)**2 = 3 + 2*sqrt(2)
 1 + sqrt(2)                             = sqrt(3 + 2*sqrt(2))
 1 =                                       sqrt(3 + 2*sqrt(2)) - sqrt(2)
The right side is brute-force calculated by interval arithmetic.


By comparison, if you use ordinary double-float arithmetic:

* (setq sqrt2 (sqrt 2d0))

1.4142135623730951d0
* (setq sqrt8 (* 2 sqrt2))

2.8284271247461903d0
* (setq surdsum (+ 3 sqrt8))

5.82842712474619d0
* (setq sqrtsurd (sqrt surdsum))

2.414213562373095d0
* (setq result (- sqrtsurd sqrt2))

0.9999999999999998d0

you get a result is not exactly 1, and you don't know
if that's the result because of roundoff error or the
result really isn't exactly 1, and there's no obvious
way to get a more accurate answer to be more confident
it might really be exactly 1. Furthermore it's easy to
get an answer that looks like exactly 1 when
mathematically it's not exactly 1:

* (sqrt 1.0000000000000002d0)

1.0d0

With unlimited-precision interval arithmetic, as I've
implemented, you can see it's definitely not 1:

* (setq flt 1.0000000000000002d0)

1.0000000000000002d0
* (multiple-value-setq (mant exp sign) (integer-decode-float flt)) exp sign

4503599627370497
*
-52
*
1
* (setq mantlo (- mant 1/2))

9007199254740993/2
* (setq manthi (+ mant 1/2))

9007199254740995/2
* (setq denom (expt 2 (- exp)))

4503599627370496
* (setq int (intervals-divide (list mantlo manthi) denom))

(9007199254740993/9007199254740992 9007199254740995/9007199254740992)
* (show-interval int)
1.0d0
1.0000000000000004d0
1.000000000000000[1..4]
NIL
*       (setq guess (interval-to-sqrt-firstguess+bounds int))

6/5
*       (setq sqrtint (interval+guess-to-sqrt-bounds int guess))

(3/4 6/5)
*       (show-interval sqrtint)
0.75d0
1.2d0
0.[7..C]
NIL
* (setq guessres (ratio-interval-to-simplest-ratio-within sqrtint))

1
* (setq sqrtint (interval-newton-sqrt int sqrtint))
0.75d0
1.2d0
0.[7..C]
[-----------------------------------!!!!!!!!!!!!!!!!!!!!!!!!!!!--------]
0.975d0
1.1428571428571428d0
0.[9..C]
[--------------------------------------!-------------------------------]
1.0d0
1.0000000000000004d0
1.000000000000000[0..4]
(1 9007199254740995/9007199254740992)
* (setq guessres (ratio-interval-to-simplest-ratio-within sqrtint))

1
* (setq sqrtint (interval-newton-sqrt int sqrtint))
1.0d0
1.0000000000000004d0
1.000000000000000[0..4]

** WARNING: New interval not included within old interval,
**  so I'll have to intersect the two to get a valid interval refinement.
[!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!-----------------------------------]
1.0d0
1.0d0
1.0000000000000000[0..H]
[!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!]
1.0d0
1.0000000000000004d0
1.000000000000000[0..4]
Interval didn't shrink enough:
Old: 1.000000000000000[0..4]
New: 1.000000000000000[0..4]
Refinement (i.e. shrinkage) factor = 1.0
** RETURN to go on:

(1 9007199254740995/9007199254740992)
* (setq guessres (ratio-interval-to-simplest-ratio-within sqrtint))

1
* (setq sqrtint (interval-newton-sqrt int sqrtint :simplify nil))
1.0d0
1.0000000000000004d0
1.000000000000000[0..4]

** WARNING: New interval not included within old interval,
**  so I'll have to intersect the two to get a valid interval refinement.
[!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!-----------------------------------]
1.0d0
1.0d0
1.0000000000000000[0..H]
Interval didn't shrink enough:
Old: 1.000000000000000[0..4]
New: 1.0000000000000000[0..H]
Refinement (i.e. shrinkage) factor = 2.0
** RETURN to go on:

(1 18014398509481987/18014398509481984)
* (setq guessres (ratio-interval-to-simplest-ratio-within sqrtint))

1
* (setq sqrtint (interval-newton-sqrt int sqrtint :simplify nil))
1.0d0
1.0d0
1.0000000000000000[0..H]

** WARNING: New interval not included within old interval,
**  so I'll have to intersect the two to get a valid interval refinement.
[-----------!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!]
1.0d0
1.0d0
1.0000000000000000[2..H]
Interval didn't shrink enough:
Old: 1.0000000000000000[0..H]
New: 1.0000000000000000[2..H]
Refinement (i.e. shrinkage) factor = 1.2
** RETURN to go on:

(36028797018963972/36028797018963971 18014398509481987/18014398509481984)
* (setq guessres (ratio-interval-to-simplest-ratio-within sqrtint))

6004799503160663/6004799503160662

Note that last line, whereby the simplest ratio within the interval
is no longer 1, but rather something much more complicated, which implies 1
is no longer in the interval.

Or note the line slightly earlier that says:
1.0000000000000000[2..H]
which clearly bounds the value away from exactly 1.

Note also the corresponding picture line:
[-----------!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!]
which shows that particular interval refinement where
the left end (lower bound) was exactly 1 but now is larger than 1:


Would you like me to set up this software as a
Web-server-side application so that you-all can play
with it?

Urgent note: I'm in desperate financial circumstances.
I desperately need a source of income to avoid becoming
homeless. Please refer me to somebody who has money to
hire me for writing software.