Tom Switzer
January 30, 2015
250

# Speed, Correctness, or Simplicity: Choose 3

This talk introduces the floating point filter implementation in Spire (spire.math.FpFilter).

January 30, 2015

## Transcript

1. Speed,  Correctness,  or  Simplicity:
Choose  3
Tom  Switzer
@9xxit

h;ps://github.com/9xxit/fpﬁlter-­‐talk

2. Overview

Floa9ng  point  is  “good  enough”…

most  of  the  9me.

3. Op9ons
Use  Double,  live  with  the  errors.

Use  higher  precision  type,
live  with  performance  loss.

But,  there  is  a  3rd  op9on…

4. Floa9ng  Point  Filters

Use  ﬂoa9ng  point  when  you  can.

Use  higher  precision  when  you  can’t.

5. Err…  Not  So  Simple
Solve  problem  using  ﬂoa9ng  point
approxima9on…

Maintain  an  error  bound  on  approxima9on.

Re-­‐evaluate  with  exact  type  if  error  too  large.

6. The  Catch

7.

What  is  the  determinant  of  my  matrix?

8. Not  Good  For:  Minimizing  Errors
in  Floa9ng  Point  Arithme9c

9.

What  is  the  sign  of  the  determinant  of  my
matrix?

10. Good  For:  Making  a  Decision

11. FpFilter[A]
Simple  wrapper:  FpFilter[Rational]
Standard  Opera2ons
+,  -­‐,  *,  /,  .sqrt,  etc

Fast  predictes
signum,  compare,  isWhole,  etc.

12. FpFilter[A]
class FpFilter[A](
apx: Double,
mes: Double,
ind: Int,
exact: => A
) {

}
ﬂoa9ng  point  approxima9on
error  bounds

13. FpFilter[A]
class FpFilter[A](
apx: Double,
mes: Double,
ind: Int,
exact: => A
) {

}
error  bounds
“Exact  Geometric  Computa2on
Burnikel,  Funke  &  Seel.

14. FpFilter[A]
class FpFilter[A](
apx: Double,
mes: Double,
ind: Int,
exact: => A
) {

}
error  bounds
thunk  for  higher  precision
Welcome  to  …

15. …  Macro  City
def abs(implicit ev: Signed[A]): FpFilter[A] = macro FpFilter.absImpl[A]
def unary_- (implicit ev: Rng[A]) : FpFilter[A] = macro FpFilter.negateImpl[A]
def +(rhs: FpFilter[A])(implicit ev: Semiring[A]): FpFilter[A] = macro FpFilter.plusImpl[A]
def -(rhs: FpFilter[A])(implicit ev: Rng[A]): FpFilter[A] = macro FpFilter.minusImpl[A]
def *(rhs: FpFilter[A])(implicit ev: Semiring[A]): FpFilter[A] = macro FpFilter.timesImpl[A]
def /(rhs: FpFilter[A])(implicit ev: Field[A]): FpFilter[A] = macro FpFilter.divideImpl[A]
def sqrt(implicit ev: NRoot[A]): FpFilter[A] = macro FpFilter.sqrtImpl[A]
def <(rhs: FpFilter[A])(implicit ev0: Signed[A], ev1: Rng[A]): Boolean = macro FpFilter.ltImpl[A]
def >(rhs: FpFilter[A])(implicit ev0: Signed[A], ev1: Rng[A]): Boolean = macro FpFilter.gtImpl[A]
def <=(rhs: FpFilter[A])(implicit ev0: Signed[A], ev1: Rng[A]): Boolean = macro FpFilter.ltEqImpl[A]
def >=(rhs: FpFilter[A])(implicit ev0: Signed[A], ev1: Rng[A]): Boolean = macro FpFilter.gtEqImpl[A]
def ===(rhs: FpFilter[A])(implicit ev0: Signed[A], ev1: Rng[A]): Boolean = macro FpFilter.eqImpl[A]
def signum(implicit ev: Signed[A]): Int = macro FpFilter.signImpl[A]

16. …  Macro  City
•  Operator  fusion
– No  intermediate  alloca9ons
•  In-­‐line  approxima9on  +  error  bounds
– Fast,  Double  arithme9c
•  Thunk  becomes  inner  defs
– Compile  down  to  private  methods

17. Turn  this…
(x + y).signum

18. …  into  this.
val fpf\$tmp\$macro\$38 = x.value;
val fpf\$apx\$macro\$39 = fpf\$tmp\$macro\$38;
val fpf\$mes\$macro\$40 = java.lang.Math.abs(fpf\$tmp\$macro\$38);
def fpf\$exact\$macro\$42 = spire.algebra.Field.apply[spire.math.Algebraic](Algebraic.AlgebraicAlgebra).fromDouble(fpf
\$tmp\$macro\$38);
val fpf\$tmp\$macro\$43 = y.value;
val fpf\$apx\$macro\$44 = fpf\$tmp\$macro\$43;
val fpf\$mes\$macro\$45 = java.lang.Math.abs(fpf\$tmp\$macro\$43);
def fpf\$exact\$macro\$47 = spire.algebra.Field.apply[Algebraic](Algebraic.AlgebraicAlgebra).fromDouble(fpf\$tmp\$macro
\$43);
val fpf\$apx\$macro\$48 = fpf\$apx\$macro\$39.+(fpf\$apx\$macro\$44);
val fpf\$mes\$macro\$49 = fpf\$mes\$macro\$40.+(fpf\$mes\$macro\$45);
def fpf\$exact\$macro\$51 = Algebraic.AlgebraicAlgebra.plus(
fpf\$exact\$macro\$42, fpf\$exact\$macro\$47);
val fpf\$err\$macro\$52 = fpf\$mes\$macro\$49.\$times(1).\$times(2.220446049250313E-16);
if (fpf\$apx\$macro\$48 > fpf\$err\$macro\$52 && fpf\$apx\$macro\$48 < Double.POSITIVE_INFINITY)
1
else
if (fpf\$apx\$macro\$48 < fpf\$err\$macro\$52.unary_\$minus && fpf\$apx\$macro\$48 > Double.NEGATIVE_INFINITY)
-1
else
if (fpf\$err\$macro\$52 == 0.0)
0
else
Algebraic.AlgebraicAlgebra.signum(fpf\$exact\$macro\$51)

