Skip to content

Commit

Permalink
Improves accuracy of vector division approximation example.
Browse files Browse the repository at this point in the history
The previous implementation was unnecessarily inaccurate:

  Measured max error of 2.465257280 ulp @ (1.2449559e-01, 5.2960633e+36)
        func (0x1.fdef16p-4, 0x1.fdfe0ap+121) = 0x1.fff0f8p-126
    func_ref (0x1.fdef16p-4, 0x1.fdfe0ap+121) = 0x1.fff0fcee3633cp-126
  4'550'824'890 tested in [0 ; inf] x [0x1.000008p-128 ; 0x1p126]
   68.8976110218998059% correctly rounded  (  3135409631)
   92.7973265084255985% faithfully rounded (+ 1087634201)
    7.2026734915744033% >= 1.0 ulp         (+  327781058)

Which is 21.7 bits of precision rather than the claimed 23.

The new implementation has lower maximum error and higher correctly-
and faithfully-rounded rates:

  Measured max error of 1.493816005 ulp @ (5.0745111e+02, 4.9807361e-01)
        func (0x1.fb737cp+8, 0x1.fe0702p-2) = 0x1.fd69eap+9
    func_ref (0x1.fb737cp+8, 0x1.fe0702p-2) = 0x1.fd69ecfcd5739p+9
  4'550'824'890 tested in [0 ; inf] x [0x1.000008p-128 ; 0x1p126]
   79.6106155822664476% correctly rounded  (  3622939709)
   99.0013070355691127% faithfully rounded (+  882436413)
    0.9986929644308947% >= 1.0 ulp         (+   45448768)

Which is 22.4 bits of precision.

Max error is greater when the denominator's reciprocal underflows:
previously about 5.26 ulp and now about 4.45 ulp.

The algorithm does not handle infinite or zero denominators (returns
NaN) or subnormal denominators whose reciprocal overflows (returns
±∞).  This commit does not document or change this behavior.

* src/vector-examples.adoc (Division approximation example): Use
`r + r (1 - y r)` formulation for 1/y estimate refinement.
  • Loading branch information
ebavier committed Dec 6, 2024
1 parent a97e505 commit cb19724
Showing 1 changed file with 7 additions and 7 deletions.
14 changes: 7 additions & 7 deletions src/vector-examples.adoc
Original file line number Diff line number Diff line change
Expand Up @@ -82,13 +82,13 @@ include::example/sgemm.S[lines=4..-1]
# v1 = v1 / v2 to almost 23 bits of precision.
vfrec7.v v3, v2 # Estimate 1/v2
li t0, 0x40000000
vmv.v.x v4, t0 # Splat 2.0
vfnmsac.vv v4, v2, v3 # 2.0 - v2 * est(1/v2)
vfmul.vv v3, v3, v4 # Better estimate of 1/v2
vmv.v.x v4, t0 # Splat 2.0
vfnmsac.vv v4, v2, v3 # 2.0 - v2 * est(1/v2)
vfmul.vv v3, v3, v4 # Better estimate of 1/v2
li t0, 0x3f800000
vmv.v.x v4, t0 # Splat 1.0
vfnmsac.vv v4, v2, v3 # 1.0 - v2 * est(1/v2)
vfmadd.vv v3, v4, v3 # Better estimate of 1/v2
vmv.v.x v4, t0 # Splat 1.0
vfnmsac.vv v4, v2, v3 # 1.0 - v2 * est(1/v2)
vfmadd.vv v3, v4, v3 # Better estimate of 1/v2
vfmul.vv v1, v1, v3 # Estimate of v1/v2
----

Expand Down

0 comments on commit cb19724

Please sign in to comment.