Argumentreduktion

Argumentreduktion går ut på att reducera argumentet x så att det ligger i det lilla intervallet [0, pi/4]. Detta är nu inte så lätt som det låter. Antag att vi vill reducera ett argument av storleksordningen 1e100. Vi vill då bestämma heltalet n så att y = x - 2 n pi och där 0<= y < 2 pi. Pi kan vi givetvis inte lagra exakt utan vi lagrar 2 pi - f, där f är ett avrundningsfel. Vårt approximativa y blir med andra ord y = x - n (2 pi - f) så att det absoluta felet i y blir n f, där n har storleksordningen x / (2 pi). Det är egentligen det relativa felet vi vill komma åt, dvs n f / y. Nu är det inte alldeles lätt att få en uppfattning om hur stort y blir; man kan givetvis räkna på det värsta fallet, när y är det minsta maskintalet större än noll. Låt oss dock, för att vara konkreta, anta att y är ungefär lika med ett. Vi vill räkna ut y med full noggrannhet i dubbel precision och tillåter därför inte ett större relativt fel än 1e-16. Så

| n f | < 1e-16, så att | x f / (2 pi) | < 1e-16 och | f | < 1e-16  2 pi / 1e100 = 6.28.. e-116

(eftersom vi antog att x var 1e100). Vi måste alltså lagra (och räkna med) 2 pi med mer än 116 siffror! Ännu värre blir det om vi tar det största maskintalet och har oturen att y blir mycket mindre än ett.

Sun har en mycket bra argumentreduktionsrutin och kan beräkna sinus, cosinus med full noggrannhet för alla maskinrepresenterbara tal. Detta ger dock upphov till vissa överraskningar. Säg att vi vill beräkna sin(1e22), en Sundator levererar då svaret -0.85220084976719 vilket är ett korrekt avrundat svar. Testar vi sedan med sin(1e23) får vi resultatet -0.32405393764300. Detta är också ett korrekt avrundat värde (alla utskrivna siffror är korrekta). Problemet är bara att sin(1e23) = 0.701140639861078469469241830480871... Lösningen på denna paradox är att 1e23 inte kan representeras exakt i datorn, så sinusrutinen räknar alltså korrekt men på ett helt annat argument än 1e23. Om man tar reda på bitrepresentationen av det tal som lagrats som approximation av 1e23 ser man att

   teckenbiten är s = 0,
   exponenten är e = 10001001011 och
   decimaldelen är d = 0101001011010000001011000111111000010100101011110110.

Talet är då (med bias = 1023 decimalt)

2^s * 2^(e-bias) * 1.d

som är 99999999999999991611392 decimalt.

Sinus för detta tal är -0.32405393764300330176... vilket stämmer utmärkt med det värde Sun levererar. Anledningen till att 1e22 inte ger några problem är att detta är den (största) tiopotens som kan representeras exakt i dubblel precision. Sunen klarar av mycket större argument. sin(1e308) ger ett korrekt avrundat värde av sinus av

100000000000000001097906362944045541740492309677311846336810682903157585
404911491537163328978494688899061249669721172515611590283743140088328307
009198146046031271664502933027185697489699588559043338384466165001178426
897626212945177628091195786707458122783970171784415105291802893207873272
974885715430223118336

Om vi tittar på Suns sinusrutin så börjar den med  driver-rutinen s_sin.c. Denna rutin avgör om argumenreduktion är nödvändig och anropar i så fall __ieee754_rem_pio2 och k_rem_pio2.c. När väl argumentet x har reducerats till y som ligger i intervallet [0, pi/4] beräknas sin(y) med hjälp av en serieutveckling i k_sin.c. (Det är inte konstigt att det tar mycket mera tid att beräkna sinus än att göra en enkel multiplikation.)

I argumentreduktionsrutinen har 2/pi lagrats med 476 decimaler. Eftersom vi vill återföra argumentet till första kvadranten subtraherar vi multipler av pi/2 (i stället för multipler av 2 pi). Om y är det reducerade värdet så gäller y = x - n (pi / 2) eller ekvivalent (2/pi)y = (2/pi)x - n (det är här 2/pi kommer in); vi går inte in på några detaljer eftersom argumenreduktionsrutinen är tämligen knepig.

Vad gör man med andra, ickeperiodiska, funktioner då?

Låt oss ta ett exempel, nämligen den naturliga logaritmem, log(x). x lagras som 2^exponent (1 + f) där f är den binära "decimaldelen". log(x) är då exponent * log(2) + log(1 + f) och "det enda" vi behöver approximera är logaritmen på ett litet intervall. I själva verket är man lite listigare.