DEMOLONG Long numbers and intervals
Contents
- Definition of long numbers
- Conversion
- Output of long numbers
- Arithmetic operations
- Matrix operations
- Output of long intervals
- Interval and non-interval operations
- Conversion between long and double
- Long numbers with error term
- Specifying extremely small errors
- Taylor series: an example
- Ill-conditioned polynomials
- Sample programs
- Enjoy INTLAB
The purpose of the long toolbox was to compute rigorous bounds for the value of certain standard functions. Those values are needed to initialize the INTLAB system. The long toolbox is slow, but fast enough to do the job.
setround(0) % set rounding to nearest longprecision(0); % default option
Recently I decided to make a true toolbox out of it supporting matrices, and many matrix operations.
Definition of long numbers
Long numbers are stored in midpoint radius representation. The midpoint is stored in an array of precision to be specified, the error is stored in one number. Long numbers or vectors are generated by the constructor long:
x = long(7) V = long([ -1 , 3 ])
long x =
7
long V =
-1 3
Conversion
Conversion of double to long by the constructor "long" is always into the current internal precision.
The internal precision may be specified by "longprecision". The call without input parameter gives the current working precision in decimals, the call with input parameter sets working precision.
p = longprecision longprecision(50)
p =
27
ans =
55
Now, the statement "x = long(7)" generates a long number with approximately 50 decimal digits. These are approximately 50 digits because the internal representation is to some base beta, a power of 2.
Output of long numbers
Output is usually to a little more than double precision. If you want to see more digits, say k, use "display" with second parameter equal to k. To see all digits, use k=0.
longprecision x = 1/long(7) display(x,40) display(x,0)
ans =
55
long x =
0.1429
long x =
0.1429
long x =
0.1429
Output of long numbers is not rigorous. All but a few of the last digits are correct.
Arithmetic operations
Long operations +,-,*,/ and ^ are supported. For example
format longdouble
x = long(reshape(1:6,2,3))
p = x'*x
===> Display long-variables as doubles
long x =
1 3 5
2 4 6
long p =
5 11 17
11 25 39
17 39 61
The output has been changed to double precision; when printing as long numbers, always the vectorized data is seen:
format longlong
p
===> Display long-variables by long representation long p(:) = 5.000000000000 * 10^+0000 1.100000000000 * 10^+0001 1.700000000000 * 10^+0001 1.100000000000 * 10^+0001 2.500000000000 * 10^+0001 3.900000000000 * 10^+0001 1.700000000000 * 10^+0001 3.900000000000 * 10^+0001 6.100000000000 * 10^+0001
Matrix operations
The determinant of the Pascal matrix is known to be always 1:
n = 5; P = pascal(n) det(P)
P =
1 1 1 1 1
1 2 3 4 5
1 3 6 10 15
1 4 10 20 35
1 5 15 35 70
ans =
1
When increasing the dimension, inevitable rounding errors may change the result, eventually producing negative results:
N = 1:20; D = zeros(size(N)); for n=N D(n) = det(pascal(n)); end [ N ; D ]
ans =
Columns 1 through 7
1.0000 2.0000 3.0000 4.0000 5.0000 6.0000 7.0000
1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000
Columns 8 through 14
8.0000 9.0000 10.0000 11.0000 12.0000 13.0000 14.0000
1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000
Columns 15 through 20
15.0000 16.0000 17.0000 18.0000 19.0000 20.0000
0.9997 1.0131 1.1568 0.8003 3.7724 -131.2330
This is avoided when changing to long precision:
format short longprecision(40); N = 1:20; D = zeros(size(N)); for n=N P = long(pascal(n)); [L,U,p] = luwpp(P); D(n) = double(parity(p)*prod(diag(U))); end [ N ; D ]
ans =
Columns 1 through 7
1.0000 2.0000 3.0000 4.0000 5.0000 6.0000 7.0000
1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000
Columns 8 through 14
8.0000 9.0000 10.0000 11.0000 12.0000 13.0000 14.0000
1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000
Columns 15 through 20
15.0000 16.0000 17.0000 18.0000 19.0000 20.0000
1.0000 1.0000 1.0000 1.0000 1.0000 1.0000
Still some rounding errors can be seen; that is avoided by increasing the precision:
longprecision(50); N = 1:20; D = zeros(size(N)); for n=N P = long(pascal(n)); [L,U,p] = luwpp(P); D(n) = double(parity(p)*prod(diag(U))); end [ N ; D ]
ans =
Columns 1 through 13
1 2 3 4 5 6 7 8 9 10 11 12 13
1 1 1 1 1 1 1 1 1 1 1 1 1
Columns 14 through 20
14 15 16 17 18 19 20
1 1 1 1 1 1 1
We use the universal routine "luwpp", LU-decomposition with partial pivoting, which can be used with a variety of data types, among them intval, affari, gradient, hessian, fl, or gfp. The latter is another possibility to produce the correct determinant - at least modulo the given prime:
gfpinit(67108859); N = 1:20; D = zeros(size(N)); for n=N P = gfp(pascal(n)); [L,U,p] = luwpp(P); D(n) = double(parity(p)*prod(diag(U))); end [ N ; D ]
Galois field prime set to 67108859
ans =
Columns 1 through 13
1 2 3 4 5 6 7 8 9 10 11 12 13
1 1 1 1 1 1 1 1 1 1 1 1 1
Columns 14 through 20
14 15 16 17 18 19 20
1 1 1 1 1 1 1
Output of long intervals
The display routine takes uncertainties into account. Only the correct digits plus some extra are displayed.
longprecision(50); x = long(1)/37; display(x,0) for i=1:100 x=x*x; x=x*37; end display(x,0)
long x(:) = 2.702702702702702702702702702702702702702702702702480790 * 10^-0002 long x(:) = 2.702702702702702702702607 * 10^-0002
Interval and non-interval operations
Computing with uncertainties may be switched off by
longinit('WithoutErrorTerm') longprecision(50); x = long(1)/37; display(x,0) for i=1:100 x=x*x; x=x*37; end display(x,0)
===> Long arithmetic computations without error terms long x(:) = 2.702702702702702702702702702702702702702702702702480790 * 10^-0002 long x(:) = 2.702702702702702702702607059944376258877795697765239694 * 10^-0002
In this case all digits including incorrect ones are displayed. Computing without error term is a usual long precision arithmetic with specified precision. Note that scalar operations suffer from quite some interpretation overhead.
Conversion between long and double
Conversion from long to double is approximately to nearest, conversion to interval is rigorous.
For example, in the following the function "longpi" calculates "pi" to the specified longprecision, "IntPi" is a true inclusion of the transcendental number "pi".
longinit('WithErrorTerm'); longprecision(100); Pi = longpi; display(Pi,0) flptPi = long2dble(Pi) IntPi = long2intval(Pi) format long infsup(IntPi)
===> Long arithmetic computations with error terms, i.e. valid long error bounds
long Pi(:) =
3.141592653589793238462643383279502884197169399375105820974944592307816406 \
286208998628034825342116 * 10^+0000
flptPi =
3.1416
intval IntPi =
3.1415
intval IntPi =
[ 3.14159265358979, 3.14159265358980]
Long numbers with error term
Long numbers may be specified with an explicit error term. For example,
longprecision(50); x = long(-1.5) display(x,0) x = addlongerror(x,1e-40) display(x,0)
long x(:) = -1.500000000000 * 10^+0000 long x(:) = -1.500000000000000000000000000000000000000000000000000000 * 10^+0000 long x(:) = -1.500000000000 * 10^+0000 long x(:) = -1.500000000000000000000000000000000000000000 * 10^+0000
defines x to be an interval with midpoint -1.5 and radius approximately 10^(-40). Only meaningful digits are stored and displayed.
Specifying extremely small errors
For very small errors leaving the range double precision floating point numbers, the error may be specified by the mantissa and the exponent of the error:
longprecision(50); x = long(2^-1000)^2; x = addlongerror(x,1,-620)
long x(:) = 8.709809816219 * 10^-0603
The final x, which is 2^(-2000), is afflicted with an error of 10^(-620).
Taylor series: an example
As an example, the following code computes the value of E = exp(x) by a Taylor series:
p = 100; longprecision(p); x = -30; t = 1; T = long(1); E = T; k = 0; while abs(t)>10^(-p) k = k+1; t = t*x/k; T = T*x/k; E = E + T; end k exp(x) display(E,0)
k =
225
ans =
9.357622968840175e-14
long E(:) =
9.357622968840174604915832223378706744958322688935880416413318619960842833 * 10^-0014
Note that for large negative values of x there quite some cancellation. This can be seen by
x = 30; t = 1; T = long(1); E = T; k = 0; while abs(t)>10^(-p) k = k+1; t = t*x/k; T = T*x/k; E = E + T; end k 1/exp(x) display(1/E,0)
k =
225
ans =
9.357622968840175e-14
long ans(:) =
9.357622968840174604915832223378706744958322688935880416413318619960842833 \
767616873665639476198700545791 * 10^-0014
Ill-conditioned polynomials
Consider the following polynomial:
P = inline(' 4999*x.^6 - 200*x.^5 + 102*x.^4 - 2*x.^3 - 2500*x.^2 + 100*x - 1 ')
P =
Inline function:
P(x) = 4999*x.^6 - 200*x.^5 + 102*x.^4 - 2*x.^3 - 2500*x.^2 + 100*x - 1
This is an example of a Bugeaud-Mignotte polynomial. The general form is
( X^n - aX + 1 )^k - 2X^(nk-k)(aX-1)^k
where a>=10, n>=3 and k>=2.
Those polynomials are constructed to have a pair of very close real roots near c=1/a+1/a^(n+1). A graph using floating-point arithmetic near c looks as follows:
e = 3e-8; c = 1/50+1/50^4; x = c*(1+linspace(-e,e)); close plot(x,P(x),x,0*x)
From the graph it is not clear whether the polynomial has no, a double or two real roots in the interval c*[1-e,1+e]. An evaluation using the long package yields the following:
y = long2dble(P(long(x))); close plot(x,y,x,0*x)
Sample programs
For sample programs using long numbers, see for example the source codes of long\longpi.m or long\@long\exp.m .
Enjoy INTLAB
INTLAB was designed and written by S.M. Rump, head of the Institute for Reliable Computing, Hamburg University of Technology. Suggestions are always welcome to rump (at) tuhh.de