Higham의 수치 알고리즘 의 정확성과 안정성은 이러한 유형의 문제를 분석 할 수있는 방법을 다루고 있다고 생각 합니다. 2 장, 특히 운동 2.8을 참조하십시오.

이 답변에서 나는 Higham의 책에서 실제로 다루지 않은 것을 지적하고 싶습니다 (그 문제에 대해서는 널리 알려지지 않은 것으로 보입니다). 이와 같은 간단한 수치 알고리즘의 속성을 증명 하려는 경우 Haskell의 sbv 와 같은 패키지를 사용하여 z3 과 같은 최신 SMT 솔버 ( 만족도 모듈로 이론 ) 의 기능을 사용할 수 있습니다 . 이것은 연필과 종이를 사용하는 것보다 다소 쉽습니다.

주어졌고 z = ( x + y ) / 2 가 x ≤ z ≤ y를 만족 하는지 알고 싶습니다 . 다음 하스켈 코드0≤x≤y

z=(x+y)/2

x≤z≤y

import Data.SBV

test1 :: (SFloat -> SFloat -> SFloat) -> Symbolic SBool
test1 fun =
  do [x, y] <- sFloats ["x", "y"]
     constrain $ bnot (isInfiniteFP x) &&& bnot (isInfiniteFP y)
     constrain $ 0 .<= x &&& x .<= y
     let z = fun x y
     return $ x .<= z &&& z .<= y

test2 :: (SFloat -> SFloat -> SFloat) -> Symbolic SBool
test2 fun =
  do [x, y] <- sFloats ["x", "y"]
     constrain $ bnot (isInfiniteFP x) &&& bnot (isInfiniteFP y)
     constrain $ x .<= y
     let z = fun x y
     return $ x .<= z &&& z .<= y

이 작업을 자동으로 수행하겠습니다 . 다음 test1 fun은 IS 명제는 그 모든 유한 수레 X , Y 와 0 ≤ X ≤ Y .x≤fun(x,y)≤y

x,y

0≤x≤y

λ> prove $ test1 (\x y -> (x + y) / 2)
Falsifiable. Counter-example:
  x = 2.3089316e36 :: Float
  y = 3.379786e38 :: Float

넘친다. 이제 다른 공식을 취한다고 가정 해보십시오 z=x/2+y/2

λ> prove $ test1 (\x y -> x/2 + y/2)
Falsifiable. Counter-example:
  x = 2.3509886e-38 :: Float
  y = 2.3509886e-38 :: Float

점진적 언더 플로로 인해 작동하지 않습니다 ( 모든 산술이 기본 2 이므로 직관적이지 않을 수 있음).(x/2)×2≠x

이제 :z=x+(y−x)/2

λ> prove $ test1 (\x y -> x + (y-x)/2)
Q.E.D.

공장! 는 Q.E.D.A는 증거 것을 test1위에서 정의 된 속성은 모든 수레를 위해 보유하고 있습니다.

동일하지만 ( 0 ≤ x ≤ y 대신)로 제한되는 것은 어떻습니까?x≤y

0≤x≤y

λ> prove $ test2 (\x y -> x + (y-x)/2)
Falsifiable. Counter-example:
  x = -3.1300826e34 :: Float
  y = 3.402721e38 :: Float

자, 넘치면 z = x + ( y / 2 – x / 2 )는 어떻습니까?y−x

z=x+(y/2−x/2)

λ> prove $ test2 (\x y -> x + (y/2 - x/2))
Q.E.D.

그래서 내가 여기에서 시도한 수식 중 가 작동하는 것처럼 보입니다 (증거와 함께). SMT 솔버 접근법은 연필과 종이로 부동 소수점 오류 분석을 수행하는 것보다 간단한 부동 소수점 수식에 대한 의심에 훨씬 빠른 방법으로 보입니다.x+(y/2−x/2)

마지막으로 정확성과 안정성의 목표는 종종 성능의 목표와 상충됩니다. 성능을 위해, 나는 컴파일러가 보다 나은 방법을 실제로 보지 못합니다 . 특히 컴파일러는 이것을 기계 명령어로 번역하는 데 많은 노력을 기울이고 있기 때문입니다.(x+y)/2

PS 이것은 모두 단 정밀도 IEEE754 부동 소수점 산술입니다. 나는 검사 배정 밀도 연산 (교체와를 함께 ), 그것도 작동합니다.x≤x+(y/2−x/2)≤y

SFloatSDouble

-ffast-math(x+y)/2

PPPS 조건부없이 간단한 대수적 표현 만 살펴 보았습니다 . 돈 해치 의 공식 은 엄청나게 좋습니다.

답변

if max(abs(x),abs(y)) >= 1.:
    return x/2. + y/2.
else:
    return (x+y)/2.

double average (double x, double y)
{
    const double C = 1; /* 0x1p-967 <= C <= 0x1p970 */
    return (C <= fabs (x)) ? (x / 2 + y / 2) : ((x + y) / 2);
}

float difference = max(x, y) - min(x, y);
return min(x, y) + (difference / 2.0);