Calcolo 4^x mod 2π per grande x

Question

Ho bisogno di calcolare sin(4^x) con x> 1000 in Matlab, con è fondamentalmente sin(4^x mod 2π) Poiché i valori all'interno della funzione sin diventano molto grandi, Matlab restituisce infinito per 4^1000. Come posso calcolare in modo efficiente questo? Preferisco evitare i grandi tipi di dati.

Penso che una trasformazione in qualcosa come sin(n*π+z) potrebbe essere una possibile soluzione.

Answer 1

È necessario fare attenzione, poiché ci sarà una perdita di precisione. La funzione sin è periodica, ma 4^1000 è un grande numero. Così efficacemente, sottraiamo un multiplo di 2 * pi per spostare l'argomento nell'intervallo [0,2 * pi).

4^1000 è circa 1e600, un numero veramente grande. Quindi farò i miei calcoli usando il mio high precision floating point tool in MATLAB. (In effetti, uno dei miei obiettivi espliciti quando ho scritto HPF era di essere in grado di calcolare un numero come sin (1e400). Anche se stai facendo qualcosa per il gusto di farlo, farlo correttamente ha ancora senso.) In questo caso , poiché so che il potere a cui siamo interessati è all'incirca 1e600, quindi eseguirò i miei calcoli in più di 600 cifre di precisione, aspettandomi che perderò 600 cifre per la cancellazione sottrattiva. Questo è un enorme problema di cancellazione sottrattiva. Pensaci. Quella operazione modulo è effettivamente una differenza tra due numeri che saranno identici per le prime 600 cifre o giù di lì!

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

Qual è il multiplo più vicino di 2 * pi che non superi questo numero? Possiamo ottenerlo con una semplice operazione.

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

Come si può vedere, le prime 600 cifre erano uguali. Ora, quando sottraiamo i due numeri,

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

Questo è il motivo per cui ho fatto riferimento ad esso come un enorme problema di cancellazione sottrattiva. I due numeri erano identici per molte cifre. Persino trasportando 1000 cifre di precisione, abbiamo perso molte cifre. Quando si sottraggono i due numeri, anche se stiamo trasportando un risultato con 1000 cifre, solo le cifre 400 dell'ordine più alto sono ora significative.

HPF è in grado di calcolare naturalmente la funzione trigonometrica. Ma come abbiamo mostrato sopra, dovremmo fidarci solo approssimativamente delle prime 400 cifre del risultato. (In alcuni problemi, la forma locale della funzione sin potrebbe farci perdere più cifre di quella.)

sin(X^1000) 
ans = 
-0.19033458127208318385994396068455455709388374041098639172943768418947125138650234240955423917696880832346734715448603532912993423621761996537053192685449334064870714463489747336279464911185192423229252660143128976923388511299599457104070322693060218958487584842139143972048735807765826659851362293280012583640059277583434162223469640779539703355744143419935430600390820454055891750089781440474478225522286222463738277009002753247363724815609283394633443329778920087022201603354152914210817007440447838392869577354385645124650950464218066771029610934877080889086985319804240164585346291661088530125354930225403524397401167317843031900829546691402971929428720760150282604082313216048252703439459284455892236101855653841958635139010896628829034919565066139672417258772760228631878006327065033172

Quindi ho ragione, e non possiamo fidarci di tutte queste cifre? Farò lo stesso calcolo, una volta in 1000 cifre di precisione, quindi una seconda volta in 2000 cifre. Calcola la differenza assoluta, quindi prendi il log10. Il risultato di 2000 cifre sarà il nostro riferimento come essenzialmente esatto rispetto al risultato di 1000 cifre.

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

Ah. Quindi, con le 1000 cifre di precisione con cui abbiamo iniziato, abbiamo perso 602 cifre. Le ultime 602 cifre nel risultato sono diverse da zero, ma continuano a completare la spazzatura. Questo era come mi aspettavo. Solo perché il tuo computer riporta alta precisione, è necessario sapere quando non fidarsi di esso.

