DEMOLONG Long numbers and intervals

Contents

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
format long

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.000000000000 * 10^+0000
long V = 
 -1.000000000000 * 10^+0000
  3.000000000000 * 10^+0000

Long vectors are always column vectors (or forced to column vectors).

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 = 
  1.428571428571 * 10^-0001
long x = 
  1.428571428571428571428571428571428571428571 * 10^-0001
long x = 
  1.428571428571428571428571428571428571428571428571184204 * 10^-0001

Output of long numbers is not rigorous. All but a few of the last digits are correct.

Arithmetic operations

Long operations +,-,*,/ and ^ are supported. Note that operations on vectors are always performed elementwise.

x = [ long(3) ; -7 ]
x*x
long x = 
  3.000000000000 * 10^+0000
 -7.000000000000 * 10^+0000
long ans = 
  9.000000000000 * 10^+0000
  4.900000000000 * 10^+0001

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.141592653589793
intval IntPi = 
   3.14159265358979
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