Sư đoàn Newton-Raphson với số nguyên lớn

Question

Tôi đang làm một lớp BigInt như một bài tập lập trình. Nó sử dụng một vectơ của phần bổ sung 2 ký tự int trong cơ sở-65536 (để các phép nhân 32 bit không tràn. Tôi sẽ tăng cơ sở khi tôi làm cho nó hoạt động hoàn toàn).

Tất cả các hoạt động toán học cơ bản được mã hóa, với một vấn đề: phân chia là đau đớn chậm với thuật toán cơ bản mà tôi có thể tạo. (Nó hoạt động như phân chia nhị phân cho mỗi chữ số của thương ... Tôi sẽ không đăng nó trừ khi ai đó muốn xem nó ....)

Thay vì thuật toán chậm, tôi muốn sử dụng Newton-Raphson để tìm (biến đổi) đối ứng và sau đó nhân (và thay đổi). Tôi nghĩ rằng tôi có đầu của tôi xung quanh những điều cơ bản: bạn đưa ra công thức (x1 = x0 (2 - x0 * ước)) đoán ban đầu tốt, và sau đó sau một số lần lặp lại, x hội tụ với nghịch đảo. Phần này có vẻ dễ dàng đủ ... nhưng tôi đang chạy vào một số vấn đề khi cố gắng áp dụng công thức này để số nguyên lớn:

Vấn đề 1:

Bởi vì tôi đang làm việc với số nguyên ... cũng .. Tôi không thể sử dụng phân số. Điều này dường như gây ra x để luôn phân kỳ (x0 * ước số phải là < 2 có vẻ như?). Trực giác của tôi nói với tôi rằng nên có một số sửa đổi phương trình mà sẽ cho phép nó để làm việc số nguyên (một số độ chính xác) nhưng tôi thực sự đấu tranh để tìm hiểu nó là gì. (Tôi thiếu kỹ năng toán học đang đánh tôi ở đây ....) Tôi nghĩ rằng tôi cần phải tìm một số phương trình tương đương, thay vì d có d * [base^somePower]? Có thể có một số phương trình như (x1 = x0 (2 - x0 * d)) hoạt động với các số nguyên không?

Vấn đề 2:

Khi tôi sử dụng công thức Newton để tìm nghịch đảo của một số con số, kết quả kết thúc lên được chỉ là một phe nhỏ dưới đây những gì câu trả lời nên ... cũ. khi cố gắng tìm đối ứng của 4 (trong hệ thập phân):

x0 = 0.3 
x1 = 0.24 
x2 = 0.2496 
x3 = 0.24999936 
x4 = 0.2499999999983616 
x5 = 0.24999999999999999999998926258176 

Nếu tôi được đại diện cho số trong cơ số 10, tôi sẽ muốn có một kết quả của 25 (và nhớ đến sản phẩm dịch phải bằng 2). Với một số đối ứng như 1/3, bạn có thể chỉ cần cắt ngắn kết quả sau khi bạn biết mình có đủ độ chính xác. Nhưng làm thế nào tôi có thể rút ra các đối ứng chính xác từ kết quả trên?

Xin lỗi nếu điều này quá mơ hồ hoặc nếu tôi yêu cầu quá nhiều. Tôi đã xem qua Wikipedia và tất cả các tài liệu nghiên cứu tôi có thể tìm thấy trên Google, nhưng tôi cảm thấy như tôi đang đập đầu vào tường. Tôi đánh giá cao sự giúp đỡ của bất cứ ai có thể cho tôi!

...

Chỉnh sửa: Got làm việc thuật toán, mặc dù nó là chậm hơn nhiều so với mong đợi của tôi. Tôi thực sự mất rất nhiều tốc độ so với thuật toán cũ của tôi, thậm chí trên con số với hàng nghìn chữ số ... Tôi vẫn còn thiếu cái gì đó. Nó không phải là một vấn đề với phép nhân, mà là rất nhanh. (Tôi thực sự sử dụng thuật toán của Karatsuba).

Đối với bất cứ ai quan tâm, đây là lần lặp hiện tại của tôi của thuật toán Newton-Raphson:

bigint operator/(const bigint& lhs, const bigint& rhs) { 
    if (rhs == 0) throw overflow_error("Divide by zero exception"); 
    bigint dividend = lhs; 
    bigint divisor = rhs; 

    bool negative = 0; 
    if (dividend < 0) { 
     negative = !negative; 
     dividend.invert(); 
    } 
    if (divisor < 0) { 
     negative = !negative; 
     divisor.invert(); 
    } 

    int k = dividend.numBits() + divisor.numBits(); 
    bigint pow2 = 1; 
    pow2 <<= k + 1; 

    bigint x = dividend - divisor; 
    bigint lastx = 0; 
    bigint lastlastx = 0; 
    while (1) { 
     x = (x * (pow2 - x * divisor)) >> k; 
     if (x == lastx || x == lastlastx) break; 
     lastlastx = lastx; 
     lastx = x; 
    } 
    bigint quotient = dividend * x >> k; 
    if (dividend - (quotient * divisor) >= divisor) quotient++; 
    if (negative)quotient.invert(); 
    return quotient; 
} 

Và đây là (thực sự xấu xí) thuật toán cũ của tôi đó là nhanh hơn:

