Thuật toán nhanh nhất để xác định x nhỏ nhất và lớn nhất tạo phương trình chính xác kép x + a == b true

Question

Trong ngữ cảnh phân tích tĩnh, tôi quan tâm đến việc xác định giá trị của x trong nhánh sau của điều kiện dưới đây:

double x; 
x = …; 
if (x + a == b) 
{ 
    … 

a và b có thể được giả định là hằng đúp chính xác (khái quát hóa để biểu thức tùy ý là phần dễ nhất của vấn đề), và trình biên dịch có thể được giả định theo chuẩn IEEE 754 nghiêm (FLT_EVAL_METHOD là 0). Chế độ làm tròn tại thời gian chạy có thể được giả định là gần nhất.

Nếu tính toán với lý trí rẻ, nó sẽ đơn giản: các giá trị cho x sẽ là số chính xác kép có trong khoảng thời gian hợp lý (b - a - 0.5 * ulp1 (b)… b - a + 0.5 * ulp2 (b)). Các giới hạn nên được bao gồm nếu b thậm chí bị loại trừ nếu b là lẻ, và ulp1 và ulp2 là hai định nghĩa hơi khác nhau của "ULP" có thể được lấy giống hệt nhau nếu không nhớ mất một chút chính xác về sức mạnh của hai.

Thật không may, tính toán với hợp lý có thể tốn kém. Xem xét khả năng khác là thu được từng giới hạn bằng cách phân đôi, trong 64 phép cộng hai lần chính xác (mỗi thao tác quyết định một bit kết quả). 128 bổ sung dấu phẩy động để có được giới hạn dưới và trên có thể nhanh hơn bất kỳ giải pháp nào dựa trên toán học.

Tôi tự hỏi liệu có cách nào cải thiện ý tưởng “bổ sung 128 điểm nổi” hay không. Trên thực tế, tôi có giải pháp riêng của mình liên quan đến các thay đổi về chế độ làm tròn và các cuộc gọi nextafter nhưng tôi không muốn làm hỏng phong cách của bất kỳ ai và khiến họ bỏ lỡ một giải pháp thanh lịch hơn so với giải pháp hiện tại của tôi. Ngoài ra tôi không chắc chắn rằng việc thay đổi chế độ làm tròn hai lần thực sự rẻ hơn 64 điểm bổ sung nổi.

Answer 1

Bạn đã đã đưa ra một giải pháp tốt đẹp và thanh lịch trong câu hỏi của bạn:

Nếu máy tính với rationals là rẻ, nó sẽ là rất đơn giản: các giá trị cho x sẽ là những con số tăng gấp đôi độ chính xác chứa trong hợp lý khoảng thời gian (b - a - 0,5 * ulp1 (b)… b - a + 0,5 * ulp2 (b)). Các giới hạn nên được bao gồm nếu b là chẵn, nếu b là lẻ, và ulp1 và ulp2 là hai định nghĩa hơi khác nhau của “ULP” có thể được chụp giống hệt nhau nếu không nhớ mất một chút chính xác về quyền hạn của hai.

Điều gì sau đây là bản phác thảo một nửa lý thuyết của giải pháp từng phần cho vấn đề dựa trên đoạn này. Hy vọng rằng tôi sẽ có cơ hội để xác thịt nó ra sớm. Để có được một giải pháp thực sự, bạn sẽ phải xử lý các subnormals, zeroes, NaNs và tất cả những thứ thú vị khác. Tôi sẽ giả định rằng a và b là, chẳng hạn như 1e-300 < |a| < 1e300 và 1e-300 < |b| < 1e300 để không xảy ra bất kỳ lúc nào.

Lỗi tràn và tràn, bạn có thể nhận được ulp1(b) từ b - nextafter(b, -1.0/0.0). Bạn có thể nhận được ulp2(b) từ nextafter(b, 1.0/0.0) - b.

Nếu b/2 <= a <= 2b, thì định lý của Sterbenz cho bạn biết rằng b - a là chính xác. Vì vậy, (b - a) - ulp1/2 sẽ là double gần nhất với giới hạn dưới và (b - a) + ulp2/2 sẽ là double gần nhất với giới hạn trên.Hãy thử các giá trị này và các giá trị ngay trước và sau và chọn khoảng thời gian rộng nhất hoạt động.

Nếu b > 2a, b - a > b/2. Giá trị tính toán của b - a bị tắt tối đa bằng một nửa ulp. Một ulp1 có tối đa hai ulp, giống như một số ulp2, vì vậy khoảng thời gian hợp lý bạn đã cho tối đa là hai vạch ulp. Tìm ra trong số năm giá trị gần nhất với công việc b-a.

Nếu a > 2b, một ulp của b-a ít nhất cũng lớn bằng một ulp b; nếu bất cứ điều gì làm việc, tôi đặt cược nó sẽ có được trong số ba giá trị gần nhất để b-a. Tôi hình dung trường hợp a và b có các dấu hiệu khác nhau hoạt động tương tự.

Tôi đã viết một chồng nhỏ mã C++ triển khai ý tưởng này. Nó không thất bại trong việc kiểm tra lông tơ ngẫu nhiên (trong một vài phạm vi khác nhau) trước khi tôi cảm thấy chán chờ đợi. Dưới đây là:

void addeq_range(double a, double b, double &xlo, double &xhi) { 
    if (a != a) return; // empty interval 
    if (b != b) { 
    if (a-a != 0) { xlo = xhi = -a; return; } 
    else return; // empty interval 
    } 
    if (b-b != 0) { 
    // TODO: handle me. 
    } 

    // b is now guaranteed to be finite. 
    if (a-a != 0) return; // empty interval 

    if (b < 0) { 
    addeq_range(-a, -b, xlo, xhi); 
    xlo = -xlo; 
    xhi = -xhi; 
    return; 
    } 

    // b is now guaranteed to be zero or positive finite and a is finite. 
    if (a >= b/2 && a <= 2*b) { 
    double upulp = nextafter(b, 1.0/0.0) - b; 
    double downulp = b - nextafter(b, -1.0/0.0); 
    xlo = (b-a) - downulp/2; 
    xhi = (b-a) + upulp/2; 
    if (xlo + a == b) { 
     xlo = nextafter(xlo, -1.0/0.0); 
     if (xlo + a != b) xlo = nextafter(xlo, 1.0/0.0); 
    } else xlo = nextafter(xlo, 1.0/0.0); 
    if (xhi + a == b) { 
     xhi = nextafter(xhi, 1.0/0.0); 
     if (xhi + a != b) xhi = nextafter(xhi, -1.0/0.0); 
    } else xhi = nextafter(xhi, -1.0/0.0); 
    } else { 
    double xmid = b-a; 
    if (xmid + a < b) { 
     xhi = xlo = nextafter(xmid, 1.0/0.0); 
     if (xhi + a != b) xhi = xmid; 
    } else if (xmid + a == b) { 
     xlo = nextafter(xmid, -1.0/0.0); 
     xhi = nextafter(xmid, 1.0/0.0); 
     if (xlo + a != b) xlo = xmid; 
     if (xhi + a != b) xhi = xmid; 
    } else { 
     xlo = xhi = nextafter(xmid, -1.0/0.0); 
     if (xlo + a != b) xlo = xmid; 
    } 
    } 
}