.FP palatino.de ]C\&[\\$1]\\$2...if \nZ=0 .so real.ref.EQdelim $$.EN.TLReal Inferno.AU.I "Eric Grosse".AI.I "Lucent Technologies, Bell Labs".I "Murray Hill NJ 07974 USA".I "ehg@bell-labs.com".\"date{19 Aug 1996, minor revisions 7 Jan 1998}.FSPreviously appeared in R.F. Boisvert (editor),.ft IThe Quality of Numerical Software: Assessment andEnhancement: Proceedings of the IFIP TC2/WG2.5Working Conference on the Quality of NumericalSoftware, Oxford,United Kingdom, 8-12 July 1996,.ft RChapman Hall,London,1997 (pp. 270-279)..FE.ABInferno is an operating system well suited to applications that need to beportable, graphical, and networked. This paper describes the fundamentalfloating point facilities of the system, including: tight rules onexpression evaluation, binary/decimal conversion, exceptions and rounding,and the elementary function library..AE.PPAlthough the focus of Inferno is interactive media, its portability acrosshardware and operating platforms, its relative simplicity, and its strengthin distributed computing make it attractive for advanced scientificcomputing as well. Since the appearance of a new operating system is arelatively uncommon event, this is a special opportunity for numericalanalysts to voice their opinion about what fundamental facilities they need.The purpose of this short paper is to describe numerical aspects of theinitial release of Inferno, and to invite comment before the tyranny ofbackward compatibility makes changes impossible..PPAn overview of Inferno is given by Dorward et. al..]C "Inferno" ,but for our immediate purposes it may suffice to say that Inferno plays therole of a traditional operating system (with compilers, process control,networking, graphics, and so on) but can run either on bare hardware or ontop of another operating system like Windows95 or Unix. Programs for.I "Inferno"are written in the language.I "Limbo"and compiled tomachine-independent object files for the.I "Dis"virtualmachine, which is then implemented with runtime compilation for bestperformance. Files are accessible over networks using the.I "Styx"protocol; together with the presentation of most system resources as filesand the manipulation of file namespaces, this permits integration of acollection of machines into a team. Limbo looks somewhat like a mixture of Cand Pascal, augmented by modules (to cope with the namespace and dynamicloading needs of large programs) and by a channel facility for convenient(coarse-grain) parallel programing. Array references are bounds-checked andmemory is garbage collected..PPThe rest of this paper covers the fundamental floating point environmentprovided by the Limbo compiler and.I "math"module, the ``elementaryfunctions,'' and finally some comments on why particular definitions werechosen or why certain facilities were included or excluded. This discussionassumes the reader is familiar with scientific computing in general and theIEEE floating point standard in particular..NH 1Floating point.PPIn Limbo, arithmetic on literal and named constants is evaluated at compiletime with all exceptions ignored. Arithmetic on variables is left by thecompiler to runtime, even if data path analysis shows the value to be acompile time constant. This implies that tools generating Limbo source mustdo their own simplification, and not expect the compiler to change $x/x$into $1$, or $-(y-x)$ into $x-y$, or even $x-0$ into $x$. Negation $-x$changes the sign of $x$; note that this not the same as $0-x$ if $x=0$..PPThe compiler may perform subexpression elimination and other forms of codemotion, but not across calls to the mode and status functions. It respectsparentheses. The evaluation order of $a+b+c$ follows the parse tree and istherefore the same as for $(a+b)+c$. These are the same rules as for Fortranand C..PPContracted multiply-add instructions (with a single rounding) are notgenerated by the compiler, though they may be used in the native.SM BLASlibraries. All arithmetic follows the IEEE floating point standard.]C "IEEEfp" ,except that denormalized numbers may not be supported; see thediscussion in section 3..PPThe most important numerical development at the language level recently hasbeen accurate binary/decimal conversion.]C "Clinger".]C "Gay".]C "SteeleWhite" .Thus printing a real using.CW "\%g"and readingit on a different machine guarantees recovering identical bits. (Limbo usesthe familiar.I "printf"syntax of C, but checks argument types againstthe format string at compile time, in keeping with its attempt to help theprogrammer by stringent type checking.) A good.I "strtod/dtoa"is,unfortunately, 1700 lines of source (15kbytes compiled), though with modestaverage runtime penalty. This code must be used in the compiler so thatcoefficients are accurately transferred to bytecodes. Smaller, faster, butsloppier, runtimes will also be provided when mandated by limited memory andspecialized use. However, programmers may assume the features described inthis paper are present in all Inferno systems intended for general computing..PPEach thread has a floating point control word (governing rounding mode andwhether a particular floating point exception causes a trap) and a floatingpoint status word (storing accumulated exception flags). Functions.I "FPcontrol"and.I "FPstatus"copy bits to the control or status word, inpositions specified by a mask, returning previous values of the bits. .I "getFPcontrol"and.I "getFPstatus"return the words unchanged..PPThe constants.I "INVAL, ZDIV, OVFL, UNFL, INEX"are non-overlappingsingle-bit masks used to compose arguments or return values. They stand forthe five IEEE exceptions:.IP •``invalid operation'' ($0/0$,$0 * infinity $,$ infinity - infinity $,$sqrt{-1}$).IP •``division by zero'' ($1/0$),.IP •``overflow'' ($1.8e308$).IP •``underflow'' ($1.1e-308$).IP •``inexact'' ($.3*.3$)..PPThe constants.I "RND_NR, RND_NINF, RND_PINF, RND_Z"are distinctbit patterns for ``round to nearest even'', ``round toward $-{infinity} $'',``round toward $+{infinity} $'', ``round toward $0$'', any of which can be setor extracted from the floating point control word using.I "RND_MASK" .For example,.IP •to arrange for the program to tolerate underflow,.I "FPcontrol(0,UNFL).".IP •to check and clear the inexact flag,.I "FPstatus(0,INEX).".IP •to set directed rounding,.I "FPcontrol(RND_PINF,RND_MASK).".PPBy default,.I "INEX"is quiet and.I "OVFL, UNFL, ZDIV,"and.I "INVAL"are fatal. By default, rounding is to nearest even, and libraryfunctions are entitled to assume this. Functions that wish to use quietoverflow, underflow, or zero-divide should either set and restore thecontrol register themselves or clearly document that the caller must do so.The ``default'' mentioned here is what a Limbo program gets if started in afresh environment. Threads inherit floating point control and status fromtheir parent at the time of spawning and therefore one can spawn a ``roundtoward 0'' shell and re-run a program to effortlessly look for roundinginstabilities in a program..NH 1Elementary functions.PPThe constants.I "Infinity, NaN, MachEps, Pi, Degree"are defined. SinceInferno has thorough support of Unicode, it was tempting to name these $infinity $, $ε $, $π $, and °, but people (or rather, theirfavorite text editing tools) may not be ready yet for non-\s-2ASCII\s0source text..I "Infinity"and.I "NaN"are the positive infinityand quiet not-a-number of the IEEE standard, double precision..I MachEpsis $2 sup {-52}$, the unit in the last place of the mantissa $1.0$.The value of.I "Pi"is the nearest machine number to themathematical value $π $..I "Degree"is$"Pi" / 180$..PPThree useful functions.I "fdim, fmax, fmin"are adopted from theNumerical C extensions.]C "NumerC" .The unusual one of these, oftendenoted $(x-y) sub {+}$, is defined by $roman "fdim" ( x , y )=x-y$ if $x > y$, else $0$. The compiler may turn these into efficient machine instruction sequences,possibly even branch-free, rather than function calls. There are two almostredundant mod functions:.I "remainder(x,y)"is as defined by the IEEEstandard (with result $roman "|" r roman "|" <= y/2$);.I "fmod(x,y)"is $x roman "mod" y$,computed in exact arithmetic with $0<= r<y$. Limbo has a ``tuple'' type,which is the natural return value in the call $(i,f)= roman "modf" ( x )$ tobreak $x$ into integer and fractional parts. The function.I "rint"rounds to an integer, following the rounding mode specified in the floatingpoint control word..PPFor a good-quality, freely-available elementary function library,.I "math"uses the IEEE subset of.I "fdlibm".]C "fdlibm" .Of course, aconforming implementation may use entirely different source, but must takecare with accuracy and with special arguments. There are the customarypower, trigonometric, Bessel, and erf functions, and specialized versions $roman "expm1"( x )=e sup x - 1$, $roman "log1p" ( x )=log ( 1 + x )$. An additional function $roman "pow10"( n ) = 10 sup n$ is defined; in the default implementation this isjust fdlibm's $roman "pow" ( 10. , n )$ but it is provided so that separatetrade-offs of accuracy and simplicity can be made.]C "MacIlroy" ..I "fdlibm"uses extra precise argument reduction, so the computed $sin (n*Pi)$is small but nonzero. If demands warrant, degree versions of thetrigonometric functions will be added, but for now the style $sin (45*Degree)$ is used.The library also provides IEEE functions.I "ilogb, scalbn, copysign, finite, isnan,"and.I "nextafter" ..PPThe functions.I "dot, norm1, norm2, iamax, gemm"are adopted from the.SM BLAS.]C "blas"to get tuned linear algebra kernels foreach architecture, possibly using extra-precise accumulators. These aredefined by $sum {{x sub i}{y sub i}}$, $sum roman | {x sub i} roman | $, $ sqrt{sum { x sub {i sup 2}}} $, $i$ suchthat $roman | {x sub i} roman | = roman max $, and $C= alpha AB + beta C$ with optional transposes on $A$and $B$. Since Limbo has only one floating-point type, there is no need herefor a precision prefix. Limbo array slices permit the calling sequences tobe more readable than in Fortran77 or C, though restricted to unit stride.This encourages better cache performance anyway. The matrix multiplyfunction.I "gemm"remains general stride (and is the foundation forother operations.]C "Kagstrom" )..PPLimbo is like C in providing singly-subscripted arrays with indexingstarting at 0. Although Limbo offers arrays of arrays, as in C, forscientific work a better choice is to adopt the style of linearizingsubscripts using Fortran storage order. This promotes easier exchange ofdata with other applications and reuses effort in organizing loops toachieve good locality. In previous language work.]C "pine" ,we implementeda C preprocessor that allowed the programmer to choose a convenient origin(such as 1) and have it compiled into 0 for the base language; because wepassed arrays as dope vectors, we were even able to allow different originsfor the same array in calling and called functions. The main lesson welearned from that experience, however, was that permutations become anightmare when there is anything but dogmatic adherence to a single origin.So for an $m$ by $n$ matrix $A$, the programmer should use loops with $0<=i<m$ and $0<= j<n$ and access $A[i+m*j]$..PPFor interoperability with foreign file formats and for saving main memory inselected applications, functions are provided for copying bits between andreals and 32-bit or 64-bit IEEE-format values..PPFinally,.I "math"provides a tuned quicksort function.I "sort(x,p)"where.I "x"is a real array and.I "p"is an int array representinga 0-origin permutation. This function leaves the contents of.I "x"untouched and rearranges.I "p"so that $x[{p sub i}]<= x[p sub {i+1}]$. This isusually what one wants to do: sort an array of abstract data types based onsome key, but without the need to actually swap large chunks of memory..NH 1Rationale.PPThis section discusses why certain numerical features were included or not..NH 2Rounding modes and accumulated exceptions.PPDirected rounding is only needed in a very few places in scientificcomputing, but in those places it is indispensable. Accumulated floatingpoint exceptions are even more useful. User trap handling is a harderproblem, and may be worth leaving for later, possibly with a default``retrospective diagnostics'' log.]C "Kahan" ..PPNote that the exception masks must be architecture independent, since theyreside in the Limbo bytecodes, and therefore the implementation involves asmall amount of bit fiddling. Still, it is efficient enough to encourageuse. It would be difficult to port to a processor that only had staticrounding modes in instruction opcodes rather than the dynamic modelspecified in section 2 of the IEEE standard. Fortunately, the Alphadoes provide both models..NH 2Sudden underflow.PPSome processor vendors make supporting gradual underflow just too hard. (Onemust struggle upon the system trap to reconstruct exactly which instructionwas executing and what the state of the registers was. On the MIPS, it issaid to be 30 pages of assembler.) So Inferno supports denormalized numbersonly if the hardware makes this easy. Providing underflow that is correctbut very slow, as some systems do, is not necessarily doing the user a favor..PPTo determine portably if a particular system offers gradual underflow, maskoff UNFL and do trial arithmetic..NH 2Speed.PPComputers with slow (software) gradual underflow usually provide a fastflush-to-0 alternative. This often suffices, though there are importantexamples where it forces an uglier and slower coding style. A worsesituation is if the hardware uses system traps for Infinity and NaNarithmetic. The resulting slowdown will make otherwise excellent and naturalalgorithms run slowly.]C "Demmel" .Sadly, even some x86 implementationsthat do non-finite arithmetic in hardware, do it relatively slowly..PPWe considered providing syntax to declare a certain program scope withinwhich precise IEEE behavior was required, and relaxing the rules outsidesuch scopes.(The numerical C extensions.]C "NumerC" use pragmafor this purpose.)These scope declarations would need to be in thebytecodes, since significant optimization may be attempted by the runtimecompiler. After some discussion, and with some trepidation, it was agreedthat instead all compilers would be required to preserve the same result andstatus as for an unoptimized version..NH 2Comparison.PPThe standard C operators.CW ==.CW !=.CW "<".CW "<=".CW ">".CW ">="are the only comparisons provided, and they behave exactlylike the ``math'' part of Table 4 of the IEEE standard. Programs interestedin handling NaN data should test explicitly. This seems to be the way mostpeople program and leads to code more understandable to nonexperts. It istrue that with more operators one can correctly write code that propagatesNaNs to a successful conclusion\-but that support has been left for later.NaN(''tag'') should be added at that same time..NH 2Precision.PPAll implementations run exclusively in IEEE double precision. If thehardware has extra-precise accumulators, the round-to-double mode is setautomatically and not changeable, in keeping with Limbo's design to haveonly one floating point type. Extended precision hardware, if available, maybe used by the built-in elementary function and.SM BLASlibraries.Also, we contemplate adding a dotsharp function that would use a very longaccumulator for very precise inner products, independent of the order ofvector elements.]C "kulisch" .But reference implementations that use onlydouble precision, avoid contracted multiply-add, and evaluate in the order 1up to n will always be available for strict portability..PPAt the time the decision was made to restrict the system to 64-bit floatingpoint, Limbo integers were almost exclusively 32-bit and the consistencyargument to have a single real type was compelling. Now that Limbo has moreinteger types the decision might be reconsidered. But so many engineersneedlessly struggle with programs run in short precision, that offering itmay do as much harm as good. On most modern computers used for generalpurpose scientific computing, 64-bit floating point arithmetic is as fast as32-bit, except for the memory traffic. In cases where the shorter precisionwould suffice and memory is a crucial concern, the programmer shouldconsider carefully scaled fixed point or specialized compression. Toefficiently interoperate with data files that use the short format,programmers may use the provided realbits32 function. While there are surelyappropriate uses for a first-class 32-bit real type, for now we followKahan's sarcastic motto ``why use lead when gold will do?''.NH 2BLAS.PPThe few.SM BLASin the core library were chosen for readability and,in case of gemm, for optimization beyond what a reasonable compiler wouldattempt. We expect that compilers will (soon) be good enough that thedifference between compiling $y+=a*x$ and calling daxpy is small. Also, asmentioned above, dot and gemm might reasonably use combined multiply-add ora long accumulator in some optional implementations..NH 2$GAMMA ( x )$.PPTo avoid confusion with the C math library, which defined.I "gamma"as $ln GAMMA $, we offer only.I "lgamma"for now. This function and .I "modf"return an (int,real) tuple rather than assigning through aninteger pointer, in keeping with Limbo's design. The opportunity has beentaken to drop some obsolete functions like.I "frexp" .Other functionsare unchanged from the C math library..NH 2Future.PPA prototype preprocessor has been written to allow the scientific programmerto write $A[i,j]$ for an $A$ that was created as a $Matrix(m,n)$ and to havethe subscript linearization done automatically. Here $Matrix$ is an Limboabstract data type containing a real array and integers $m$, $n$, and columnstride $lda$ used as in typical Fortran calling sequences..PPThe Limbo compiler is soon expected to implement the type.I "complex" ..PPHigher level numerical libraries will also be provided, and although thattopic is beyond the scope of this paper, opinions about what should comefirst would be welcome..PPDistributed computing has not been mentioned here because it involvesrelatively few considerations specific to floating point computation.However, it may be worth noting that in the default environment (withunderflow trapped, so that presence or absence of denormalized numbers isnot significant) programs run independently on heterogeneous machinesnevertheless get precisely identical results, even with respect to threadscheduling. This implies that certain communication steps can be avoided,and that regression testing is considerably simplified..PPPlease direct comments on these numerical aspects of Inferno to Eric Grosse.More general technical comments can be directed to Vita Nuova.CW comments@vitanuova.com ). (I am grateful to Joe Darcy, Berkeley,to David Gay, Bell Labs, to David Hook, University of Melbourne,and to participants of the IFIP WG2.5 WorkingConference on Quality of Numerical Software for insightful comments on afirst draft of this paper..\"the principal developers of Inferno: Sean Dorward, Rob Pike, Dave Presotto, Howard Trickey, and Phil Winterbottom..SHTrademarks.LPInferno, Limbo, and Dis are trademarks of Vita Nuova Holdings Limited.Unix is a trademark of Unix Systems Laboratories.Windows95 is a trademark of Microsoft..EQdelim off.EN.SHReferences.nr PS -1.nr VS -1.LP.nr [N 0 1.de ]N.IP \\n+([N..if \nZ>0 .tm \\$1 \\n([N......Lines 1–500 of 611. Load more
.FP palatino .de ]C \&[\\$1]\\$2 .. .if \nZ=0 .so real.ref .EQ delim $$ .EN .TL Real Inferno .AU .I "Eric Grosse" .AI .I "Lucent Technologies, Bell Labs" .I "Murray Hill NJ 07974 USA" .I "ehg@bell-labs.com" .\"date{19 Aug 1996, minor revisions 7 Jan 1998} .FS Previously appeared in R.F. Boisvert (editor), .ft I The Quality of Numerical Software: Assessment and Enhancement: Proceedings of the IFIP TC2/WG2.5 Working Conference on the Quality of Numerical Software, Oxford, United Kingdom, 8-12 July 1996, .ft R Chapman Hall, London, 1997 (pp. 270-279). .FE .AB Inferno is an operating system well suited to applications that need to be portable, graphical, and networked. This paper describes the fundamental floating point facilities of the system, including: tight rules on expression evaluation, binary/decimal conversion, exceptions and rounding, and the elementary function library. .AE .PP Although the focus of Inferno is interactive media, its portability across hardware and operating platforms, its relative simplicity, and its strength in distributed computing make it attractive for advanced scientific computing as well. Since the appearance of a new operating system is a relatively uncommon event, this is a special opportunity for numerical analysts to voice their opinion about what fundamental facilities they need. The purpose of this short paper is to describe numerical aspects of the initial release of Inferno, and to invite comment before the tyranny of backward compatibility makes changes impossible. .PP An overview of Inferno is given by Dorward et. al. .]C "Inferno" , but for our immediate purposes it may suffice to say that Inferno plays the role of a traditional operating system (with compilers, process control, networking, graphics, and so on) but can run either on bare hardware or on top of another operating system like Windows95 or Unix. Programs for .I "Inferno" are written in the language .I "Limbo" and compiled to machine-independent object files for the .I "Dis" virtual machine, which is then implemented with runtime compilation for best performance. Files are accessible over networks using the .I "Styx" protocol; together with the presentation of most system resources as files and the manipulation of file namespaces, this permits integration of a collection of machines into a team. Limbo looks somewhat like a mixture of C and Pascal, augmented by modules (to cope with the namespace and dynamic loading needs of large programs) and by a channel facility for convenient (coarse-grain) parallel programing. Array references are bounds-checked and memory is garbage collected. .PP The rest of this paper covers the fundamental floating point environment provided by the Limbo compiler and .I "math" module, the ``elementary functions,'' and finally some comments on why particular definitions were chosen or why certain facilities were included or excluded. This discussion assumes the reader is familiar with scientific computing in general and the IEEE floating point standard in particular. .NH 1 Floating point .PP In Limbo, arithmetic on literal and named constants is evaluated at compile time with all exceptions ignored. Arithmetic on variables is left by the compiler to runtime, even if data path analysis shows the value to be a compile time constant. This implies that tools generating Limbo source must do their own simplification, and not expect the compiler to change $x/x$ into $1$, or $-(y-x)$ into $x-y$, or even $x-0$ into $x$. Negation $-x$ changes the sign of $x$; note that this not the same as $0-x$ if $x=0$. .PP The compiler may perform subexpression elimination and other forms of code motion, but not across calls to the mode and status functions. It respects parentheses. The evaluation order of $a+b+c$ follows the parse tree and is therefore the same as for $(a+b)+c$. These are the same rules as for Fortran and C. .PP Contracted multiply-add instructions (with a single rounding) are not generated by the compiler, though they may be used in the native .SM BLAS libraries. All arithmetic follows the IEEE floating point standard .]C "IEEEfp" , except that denormalized numbers may not be supported; see the discussion in section 3. .PP The most important numerical development at the language level recently has been accurate binary/decimal conversion .]C "Clinger" .]C "Gay" .]C "SteeleWhite" . Thus printing a real using .CW "\%g" and reading it on a different machine guarantees recovering identical bits. (Limbo uses the familiar .I "printf" syntax of C, but checks argument types against the format string at compile time, in keeping with its attempt to help the programmer by stringent type checking.) A good .I "strtod/dtoa" is, unfortunately, 1700 lines of source (15kbytes compiled), though with modest average runtime penalty. This code must be used in the compiler so that coefficients are accurately transferred to bytecodes. Smaller, faster, but sloppier, runtimes will also be provided when mandated by limited memory and specialized use. However, programmers may assume the features described in this paper are present in all Inferno systems intended for general computing. .PP Each thread has a floating point control word (governing rounding mode and whether a particular floating point exception causes a trap) and a floating point status word (storing accumulated exception flags). Functions .I "FPcontrol" and .I "FPstatus" copy bits to the control or status word, in positions specified by a mask, returning previous values of the bits. .I "getFPcontrol" and .I "getFPstatus" return the words unchanged. .PP The constants .I "INVAL, ZDIV, OVFL, UNFL, INEX" are non-overlapping single-bit masks used to compose arguments or return values. They stand for the five IEEE exceptions: .IP • ``invalid operation'' ($0/0$,$0 * infinity $,$ infinity - infinity $,$sqrt{-1}$) .IP • ``division by zero'' ($1/0$), .IP • ``overflow'' ($1.8e308$) .IP • ``underflow'' ($1.1e-308$) .IP • ``inexact'' ($.3*.3$). .PP The constants .I "RND_NR, RND_NINF, RND_PINF, RND_Z" are distinct bit patterns for ``round to nearest even'', ``round toward $-{infinity} $'', ``round toward $+{infinity} $'', ``round toward $0$'', any of which can be set or extracted from the floating point control word using .I "RND_MASK" . For example, .IP • to arrange for the program to tolerate underflow, .I "FPcontrol(0,UNFL)." .IP • to check and clear the inexact flag, .I "FPstatus(0,INEX)." .IP • to set directed rounding, .I "FPcontrol(RND_PINF,RND_MASK)." .PP By default, .I "INEX" is quiet and .I "OVFL, UNFL, ZDIV," and .I "INVAL" are fatal. By default, rounding is to nearest even, and library functions are entitled to assume this. Functions that wish to use quiet overflow, underflow, or zero-divide should either set and restore the control register themselves or clearly document that the caller must do so. The ``default'' mentioned here is what a Limbo program gets if started in a fresh environment. Threads inherit floating point control and status from their parent at the time of spawning and therefore one can spawn a ``round toward 0'' shell and re-run a program to effortlessly look for rounding instabilities in a program. .NH 1 Elementary functions .PP The constants .I "Infinity, NaN, MachEps, Pi, Degree" are defined. Since Inferno has thorough support of Unicode, it was tempting to name these $infinity $, $ε $, $π $, and °, but people (or rather, their favorite text editing tools) may not be ready yet for non-\s-2ASCII\s0 source text. .I "Infinity" and .I "NaN" are the positive infinity and quiet not-a-number of the IEEE standard, double precision. .I MachEps is $2 sup {-52}$, the unit in the last place of the mantissa $1.0$. The value of .I "Pi" is the nearest machine number to the mathematical value $π $. .I "Degree" is $"Pi" / 180$. .PP Three useful functions .I "fdim, fmax, fmin" are adopted from the Numerical C extensions .]C "NumerC" . The unusual one of these, often denoted $(x-y) sub {+}$, is defined by $roman "fdim" ( x , y )=x-y$ if $x > y$, else $0$. The compiler may turn these into efficient machine instruction sequences, possibly even branch-free, rather than function calls. There are two almost redundant mod functions: .I "remainder(x,y)" is as defined by the IEEE standard (with result $roman "|" r roman "|" <= y/2$); .I "fmod(x,y)" is $x roman "mod" y$, computed in exact arithmetic with $0<= r<y$. Limbo has a ``tuple'' type, which is the natural return value in the call $(i,f)= roman "modf" ( x )$ to break $x$ into integer and fractional parts. The function .I "rint" rounds to an integer, following the rounding mode specified in the floating point control word. .PP For a good-quality, freely-available elementary function library, .I "math" uses the IEEE subset of .I "fdlibm" .]C "fdlibm" . Of course, a conforming implementation may use entirely different source, but must take care with accuracy and with special arguments. There are the customary power, trigonometric, Bessel, and erf functions, and specialized versions $roman "expm1"( x )=e sup x - 1$, $roman "log1p" ( x )=log ( 1 + x )$. An additional function $roman "pow10"( n ) = 10 sup n$ is defined; in the default implementation this is just fdlibm's $roman "pow" ( 10. , n )$ but it is provided so that separate trade-offs of accuracy and simplicity can be made .]C "MacIlroy" . .I "fdlibm" uses extra precise argument reduction, so the computed $sin (n*Pi)$ is small but nonzero. If demands warrant, degree versions of the trigonometric functions will be added, but for now the style $sin (45*Degree)$ is used. The library also provides IEEE functions .I "ilogb, scalbn, copysign, finite, isnan," and .I "nextafter" . .PP The functions .I "dot, norm1, norm2, iamax, gemm" are adopted from the .SM BLAS .]C "blas" to get tuned linear algebra kernels for each architecture, possibly using extra-precise accumulators. These are defined by $sum {{x sub i}{y sub i}}$, $sum roman | {x sub i} roman | $, $ sqrt{sum { x sub {i sup 2}}} $, $i$ such that $roman | {x sub i} roman | = roman max $, and $C= alpha AB + beta C$ with optional transposes on $A$ and $B$. Since Limbo has only one floating-point type, there is no need here for a precision prefix. Limbo array slices permit the calling sequences to be more readable than in Fortran77 or C, though restricted to unit stride. This encourages better cache performance anyway. The matrix multiply function .I "gemm" remains general stride (and is the foundation for other operations .]C "Kagstrom" ). .PP Limbo is like C in providing singly-subscripted arrays with indexing starting at 0. Although Limbo offers arrays of arrays, as in C, for scientific work a better choice is to adopt the style of linearizing subscripts using Fortran storage order. This promotes easier exchange of data with other applications and reuses effort in organizing loops to achieve good locality. In previous language work .]C "pine" , we implemented a C preprocessor that allowed the programmer to choose a convenient origin (such as 1) and have it compiled into 0 for the base language; because we passed arrays as dope vectors, we were even able to allow different origins for the same array in calling and called functions. The main lesson we learned from that experience, however, was that permutations become a nightmare when there is anything but dogmatic adherence to a single origin. So for an $m$ by $n$ matrix $A$, the programmer should use loops with $0<= i<m$ and $0<= j<n$ and access $A[i+m*j]$. .PP For interoperability with foreign file formats and for saving main memory in selected applications, functions are provided for copying bits between and reals and 32-bit or 64-bit IEEE-format values. .PP Finally, .I "math" provides a tuned quicksort function .I "sort(x,p)" where .I "x" is a real array and .I "p" is an int array representing a 0-origin permutation. This function leaves the contents of .I "x" untouched and rearranges .I "p" so that $x[{p sub i}]<= x[p sub {i+1}]$. This is usually what one wants to do: sort an array of abstract data types based on some key, but without the need to actually swap large chunks of memory. .NH 1 Rationale .PP This section discusses why certain numerical features were included or not. .NH 2 Rounding modes and accumulated exceptions .PP Directed rounding is only needed in a very few places in scientific computing, but in those places it is indispensable. Accumulated floating point exceptions are even more useful. User trap handling is a harder problem, and may be worth leaving for later, possibly with a default ``retrospective diagnostics'' log .]C "Kahan" . .PP Note that the exception masks must be architecture independent, since they reside in the Limbo bytecodes, and therefore the implementation involves a small amount of bit fiddling. Still, it is efficient enough to encourage use. It would be difficult to port to a processor that only had static rounding modes in instruction opcodes rather than the dynamic model specified in section 2 of the IEEE standard. Fortunately, the Alpha does provide both models. .NH 2 Sudden underflow .PP Some processor vendors make supporting gradual underflow just too hard. (One must struggle upon the system trap to reconstruct exactly which instruction was executing and what the state of the registers was. On the MIPS, it is said to be 30 pages of assembler.) So Inferno supports denormalized numbers only if the hardware makes this easy. Providing underflow that is correct but very slow, as some systems do, is not necessarily doing the user a favor. .PP To determine portably if a particular system offers gradual underflow, mask off UNFL and do trial arithmetic. .NH 2 Speed .PP Computers with slow (software) gradual underflow usually provide a fast flush-to-0 alternative. This often suffices, though there are important examples where it forces an uglier and slower coding style. A worse situation is if the hardware uses system traps for Infinity and NaN arithmetic. The resulting slowdown will make otherwise excellent and natural algorithms run slowly .]C "Demmel" . Sadly, even some x86 implementations that do non-finite arithmetic in hardware, do it relatively slowly. .PP We considered providing syntax to declare a certain program scope within which precise IEEE behavior was required, and relaxing the rules outside such scopes. (The numerical C extensions .]C "NumerC" use pragma for this purpose.) These scope declarations would need to be in the bytecodes, since significant optimization may be attempted by the runtime compiler. After some discussion, and with some trepidation, it was agreed that instead all compilers would be required to preserve the same result and status as for an unoptimized version. .NH 2 Comparison .PP The standard C operators .CW == .CW != .CW "<" .CW "<=" .CW ">" .CW ">=" are the only comparisons provided, and they behave exactly like the ``math'' part of Table 4 of the IEEE standard. Programs interested in handling NaN data should test explicitly. This seems to be the way most people program and leads to code more understandable to nonexperts. It is true that with more operators one can correctly write code that propagates NaNs to a successful conclusion\-but that support has been left for later. NaN(''tag'') should be added at that same time. .NH 2 Precision .PP All implementations run exclusively in IEEE double precision. If the hardware has extra-precise accumulators, the round-to-double mode is set automatically and not changeable, in keeping with Limbo's design to have only one floating point type. Extended precision hardware, if available, may be used by the built-in elementary function and .SM BLAS libraries. Also, we contemplate adding a dotsharp function that would use a very long accumulator for very precise inner products, independent of the order of vector elements .]C "kulisch" . But reference implementations that use only double precision, avoid contracted multiply-add, and evaluate in the order 1 up to n will always be available for strict portability. .PP At the time the decision was made to restrict the system to 64-bit floating point, Limbo integers were almost exclusively 32-bit and the consistency argument to have a single real type was compelling. Now that Limbo has more integer types the decision might be reconsidered. But so many engineers needlessly struggle with programs run in short precision, that offering it may do as much harm as good. On most modern computers used for general purpose scientific computing, 64-bit floating point arithmetic is as fast as 32-bit, except for the memory traffic. In cases where the shorter precision would suffice and memory is a crucial concern, the programmer should consider carefully scaled fixed point or specialized compression. To efficiently interoperate with data files that use the short format, programmers may use the provided realbits32 function. While there are surely appropriate uses for a first-class 32-bit real type, for now we follow Kahan's sarcastic motto ``why use lead when gold will do?'' .NH 2 BLAS .PP The few .SM BLAS in the core library were chosen for readability and, in case of gemm, for optimization beyond what a reasonable compiler would attempt. We expect that compilers will (soon) be good enough that the difference between compiling $y+=a*x$ and calling daxpy is small. Also, as mentioned above, dot and gemm might reasonably use combined multiply-add or a long accumulator in some optional implementations. .NH 2 $GAMMA ( x )$ .PP To avoid confusion with the C math library, which defined .I "gamma" as $ln GAMMA $, we offer only .I "lgamma" for now. This function and .I "modf" return an (int,real) tuple rather than assigning through an integer pointer, in keeping with Limbo's design. The opportunity has been taken to drop some obsolete functions like .I "frexp" . Other functions are unchanged from the C math library. .NH 2 Future .PP A prototype preprocessor has been written to allow the scientific programmer to write $A[i,j]$ for an $A$ that was created as a $Matrix(m,n)$ and to have the subscript linearization done automatically. Here $Matrix$ is an Limbo abstract data type containing a real array and integers $m$, $n$, and column stride $lda$ used as in typical Fortran calling sequences. .PP The Limbo compiler is soon expected to implement the type .I "complex" . .PP Higher level numerical libraries will also be provided, and although that topic is beyond the scope of this paper, opinions about what should come first would be welcome. .PP Distributed computing has not been mentioned here because it involves relatively few considerations specific to floating point computation. However, it may be worth noting that in the default environment (with underflow trapped, so that presence or absence of denormalized numbers is not significant) programs run independently on heterogeneous machines nevertheless get precisely identical results, even with respect to thread scheduling. This implies that certain communication steps can be avoided, and that regression testing is considerably simplified. .PP Please direct comments on these numerical aspects of Inferno to Eric Grosse. More general technical comments can be directed to Vita Nuova .CW comments@vitanuova.com ). ( I am grateful to Joe Darcy, Berkeley, to David Gay, Bell Labs, to David Hook, University of Melbourne, and to participants of the IFIP WG2.5 Working Conference on Quality of Numerical Software for insightful comments on a first draft of this paper. .\"the principal developers of Inferno: Sean Dorward, Rob Pike, Dave Presotto, Howard Trickey, and Phil Winterbottom. .SH Trademarks .LP Inferno, Limbo, and Dis are trademarks of Vita Nuova Holdings Limited. Unix is a trademark of Unix Systems Laboratories. Windows95 is a trademark of Microsoft. .EQ delim off .EN .SH References .nr PS -1 .nr VS -1 .LP .nr [N 0 1 .de ]N .IP \\n+([N. .if \nZ>0 .tm \\$1 \\n([N .. ....