19. Examples

20. 2D  Orienta2on

21. p
q
r

22. p
q
r

23. p
q
r
RIGHT

24. p
r
q

25. p
r
q
LEFT

26. p
r
q

27. p
r
q
NO  TURN

28. trait Turn[@spec A] {
def apply(
px: A, py: A,
qx: A, qy: A,
rx: A, ry: A
): Int
}

29. object FastTurn extends
Turn[Double] {
def apply(
px: Double, py: Double,
qx: Double, qy: Double,
rx: Double, ry: Double
): Int = signum {
(qx - px) * (ry - py) -
(rx - px) * (qy - py)
}
}

30. Accuracy  of  Fast  Turn

31. object ExactTurn extends Turn[Double] {
def apply(
px: Double, py: Double,
qx: Double, qy: Double,
rx: Double, ry: Double
): Int = {
val pxa = Algebraic(px)
val pya = Algebraic(py)
val qxa = Algebraic(qx)
val qya = Algebraic(qy)
val rxa = Algebraic(rx)
val rya = Algebraic(ry)
((qxa - pxa) * (rya - pya) –
(rxa - pxa) * (qya - pya)).signum
}
}

32. 10,000x  Slower!

33. Let’s  try  again…

34. object FilteredTurn extends Turn[Double] {
def apply(
px: Double, py: Double,
qx: Double, qy: Double,
rx: Double, ry: Double
): Int = {
val pxf = FpFilter.exact[Algebraic](px)
val pyf = FpFilter.exact[Algebraic](py)
val qxf = FpFilter.exact[Algebraic](qx)
val qyf = FpFilter.exact[Algebraic](qy)
val rxf = FpFilter.exact[Algebraic](rx)
val ryf = FpFilter.exact[Algebraic](ry)
((qxf - pxf) * (ryf - pyf) –
(rxf - pxf) * (qyf - pyf)).signum
}
}

35. FilteredTurn  Speed
Rela9ve  to  FastTurn

36. Polynomial  Root  Finding

37. Polynomial[A]

38. Interval[A]
Root

Real  Roots.”  John  AbboT.
QIR  for  short.

40. QIR  (N=8)

41. QIR  (N=8)

42. QIR  (N=8)

43. QIR  (N=8)

44. QIR  (N=8)

45. QIR  (N=8)

46. QIR  (N=8)

47. QIR  (N=8)

48. QIR  (N=16)

49. QIR  (N=16)

50. QIR  (N=16)

51. QIR
•  Requires  2  polynomial  evalua9ons
– High  precision  generally  required
– Under  some  assump9ons
•  Occasionally  fails  when  assump9ons  not  met
– Fallsback  to  bisec9on!

52. QIR  (N=8)

53. QIR  (N=8)

54. QIR  (N=8)
FAILED!

55. Falls  back  to  bisec9on…

56. Bisec9on

57. Bisec9on

58. Bisec9on

59. Bisec9on

60. Bisec9on

61. Bisec9on

62. Bisec9on

63. Bisec9on

Requires  only  sign  tests

Converges  slowly,  1  bit  at  a  9me!

64. trait SignTest[A] {
def apply(
poly: Polynomial[A],
x: A
): Sign
}

65. final class FilteredSignTest[@sp A: Semiring](
implicit A: IsAlgebraic[A]
) extends SignTest[A] {
def apply(poly: Polynomial[A], x: A): Sign = {
val x0 = FpFilter(A.toDouble(x), A.toAlgebraic(x))
@tailrec
def loop(acc: FpFilter[Algebraic], i: Int): Sign =
if (i >= 0) {
val c = poly.nth(i)
val cftr = FpFilter(A.toDouble(c),
A.toAlgebraic(c))
loop(cftr + acc * x0, i - 1)
} else {
Sign(acc.signum)
}
loop(FpFilter.approx(Algebraic.Zero), poly.degree)
}
}

66. Accuracy  Using  Double

67. Speed  Up  from  Exact  Sign  Test
Fast  (d=8)
Fast  (d=16)
Fast  (d=32)
Filtered  (d=8)
Filtered  (d=16)
Filtered  (d=16)

68. Summary
•  Works  like  any  other  number  type
– Operator  fusion  +  inlining  within  expressions
•  Speeds  up  predicates
– Sign  tests,  comparisons,  etc.
•  Near-­‐Double  performance
– 2-­‐4x  in  most  cases
h;p://github.com/9xxit/fpﬁlter-­‐talk

69. Thanks!
h;p://github.com/non/spire

Tom  Switzer  @9xxit