Safe addition: Difference between revisions

Content added Content deleted
(Added Hare)
m (syntax highlighting fixup automation)
Line 28:
=={{header|Ada}}==
An interval type based on Float:
<langsyntaxhighlight Adalang="ada">type Interval is record
Lower : Float;
Upper : Float;
end record;</langsyntaxhighlight>
Implementation of safe addition:
<langsyntaxhighlight Adalang="ada">function "+" (A, B : Float) return Interval is
Result : constant Float := A + B;
begin
Line 51:
return (Float'Adjacent (0.0, Float'First), Float'Adjacent (0.0, Float'Last));
end if;
end "+";</langsyntaxhighlight>
The implementation uses the attribute T'Machine_Rounds to determine if rounding is performed on inexact results. If the machine rounds a symmetric interval around the result is used. When the machine does not round, it truncates. Truncating is rounding towards zero. In this case the implementation is twice better (in terms of the result interval width), because depending on its sign, the outcome of addition can be taken for either of the bounds. Unfortunately most of modern processors are rounding.
 
Test program:
<langsyntaxhighlight Adalang="ada">with Ada.Text_IO; use Ada.Text_IO;
procedure Test_Interval_Addition is
-- Definitions from above
Line 64:
begin
Put (1.14 + 2000.0);
end Test_Interval_Addition;</langsyntaxhighlight>
Sample output:
<pre> 2.00113989257813E+03, 2.00114013671875E+03</pre>
 
=={{header|AutoHotkey}}==
<langsyntaxhighlight lang="ahk">Msgbox % IntervalAdd(1,2) ; [2.999999,3.000001]
 
SetFormat, FloatFast, 0.20
Line 82:
err:=0.1**(SubStr(A_FormatFloat,3) > 15 ? 15 : SubStr(A_FormatFloat,3))
Return "[" a+b-err ","a+b+err "]"
}</langsyntaxhighlight>
 
=={{header|C}}==
Line 97:
 
=== C99 with fesetround() ===
<langsyntaxhighlight lang="c">#include <fenv.h> /* fegetround(), fesetround() */
#include <stdio.h> /* printf() */
Line 146:
}
return 0;
}</langsyntaxhighlight>
Output:
<pre>
Line 169:
{{works with|MinGW}}
 
<langsyntaxhighlight lang="c">#include <float.h> /* _controlfp() */
#include <stdio.h> /* printf() */
 
Line 217:
}
return 0;
}</langsyntaxhighlight>
 
=== fpsetround() ===
{{works with|OpenBSD|4.8/amd64}}
 
<langsyntaxhighlight lang="c">#include <ieeefp.h> /* fpsetround() */
#include <stdio.h> /* printf() */
 
Line 270:
}
return 0;
}</langsyntaxhighlight>
 
Output from OpenBSD: <pre>1 + 2 =
Line 291:
=={{header|C sharp|C#}}==
{{trans|Java}}
<langsyntaxhighlight lang="csharp">using System;
 
namespace SafeAddition {
Line 330:
}
}
}</langsyntaxhighlight>
{{out}}
<pre>(1.2 + 0.03) is in the range (1.23, 1.23)</pre>
Line 336:
=={{header|C++}}==
{{trans|C#}}
<langsyntaxhighlight lang="cpp">#include <iostream>
#include <tuple>
 
Line 378:
 
return 0;
}</langsyntaxhighlight>
{{out}}
<pre>(1.200000 + 0.030000) is in the range (1.2299998998641968, 1.2300001382827759)</pre>
Line 384:
=={{header|D}}==
{{trans|Kotlin}}
<langsyntaxhighlight Dlang="d">import std.traits;
auto safeAdd(T)(T a, T b)
if (isFloatingPoint!T) {
Line 399:
auto r = safeAdd(a, b);
writefln("(%s + %s) is in the range %0.16f .. %0.16f", a, b, r.d, r.u);
}</langsyntaxhighlight>
 
{{out}}
Line 416:
E currently inherits [[Java]]'s choices of [[IEEE]] floating point behavior, including round to nearest. Therefore, as in the Ada example, given ignorance of the actual direction of rounding, we must take one step away in ''both'' directions to get a correct interval.
 
<langsyntaxhighlight lang="e">def makeInterval(a :float64, b :float64) {
require(a <= b)
def interval {
Line 432:
}
return interval
}</langsyntaxhighlight>
 
<langsyntaxhighlight lang="e">? makeInterval(1.14, 1.14) + makeInterval(2000.0, 2000.0)
# value: [2001.1399999999999, 2001.1400000000003]</langsyntaxhighlight>
 
E provides only 64-bit "double precision" floats, and always prints them with sufficient decimal places to reproduce the original floating point number exactly.
Line 441:
=={{header|Forth}}==
{{trans|Tcl}}
<langsyntaxhighlight lang="forth">c-library m
s" m" add-lib
\c #include <math.h>
Line 455:
 
: savef+ ( F: r1 r2 -- r3 r4 ) \ r4 <= r1+r2 <= r3
f+ fdup fstepup fswap fstepdown ;</langsyntaxhighlight>
{{output}}
<pre>
Line 465:
=={{header|Go}}==
{{trans|Tcl}}
<langsyntaxhighlight lang="go">package main
 
import (
Line 494:
a, b := 1.2, .03
fmt.Println(a, b, safeAdd(a, b))
}</langsyntaxhighlight>
Output:
<pre>
Line 502:
=={{header|Hare}}==
{{trans|C}}
<langsyntaxhighlight lang="hare">use fmt;
use math;
 
Line 528:
math::setround(orig);
return interval{a = r0, b = r1};
};</langsyntaxhighlight>
{{out}}
<pre>
Line 539:
=={{header|J}}==
J uses 64 bit IEEE floating points, providing 53 binary digits of accuracy.
<langsyntaxhighlight lang="j"> err =. 2^ 53-~ 2 <.@^. | NB. get the size of one-half unit in the last place
safeadd =. + (-,+) +&err
0j15": 1.14 safeadd 2000.0 NB. print with 15 digits after the decimal
2001.139999999999873 2001.140000000000327</langsyntaxhighlight>
 
