1

I'll assume this is not considered off-topic since simd is in the swift repo and about the only mention I could find of simd_orient, except its very sparse documentation, was @scanon introducing it at WWDC 2016.

If there are more thorough documentation of simd_orient (and other barely documented simd functions), I'd appreciate pointers to such.

For example, I assume that the orientation() function below is making proper use of it, but I'd like to know if using it for calculating the area and testing if the triangle contains a point is recommended or not:

struct Triangle2D {
    let a, b, c: float2
    
    init(_ a: float2, _ b: float2, _ c: float2) {
        (self.a, self.b, self.c) = (a, b, c)
    }
    
    /// Returns a positive value if the triangle is positively oriented, a
    /// negative value if it is negatively oriented and zero if it is
    /// degenerate (the three points are colinear).
    func orientation() -> Float {
        return simd_orient(a, b, c)
    }

    /// Returns the area of the triangle. The area is positive if the triangle
    /// is positively oriented, negative if it's negatively oriented and zero
    /// if it's degenerate.
    func area() -> Float {
        return orientation() / 2
    }
    
    /// Returns true if the point is within or on the edge of the triangle.
    /// The triangle has to be positively oriented.
    func contains(point p: float2) -> Bool {
        if simd_orient(a, b, p) < 0 { return false }
        if simd_orient(b, c, p) < 0 { return false }
        if simd_orient(c, a, p) < 0 { return false }
        return true
    }
    
    /* ... */
}

Also, I'm not entirely sure how the version of simd_orient that takes two rather than three float2s works. For example, is it always exactly like the three-argument one with the last argument being float2(0, 0):

$ swift
Welcome to Apple Swift version 4.2 (swiftlang-1000.0.36 clang-1000.0.4). Type :help for assistance.
  1> import simd
  2> simd_orient(float2(0.1, 2.3), float2(4.5, 6.7)) 
$R0: Float = -9.67999935
  3> simd_orient(float2(0, 0), float2(0.1, 2.3), float2(4.5, 6.7))
$R1: Float = -9.6800003
  4> simd_orient(float2(0.1, 2.3), float2(4.5, 6.7), float2(0, 0))
$R2: Float = -9.67999935

?

Is there any where I can read up on these things?

2

Probably off-topic here, but I'll point out that there's very extensive documentation of these operations in the header where they're defined ($SDKROOT/usr/include/simd/geometry.h).

Your area function will not work. You should not attempt to interpret the result of simd_orient beyond the sign. It happens to give you the area in simple cases, but in hard-to-round cases it will not. Your other functions are correct as written.

Because the sign is all that matters, -9.67999935 and -9.6800003 are the same result as far as the semantics of the simd_orient function are concerned. We considered having it return (-1, 0 +1) to avoid this confusion, but that would involve doing some extra work for no real benefit.

2 Likes
3

What would you suggest as a working way to compute the area?
I mean would this:

    func area() -> Float {
        return simd_length(simd_cross(b - a, c - a)) / 2
    }

or perhaps rather this:

    func area() -> Float {
        return ((b.x - a.x) * (c.y - a.y) - (c.x - a.x) * (b.y - a.y)) / 2
    }

be any better?

Thank's a lot for the pointer to the header file (whose contents are great)!

4

Miscalculating Area and Angles of a Needle-like Triangle gives a numerically stable algorithm for the area of a triangle in §2:

First sort a, b, c so that a ≥ b ≥ c ; this can be done at the cost of at most three comparisons. If c-(a-b) < 0 then the data are not side-lengths of a real triangle; otherwise compute its area

∆(a, b, c) : = (1/4)√( (a+(b+c)) (c-(a-b)) (c+(a-b)) (a+(b-c)) ) .

Do not remove parentheses from this formula! It cannot give rise to √(< 0) .

1 Like
5

But in that algorithm a, b, c are the side lengths, and assumes they are given and exact, and it is compared to Heron's formula. I guess I could calculate the side lengths using simd_precise_distance, but I'm not sure how this method would compare to the one I gave above (where a b c are vertices), repeated here:

    func area() -> Float {
        return ((b.x - a.x) * (c.y - a.y) - (c.x - a.x) * (b.y - a.y)) / 2
    }

EDIT: Here is a program with an example triangle for which the two methods show quite different areas:

import simd

struct Triangle2D {
    let a, b, c: float2
    
