Skip to content

Commit

Permalink
Improves accuracy of vector numerical examples (#1773)
Browse files Browse the repository at this point in the history
* Improves accuracy of vector division approximation example.

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.

* Improves accuracy of vector square root approximation example.

The previous implementation was unnecessarily inaccurate:

  Measured max error of 2.647303869947791 ulp @ 2.0831878e-19
        func (0x1.ebe0fap-63) = 0x1.f5d6b2p-32
    func_ref (0x1.ebe0fap-63) = 0x1.f5d6acb494965p-32
  2'139'095'041 tested in [0 , inf]
   57.2114955877736548% correctly rounded  (  1223808265)
   88.4557938629712339% faithfully rounded (+  668345235)
   11.5442061370287661% >= 1.0 ulp         (+  246941541)

Which is 21.6 bits of precision rater than the claimed 23.

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

  Measured max error of 0.800413005053997 ulp @ 3.0387554e-29
        func (0x1.342a9ap-95) = 0x1.8d378p-48
    func_ref (0x1.342a9ap-95) = 0x1.8d378199cfbbcp-48
  2'139'095'041 tested in [0 , inf]
   87.4400844352198163% correctly rounded  (  1870426510)
  100.0000000000000000% faithfully rounded (+  268668531)

With is 23.3 bits of precision.

* src/vector-examples.adoc (Square root approximation example): Use
`r - (0.5 r) ((x r) r - 1)` formulation for 1/sqrt(x) estimate
refinement.
  • Loading branch information
ebavier authored Dec 14, 2024
1 parent 89e3947 commit 44e4965
Showing 1 changed file with 24 additions and 24 deletions.
48 changes: 24 additions & 24 deletions src/vector-examples.adoc
Original file line number Diff line number Diff line change
Expand Up @@ -82,38 +82,38 @@ 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
----

=== Square root approximation example

----
# v1 = sqrt(v1) to almost 23 bits of precision.
# v1 = sqrt(v1) to more than 23 bits of precision.
fmv.w.x ft0, x0 # Mask off zero inputs
vmfne.vf v0, v1, ft0 # to avoid div by zero
vfrsqrt7.v v2, v1, v0.t # Estimate 1/sqrt(x)
vmfne.vf v0, v2, ft0, v0.t # Additionally mask off +inf inputs
li t0, 0x40400000
vmv.v.x v4, t0 # Splat 3.0
vfmul.vv v3, v1, v2, v0.t # x * est
vfnmsub.vv v3, v2, v4, v0.t # - x * est * est + 3
vfmul.vv v3, v3, v2, v0.t # est * (-x * est * est + 3)
li t0, 0x3f000000
fmv.w.x ft0, t0 # 0.5
vfmul.vf v2, v3, ft0, v0.t # Estimate to 14 bits
vfmul.vv v3, v1, v2, v0.t # x * est
vfnmsub.vv v3, v2, v4, v0.t # - x * est * est + 3
vfmul.vv v3, v3, v2, v0.t # est * (-x * est * est + 3)
vfmul.vf v2, v3, ft0, v0.t # Estimate to 23 bits
vfmul.vv v1, v2, v1, v0.t # x * 1/sqrt(x)
vmfne.vf v0, v1, ft0 # to avoid DZ exception
vfrsqrt7.v v2, v1, v0.t # Estimate r ~= 1/sqrt(v1)
vmfne.vf v0, v2, ft0, v0.t # Mask off +inf to avoid NV
li t0, 0x3f800000
fli.s ft0, 0.5
vmv.v.x v5, t0 # Splat 1.0
vfmul.vv v3, v1, v2, v0.t # t = v1 r
vfmul.vf v4, v2, ft0, v0.t # 0.5 r
vfmsub.vv v3, v2, v5, v0.t # t r - 1
vfnmsac.vv v2, v3, v4, v0.t # r - (0.5 r) (t r - 1)
# Better estimate of 1/sqrt(v1)
vfmul.vv v1, v1, v2, v0.t # t = v1 r
vfmsub.vv v2, v1, v5, v0.t # t r - 1
vfmul.vf v3, v1, ft0, v0.t # 0.5 t
vfnmsac.vv v1, v2, v3, v0.t # t - (0.5 t) (t r - 1)
# ~ sqrt(v1) to about 23.3 bits
----

=== C standard library strcmp example
Expand Down

0 comments on commit 44e4965

Please sign in to comment.