Säkra summor i Fortran 95

För (rätt många) år sedan var det populärt med så kallad intervallaritmetik. Ett flyttal representeras då som ett intervall med hjälp av två flyttal. Genom att systematiskt avrunda nedåt respektive uppåt kan vi se till att ett beräknat resultat hela tiden ligger i intervallet. IEEE-standarden har stöd för sådana avrundningar vilket göra att det inte är så svårt att implementera intervallaritmetik. Det finns flera nackdelar, som förklarar varför intevallaritmetik aldrig slog igenom. Det krävs ju dubbelt så mycket minne och dubbelt så många räkneoperationer. Om man inte har speciella algoritmer så kan intervallen bli väldigt stora. Det är t.ex. inte så intressant att få veta att ett resultat som borde vara omkring 10 ligger i intervallet [-10000, 20000] till exempel.
 
Här följer några exempel. Jag har utnyttjat Suns Fortran95-kompilator som har stöd för intervallaritmetik (det är dock inte en del av standarden för Fortran95). 

% cat interval1.f95
program interval_ex
  interval :: a, b   ! a och b är två variabler av typen interval

  a = [2.0, 3.0]     ! två intervallkonstanter
  b = [-1.0, 3.0]

  print*, a + b, a - b, a * b, a / b, b / a  ! räkna lite

  a = [0.1]  ! kort för [0.1, 0.1]

  print*, a
  print*, a * 10

end program interval_ex

% f95 -xia interval1.f95      xia slår på intervallaritmetiken

% a.out
 [1.0,6.0] [-1.0,4.0] [-3.0,9.0] [-Inf,Inf] [-0.5,1.5]
 [0.099999999999999991,0.10000000000000001]            OBS: avrundat åt två håll
 [0.99999999999999988,1.0000000000000003]

Man kan notera att a + b etc. fungerar som man tror. Notera att a / b ger division med noll eftersom nollan tillhör intervallet b. I dessa enkla exempel räknar vi exakt, men det är inte fallet med a = [0.1]. Vi vet att  0.1 inte kan lagras exakt i enkel- eller dubbel precision så efter tilldelningen kommer a att innehålla 0.1 avrundat nedåt och uppåt.

Här följer ett summationsprogram:
  

% cat summor.f95
program summor
  implicit none

  integer             :: k
  interval (kind = 4) :: summa  ! summa är en intervallvariabel i enkel prec.

  summa = [0.0]
  do k = 1, 1000000
    summa = summa + [1.0] / k
  end do

  print*, "1 + 1/2 + ... + 1/1000000 = ", summa, wid(summa)  ! wid = intervallets längd

  summa = [0.0]
  do k = 1000000, 1, -1
    summa = summa + [1.0] / k
  end do
  
  print*, "1/1000000 + ... + 1/2 + 1 = ", summa, wid(summa)

end program summor

% f95 -xia summor.f95            

% a.out
 1 + 1/2 + ... + 1/1000000 =  [13.971489,14.924348] 0.952858
 1/1000000 + ... + 1/2 + 1 =  [14.350338,14.436089] 0.085749626

Notera bredden på intervallen! Det finns några andra intervallfunktioner förutom wid. inf ger intervallets minsta värde och sup dess största. mid ger mittpunkten.

För att se hur väl intervallaritmetik fungerar för en riktig tillämpning så konverterade jag ett program för lösandet av linjära ekvationssystem. Det är en rätt enkel process. Alla typdeklarationer för flyttalstyperna får konverteras till motsvarande intervalltyp. En del flyttalskonstanter får också skrivas om. Det enda problem jag hade var att konvertera jämförelseoperationen vid den partiella pivoteringen.

En if-sats av typen if ( a < b ) then ... kommer att ge kompileringsfel om a och b är intervall. Detta är ju inte så konstigt ty vad menar vi med "mindre än" om t.ex. t.ex. a = [1, 3] och b = [2, 4]?

Sun har implementerat tre varianter av "mindre än":

Relationsoperatorerna har namnen .clt. (certainly less than), .plt och .slt..

I en test med A x = b där A hade dimension 500 och konditionstal 172 så varierade vidden på x-intervallen mellan 10-5 och 1.69 fastän maxfelet i lösningen var så litet som 6 10-15. Beräkningen gjordes i dubbel precision.  Även om det kan vara enkelt att konvertera programmen så är således resultatet inte alltid så användbart.

Om man inte har tillgång till Suns Fortran 95 så kan man implementera intervallaritmetik på egen hand. Det förenklas givetvis om man har språk som stödjer operatoröverlagring, som C++ och Fortran (i viss utsträckning), eftersom syntaxen blir så mycket trevligare och då det gör konverterandet av gamla program rätt så enkel. Java stödjer inte överlagring, vilket är en stor brist i fall som detta.

Här följer ett exempel där jag har beräknat våra två summor i enkel precision och med de fyra avrundningsmoderna. Dessa är avrundning till närmaste, avrundning mot +oändligheten (alltid uppåt), avrundning mot minus oändligheten samt avrundning mot noll (nedåt för positiva tal och uppåt för negativa). I tabellen nedan har vi summorna s1 och s2 (som i labben) och avrundningsriktningarna (positive = +oändligheten, negative = -oändligheten). Programmet, round.f90, är skrivet i Fortran 90 och utnyttjar inte överlagring. error är skillnaden mellan enkelprecisionsvärdena och det exakta (avrundat till dubbel precision). 

type            s1        error      s2        error
nearest  :  14.3573580  -0.0354  14.3926516  -0.0001
tozero   :  13.9714899  -0.4212  14.3503389  -0.0423
positive :  14.9243479   0.5317  14.4360886   0.0434
negative :  13.9714899  -0.4213  14.3503389  -0.0424

Notera att avrundning till närmaste ger ett mycket bättre resultat. Jämför också med summationsprogrammet ovan.

Slutligen ett litet problem. Kan Du förklara följande? 

% cat interval3.f95 
program sin_ex

  print*, wid([1e20]), sin([1e20]), wid(sin([1e20]))
  print*, wid([1e21]), sin([1e21]), wid(sin([1e21]))
  print*, wid([1e22]), sin([1e22]), wid(sin([1e22]))
  print*, wid([1e23]), sin([1e23]), wid(sin([1e23]))

end program sin_ex

% f95 -xia interval3.f95
% a.out
 0.0E+0 [-0.64525128526578091,-0.64525128526578079] 1.1102230246251565E-16
 0.0E+0 [-0.66712017707180494,-0.66712017707180482] 1.1102230246251565E-16
 0.0E+0 [-0.85220084976718891,-0.85220084976718879] 1.1102230246251565E-16
 16777216.0 [-1.0,1.0] 2.0

Back