    init(_ a: float2, _ b: float2, _ c: float2) {
        (self.a, self.b, self.c) = (a, b, c)
    }
    
    /// Returns a positive value if the triangle is positively oriented, a
    /// negative value if it is negatively oriented and zero if it is
    /// degenerate (the three points are colinear).
    func orientation() -> Float {
        return simd_orient(a, b, c)
    }

    /// Returns the area of the triangle. The area is positive if the triangle
    /// is positively oriented, negative if it's negatively oriented and zero
    /// if it's degenerate.
    func areaFast() -> Float {
        // See: https://forums.swift.org/t/15725/2
        //return orientation() / 2 // <-- Not correct according to above link.
        //return simd_length(simd_cross(b - a, c - a)) / 2 // <- not signed
        return ((b.x - a.x) * (c.y - a.y) - (c.x - a.x) * (b.y - a.y)) / 2
    }
    
    /// Returns the area of the triangle, no matter the orientation, or -1 if
    /// the computed side length of the triangle cannot be a real triangle.
    func areaPrecise() -> Float {
        // https://people.eecs.berkeley.edu/~wkahan/Triangle.pdf
        // First sort sa, sb, sc such that sa >= sb >= sc:
        var sa = simd_precise_distance(a, b)
        var sb = simd_precise_distance(b, c)
        var sc = simd_precise_distance(c, a)
        if (sa < sc) { swap(&sa, &sc) }
        if (sa < sb) {
            swap(&sa, &sb)
        } else if (sb < sc) {
            swap(&sb, &sc)
        }
        // If sc-(sa-sb) < 0, the data are not side-lengths of a real triangle:
        if sc - (sa - sb) < 0 { return -1 }
        // Area is:
        // (1/4)√( (sa+(sb+sc))*(sc-(sa-sb))*(sc+(sa-sb))*(sa+(sb-sc)) )
        return sqrt( (sa+(sb+sc))*(sc-(sa-sb))*(sc+(sa-sb))*(sa+(sb-sc)) ) / 4
    }
}

let triangle = Triangle2D(float2(0,0), float2(15_000, 0), float2(5_000.001, 0.001))
print("areaFast:   ", triangle.areaFast())
print("areaPrecise:", triangle.areaPrecise())

It will print:

areaFast:    7.5000005
areaPrecise: 0.0

I tried to ask Wolfram Alpha for the true area but for some reason it doesn't think those vertices form a possible triangle, but AFAICS it's no less possible than eg this:

let triangle = Triangle2D(float2(0, 0), float2(15_000, 0), float2(5_000, 10))
print("areaFast:   ", triangle.areaFast())
print("areaPrecise:", triangle.areaPrecise())

for which Wolfram Alpha agrees with areaFast:

areaFast:    75000.0
areaPrecise: 74115.914

And here is another example for which areaFast is clearly right and areaPrecise is wrong:

let triangle = Triangle2D(float2(1, 0), float2(0, 1000000000), float2(0, -1000000000))
print("areaFast:   ", triangle.areaFast())    // 1e+09
print("areaPrecise:", triangle.areaPrecise()) // 0.0

(Again, sorry if this is considered off-topic, I'll trust a moderator to intervene if that's the case.)

6

I think your areaFast is correct, as it is just half the cross? product of (the vector from a to b) and (the vector from a to c), which is itself a way to calculate signed area.

7

Well, from a strictly mathematical standpoint, the area of that triangle is exactly 7.5, since it has a base of 15,000 and a height of 0.001.

I haven’t looked into numerically stable algorithms for calculating the area of a triangle specified by floating-point coordinates, so I don’t know if there’s a good way to reliably produce the correct result, but that result ought to be 7.5 in this example.

8

Actually, the result should be 7.50000048.

How I got that number is:

  • Get all relevant number's debugDescriptions,
    -- as those are the values the resulting floats are actually being interpreted as.
  • Put those into Wolfram Alpha, and extract it's result.
    -- namely "7.500000375"
  • convert that into a float via Float("7.500000375")
  • Get it's debugDescription
    -- which is 7.50000048
1 Like
9

How exactly did you get eg Float(0.001) to print as 0.00100000005? The only way I can see it is in the resulting values reported by the REPL:

$ swift
Welcome to Apple Swift version 4.2 (swiftlang-1000.0.36 clang-1000.0.4). Type :help for assistance.
  1> let x: Float = 0.001
