[PDF] How do you compute the midpoint of an interval?

View - Download How do you compute the midpoint of an interval?

Previous PDF

Next PDF

>G A/, ?H@yy8dee9R ?iiTb,ff?HXb+B2M+2f?H@yy8dee9Rpk am#KBii2/ QM Rd T` kyR9 >GBb KmHiB@/Bb+BTHBM`v QT2M ++2bb `+?Bp2 7Q` i?2 /2TQbBi M/ /Bbb2KBMiBQM Q7 b+B@

2MiB}+ `2b2`+? /Q+mK2Mib- r?2i?2` i?2v `2 Tm#@

HBb?2/ Q` MQiX h?2 /Q+mK2Mib Kv +QK2 7`QK

i2+?BM; M/ `2b2`+? BMbiBimiBQMb BM 6`M+2 Q` #`Q/- Q` 7`QK Tm#HB+ Q` T`Bpi2 `2b2`+? +2Mi2`bX /2biBMû2 m /ûT¬i 2i ¨ H /BzmbBQM /2 /Q+mK2Mib b+B2MiB}[m2b /2 MBp2m `2+?2`+?2- Tm#HBûb Qm MQM-

Tm#HB+b Qm T`BpûbX

>Qr /Q vQm +QKTmi2 i?2 KB/TQBMi Q7 M BMi2`pH\ ?6?`?û?/?û?`?B?+ ?:?Q?m??H??`?/ hQ +Bi2 i?Bb p2`bBQM,

6`û/û`B+ :QmH`/X >Qr /Q vQm +QKTmi2 i?2 KB/TQBMi Q7 M BMi2`pH\X *J h`Mb+iBQMb QM Ji?2@

