Speed, Correctness, or Simplicity: Choose 3

Speed, Correctness, or Simplicity: Choose 3

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

B2ecd194174398696f1d9a80bf6a14fa?s=128

Tom Switzer

January 30, 2015
Tweet

Transcript

  1. 1.

    Speed,  Correctness,  or  Simplicity:     Choose  3   Tom

     Switzer   @9xxit     h;ps://github.com/9xxit/fpfilter-­‐talk  
  2. 3.

    Op9ons   Use  Double,  live  with  the  errors.    

    Use  higher  precision  type,   live  with  performance  loss.     But,  there  is  a  3rd  op9on…  
  3. 4.

    Floa9ng  Point  Filters     Use  floa9ng  point  when  you

     can.     Use  higher  precision  when  you  can’t.  
  4. 5.

    Err…  Not  So  Simple   Solve  problem  using  floa9ng  point

      approxima9on…     Maintain  an  error  bound  on  approxima9on.     Re-­‐evaluate  with  exact  type  if  error  too  large.  
  5. 11.

    FpFilter[A]   Simple  wrapper:  FpFilter[Rational] Standard  Opera2ons   +,  -­‐,

     *,  /,  .sqrt,  etc     Fast  predictes   signum,  compare,  isWhole,  etc.  
  6. 12.

    FpFilter[A]   class FpFilter[A]( apx: Double, mes: Double, ind: Int,

    exact: => A ) { … } floa9ng  point  approxima9on   error  bounds  
  7. 13.

    FpFilter[A]   class FpFilter[A]( apx: Double, mes: Double, ind: Int,

    exact: => A ) { … } error  bounds   “Exact  Geometric  Computa2on   Using  Cascading.”   Burnikel,  Funke  &  Seel.  
  8. 14.

    FpFilter[A]   class FpFilter[A]( apx: Double, mes: Double, ind: Int,

    exact: => A ) { … } error  bounds   thunk  for  higher  precision   Welcome  to  …  
  9. 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]
  10. 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  
  11. 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)
  12. 28.

    trait Turn[@spec A] { def apply( px: A, py: A,

    qx: A, qy: A, rx: A, ry: A ): Int }
  13. 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) } }
  14. 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 } }
  15. 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 } }
  16. 51.

    QIR   •  Requires  2  polynomial  evalua9ons   – High  precision

     generally  required   •  Very  fast  convergence  (quadra9c)   – Under  some  assump9ons   •  Occasionally  fails  when  assump9ons  not  met   – Fallsback  to  bisec9on!  
  17. 63.

    Bisec9on     Requires  only  sign  tests     Converges

     slowly,  1  bit  at  a  9me!  
  18. 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) } }
  19. 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)  
  20. 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/fpfilter-­‐talk