x: Float = 0.00100000005
  2> print(x.debugDescription)
0.001
  3> dump(x)
- 0.001
$R0: Float = 0.00100000005

But note that its debugDescription (and dump) is still "0.001".

And this command line app:

let x: Float = 0.001
print(x)
print(x.debugDescription)
dump(x)
print(String(describing: x))

will print:

0.001
0.001
- 0.001
0.001

EDIT: Aha, turns out that it is a matter of compiler versions, the above is the default toolchain of Xcode 10 beta 6, that is:

$ swiftc --version
Apple Swift version 4.2 (swiftlang-1000.0.36 clang-1000.10.44)
Target: x86_64-apple-darwin17.7.0

But if I instead compile the above program with:

$ swiftc --version
Apple Swift version 4.1.2 (swiftlang-902.0.54 clang-902.0.39.2)
Target: x86_64-apple-darwin17.7.0

It prints:

0.001
0.00100000005
- 0.00100000005
0.001

So, Swift 4.2 changed the way it prints float values, and I now remember noticing a similar change before and mentioning in another thread, the answer being:

The REPL seems to be using the old algorithm still though, printing a non-minimal-length decimal value instead of the minimal-length decimal value that parses back into the original floating-point value.

10

The REPL uses lldb’s data formatters, which probably use the libc float printer. It’s worth filing a bug against lldb to use the better algorithm as well.

2 Likes
11

I might be missing something obvious here now but: How (in Swift) do I get the exact decimal representation for a specific (let's say single precision) floating point value?

Ie, 0.001 is just the minimal-length decimal value needed to parse back into the same Float value, but I want the exact decimal representation (unrelated to Float) for the value that happens to be the closest value to 0.001 representable as a Float (that value is not 0.001).

12

I don’t think there’s an off the shelf mechanism for that, since it’s not particularly useful. You can use the %a format specifier with printf-derived APIs to print the exact value using hex float notation. Otherwise, you could extract the exponent and significand fields yourself and compute significand * 2**(exponent-bias) yourself.

1 Like
13

I'm not sure why Decimal doesn't have an initializer that takes a Float, but since every Float value is exactly representable as a Double, then perhaps this is correct:

Decimal(Double((Float(0.001)))) // 0.0010000000474974513152

If so, then the reasoning in @Dante-Broggi's above post (which afaics is mistaking Swift's (old, inaccurate) minimal-decimal-values-that-parses-back-to-the-same-float-values for the exact decimal representations of the float values) could be properly done like this instead:

    let correctFloatAreaAsDecimal =
        (Decimal(15_000) * Decimal(Double(Float(0.001)))) / 2.0
    print(correctFloatAreaAsDecimal)// 7.500000356230884864
    print(Float("7.500000356230884864")!) // 7.5000005
    // And of course, 7.5000005 is just the minimal-length decimal value that
    // parses back to the same closest value representable as a Float, a value
    // whose exact decimal representation is (supposedly):
    print(Decimal(Double(Float(7.5000005)))) // 7.500000476837158912

Which means that areaFast() did actually produce the best result it could (at least for this particular triangle) given that 0.001 isn't exactly representable as a Float value.

But according to this online converter, the exact decimal representation for the value Float(0.001) is:
0.001000000047497451305389404296875
and not, as above:
0.0010000000474974513152

14 
1 > import Foundation
2 > String(format: "%.64f", 0.001 as Float)
String = "0.0010000000474974513053894042968750000000000000000000000000000000"
1 Like
15

Yes, and just as a further note on the problems that are still there when using Decimal as i did above:
Having the exact (decimal) representation, like:
0.001000000047497451305389404296875 for Float(0.001),
is still not enough, since we'd need an arbitrary-precision (decimal) number type to do exact calculations on them (Decimal has 38 digit mantissa).

But anyway, I guess the take away is that these numbers:
7.5000003562308847904205322265625
7.500000356230884864
7.500000375
7.50000048
7.5000005
are all mapped to one and the same Float value, a value which when exactly represented in decimal is:
7.500000476837158203125
and this is why the new float printing algorithm in Swift 4.2 now prints it as
7.5000005 rather than
7.50000048

And among the methods mentioned here to compute the area of a triangle, areaFast() seems to be most correct / least problematic.

1 Like