KiB+H aQ7ir`2- kyR9- 9y UkV- RyXRR98fk9Nj33kX ?H@yy8dee9Rpk A

How do you compute the midpoint of an interval?

FR

EDERIC GOUALARD, CNRS, LINA, UMR 6241

The algorithm that computes the midpoint of an interval with oating-point bounds requires some careful

devising to handle all possible inputs correctly. We review several implementations from prominent C/C++

interval arithmetic packages and analyze their potential failure to deliver the expected results. We then

show how to amend them to avoid common pitfalls. The results presented are also relevant to non-interval

arithmetic computation such as the implementation of bisection methods. Enough background on IEEE 754

oating-point arithmetic is provided for this paper to serve as a practical introduction to the analysis of

oating-point computation.

Categories and Subject Descriptors: G.1.0 [General]: Computer arithmetic; Error analysis; Interval arith-

metic General Terms: Algorithms, Reliability, Experimentation

Additional Key Words and Phrases:

oating-point number, IEEE 754 standard, interval arithmetic, mid- point, rounding error

1. INTRODUCTION

In his 1966 report [

Forsythe 1966

] \How do you solve a quadratic equation?", Forsythe con- siders the seemingly simple problem of reliably solving a quadratic equation on a computer using oating-point arithmetic. Forsythe's goal is both to warn a large audience away from unstable classical textbook formulae as well as get them acquainted with the characteristics of pre-IEEE 754 standard oating-point arithmetic, a dual objective shared by his later paper \Pitfalls in Computation, or Why a Math Book isn't Enough"[Forsythe 1970]. Following Forsythe's track, we consider here the problem of computing a good approxi- mation to the midpoint between two oating-point numbers. We strive to provide both a reliable algorithm for midpoint computation and an introduction to oating-point compu- tation according to the IEEE 754 standard [

IEEE 1985

1. Given two real numbersaandbfromR, witha6b, the midpointm(I) of the closed intervalI= [a;b] =fx2Rja6x6bgcan easily be obtained with the straightforward formula: m ([a;b]) =a+b2 :(1)1 A new version of the original standard from 1985 was approved and released in 2008 [

IEEE 2008

]. Most prominently, it considers both binary and decimal representations of oating-point numbers, while the 1985

version only considered a binary representation. In this paper, we only consider the features standardized

by the 1985 version|by far the most used to this day|for the sake of simplicity.Part of the work presented here was funded byIF CPAR/CEFIPRAPro ject4502-1.

Author's address: F. Goualard, Universite de Nantes. Nantes Atlantique Universite. CNRS, LINA, UMR

6241, 2 rue de la Houssiniere, BP 92208, F-44322 NANTES CEDEX 3.

Permission to make digital or hard copies of part or all of this work for personal or classroom use is

granted without fee provided that copies are not made or distributed for prot or commercial advantage

and that copies show this notice on the rst page or initial screen of a display along with the full citation.

Copyrights for components of this work owned by others than ACM must be honored. Abstracting with

credit is permitted. To copy otherwise, to republish, to post on servers, to redistribute to lists, or to use any

component of this work in other works requires prior specic permission and/or a fee. Permissions may be

requested from Publications Dept., ACM, Inc., 2 Penn Plaza, Suite 701, New York, NY 10121-0701 USA, fax +1 (212) 869-0481, or permissions@acm.org.

©YYYY ACM 0098-3500/YYYY/01-ARTA$10.00

DOI 10.1145/0000000.0000000 http://doi.acm.org/10.1145/0000000.0000000 ACM Transactions on Mathematical Software, Vol. V, No. N, Article A, Publication date: January YYYY.

A:2F. Goualard

But, ifaandbare two

oating-point numbersfrom a setFofnite oating-point numbers, the suma+bmay not be a oating-point number itself, and we therefore have to take care of rounding it correctly to ensure that our oating-point implementation of m(I) does not violate the fundamental property: m (I)2I;(2) viz., that the oating-point midpoint of an intervalIshould belong toI. This property is a prerequisite to use the midpoint operator in the dichotomic search process of some numerical algorithm, even though it is not sucient to ensure that the search space is separated into sub-spaces of the same size

2. There are, however, problems

that require more from a midpoint operator,viz., that it computes the oating-point number closest to the real midpoint ofI3:

8v2F:jcvj>jcm(I)j;withc=a+b2

;c2R:(3)

In interval arithmetic [

Moore 1966

Alefeld and Herzb erger1983

Ha yes2003

], some ex- amples of such problems are: Ona verage,the least o verestimationin dening a centered form[Moore 1966;Ratsc hek 1980
] of a rational function is obtained by choosing a value closest to the midpoint as developing point [

Ratschek and Rokne 1984

, p. 41];

Forthe Kra wczykop erator[

Neumaier 1990

, p. 177], the midpoint is the value that gives the tightest result over all other points of the interval under consideration [

Neumaier 1990

Theo. 5.1.9, p. 178].

Inclusion and precision are not, however, the only possible problems for oating-point im- plementations of the midpoint operator: what should be, for example, the midpoint of an empty interval? As we shall see in Section 2 , the IEEE 754 standard comes to the rescue here by dening a specialNot-a-Numbervalue (NaN) that should be used for such situation; another diculty lies in the fact that, as we will see in the next section, one or both bounds of an interval with oating-point bounds may be an innite value \1". What should then be the midpoint of an interval of the form [a;+1], [1;b], or even [1;+1]? Though the midpoint is mathematically undened, many applications (e.g., constraint solvers such as IBEX [

Chabert et al. 2012

] or Realpaver [

Granvilliers and Benhamou 2006

]) typically expect nite values and would behave incorrectly, should the midpoint operator return a NaN for these intervals. As a consequence, the current draft of the future IEEE P1788 standard on interval arithmetic denes the midpoint operator as follows:8>>>< >>:m (?) is NaN; m ([1;b]) =realmax; m ([a;+1]) =realmax; m ([1;+1]) = 0; m ([a;b]) =rndnrh(a+b)=2i;for (a;b)2F2;(4) whererndnrhxiis dened in the next section as returning the oating-point value closest to the real numberx, andrealmaxis the largest positive element ofF. For non-empty intervals with nite bounds, Equation ( 3 ) implies Equation ( 2 ). However, some implementations may choose to ensure only the latter and a relaxation of the former. Due to its many uses (centered forms, Krawczyk operator, Newton operator [

Moore 1966

computation of preconditioning matrices [

Kearfott 1990

], representation of intervals in the2

Where the sizew(I) of an intervalI= [a;b] with

oating-point bounds is dened asw(I) =ba.

3If the real midpoint is midway between two

oating-point numbers, there exist two possible values for m (I) that satisfy Equation (3). ACM Transactions on Mathematical Software, Vol. V, No. N, Article A, Publication date: January YYYY. How do you compute the midpoint of an interval? A:3 midpoint/radius format [

Rump 1999a

], ...), the midpoint operator is a staple of interval arithmetic libraries. It is, therefore, paramount that its oating-point implementation at least satises Equation ( 2 ). Accuracy, as stipulated by Equation ( 3 ), is also desirable, as seen in the examples above. Nevertheless, we will see in Section 3 that s omeform ulae implemented in popular C/C++ interval libraries may not ensure even the containment requirement for some inputs. However, this study should not be considered as a report on the quality of these libraries since dierent implementations of the midpoint operator might suit specic needs, and some libraries were dened with no provision to support intervals with innite bounds in the rst place. We are only concerned here with the properties of the formulae these libraries implement.

In Section

3 , we analyze the various formulae both theoretically and practically; contrary to most expositions, we consider the impact of both over ow and under ow on the accuracy and correctness of the formulae. The error analysis conducted in this paper requires slightly more than a simple working knowledge of oating-point arithmetic as dened by the IEEE 754 standard. As a conse- quence, the basic facts on oating-point arithmetic required in Section 3 are presen tedin

Section

2 for the sak eof self-con tainedness. It turns out that the study of the correct implementation of a oating-point midpoint operator may serve as a nice introduction to many important aspects of oating-point computation at large: the formulae studied are simple enough for their analysis to be easily understandable, while the set of problems raised is suciently broad in scope to be of general interest. We then hope that this paper will be valuable as both a technical presentation of reliable, accurate, and fast methods to compute a midpoint as well as an introduction to error analysis of oating-point formulae.

2. FLOATING-POINT ARITHMETIC IN A NUTSHELL

According to the IEEE 754 standard [

IEEE 1985

], a oating-point number'is represented by a sign bits, a signicandm(wheremis a bit string of the form \0:f" or \1:f", withf thefractional part) and an integral exponentE: '= (1)sm2E:(5) The IEEE 754 standard denes several formats varying in the number of bitsl(f) and l(E) allotted to the representation offandE, the most prominent ones beingsingle preci- sion|l(E);l(f)=8;23|anddouble precision|l(E);l(f)=11;52. We will also use for pedagogical purposes anad hocIEEE 754 standard-complianttinyprecision format|l(E);l(f)=3;3. Wherever possible, the signicand must be of the form \1:f" since it is the form that stores the largest number of signicant gures for a given size ofm: '=0:0110120 =0:110121 =1:10122: Floating-point numbers with such a signicand are callednormal numbers. Such prevalence is given to normal numbers that the leading \1" is left implicit in the representation of an IEEE 754 number, and only the fractional partfis stored (see Figure1 ). The exponentEis a signed integer stored as a biased exponente=E+ bias, with bias = 2 l(E)11. The biased exponenteis a non-negative integer that ranges from e min= 0 toemax= 2l(E)1. However, for the representation of normal numbers,Eonly ranges fromEmin= (eminbias) + 1 toEmax=emaxbias1 because the smallest and largest values are reserved for special purposes (see below). As an example of what ACM Transactions on Mathematical Software, Vol. V, No. N, Article A, Publication date: January YYYY.

A:4F. Goualard

precedes, the bias for thetinyformat is equal to 3,eranges from 0 to 7, andEranges from

2 to +3.b6

s b5b4b3 e b2b1b0 fFig. 1. Binary representation as a seven bit string of atiny oating-point number. Consider the binary number= 1:0011. It cannot be represented as atiny oating-point number since its fractional part has four bits, and thetinyformat has room for only three.

It therefore has to be rounded to the

oating-point number hiaccording to one of four rounding directions

4(see Figure2 ):

Roundingto ward0:

hi=rndzrhi;

Roundingto neares t-even:

hi=rndnrhi;

Roundingto ward1:

hi=rnddnhi;

Roundingto ward+ 1:

hi=rnduphi.0+1 rndzrhirndnrhirnddnhirnduphi Fig. 2. Rounding a real number according to the IEEE 754 standard. Note the use of angles \hi" instead of more conventional parentheses for the rounding operators. They are used to express the fact that each individual value and/or operation comprising the expression is individually rounded according to the leading operator. For example: h1+2i (1) + (2));8(1;2)2R2: When rounding to nearest, ifis equidistant from two consecutive oating-point numbers, it is rounded to the one whose rightmost bit of the fractional part is zero (the \even" number). It is not possible to represent 0 as a normal number. Additionally, consider the number '= 0:0001120:

To store it as atiny

oating-point normal number requires shifting the leftmost \1" of the fractional part to the left of the radix point: '= 1:10024:4

The actual direction chosen may depend on the settings of the Floating-Point Unit at the time, or alter-

natively, on the machine instruction used, for some architectures. ACM Transactions on Mathematical Software, Vol. V, No. N, Article A, Publication date: January YYYY. How do you compute the midpoint of an interval? A:5 However, doing so requires an exponent smaller thanEmin. It is nevertheless possible to represent', provided we accept to store it with a \0" to the left of the radix point: '= 0:01122: Numbers with a signicand of the form \0:f" (withf6= 0) are calledsubnormal numbers. Their introduction is necessary to reduce the large gap that would otherwise occur around

0 (compare Figure

quotesdbs_dbs15.pdfusesText_21

[PDF] exercice graphique 5eme

[PDF] test statistique seconde

[PDF] qu'est ce que le contrôle stratégique

[PDF] le système d'information contribue-t-il ? l'efficacité de la prise de décision

[PDF] maths 1ere es suites exercices corrigés

[PDF] une personne loue une maison ? partir du 1er janvier 1991

[PDF] controle sur les suites numériques 1ere s

[PDF] ds suite arithmétique 1s

[PDF] ds suite arithmétique et géométrique

[PDF] ds recurrence

[PDF] devoir raisonnement par recurrence

[PDF] controle recurrence ts

[PDF] dm de maths terminale s recurrence

[PDF] calculer u1 et u2 la suite un est elle arithmétique géométrique

[PDF] ds suites arithmétiques et géométriques 1ere s