Possiamo eseguire il calcolo senza ricorrere a uno strumento di alta precisione? Stai attento. Ad esempio, supponiamo di usare un tipo di calcolo powermod? Quindi, calcola la potenza desiderata, mentre prendi il modulo ad ogni passo. Così, fatto in doppia precisione:

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

Ah, ma ricordate che la vera risposta è stata -0,19033458127208318385994396068455455709388 ...

Quindi non v'è essenzialmente nulla di significativo rimanente. Abbiamo perso tutte le nostre informazioni in quel calcolo. Come ho detto, è importante stare attenti.

Quello che accadeva era dopo ogni fase di quel ciclo, abbiamo subito una piccola perdita nel calcolo del modulo. Ma poi abbiamo moltiplicato la risposta per 4, che ha fatto crescere l'errore di un fattore 4, e poi un altro fattore di 4, ecc. E, naturalmente, dopo ogni passo, il risultato perde un pochino alla fine del numero . Il risultato finale era completo crapola.

Osserviamo l'operazione per una potenza inferiore, solo per convincerci di quello che è successo. Qui per esempio, prova la ventesima potenza. Usando la doppia precisione,

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

Ora, utilizzare un ciclo in un calcolo powermod, prendendo il mod dopo ogni passaggio. In sostanza, questo scarta multipli di 2 * pi dopo ogni passo.

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

Ma è il valore corretto? Ancora una volta, userò hpf per calcolare il valore corretto, mostrando le prime 20 cifre di quel numero. (Da quando ho fatto il calcolo di 50 cifre totali, sarò assolutamente fiducia il primo 20 di loro.)

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

Infatti, mentre i risultati in doppia precisione si impegnano a l'ultima cifra indicata, quelli doppio i risultati erano entrambi errati dopo la 5a cifra significativa. A quanto pare, ANCORA abbiamo bisogno di trasportare più di 600 cifre di precisione per questo loop per produrre un risultato di qualsiasi significato.

Infine, per uccidere completamente questo cavallo morto, potremmo chiedere se è possibile eseguire un calcolo di powermod migliore. Cioè, sappiamo che 1000 può essere scomposto in una forma binaria (uso DEC2BIN) come:

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

Possiamo usare uno schema di quadratura ripetuto per espandere la grande potenza con un minor numero di moltiplicazioni, e quindi causare errori accumulati? In sostanza, potremmo provare a calcolare

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

Tuttavia, fare questo ripetutamente quadrando 4, quindi prendendo il mod dopo ogni operazione.Ciò tuttavia non riesce, poiché l'operazione modulo rimuoverà solo i multipli interi di 2 * pi. Dopotutto, la mod è davvero progettata per funzionare su interi. Quindi guarda cosa succede. Possiamo esprimere 4^2 come:

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

Possiamo semplicemente quadrare il resto, quindi riprendere il mod? NO!

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

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

Siamo in grado di capire che cosa è successo quando espandiamo questa forma:

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

Che cosa si ottiene quando si rimuovono multipli interi di 2 * pi? Hai bisogno di capire perché il ciclo diretto mi ha permesso di rimuovere multipli interi di 2 * pi, ma l'operazione di squadratura sopra non lo fa. Ovviamente anche il ciclo diretto non è riuscito a causa di problemi numerici.

Answer 2

Vorrei prima ridefinire la domanda come segue: calcolare 4^1000 modulo 2pi. Quindi abbiamo diviso il problema in due.

utilizzare alcuni matematica inganno:

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

Quindi, si può solo moltiplicare 4 tante volte da solo e informatica 2pi mod ogni fase. La vera domanda da porsi, ovviamente, è qual è la precisione di questa cosa. Ciò richiede un'attenta analisi matematica. Potrebbe essere o non essere una schifezza totale.

Answer 3

In seguito al suggerimento di Pavel con mod ho trovato una funzione mod per le alte potenze su mathwors.com. bigmod(number,power,modulo) NON può calcolare 4^4000 mod 2π. Perché funziona solo con numeri interi come modulo e non con decimali.

Questa dichiarazione non è più corretta: sin(4^x) è sin(bidmod(4,x,2*pi)).