How to implement rounding in an all-purpose stack language using different types?

Question

Disclaimer: If you are not terribly interested in numerics and mathematical processes, this is most likely nothing for you.

I am currently a bit stuck in a development process of a private project I followed for a long time and I think I need some input because I cannot decide how to do it best. Unfortunately I must explain the problem a bit.

EDIT: This is NOT about rounding implementation. This is already done, both for float/doubles and arbitrary precision. You can assume I have enough knowledge about numerics, rounding modes and floating-point problems, the problem is how to design rounding operators in a stack-based program.

I am interested in the reproducibility problem in numerical computations. Due to the immense power (billions of operations per second) and certain problems by the languages used for number crunching (C and FORTRAN), people cannot check what the machine is actually doing in contrast what the machine is supposedly doing (Deteoriation of language standards like indifference to qNaNs and sNaNs, allow to "optimize" a+(b+c) == (a+b)+c, x86 vs MME floating-point, silent FMA support etc. etc.). This view is shared by William Kahan, the designer of the x86 floating-point architecture.

So I worked on a stack-based language like FORTH which allows reproducible results. You have types: aggregates (Vector,Matrix,Tensor) and numbers (complex, reals). You put two space instances on the stack, execute an ADD, and you get the resulting sum. If the spaces were two floats, you get a float. If the spaces were two doubles, you get a double. If the two spaces were arbitrary precision numbers, you get an arbitrary precision number. The same with vectors and matrices.

My aching problem is rounding. If we exclude exact rationals, after division every higher operation (algebraic and transcendental operations) requires rounding. There are two rounding operations which are already implemented: place rounding, setting the exact position of the rounding (rounding down: 2.345 using 0 => 2 / 2.345 using -1 => 2.3) and significant digits rounding by setting the length of the result (rounding down: 2.345 using 1 = 2 / 2.345 using 2 = 2.3).

The approaches:

• One ring to rule them all: Only one rounding setting for all available types. Problem: Some operations are exact and should be executed as such; arbitrary precision numbers offer exact addition and exact multiplication. If I define e.g. ADD_EXACT I include a keyword which is not implemented by most datatypes and will be therefore a normal ADD and the danger is that I (or others) forget to use ADD_EXACT when necessary. Another problem is that it is makes sense to use different rounding modes for some operations: place rounding for addition and significant digit rounding for multiplication. I fear to swamp the program with unnecessary switches to set the rounding mode.

• Perhaps some more rings....: Several settable rounding modes. Allows adding and multiplication independently. Problem: Much more possibilities (set and retrieve rounding mode, rounding operation and rounding parameter). How many choices do I offer ? Only two and risking that if e.g. FMA cannot be supported or many more and risking to overwhelm the program with much too much rounding bureaucracy ?

Rounding is a crucial operation and once implemented, it cannot be changed easily. This must be solved as good as possible.

What I need is a simple and elegant model and I fear that I have worked so long on it that I simply do not see it and/or need a plan to find out the right solution. Help is very appreciated and I will update this post to clarify some aspects I may have forgotten.

EDIT: Performance: Is not the design goal. It does not help you to get an incredibly fast wrong answer. It does not help you if two programmers on different machines get a different or (even similar) result if noone knows what is the right answer (or that even both are wrong). The goal of the program is exactly that: Find out when IEEE 754R precision will fail and find indicators and methods to validate numeric results. If you know that double precision and the tool will get you to the solution, by all means use C and FORTRAN.

Answer 1

I couple of notes:

• I don't understand your ADD_EXACT problem. If the operation is exact, don't round. Alternatively, if there has to be a rounding step (this includes most FPU operations), it must not introduce error.
• FMA should probably be a completely separate operation from both addition and multiplication, since its rounding behavior is fundamentally different: It rounds only once. The two rounding steps of separate-multiply-add would have to be prescient to match FMA rounding, regardless of the exact rounding strategy.
• The terms you use to talk about rounding worry me. It sounds like you're unaware of, or departing from, binary floating point and how rounding can actually vary in IEE-754 and in general (round to nearest even, round towards zero, round towards positive infinity, etc.).
• You seem inexplicably focused on the aesthetic but largely inconsequential concern of how to name the digit to which one rounds. No mention of the actual rounding strategy (see above).
• Design-wise, there should be a default rounding mode, with the option to explicitly use others. Whether to do this choice on an operation-by-operation basis, with a lexically scoped with_rounding mode { ... code ... } block, or dynamically scoped (i.e. set a global or thread-local variable that controls all rounding everywhere) is a very hard question. All options have severe drawbacks and I don't feel qualified to have an opinion on this.
• You pitch your design as a potential alternative to FORTRAN, C, etc. for number crunching, but making performance competitive with (occasionally inaccurate) FPUs built in silicon is hard to impossible, and you don't seem worried about it. You will have to deal with binary, fixed-width floating point (i.e., IEEE-754) and the idiosyncrasies of its hardware implementations. (The essentially dynamic typing won't help either.)

Answer 2

Due to the inexplicable power of rubber duck debugging I have now a solution I find better than all the other approaches I thought of.

I will implement a rounding stack. When a subroutine is entered, the last entry will be duplicated and pushed as first element on the stack (It will also be used entry if the stack is completely cleared). The stack is strictly local, access to rounding values of superroutines is not possible. When the subroutine is left, the stack will be cleared; no survivors.

The stack has the following operations: ROUND_DUP (to save the original state if necessary), ROUND_PUSH , ROUND_POP, ROUND_CYCLE (see later) and accessor functions for the current rounding.

ROUND_CYCLE cycles through the rounding modes. If we have 3 rounding modes on the stack, this means rounding mode 1,2,3,1,2,3...If we have only one rounding mode on the stack, the instruction does nothing.

Augmented are test functions: EXACT_ADD_SUPPORTED EXACT_MUL_SUPPORTED (It is not advisable to use only IS_ARBITRARY_PRECISION because it is possible that the system may have an long accumulator for adding (even in hardware) and does not support exact multiplication).

So I can write general functions for both arbitrary precision and fixed precision.

EXACT_ADD_SUPPORTED IF ROUND_PUSH EXACT ENDIF DO .... ROUND_CYCLE ADD ROUND_CYCLE LOOP 

If this is arbitrary precision, the two rounding values are switched and I have exact addition. If this is fixed precision, only one rounding value remains and is used; ROUND_CYCLE does nothing. So I can set the necessary roundings at the beginning of the program, making it easy to find and understand. Given that CYCLE does only switches, I must think before how the rounding modes must be set during the execution, avoiding superfluous and unstructured approaches. I have only one current rounding mode, avoiding overhead.

If you have improvements and criticism, be invited to share.