=={{header|Java}}==
{{trans|Kotlin}}
<langsyntaxhighlight Javalang="java">public class SafeAddition {
private static double stepDown(double d) {
return Math.nextAfter(d, Double.NEGATIVE_INFINITY);
Line 565:
System.out.printf("(%.2f + %.2f) is in the range %.16f..%.16f", a, b, result[0], result[1]);
}
}</langsyntaxhighlight>
{{out}}
<pre>(1.20 + 0.03) is in the range 1.2299999999999998..1.2300000000000002</pre>
Line 571:
=={{header|Julia}}==
Julia has the IntervalArithmetic module, which provides arithmetic with defined precision along with the option of simply computing with actual two-number intervals:
<syntaxhighlight lang="julia">
<lang Julia>
julia> using IntervalArithmetic
 
Line 591:
julia> a + b
[0.399999, 0.900001]
</syntaxhighlight>
</lang>
 
=={{header|Kotlin}}==
{{trans|Tcl}}
<langsyntaxhighlight lang="scala">// version 1.1.2
 
fun stepDown(d: Double) = Math.nextAfter(d, Double.NEGATIVE_INFINITY)
Line 607:
val b = 0.03
println("($a + $b) is in the range ${safeAdd(a, b)}")
}</langsyntaxhighlight>
 
{{out}}
Line 619:
=={{header|Nim}}==
{{trans|C}}
<langsyntaxhighlight lang="nim">import fenv, strutils
 
proc `++`(a, b: float): tuple[lower, upper: float] =
Line 638:
echo x.ff, " + ", y.ff, " ="
echo " [", d.ff, ", ", u.ff, "]"
echo " size ", (u - d).ff, "\n"</langsyntaxhighlight>
Output:
<pre>1.0000000000000000 + 2.0000000000000000 =
Line 658:
=={{header|Perl}}==
There are ways avoid this whole problem (e.g. <code>Math::Decimal</code>), but another module suffices for this task.
<langsyntaxhighlight lang="perl">use strict;
use warnings;
use Data::IEEE754::Tools <nextUp nextDown>;
Line 668:
}
 
printf "%.17f (%.17f, %.17f)\n", safe_add (1/9,1/7);</langsyntaxhighlight>
{{out}}
<pre>0.25396825396825395 (0.25396825396825390, 0.25396825396825401)</pre>
Line 680:
Not surprisingly you have to get a bit down and dirty to manage this sort of stuff in Phix.
 
<!--<langsyntaxhighlight Phixlang="phix">-->
<span style="color: #008080;">include</span> <span style="color: #000000;">builtins</span><span style="color: #0000FF;">\</span><span style="color: #000000;">VM</span><span style="color: #0000FF;">\</span><span style="color: #000000;">pFPU</span><span style="color: #0000FF;">.</span><span style="color: #000000;">e</span> <span style="color: #000080;font-style:italic;">-- :%down53 etc</span>
Line 734:
<span style="color: #7060A8;">printf</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #008000;">" size %.16g\n\n"</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">high</span> <span style="color: #0000FF;">-</span> <span style="color: #000000;">low</span><span style="color: #0000FF;">);</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
<!--</langsyntaxhighlight>-->
{{out}}
<pre>
Line 763:
Python doesn't include a module that returns an interval for safe addition, however, it does [http://docs.python.org/library/math.html#math.fsum include] a routine for performing additions of floating point numbers whilst preserving precision:
 
