2012-12-02 40 views
6

Tôi cần tính sin(4^x) với x> 1000 trong Matlab, về cơ bản là sin(4^x mod 2π) Vì các giá trị bên trong hàm sin trở nên rất lớn, Matlab trả về vô hạn cho 4^1000. Làm thế nào tôi có thể tính toán hiệu quả điều này? Tôi muốn tránh các loại dữ liệu lớn.Tính 4^x mod 2π cho lớn x

Tôi nghĩ rằng việc chuyển đổi thành một thứ như sin(n*π+z) có thể là giải pháp khả thi.

+1

Bạn cũng có thể muốn xem [mathoverflow] (http://mathoverflow.net/) cho các câu hỏi về toán học - mặc dù điều này không liên quan đến lập trình. – Jeff

+2

Tôi có linh cảm rằng mọi phương pháp sẽ dẫn đến rác do số học chính xác gấp đôi: ví dụ: tội lỗi tính toán (2 * pi * 10^i) cho i = 1: 100 cho junk trong Matlab (hoặc C cho rằng vấn đề ... nó thực sự là vấn đề số học chính xác gấp đôi), và tính toán tương tự mod (2 * pi * 10^i + 0,0001, 2 * pi) cho junk cho i = 1: 100. – db1234

+0

Đó cũng là nỗi sợ của tôi. Do đó, tôi thường cố gắng tìm các giải pháp toán học trước khi sử dụng vòng lặp. Cảm ơn bạn. – Bene

Trả lời

13

Bạn cần phải cẩn thận vì sẽ mất độ chính xác. Hàm sin là định kỳ, nhưng 4^1000 là một số lớn. Vì vậy, có hiệu quả, chúng tôi trừ một bội số của 2 * pi để di chuyển đối số vào khoảng [0,2 * pi).

4^1000 là khoảng 1e600, một số thực sự lớn. Vì vậy, tôi sẽ làm tính toán của mình bằng cách sử dụng high precision floating point tool in MATLAB. (Trong thực tế, một trong những mục tiêu rõ ràng của tôi khi tôi viết HPF là có thể tính toán một số giống như tội lỗi (1e400). Ngay cả khi bạn đang làm một cái gì đó cho niềm vui của nó, làm nó vẫn hợp lý.) Trong trường hợp này , vì tôi biết rằng sức mạnh chúng ta quan tâm là khoảng 1e600, sau đó tôi sẽ tính toán của tôi với hơn 600 chữ số chính xác, hy vọng rằng tôi sẽ mất 600 chữ số bằng cách hủy bỏ trừ. Đây là một vấn đề hủy lớn. Hãy suy nghĩ về nó. Hoạt động mô đun đó có hiệu quả là sự khác biệt giữa hai con số sẽ giống hệt nhau cho 600 chữ số đầu tiên hoặc hơn thế!

X = hpf(4,1000); 
X^1000 
ans = 
114813069527425452423283320117768198402231770208869520047764273682576626139237031385665948631650626991844596463898746277344711896086305533142593135616665318539129989145312280000688779148240044871428926990063486244781615463646388363947317026040466353970904996558162398808944629605623311649536164221970332681344168908984458505602379484807914058900934776500429002716706625830522008132236281291761267883317206598995396418127021779858404042159853183251540889433902091920554957783589672039160081957216630582755380425583726015528348786419432054508915275783882625175435528800822842770817965453762184851149029376 

Điểm gần nhất của 2 * pi không vượt quá con số này là bao nhiêu? Chúng ta có thể làm được điều đó bằng một thao tác đơn giản.

twopi = 2*hpf('pi',1000); 
twopi*floor(X^1000/twopi) 
ans = 114813069527425452423283320117768198402231770208869520047764273682576626139237031385665948631650626991844596463898746277344711896086305533142593135616665318539129989145312280000688779148240044871428926990063486244781615463646388363947317026040466353970904996558162398808944629605623311649536164221970332681344168908984458505602379484807914058900934776500429002716706625830522008132236281291761267883317206598995396418127021779858404042159853183251540889433902091920554957783589672039160081957216630582755380425583726015528348786419432054508915275783882625175435528800822842770817965453762184851149029372.6669043995793459614134256945369645075601351114240611660953769955068077703667306957296141306508448454625087552917109594896080531977700026110164492454168360842816021326434091264082935824243423723923797225539436621445702083718252029147608535630355342037150034246754736376698525786226858661984354538762888998045417518871508690623462425811535266975472894356742618714099283198893793280003764002738670747 

Như bạn có thể thấy, 600 chữ số đầu tiên giống nhau. Bây giờ, khi chúng tôi trừ hai số,

X^1000 - twopi*floor(X^1000/twopi) 
ans = 
3.333095600420654038586574305463035492439864888575938833904623004493192229633269304270385869349155154537491244708289040510391946802229997388983550754583163915718397867356590873591706417575657627607620277446056337855429791628174797085239146436964465796284996575324526362330147421377314133801564546123711100195458248112849130937653757418846473302452710564325738128590071680110620671999623599726132925263826 

Đây là lý do tại sao tôi gọi nó là một vấn đề hủy lớn trừ. Hai số này giống hệt nhau cho nhiều chữ số. Thậm chí mang theo 1000 chữ số chính xác, chúng tôi đã mất nhiều chữ số. Khi bạn trừ hai số, mặc dù chúng tôi đang mang kết quả với 1000 chữ số, chỉ có 400 chữ số cao nhất hiện có ý nghĩa.

HPF có thể tính toán chức năng trig của khóa học. Nhưng như chúng tôi đã trình bày ở trên, chúng ta chỉ nên tin tưởng khoảng 400 chữ số đầu tiên của kết quả. (Trên một số vấn đề, hình dạng cục bộ của hàm sin có thể khiến chúng ta mất nhiều chữ số hơn số đó.)

sin(X^1000) 
ans = 
-0.19033458127208318385994396068455455709388374041098639172943768418947125138650234240955423917696880832346734715448603532912993423621761996537053192685449334064870714463489747336279464911185192423229252660143128976923388511299599457104070322693060218958487584842139143972048735807765826659851362293280012583640059277583434162223469640779539703355744143419935430600390820454055891750089781440474478225522286222463738277009002753247363724815609283394633443329778920087022201603354152914210817007440447838392869577354385645124650950464218066771029610934877080889086985319804240164585346291661088530125354930225403524397401167317843031900829546691402971929428720760150282604082313216048252703439459284455892236101855653841958635139010896628829034919565066139672417258772760228631878006327065033172

Tôi có quyền không và chúng tôi không thể tin tưởng tất cả các chữ số này? Tôi sẽ làm cùng một tính toán, một lần trong 1000 chữ số của độ chính xác, sau đó một lần thứ hai trong 2000 chữ số. Tính toán sự khác biệt tuyệt đối, sau đó lấy log10. Kết quả 2000 chữ số sẽ là tham chiếu của chúng tôi về cơ bản chính xác so với kết quả 1000 chữ số.

double(log10(abs(sin(hpf(4,[1000 0])^1000) - sin(hpf(4,[2000 0])^1000)))) 
ans = 
     -397.45 

Ah. Vì vậy, trong số 1000 chữ số chính xác mà chúng tôi đã bắt đầu, chúng tôi mất 602 chữ số. 602 chữ số cuối cùng trong kết quả là khác 0 nhưng vẫn hoàn thành rác. Điều này đúng như tôi mong đợi. Chỉ vì máy tính của bạn báo cáo độ chính xác cao, bạn cần phải biết khi nào không tin tưởng nó.

Chúng tôi có thể thực hiện tính toán mà không nhờ đến một công cụ có độ chính xác cao không? Hãy cẩn thận. Ví dụ, giả sử chúng ta sử dụng một loại tính toán powermod? Do đó, tính toán công suất mong muốn, trong khi lấy mô đun ở mọi bước. Vì vậy, thực hiện trong độ chính xác kép:

X = 1; 
for i = 1:1000 
    X = mod(X*4,2*pi); 
end 
sin(X) 
ans = 
     0.955296299215251 

Ah, nhưng hãy nhớ rằng câu trả lời đúng là -0,19033458127208318385994396068455455709388 ...

Vì vậy, về cơ bản là không có gì ý nghĩa còn lại. Chúng tôi đã mất tất cả thông tin trong tính toán đó. Như tôi đã nói, điều quan trọng là phải cẩn thận.

Điều gì đã xảy ra sau mỗi bước trong vòng lặp đó, chúng tôi phát sinh một sự mất mát nhỏ trong tính toán mô đun. Nhưng sau đó chúng tôi nhân câu trả lời bằng 4, điều này làm cho lỗi tăng thêm 4 lần, và sau đó là một hệ số 4, v.v. Và tất nhiên, sau mỗi bước, kết quả sẽ mất một chút nhỏ ở cuối số . Kết quả cuối cùng là crapola hoàn chỉnh.

Cho phép xem xét hoạt động cho một sức mạnh nhỏ hơn, chỉ để thuyết phục chính mình những gì đã xảy ra. Ở đây ví dụ, hãy thử sức mạnh thứ 20. Sử dụng độ chính xác gấp đôi,

mod(4^20,2*pi) 
ans = 
      3.55938555711037 

Bây giờ, hãy sử dụng vòng lặp trong tính toán powermod, lấy mod sau mỗi bước. Về cơ bản, điều này loại bỏ bội số của 2 * pi sau mỗi bước.

X = 1; 
for i = 1:20 
    X = mod(X*4,2*pi); 
end 
X 
X = 
      3.55938555711037 

Nhưng đó có phải là giá trị chính xác không? Một lần nữa, tôi sẽ sử dụng hpf để tính toán giá trị chính xác, hiển thị 20 chữ số đầu tiên của số đó. (Kể từ khi tôi đã thực hiện việc tính toán trong tổng số 50 chữ số, tôi hoàn toàn tin tưởng sẽ 20 đầu tiên của họ.)

mod(hpf(4,[20,30])^20,2*hpf('pi',[20,30])) 
ans = 
      3.5593426962577983146 

Trong thực tế, trong khi các kết quả trong độ chính xác kép đồng ý với các chữ số cuối cùng cho thấy, những đôi kết quả cả hai đều thực sự sai lầm qua chữ số thứ 5 đáng kể. Khi nó quay ra, chúng tôi VẪN cần phải mang hơn 600 chữ số chính xác cho vòng lặp này để tạo ra kết quả của bất kỳ ý nghĩa nào.

Cuối cùng, để tiêu diệt hoàn toàn con ngựa chết này, chúng tôi có thể hỏi liệu có thể thực hiện tính toán powermod tốt hơn không. Nghĩa là, chúng ta biết rằng 1000 có thể được chia ra thành một dạng nhị phân (sử dụng DEC2BIN) như:

512 + 256 + 128 + 64 + 32 + 8 
ans = 
     1000 

Chúng ta có thể sử dụng một chương trình bình phương lặp đi lặp lại để mở rộng quyền lực lớn với phép nhân ít hơn, và do đó gây ra lỗi ít tích lũy? Về cơ bản, chúng tôi có thể cố gắng tính

4^1000 = 4^8 * 4^32 * 4^64 * 4^128 * 4^256 * 4^512 

Tuy nhiên, làm điều này bằng cách lặp lại bình phương 4, sau đó lấy mod sau mỗi thao tác.Tuy nhiên, điều này thất bại, vì thao tác modulo sẽ chỉ loại bỏ bội số nguyên của 2 * pi. Sau khi tất cả, mod thực sự được thiết kế để làm việc trên các số nguyên. Vì vậy, nhìn vào những gì sẽ xảy ra. Chúng tôi có thể diễn tả 4^2 là:

4^2 = 16 = 3.43362938564083 + 2*(2*pi) 

Chúng ta chỉ có thể bỏ phần còn lại, sau đó lấy lại mod không? KHÔNG!

mod(3.43362938564083^2,2*pi) 
ans = 
      5.50662545075664 

mod(4^4,2*pi) 
ans = 
      4.67258771281655 

Chúng ta có thể hiểu những gì đã xảy ra khi chúng tôi mở rộng hình thức này:

4^4 = (4^2)^2 = (3.43362938564083 + 2*(2*pi))^2 

bạn sẽ nhận được gì khi bạn loại bỏ bội INTEGER của 2 * pi? Bạn cần phải hiểu lý do tại sao vòng lặp trực tiếp cho phép tôi loại bỏ bội số nguyên của 2 * pi, nhưng hoạt động bình phương trên không. Tất nhiên, vòng lặp trực tiếp thất bại quá vì các vấn đề số.

+0

Nó không phải là một giải pháp tuyệt vời. Tuy nhiên nó hoạt động và hiệu suất là ok. Cảm ơn bạn. – Bene

+0

Vấn đề là, tính toán mod cho một hàm tuần hoàn là một kẻ giết người ở đây. Mods cho số nguyên là dễ dàng, vì không có mất mát. Nhưng mod đối với 2pi không phải là dễ dàng như vậy. Một trong những mục tiêu của tôi khi tôi viết HPF là tính toán tội lỗi (1e400). double (sin (hpf ('1e400', 1000))) = - 0.998538231983098. Giá trị đó là chính xác. –

+0

Có vẻ như tôi đã thực sự tính toán 'sin (4^4) ~ -0.9992'. 'tội lỗi (4^1000) ~ -0.19033', tôi nghĩ vậy. – DSM

5

Lần đầu tiên tôi xác định lại câu hỏi như sau: tính 4^1000 modulo 2pi. Vì vậy, chúng tôi đã phân chia vấn đề thành hai.

Sử dụng một số thủ đoạn gian trá toán:

(a+2pi*K)*(b+2piL) = ab + 2pi*(garbage)

Do đó, bạn có thể chỉ nhân 4 nhiều lần bởi chính nó và tính toán mod 2pi mọi giai đoạn. Tất nhiên, câu hỏi thực sự đặt ra là độ chính xác của thứ này là gì. Điều này cần phân tích toán học cẩn thận. Nó có thể hoặc không thể là một crap tổng.

+0

Tôi có linh cảm rằng phương pháp này, hoặc bất kỳ phương pháp nào sẽ dẫn đến rác do phép tính số học chính xác gấp đôi: ví dụ: tính toán sin (2 * pi * 10^i) cho i = 1: 100 cho junk, và tương tự với mod máy tính (2 * pi * 10^i + 0,0001, 2 * pi) – db1234

+0

có .. có thể. Người ta cần phải cẩn thận với arithmetics máy. –

+1

Tôi tìm thấy một chức năng mod cho quyền hạn và số lượng lớn: http://www.mathworks.com/matlabcentral/fileexchange/7908-big-modulo-function/content/bigmod.m Đây có thể là giải pháp làm việc. – Bene

3

Theo gợi ý của Pavel với mod tôi đã tìm thấy một hàm mod cho các cường độ cao trên mathwors.com. bigmod(number,power,modulo) KHÔNG thể tính 4^4000 mod 2π. Bởi vì nó chỉ làm việc với các số nguyên như modulo chứ không phải với số thập phân.

Tuyên bố này không chính xác nữa: sin(4^x)sin(bidmod(4,x,2*pi)).

+0

Tôi không tin điều này. Trên Octave (phiên bản mã nguồn mở của Matlab), tôi cố gắng sử dụng thường trình bigmod.m bằng cách gõ 'sin (2 * pi * bigmod (4, 100, 2 * pi))' và kết quả không bằng không. Tôi khá chắc chắn rằng đây cũng là trường hợp của Matlab. – db1234

+1

@dblazevski Trừ khi tôi bị nhầm lẫn nặng, nó không phải là số không ... '4^100 mod 2pi' không phải là số nguyên. – Dougal

+0

Opps ... sai lầm ngu ngốc, nhờ @Dougal, tôi đoán một thử nghiệm sẽ là 'tội lỗi (bigmod (4 * 2 * pi, 100, 2 * pi))', và điều đó cho số không vì nó là nghĩa vụ phải làm – db1234

Các vấn đề liên quan