bigint operator/(const bigint& lhs, const bigint & rhs) { 
    if (rhs == 0) throw overflow_error("Divide by zero exception"); 
    bigint dividend = lhs; 
    bigint divisor = rhs; 

    bool negative = 0; 
    if (dividend < 0) { 
     negative = !negative; 
     dividend.invert(); 
    } 
    if (divisor < 0) { 
     negative = !negative; 
     divisor.invert(); 
    } 

    bigint remainder = 0; 
    bigint quotient = 0; 
    while (dividend.value.size() > 0) { 
     remainder.value.insert(remainder.value.begin(), dividend.value.at(dividend.value.size() - 1)); 
     remainder.value.push_back(0); 
     remainder.unPad(); 
     dividend.value.pop_back(); 

     if (divisor > remainder) { 
      quotient.value.push_back(0); 
     } else { 
      int count = 0; 
      int i = MSB; 
      bigint value = 0; 
      while (i > 0) { 
       bigint increase = divisor * i; 
       bigint next = value + increase; 
       if (next <= remainder) { 
        value = next; 
        count += i; 
       } 
       i >>= 1; 
      } 
      quotient.value.push_back(count); 
      remainder -= value; 
     } 
    } 

    for (int i = 0; i < quotient.value.size()/2; i++) { 
     int swap = quotient.value.at(i); 
     quotient.value.at(i) = quotient.value.at((quotient.value.size() - 1) - i); 
     quotient.value.at(quotient.value.size() - 1 - i) = swap; 
    } 

    if (negative)quotient.invert(); 
    quotient.unPad(); 
    return quotient; 
} 

Answer 1

Trước hết, bạn có thể thực hiện phân chia theo thời gian O(n^2) và với hằng số hợp lý, vì vậy không chậm hơn nhiều so với phép nhân. Tuy nhiên, nếu bạn sử dụng thuật toán giống như Karatsuba hoặc thậm chí là FFT thuật toán nhân được dựa trên, thì bạn thực sự có thể tăng tốc thuật toán phân chia của mình bằng cách sử dụng Newton-Raphson.

Lặp lại Newton-Raphson để tính toán nghịch đảo của x là q[n+1]=q[n]*(2-q[n]*x).

Giả sử chúng tôi muốn tính floor(2^k/B) trong đó B là số nguyên dương. WLOG, B≤2^k; nếu không, thương là 0. Lặp lại Newton-Raphson cho x=B/2^k sản lượng q[n+1]=q[n]*(2-q[n]*B/2^k). chúng tôi có thể sắp xếp nó như

q[n+1]=q[n]*(2^(k+1)-q[n]*B) >> k

Mỗi lần lặp của loại hình này đòi hỏi chỉ phép nhân số nguyên và thay đổi chút. Liệu nó có hội tụ với floor(2^k/B) không? Không cần thiết. Tuy nhiên, trong trường hợp xấu nhất, nó cuối cùng thay đổi giữa floor(2^k/B) và ceiling(2^k/B) (Chứng minh điều đó!). Vì vậy, bạn có thể sử dụng một số thử nghiệm không thông minh để xem bạn có đang trong trường hợp này hay không và trích xuất floor(2^k/B). ("kiểm tra không thông minh" này sẽ nhanh hơn rất nhiều so với phép nhân trong mỗi lần lặp; Tuy nhiên, sẽ tốt hơn nếu tối ưu hóa điều này).

Thật vậy, tính toán floor(2^k/B) đủ để tính floor(A/B) cho bất kỳ số nguyên dương nào A,B. Hãy k sao cho A*B≤2^k và xác minh floor(A/B)=A*ceiling(2^k/B) >> k.

Cuối cùng, tối ưu hóa đơn giản nhưng quan trọng cho phương pháp này là cắt bớt phép nhân (tức là chỉ tính các bit cao hơn của sản phẩm) trong các lần lặp đầu của phương pháp Newton-Raphson. Lý do để làm như vậy, là kết quả của các lần lặp lại đầu tiên là xa thương, và nó không quan trọng để thực hiện chúng không chính xác. (Tinh chỉnh đối số này và chỉ ra rằng nếu bạn làm điều này một cách thích hợp, bạn có thể chia hai số nguyên ≤n-bit trong thời gian O(M(2n)), giả sử bạn có thể nhân hai số nguyên ≤k-bit trong thời gian M(k) và M(x) là hàm lồi ngày càng tăng).

Answer 2

Newton- Raphson là một thuật toán xấp xỉ - không thích hợp để sử dụng trong toán học số nguyên. Bạn sẽ nhận được lỗi làm tròn mà sẽ dẫn đến các loại vấn đề bạn đang nhìn thấy. Bạn có thể thực hiện vấn đề với các số dấu phẩy động và sau đó xem nếu bạn nhận được số nguyên, chính xác đến một số chữ số được chỉ định (xem đoạn tiếp theo)

Về vấn đề thứ hai, chọn độ chính xác (số chữ số thập phân) muốn cho độ chính xác và tròn với độ chính xác đó. Nếu bạn chọn hai mươi chữ số chính xác trong bài toán, bạn sẽ làm tròn đến 0,25. Bạn chỉ cần lặp lại cho đến khi các chữ số chính xác yêu cầu của bạn ổn định. Nói chung, đại diện cho số không hợp lý trên một máy tính thường giới thiệu sự thiếu chính xác.

Answer 3

Nếu tôi thấy điều này đúng thì một cải tiến chính là chọn giá trị bắt đầu tốt cho x.Biết bao nhiêu chữ số ước có bạn biết được nơi các bit quan trọng nhất của nghịch đảo có được, như

1/x = pow(2,log2(1/x)) 
1/x = pow(2,-log2(x)) 
1/x >= pow(2,-floor(log2(x))) 

sàn (log2 (x)) chỉ đơn giản là chỉ số của tập bit quan trọng nhất.