<langsyntaxhighlight lang="python">>>> sum([.1, .1, .1, .1, .1, .1, .1, .1, .1, .1])
0.9999999999999999
>>> from math import fsum
>>> fsum([.1, .1, .1, .1, .1, .1, .1, .1, .1, .1])
1.0</langsyntaxhighlight>
 
=={{header|Racket}}==
 
<syntaxhighlight lang="racket">
<lang Racket>
#lang racket
 
Line 799:
;; literals
(bf+ (bf "1.14") (bf "2000.0"))
</syntaxhighlight>
</lang>
 
=={{header|Raku}}==
Line 831:
 
 
<syntaxhighlight lang="raku" perl6line>say "Floating points: (Nums)";
say "Error: " ~ (2**-53).Num;
 
Line 858:
say "\nRaku default stringification for 1.5**63:\n" ~ $rat; # standard stringification
say "\nRat::Precise stringification for 1.5**63:\n" ~$rat.precise; # full precision
}</langsyntaxhighlight>
{{out}}
<pre>Floating points: (Nums)
Line 887:
constrained by how much virtual memory is (or can be) allocated.
<br>Eight million digits seems about a practical high end, however.
<langsyntaxhighlight lang="rexx">numeric digits 1000 /*defines precision to be 1,000 decimal digits. */
 
y=digits() /*sets Y to existing number of decimal digits.*/
 
numeric digits y + y%10 /*increase the (numeric) decimal digits by 10%.*/</langsyntaxhighlight><br><br>
 
=={{header|Ruby}}==
Line 898:
When adding <tt>BigDecimal</tt> values, <tt>a + b</tt> is always safe. This example uses <tt>a.add(b, prec)</tt>, which is not safe because it rounds to <tt>prec</tt> digits. This example computes a safe interval by rounding to both floor and ceiling.
 
<langsyntaxhighlight lang="ruby">require 'bigdecimal'
require 'bigdecimal/util' # String#to_d
 
Line 920:
["0.1", "0.00002"],
["0.1", "-0.00002"],
].each { |a, b| puts "#{a} + #{b} = #{safe_add(a, b, 3)}" }</langsyntaxhighlight>
 
Output: <pre>1 + 2 = 0.3E1..0.3E1
Line 928:
 
=={{header|Scala}}==
<langsyntaxhighlight Scalalang="scala">object SafeAddition extends App {
val (a, b) = (1.2, 0.03)
val result = safeAdd(a, b)
Line 940:
println(f"($a%.2f + $b%.2f) is in the range ${result.head}%.16f .. ${result.last}%.16f")
 
}</langsyntaxhighlight>
 
=={{header|Swift}}==
<langsyntaxhighlight lang="swift">let a = 1.2
let b = 0.03
 
print("\(a) + \(b) is in the range \((a + b).nextDown)...\((a + b).nextUp)")</langsyntaxhighlight>
 
{{out}}
Line 955:
<br>
{{libheader|critcl}}
<langsyntaxhighlight lang="tcl">package require critcl
package provide stepaway 1.0
critcl::ccode {
Line 966:
critcl::cproc stepdown {double value} double {
return nextafter(value, -DBL_MAX);
}</langsyntaxhighlight>
With that package it's then trivial to define a "safe addition" that returns an interval as a list (lower bound, upper bound).
<langsyntaxhighlight lang="tcl">package require stepaway
proc safe+ {a b} {
set val [expr {double($a) + $b}]
return [list [stepdown $val] [stepup $val]]
}</langsyntaxhighlight>
 
 
=={{header|Transd}}==
<langsyntaxhighlight lang="scheme">#lang transd
 
MainModule : {
Line 988:
(lout "(+ " a " " b ") is in the range: " prec: 20 (safeAdd a b))
)
}</langsyntaxhighlight>{{out}}
<pre>
(+ 1.2 0.03) is in the range: [1.2299999999999997602, 1.2300000000000002043]
Line 997:
 
However, if Wren is embedded in a suitable C program, then we can ask the latter to call it for us and pass back the result.
<langsyntaxhighlight lang="ecmascript">/* safe_addition.wren */
class Interval {
construct new(lower, upper) {
Line 1,021:
var a = 1.2
var b = 0.03
System.print("(%(a) + %(b)) is in the range %(Interval.safeAdd(a,b))")</langsyntaxhighlight>
 
which we embed in the following C program and run it (using GCC 9.3.0).
 
Note that Wren's built-in print statement never displays numbers with more than 14 digit accuracy. However, if the interval bounds had been displayed directly from C with 17 digit accuracy (the maximum for the 'double' type), there would have been a marginal difference between them, namely: [1.2299999999999998, 1.2300000000000002].
<langsyntaxhighlight Clang="c">#include <stdlib.h>
#include <stdio.h>
#include <math.h>
Line 1,108:
free(script);
return 0;
}</langsyntaxhighlight>
 
{{out}}