GIER/GAIII DPA double ln() fejlsøgning

Indledning

I januar 2019 spurgte Mogens mig, om jeg måske kendte noget til Chebyshev polynomier og approksimation af matematiske funktioner. Under arbejdet med programmerne i dokumentet Doppelte Präzision Arithmetik in GIER ALGOL (DPA dokumentet) havde Mogens opdaget, at programmet til at afprøve funktionen double ln() til at beregne den naturlige logaritme i dobbelt præcision, kørt på Mogens' simulator med Mogens' indtastede HJÆLP, GAIII og DPA-programmer, ikke gav samme resultater som angivet i DPA dokumentet. Så spørgsmålet var, kunne man finde en forklaring på dette fænomen?

Og det var jo en spændende opgave: Jeg er ikke helt uden relevant baggrund: Numerisk analyse fra Aarhus Universitet 1976-77 med Ole Østerby, lærebog Ralston: "A First Course in Numerical Analysis", 1965, og studier af (utvivlsomt) Peter Naurs implementation af matematiske standard funktioner i GA4, f.eks. arctan() og exp(), baggrund bl.a. Lance: "Numerical Methods for High Speed Computers", 1960.

Undersøgelsen fulgte to spor: Dels analyser af det foreliggende materiale, koden til double ln() og de rapporterede testresultater. Dels (gen-)oparbejdelse af forståelsen af de forskellige metoder til approksimation af funktioner, hvor bl.a. de omtalte Chebyshev polynomier ofte spiller en rolle.

Det blev analysen af koden, herunder en selvstændig dublering (til en god tilnærmelse, i hvert fald) af double ln() i Haskell, der endte med at afsløre synderen: Der mangler et ciffer, 6, i een af koefficienterne til det 25.-grads polynomium, der indgår i beregningen af logaritmefunktionen: I DPA dokumentet ser linien i koden til double ln() funktionen efter den, der bærer kommentaren "; a2", således ud, her gengivet sammen med de to foregående linier, de to efterfølgende linier og en passende markør:

     279/1016/ 744/ 971
    1023/1002/ 262/  25        ; a2
     358/  69/ 715/  6                 <- Fejl i denne linie
       0/   2/ 498/ 844        ; a3
     351/ 747/ 646/  94

Erstatter man den markerede linie med

     358/  69/ 715/  66

stemmer resultaterne af testkørslerne på Mogens' simulator med resultaterne fra DPA dokumentet.

I det følgende gives en kort introduktion til logaritmer og koden til double ln() præsenteres. Derefter opsummeres det forløb, der endte med at afsløre fejlen, uden så mange dikkedarer. Endelig følger, hæmningsløst, en masse detaljer om undersøgelserne.

Og lige for en ordens skyld: Der er altså tale om en fejl i den udskrift af koden, som er gengivet i DPA dokumentet. Det er ikke Mogens, der har tastet noget forkert ind.

Logaritmer

Hvis x = a^y, hvor ^ er potensopløftning, er y logaritmen med grundtal a til x. Jeg vil skrive y = loga(x), så x altid er = a^loga(x) for positive x. For eksempel er 2^10 = 1024, så log2(1024) = 10. Et andet eksempel: log10(100) = 2 fordi 100 = 10^2. Specielt er loga(1) = 0 for alle grundtal a.

Også for negative eksponenter giver logaritmen fin mening: For eksempel er 0.01 = 10^(-2), så log10(0.01) = -2. (Jeg bruger punktum som decimal separator.)

Oprindelig er logaritmer opfundet til hjælp ved beregninger og den helt grundlæggende formel, der kan bruges til at erstatte multiplikation med addition, er:

x*y = 10^(log10(x)+log10(y))

Man kan altså regne x*y ud ved at tage log10 til x og log10 til y, addere dem og så opløfte 10 til denne sum. Det lyder jo ikke umiddelbart helt som en besparelse, men med passende tabeller til at slå log10 og 10-potensopløftninger (også kaldet antilogaritmer) op, har det faktisk hjulpet ved håndberegninger: Pierre-Simon de Laplace (1749-1827) skal have udtalt, at Logarithms, by shortening the labors, doubled the life of the astronomer (William Dunham i "Euler: The Master of Us All", med henvisning til Victor Katz, "A History of Mathematics: An Introduction", p. 420).

Formlen for x*y gælder for alle grundtal, altså:

x*y = a^(loga(x) + loga(y))

For at se det, kan vi i første omgang skrive

x*y = a^(loga(x)) * a^(loga(y))

Nu gælder det jo, at

a^p * a^q = a^(p+q)

F.eks. er både 2^3 * 2^2 = 8*4 = 32 og 2^(3+2) = 2^5 = 32. Men det vil så sige, at

x*y = a^(loga(x)) * a^(loga(y)) = a^(loga(x) + loga(y))

og så er den hjemme.

Så længe vi nøjes med hele tal og ikke mindst heltalspotenser, er det ikke så vanskeligt. Men matematikere generaliserer så gerne og hvis vi ser, at log10(100) = 2 og log10(1000) = 3, så vil man jo straks spørge, hvad er mon log10(500), f.eks.? Det er jo nok et tal mellem 2 og 3, men hvilket tal? Som svar vil jeg nøjes med at notere, at det har matematikerne selvfølgelig styr på, så man kan finde en passende logaritme til alle tal større end 0.

Logaritmer med forskellige grundtal er proportionale, således at

loga(x) = logb(x) * loga(b)

Så hvis vi har logb(x), kan vi finde loga(x) ved at gange med loga(b), som jo er en konstant. Det indser vi på denne måde:

x = b ^ (logb(x)) = (a ^ loga(b)) ^ (logb(x))

Nu gælder

(a^p)^q = a^(p*q)

F.eks. er både (2^3)^2 = 8^2 = 64 og 2^(3*2) = 2^6 = 64. Derfor kan vi gå videre med

(a ^ loga(b)) ^ (logb(x)) = a ^ (loga(b) * logb(x))

Så

x = a ^ (loga(b) * logb(x))

Så tager vi loga på begge sider

loga(x) = loga(b) * logb(x)

som (bortset fra faktorernes orden) var det, vi ledte efter.

Der gælder:

loga(b^x) = x * loga(b)

Det ser vi ved at bruge proportionaliteten:

loga(b^x) = logb(b^x) * loga(b) = x * loga(b)

Den naturlige logaritme: loge() eller ln()

Som beskrevet findes logaritmer til forskellige grundtal og to logaritmer med forskellige grundtal adskiller sig blot ved en proportionalitetskonstant. Til beregninger i det velkendte 10-talssystem har log10() som nævnt gjort god tjeneste. Imidlertid har matematikere, mindst siden Euler, kastet særlig kærlighed efter logaritmen med det specielle grundtal, kaldet e, der har værdien sådan ca. 2.718281828459045..., hvor decimalerne fortsætter i en svimlende uendelighed. Og hvorfor nu det?

Tjah, der kan nævnes flere grunde: En af dem er, at den naturlige logaritme, loge(x) eller ln(x) som den også benævnes, kan findes som arealet fra t=1 til t=x under den særligt simple hyperbel f(t) = 1/t. Noget tilsvarende gør sig selvfølgelig gældende for andre logaritmer, men så med lidt proportionalitetskonstanter til at forvirre billedet.

Det skal også nævnes, at potensopløftningsfunktionen e^x, den omvendte funktion af ln(), normalt betegnes exp(), eksponentialfunktionen.

double ln()

Koden for double ln() er fisket ud af Mogens' indtastede DPA-kode, som er den gaIIIdp.flx, man finder på:

http://datamuseum.dk/site_dk/software/scripts/showitem.pl?itemid=319

I Mogens' udgave er fejlen rettet, så jeg har introduceret fejlen igen som den fremgår af DPA dokumentet og tilføjet linienumre:

  1 [double ln(x)]
  2
  3 _b k=c60+e70,i=10,a6,b10
  4
  5 [track 1]
  6
  7 _d b=40i
  8
  9       ncn (c47)     t    2       ; _i_f number of parameters|=1
 10       ps   15      ,hh   c55     ; _t_h_e_n alarm(|<<param|>)
 11       ps  (c46)    ,pm   s1      ; M:=description(x)
 12       grn  1c37    ,hs   c22     ; UV1:=0; _g_o_t_o take formal
 13       arnf(c36)    ,gm   c47     ; RF:=x[0]; c47:=exponent(x)
 14       hv   ra3      LZ           ; _i_f x=0 _t_h_e_n _g_o_t_o error
 15       pm  (c36)     t    1 V NT  ; _i_f x>0 _t_h_e_n M:=x[1]; skip next
 16 a3:   ps   5       ,hh   c55     ; error: alarm(|<<ln|>)
 17       hs   c65     ,qq   1c60.29+1 ; take constant track
 18       sr   c40     ,tl   1       ; t:=(numerus(x)-1)*2
 19       sr   c40     ,tl   11      ;     -1
 20       gr   c37     ,gm   1c37
 21       grn  c68     ,gr   c54     ; b1:=0; _c_o_m_m_e_n_t b1: c68,c54
 22       pm   sb5     ,gs   ra2     ; b:=a[n]; _c_o_m_m_e_n_t b: R,M
 23 a1:   pp   b5      ,pp   p-1     ; p:=number of constant cells
 24                                  ; again: p:=p-1
 25       gm   c36     ,pm   c68     ; save b[1]
 26       gm   ra5     ,pm   c54     ; b2:=b1; _c_o_m_m_e_n_t b2: ra5,sb9
 27       gm   sb9     ,pm   c36
 28       gr   c68     ,gm   c54     ; b1:=b
 29       mkn  c37     ,pm   c68     ; b:=b1*t
 30       mk   1c37    ,ml   c37
 31       tl   1       ,sr   ra5     ;    *2 - b2
 32 a4:   gr   c36      X
 33 a2:   ar   p       ,sr   sb9     ;    + a[p]
 34       tl   -39     ,ar   c36
 35       bs   pb8     ,hv   ra6
 36       pp   p-1     ,ar  (ra2)    ; _i_f double cell constant
 37 a6:   bs   p       ,hh   ra1     ; _i_f p>0 _t_h_e_n _g_o_t_o again
 38       sr   ra5     ,ud   ra4     ; b:=(b-b2)/2
 39       sr   sb9     ,tl   -39
 40       ar   c36     ,ar   c47     ; double ln:=(b+exponent(x))
 41       gr   c68     ,grn  c36
 42       mkn  sb6     ,pm   c68     ;             * ln(2)
 43       mk   sb7     ,ml   sb6
 44       tl   1
 45       nl   ra      ,tl   -10     ; normalize; shift to fl.p. position
 46       gm   1c37    ,ps  (c46)
 47 a:    pm            t    8 D     ; exponent(double ln)
 48       grf  c37     ,hv   c51     ; exit
 49 a5:   0                          ; working location b2[0]
 50
 51   _b k=e90,i=0
 52     _d i=e89
 53     qq 144.9+c60.19              ; pass 6 word
 54     _tdoubleln;
 55     _d e71=e71+e68, e72=e72+e68, e89=e89+e68+e68+e68
 56   _e
 57
 58 _d c60=c60+e68+e68
 59 #072
 60 [double ln(x) track 2]
 61
 62 b:       0/ 556/ 144/ 566        ; constants*2|&(-10): a0
 63         53/ 381/ 566/ 463
 64          0/ 253/ 479/ 229        ; a1
 65        279/1016/ 744/ 971
 66       1023/1002/ 262/  25        ; a2
 67        358/  69/ 715/  6
 68          0/   2/ 498/ 844        ; a3
 69        351/ 747/ 646/  94
 70       1023/1023/ 696/ 282        ; a4
 71         30/ 184/ 831/  56
 72          0/   0/  44/1006        ; a5
 73        258/ 507/ 538/ 494
 74       1023/1023/1017/ 582        ; a6
 75         52/ 752/ 508/ 222
 76          0/   0/   0/ 968        ; a7
 77        275/ 594/ 951/ 406
 78       1023/1023/1023/ 878        ; a8
 79        305/ 664/ 256/   8
 80          0/   0/   0/  22        ; a9
 81         89/ 773/ 283/1007
 82       1023/1023/1023/1020        ; a10
 83 b10:   294/ 822/ 514/ 497
 84        273/ 466/ 405/ 750        ; a11
 85        980/1016/  70/ 409        ; a12
 86          6/ 830/ 842/ 171        ; a13
 87       1022/ 936/ 806/ 455        ; a14
 88          0/ 177/ 966/ 310        ; a15
 89       1023/ 995/ 386/ 929        ; a16
 90          0/   4/ 636/ 863        ; a17
 91       1023/1023/ 257/  87        ; a18
 92          0/   0/ 124/ 672        ; a19
 93       1023/1023/1003/ 698        ; a20
 94          0/   0/   3/ 328        ; a21
 95       1023/1023/1023/ 467        ; a22
 96          0/   0/   0/  91        ; a23
 97       1023/1023/1023/1009        ; a24
 98 b1:      0/   0/   0/   3        ; a25
 99 b2:    177/ 456/ 383/ 500        ; ln(2)
100        231/ 755/ 350/ 316
101 b4:   0                          ; working location b2[1]
102
103 _d b3=1b2, b5=b1-b, b6=b2-b, b7=b3-b, b9=b4-b, b8=b-b10
104
105 _e

Og hvordan virker den så? Nogenlunde således:

Vi skal beregne ln(x), den naturlige logaritme til x, logaritmen med grundtallet e = 2.718281828459045...
double ln() beregner ln(x) som log2(x) * ln(2), altså som logaritmen med grundtallet 2 til x, ganget med konstanten ln(2) = 0.6931471805599453 ... Her bruges proportionaliteten beskrevet tidligere.
For at beregne log2(x), udnytter double ln(), at tallet x, som er et flydende tal med dobbelt længde, er repræsenteret som y*2^n, hvor 1 <= y < 2 og n er et heltal. (x antages > 0. Hvis x <= 0, melder double ln() en fejl.) Så log2(x) = log2(y*2^n) = log2(y) + log2(2^n) = log(y) + n, i henhold til formlerne beskrevet tidligere. y er det tal, der betegnes med numerus(x) i kommentaren til linie 18 i koden.
For at beregne log2(y), hvor nu 1 <= y < 2, transformerer double ln() først dette interval mellem 1 og 2 til intervallet mellem -1 og 1. Præcist t = 2*y - 3, hvorved -1 <= t < 1 (linie 18-19 i koden).
log2(y) approksimeres nu med et 25. grads polynomium i t, koefficienter a0, a1, ..., a25 (linie 62-98) indgår, resultat til variablen b (linie 29-33). Det skal bemærkes, at a0, ..., a25 ikke direkte er koefficienterne til polynomiet, men de indgår, via noget skalering: Se kommentarerne til linierne 29-33, hvis man da ikke direkte kan gennemskue koden.
Endelig beregnes double ln(x) som (b + n [2-tals eksponenten til x]) * ln(2) (linie 40-43).

Mogens' materiale

I tilknytning til DPA dokumentet, havde Mogens transskriberet double ln() til en form, der kunne accepteres af GA4's maskinkodeoversætter (passage 9), som ikke understøtter de slip'ske /-konstanter. Derved havde han fået en udgave, der kunne køres på vores fysiske GIER, så han kunne konstatere, at forskellen ikke skyldtes en forskel i hans simulator.

Her et udsnit af Mogens' program:

_b_e_g_i_n
   _p_r_o_c_e_d_u_r_e double ln(a,a1,b,b1);
   _v_a_l_u_e b,b1;
   _r_e_a_l a,a1,b,b1;
   _b_e_g_i_n
      _c_o_m_m_e_n_t a:=ln(b);
      _b_o_o_l_e_a_n bln;
      _c_o_r_e _c_o_d_e bln,a,a1,b,b1;
      ... (transskriberet double ln() fra DPA dokumentet)
      _e;
      gier(bln)
   _e_n_d double ln;
   _p_r_o_c_e_d_u_r_e double write(n,a,a1);
   _v_a_l_u_e n,a,a1;
   _i_n_t_e_g_e_r n;
   _r_e_a_l a,a1;
   _b_e_g_i_n
      _b_o_o_l_e_a_n bwrite;
      _c_o_r_e _c_o_d_e bwrite,n,a,a1;
      ... (transskriberet double write() fra DPA dokumentet)
      _e;
      gier(bwrite)
   _e_n_d double write;
   _r_e_a_l a,a1,b,b1;
   select(16);
   b:=1;
   b1:=0;
   a:=a1:=1234;
   double ln(a,a1,b,b1);
   writecr;
   writetext(|<<Bad:  -1.45689666927739658319'- 19|>);
   writecr;
   writetext(|<<Good: -3.38813178901720135627'- 21|>);
   writecr;
   writetext(|<<This: |>);
   double write(20,a,a1);
   writecr;
_e_n_d;

Programmet beregner double ln(1) og udskriver den lakoniske besked:

Bad:  -1.45689666927739658319'- 19
Good: -3.38813178901720135627'- 21
This: -1.45689666927739658319'- 19

I denne udskrift stammer konstanten Bad: -1.45689666927739658319'- 19 fra kørslen af DPA testprogrammet med GAIII systemet på simulatoren mens konstanten Good: -3.38813178901720135627'- 21 er den rapporterede værdi fra DPA dokumentet. Så udskriften skal ses som en indikation af, at den værdi (This:) programmet har beregnet matcher den forkerte værdi af ln(1) (Bad:) og ikke den rapporterede fra DPA dokumentet (Good:).

Fejlen afdækket

En af mulighederne for at forklare den observerede forskel mellem de rapporterede testresultater og de faktiske testresultater er jo, at der er en simpel forskel af en art mellem den kode, der har produceret de rapporterede testresultater og så koden, som Mogens har tastet ind. En mulighed ville f.eks. være, at der var en forskel i een af koefficienterne til det 25. grads polynomium, der bruges til at approksimere log2(y).

Som grundlag for at efterprøve denne mulighed, havde jeg implementeret double ln() som et program i det funktionelle programmeringssprog Haskell. Ikke en fuldstændig eksakt kopi, forskelle i afrunding og lignende kunne stadig gøre sig gældende. Men en rigtig god dublet, der kunne eftergøre beregningerne meget tæt.

Og en forskel i en polynomiekoefficient, hvordan kunne den mon se ud? Tjah, lad os tage den linie, hvor vi ender med at finde fejlen, som eksempel:

 67        358/  69/ 715/  6

Disse koefficienter er jo skrevet op i henhold til slip-notation som en gruppe med 4 decimale tal, der hver ender i 10 bits af det 40-bits tal, som ønskes lagret. Så hver af de 4 tal ligger mellem 0 og 1023. Og vi regner jo i dobbelt præcision, så de store af koefficienterne har to sådanne linier.

Overvejer vi det første tal, 358, kunne nærliggende muligheder for forskelle jo være, at et enkelt af cifrene var anderledes. F.eks. 458 i stedet for 358. Og det har jeg så forsøgt på den forholdsvis grove facon, at alle muligheder, der fremkommer ved at addere eller subtrahere 1, 2, ... 9 i hver position, er blevet prøvet. Så for eksemplet 358 altså

359 360 361 362 363 364 365 366 367

ved at addere 1, 2, ... 9 i enernes position og

357 356 355 354 353 352 351 350 349

ved at subtrahere 1, 2, ... 9 i enernes position. Og videre

368 378 388 398 408 418 428 438 448

ved at addere 10, 20, ... 90 og

348 338 328 318 308 298 288 278 268

ved at subtrahere 10, 20, ... 90, altså ved at addere og subtrahere 1, 2, ... 9 i 10'ernes position. Og så fremdeles, også for hundrederne og tusinderne. Det er klart, at en hel masse af disse muligheder slet ikke svarer til simple forskelle i koden. Men de var nemme at generere og forbedringer kunne jo afvente resultaterne fra disse simple justeringer af koden.

Og for lige at gøre op, der er 37 ord med disse konstanter, hver med 4 tal, som hver skal afprøves i 18 mulige varianter i hver af de 4 cifferpositioner. I alt 10656 muligheder. Og tanken er nu at beregne double ln(x) for et eller flere passende x'er, i hver af de 10656 varianter, og så se, om man måske kommer tættere på det rapporterede resultat.

Efter forskellige indledende armbøjninger, fremkom resultatet, opsummeret i denne tabel, der afdækkede fejlen:


i	c[i]	d	c[i]+d	DPAln(1.47)-Haskell ln[i,d](1.47)	DPAln(1.93)-Haskell ln[i,d](1.93)	Sum relativ forskel	Haskell ln[i,d](1)
39	56	-60	-4	-7.81e-21	1.73e-21	2.29e-20	-2.85e-19
103	455	80	535	3.72e-21	-8.88e-21	2.32e-20	4.34e-20
103	455	90	545	-1.19e-20	-6.78e-22	3.20e-20	6.69e-20
71	8	-60	-52	4.08e-21	-1.60e-20	3.49e-20	-2.85e-19
23	6	60	66	-1.08e-20	-6.94e-21	3.87e-20	-3.52e-21
147	3	50	53	1.19e-20	5.49e-21	3.93e-20	-2.62e-19
123	672	-60	612	9.76e-22	2.86e-20	4.61e-20	-3.52e-21
87	497	60	557	1.28e-20	1.00e-20	4.84e-20	-3.52e-21
71	8	-70	-62	-1.67e-20	-6.21e-21	5.29e-20	-3.09e-19
87	497	70	567	-6.60e-21	2.41e-20	5.38e-20	2.00e-20

I tabellen er c[i] det i'te tal, regnet fra i = 0, i den liste af 148 10-bits tal, der udgør de 26 koefficienter (a0, a1, ..., a25) til det approksimerende 25.-grads polynomium. Så f.eks. c[0] = 0, c[1] = 556, c[2] = 144 og c[3] = 566 fra linie 62 i koden. Og fejlen således, at c[23] skal rettes fra 6 til 66. Tallet d er en justering af c[i], -9, -8, ..., 9 ganget med en 10-potens, og c[i]+d den justerede værdi.

Første linie i tabellen handler således om c[39], som er sidste 10-bits gruppe med værdien 56 i linie 71 i koden. Justeringen d er -60, så den justerede værdi er 56-60 = -4. Ikke nogen særlig sandsynlig forskel, selvfølgelig, men altså en af de afprøvede kandidater.

For at vægte de forskellige kandidater, bruges forskellen mellem DPAln(x) og Haskell ln[i,d](x): DPAln(x) er værdien af double ln(x), som den er angivet i DPA dokumentet. For x mellem 1 og 2, har vi rådighed over:

Test 10.127 double ln, double exp, exlog:
  DPAln(1.47) = 3.85262400790644933614 * 10^(-1)
  DPAln(1.93) = 6.57520002916794183895 * 10^(-1)
Test 10.128 double ln
  DPAln(1.00) = -3.38813178901720135627 * 10^(-21)

Og Haskell ln[i,d](x) er værdien, beregnet med double ln() metoden, med det tidligere omtalte Haskell kode, der i hvert fald forsøger at følge double ln() nøje, med c[i] justeret med d.

Så igen for første linie i tabellen, ser vi denne forskel for x = 1.47 udregnet til -7.81e-21. Og for x = 1.93 er forskellen 1.73e-21.

Til vægtningen af de forskellige kandidater til justering, er sådanne absolutte forskelle ikke så egnede. I stedet bruges relative forskelle, altså den absolutte forskel, normeret med værdien af de sammenlignede størrelser. Det kan gøres på flere måder, jeg har brugt

relativ forskel(a,b) = 2*(abs(a-b))/(abs(a)+abs(b))

Og det er summen af de to relative forskelle, for x = 1.47 og x = 1.93, der står i kolonnen "Sum relativ forskel". Denne kolonne bruges så som vægten af justeringskandidaten [i,d] og tabellen er de 10 kandidater med laveste vægt, altså de 10 mindste summer af relative forskelle.

Så puster vi lidt ud og kigger på kandidaterne: Altså, den første ændrer som sagt c[39] fra 56 til -4, ikke så sandsynlig en forskel.

De to næste handler om c[103] = 455, der prøves ændret til både 535 og 545. Her kiggede jeg mest på den anden, der jo godt kunne tænkes, altså en ombytning af to cifre. (Ikke noget, jeg oprindelig havde overvejet, må jeg indrømme.)

Fjerde mulighed, justering af c[71] fra 8 til -52: Nej.

Og så den femte, c[23] ændres fra 6 til 66. Jeg sad jo og kiggede tilbage i koden for hver kandidat og jeg må sige, at lige med det samme blev jeg nærmest fuldstændig sikker, for kig engang tilbage på linie 67, som det handler om: I netop denne linie overholdes de pæne kolonner ikke: I alle øvrige linier er tallene stillet op med sidste ciffer i samme kolonne. Så det passer præcis med det manglende ciffer, 6.

En hurtig justering af Mogens' GA4 program og kørsel på simulatoren gav den forløsende udskrift

Bad:  -1.45689666927739658319'- 19
Good: -3.38813178901720135627'- 21
This: -3.38813178901720135627'- 21

Det så jo lovende ud. Så besked til Mogens, der svarede tilbage efter en halv times tid, at med den angivne justering passede double ln() test resultaterne nu med de rapporterede.

Sidste kolonne i tabellen er blot den beregnede værdi af Haskell ln[i,d](1), som teoretisk skulle være = 0, men som blot er et meget lille tal. Denne kolonne ses også at pege på [23,60] som en lovende kandidat.

Endelig kan vi notere, at den ukorrigerede double ln() har vægten 4.48e-19, så det er ikke meget, der skiller fårene fra bukkene.

Detaljer: Multilængde aritmetik i Haskell

Dobbelt-længde på GIER er en del, 78-79 bits i fastkomma, 68-69 bits flydende, forudsat de gængse 10 bits er sat af til eksponenten: For dobbelt-længde flydende, noget der ligner 20 betydende decimale cifre. Så gængs "double" (ofte IEEE standard dit og dat) med et antal betydende decimale cifre omkring 15 kommer til kort her.

Haskell har indbygget multilængde heltal, men ikke multilængde flydende i nogen væsentlig henseende. Der findes biblioteker til multilængde flydende beregninger, jeg fandt

http://hackage.haskell.org/package/rounded
http://hackage.haskell.org/package/haskell-mpfr

begge baseret på

https://www.mpfr.org/

men denne vej så umiddelbart utiltalende indviklet ud.

Så jeg valgte en anden vej, nemlig at bruge Haskell's indbyggede rationelle tal, hvor både tæller og nævner er multilængde heltal. Haskell bruger % som brøkstreg i denne sammenhæng, så man kan skrive f.eks. 1%2 + 1%3 og få resultatet 5%6 i betydningen 1/2+1/3 = 5/6. Tanken er, at hvis man gennemfører beregninger med rationelle tal, får man jo faktisk et eksakt resultat, hvis man antager, at alle indgående størrelser også er eksakte. Det kan man selvfølgelig ikke gå ud fra, men i det mindste kan man være vis på, at der ikke tilføres afrundingsfejl under beregningerne.

Som eksempel, lad os beregne ln(2) med denne teknik. Som det måske er bekendt og i hvert fald væsentligt for det aktuelle emne, approksimation af matematiske funktioner, findes der rækkeudviklinger til den slags. Jeg kan f.eks. slå op i min gummibibel (CRC Standard Mathematical Tables) eller kigge på

https://da.wikipedia.org/wiki/Naturlig_logaritme

og finde, at

ln(x) = (x-1) - (x-1)^2/2 + (x-1)^3/3 - (x-1)^4/4 + ...

altså med x = 2, at

ln(2) = 1 - 1/2 + 1/3 - 1/4 + ...

I Haskell kan man, om kronket eller elegant vil jeg lade være usagt, skrive

[ (-1)^(k-1)%k | k <- [1..10] ]

og få de første 10 led i summen frem:

[1%1,(-1)%2,1%3,(-1)%4,1%5,(-1)%6,1%7,(-1)%8,1%9,(-1)%10]

Lægger men sammen får man først 1627 % 2520 og med divisionen udført, at ln(2) er ca. lig med 0.645634920634920... Afrundet til 5 decimaler er ln(2) = 0.69315, så det er jo ikke specielt tæt på. Summen konvergerer faktisk langsomt, man skal ud på tusindvis af led, bare for få snoldede 4 decimaler af ln(2).

(Her skal jeg lige indskyde, at hvis man vil prøve Haskell selv, findes der råd på https://www.haskell.org/. Kvik guide: Der findes flere implementationer, flagskibet ghc (Glasgow Haskell Compiler) findes til i hvert fald Fedora Linux og som FreeBSD port. Til Windows brugte jeg, sidst jeg havde behov, noget, der hed Haskell Platform, men det er længe siden. Der er også den mere beskedne Hugs, længe ikke vedligeholdt, men fuldt i stand til at håndtere den aktuelle sag. Der findes Fedora Linux installation (hugs98) og en FreeBSD port (som dog ikke virkede for mig). Til Windows og almindelig oplysning om Hugs, kig på https://www.haskell.org/hugs/.)

Videre med rækkeudviklinger for ln(x), findes

ln(x) = 2*(x-1)/(x+1) + 2*((x-1)/(x+1))^3/3 + 2*((x-1)/(x+1))^5/5 + ...

altså med x = 2

ln(2) = 2/3 + 2/3^3/3 + 2/3^5/5 + 2/3^7/7 + ...

med meget bedre konvergens: 10 led er nok til at skaffe os de første 10 decimaler:

ln(2) = 0.6931471805...

Og lige for klarhedens skyld, i alle detaljer: Vi tager de 10 led:

take 10 [ 2 % (3^k*k) | k <- [1,3..] ]

De ser således ud:

[2 % 3,2 % 81,2 % 1215,2 % 15309,2 % 177147,2 % 1948617,2 % 20726199,2 % 215233605,2 % 2195382771,2 % 22082967873]

Dem adderer vi og får denne bombastiske brøk:

1302374561632216 % 1878929321474205

For læseligheden tager vi tælleren og ganger med 10^10 og heltalsdividerer med nævneren

(1302374561632216 * 10^10) `div` 1878929321474205

så vi kan se de første 10 decimaler af vores tilnærmelse til ln(2):

6931471805

Altså, gange med en 10-potens og heltalsdividere er jo for at bevare alle cifrene: Hvis vi i stedet konverterer tæller og nævner til gængs "double" og dividerer, er vi jo straks nede på gængs "double" præcision, og det duer jo ikke i almindelighed.

De GIER'ske dobbelt-længde konstanter

Til at dublere double ln() beregningen må vi jo kunne forstå, hvordan de mange konstanter, der indgår, egentlig skal fortolkes. Altså, når der nu står f.eks.

 62 b:       0/ 556/ 144/ 566        ; constants*2|&(-10): a0
 63         53/ 381/ 566/ 463

som forventeligt skal fortolkes til at repræsentere en koefficient, a0, hvad er det så egentlig for en værdi, der er tale om? Til at befordre vores forståelse, kan vi bruge konstanten ln(2), der indgår på denne måde:

 99 b2:    177/ 456/ 383/ 500        ; ln(2)
100        231/ 755/ 350/ 316

Umiddelbart repræsenterer de to linier heltallene

177*2^30 + 456*2^20 + 383*2^10 + 500 = 190530846196
231*2^30 + 755*2^20 + 350*2^10 + 316 = 248826394940

i henhold til den slip'ske /-notation.

Det skal med det samme siges, at disse talkonstanter ikke er i den flydende repræsentation, der ellers bruges af DPA-koden: I stedet er det en fastkomma repræsentation, som den vi kender fra heltal, men, lidt uvant i nyere tid, med et underforstået binært punkt et fast sted i den binære repræsentation. Et eksempel på sådan en fastkomma repræsentation beskrives i Lærebog i kodning for GIER, bind 1, en fortolkning af de 40 bits til tal, som et GIER ord består af, hvor der står et underforstået binært punkt mellem den to mest betydende bits i ordet. Som skrevet står:

... hver celle består af 42 bits, som vi nummererer fra 0 til 41
inclusive; de første 40, bits nr. 0-39, kan bruges til tallagring ...
Ved regning med fast komma fungerer maskinens indbyggede aritmetiske
operationer, som om kommaet er placeret mellem position nr. 0 og
position nr. 1; tal med 0 i position 0 opfattes som tal i intervallet
0 <= x < 1, mens tal med 1 i position 0 opfattes som tal i intervallet
-1 <= x < 0, nemlig som talværdien af binærbrøken (et tal mellem 1 og
2) minus 2. Denne talrepræsentation kan karakteriseres ved, at
maskinen regner modulo 2 i intervallet -1 <= x < 1 med 39 binære cifre
efter kommaet.

Videre med ln(2) finder vi, at 190530846196, heltalskonstanten, der repræsenteres af den første linie konstanter, divideret med 2^38, giver 0.69314..., som fint passer med ln(2). Der er altså ikke binært punkt ("komma") mellem bit nr. 0 og 1, men mellem bit nr. 1 og 2.

Til den anden del af ln(2), 248826394940, finder vi, at den mest signifikante bit, med nr. 0, ikke skal tælles i regnestykket, så 2-potensen, der skal divideres med, bliver 38+39=77 og dermed

190530846196/2^38 + 248826394940/2^77 = 0.6931471805599453094172290...

med 25 decimaler. I god overensstemmelse med de første 25 decimaler af ln(2), beregnet med ovenstående rækkeudvikling:

ln(2) = 0.6931471805599453094172321...

Fortolkningen, hvor vi ikke bruger den mest signifikante bit af den anden del af et dobbelt tal, stemmer med konventionen for GIER dobbelt-præcision operationer med fast komma, som er beskrevet i "GIER - A Danish Computer of Medium Size", IEEE Transactions on Electronic Computers, Vol. EC-12, Number 5, December, 1963.

Denne fortolkning kan vi så forsøge at bruge på de øvrige konstanter, a'erne, der indgår i beregningen af det approksimerende 25.-grads polynomium. Så i første omgang finder vi f.eks. for a0, der er skrevet ved

 62 b:       0/ 556/ 144/ 566        ; constants*2|&(-10): a0
 63         53/ 381/ 566/ 463

at en tilsvarende fortolkning giver

a0 = 0.0021215101809657185364938...

med 25 decimaler. Nu skal vi ikke dvæle for meget ved den præcise værdi af disse konstanter: Som der står i kommentaren ("constants*2|&(-10)", der jo skal læses som "constants*2^(-10)") er der åbenbart allerede foretaget en skalering med 2^(-10). Og også i koden, efter evaluering af det approksimerende polynomium, sker der noget skalering med 2-potenser af den beregnede værdi (linierne 44-45). Så i første omgang prøver vi blot en fornuftig fortolkning, der ser ud til at passe på omstændighederne, og så ser vi, hvad der kommer ud af det.

En ting mangler vi stadig: Negative tal. F.eks. finder vi for konstanten a2

 66       1023/1002/ 262/  25        ; a2
 67        358/  69/ 715/  6

med lutter binære 1-taller i den mest signifikante ende. Her er GIER-konventionen jo at bruge 2-komplement, altså hvis den mest signifikante bit er 1, skal tallet opfattes som det negative tal, der fremkommer ved at subtrahere 2-potensen svarende til denne bit. Med denne konvention finder vi:

a2 = -0.0000829472202918561784541...

For små konstanter, der kan nøjes med et enkelt ord, finder vi f.eks. for

 84        273/ 466/ 405/ 750        ; a11

at

a11 = 0.0000000000019430179425341...

For små tal, der også er negative, bruger vi igen 2-komplement konventionen, så vi f.eks. for

 85        980/1016/  70/ 409        ; a12

finder, at

a12 = -0.0000000000003055884104830...

Men, som sagt, den præcise størrelsesorden af disse konstanter tager vi ikke så nøje i første omgang.

Justering: Beregning af log2(2)

Til endelig afgørelse af spørgsmålet om, hvad er det egentlig for en skalering, der bliver brugt, vil vi blot efter bedste evne implementere det 25.-grads polynomium, der approksimerer log2(y), og så se, hvad der kommer ud af det.

Heldigvis er kommentarerne til evalueringen udmærkede. En nærmere granskning afslører, at evalueringen gentager et grundlæggende skridt, hvor sættet med de tre størrelser (b,b1,b2), startende med (0,0,0), i hvert skridt bringes over i:

( 2*b*t - b1 + a[p], b, b1 )

Her er t argumentet til log2(y), skaleret til intervallet -1 <= t < 1 ved t = 2*y - 3 og a[p] er en koefficient fra listen a0, a1, ... Og slutresultatet fremkommer som værdien af (b-b2)/2 for det sidste sæt.

Som vi erindrer, her beregner vi jo log2(), logaritmen med grundtallet 2. Så vi kan jo prøve at beregne log2(y) med y = 2 og se, om vi kommer i nærheden af resultatet, som jo gerne skulle være = 1. Vi skal lige huske at sætte t = 2*y - 3, som jo er = 1 for y = 2.

Og, regne, regne, ud kommer tallet 0.0019531249999999999998047 og det ligner selvfølgelig noget, der kunne være en pæn brøk. I minsandten: 1/0.001953125 = 512 så vi skal bare huske at skalere alle a'erne med en faktor 512.

Således er vi i hus og kan beregne log2() med vores Haskell dublet, og dermed double ln() ved efterfølgende multiplikation med ln(2). Og gennemføre de beskrevne eksperimenter med varianter, der leder til fejlen.

Det approksimerende polynomium

Lidt af nysgerrighed, men også for eventuelt at kaste lys over, hvad det egentlig er for et approksimerende polynomium, der bliver brugt, har jeg brugt lidt Haskell-trickeri til at beregne nogle mere eksplicitte koefficienter til polynomiet. Sagen er jo, at der indgår disse konstanter a0, a1, ..., men på grund af den lidt kringlede logik med tallene b, b1 og b2, er det jo ikke lysende klart, hvad der egentlig foregår. Om så en mere eksplicit opskrivning af polynomiet vil hjælpe vil måske vise sig.

Uanset: Hvis vi optrevler (b,b1,b2)-logikken, finder vi, at det polynomium i t, der faktisk beregnes, ser således ud, koefficienterne skrevet med 25 decimaler:

p(t) =
  +0.5849625007211561815096339
  +0.4808983469629878026568034*t
  -0.0801497244938312998035730*t^2
  +0.0178110498875180551500490*t^3
  -0.0044527624718795391341261*t^4
  +0.0011874033258348311368861*t^5
  -0.0003298342571758053776848*t^6
  +0.0000942383591894804199112*t^7
  -0.0000274861881033373368909*t^8
  +0.0000081440557606494295983*t^9
  -0.0000024432166824371603919*t^10
  +0.0000007403685703850915178*t^11
  -0.0000002262239376452090944*t^12
  +0.0000000696077350481605350*t^13
  -0.0000000215446253060136427*t^14
  +0.0000000067020082261670665*t^15
  -0.0000000020956579804831676*t^16
  +0.0000000006585383349744233*t^17
  -0.0000000002055307035675468*t^18
  +0.0000000000638848973721906*t^19
  -0.0000000000219060325434838*t^20
  +0.0000000000075743855632026*t^21
  -0.0000000000013997691894473*t^22
  +0.0000000000002273736754432*t^23
  -0.0000000000004263256414560*t^24
  +0.0000000000001705302565824*t^25

Og det er jo, skal vi huske, approksimationen af log2(y) for 1 <= y < 2, hvor t = 2*y - 3, så -1 <= t < 1.

Lad os sige, at d[i] er koefficienten til t^i i p(t), så f.eks. d[0] = 0.5849625007211561815096339 og d[10] = -0.0000024432166824371603919.

Det første tal her, d[0] = 0.58496..., det er jo p(0), altså p(t) for t = 0. Og regner vi baglæns tilbage til y med y = (t + 3)/2, skulle dette jo så være log2(3/2). Og ganske rigtigt: Min regnemaskine siger

log2(3/2) = ln(3/2)/ln(2) = 0.584962500721156

hvor vi bruger log2(x) = ln(x)/ln(2), fordi ln() er mere almindeligt tilgængelig end log2(). Så

d[0] = ln(3/2)/ln(2)

Inspireret af dette, prøver vi at gange de øvrige koefficienter med ln(2). Her er de første par stykker:

d[1]*ln(2) = 0.333333333333333
d[2]*ln(2) = -0.0555555555555556
d[3]*ln(2) = 0.0123456790123457

Det ser jo påfaldende "pænt" ud, så vi reciprokker og finder:

1/(d[1]*ln(2)) = 3.0                    =  1*3^1
1/(d[2]*ln(2)) = -18.0                  = -2*3^2
1/(d[3]*ln(2)) = 81.0000000000001   ca. =  3*3^3
1/(d[4]*ln(2)) = -323.999999999998  ca. = -4*3^4
1/(d[5]*ln(2)) = 1214.9999999997    ca. =  5*3^5
1/(d[6]*ln(2)) = -4374.00000000604  ca. = -6*3^6
1/(d[7]*ln(2)) = 15309.0000006071   ca. =  7*3^7
1/(d[8]*ln(2)) = -52487.9999898492  ca. = -8*3^8
1/(d[9]*ln(2)) = 177146.999393079   ca. =  9*3^9
1/(d[10]*ln(2)) = -590490.009035893 ca. = -10*3^10
...

Så det ser jo ud til, at for k = 1, 2, ..., er

1/(d[k]*ln(2)) = (-1)^(k-1) * k * 3^k

eller med andre ord

d[k]*ln(2) = (-1)^(k-1) / k / 3^k

eller

d[k] = (-1)^(k-1) / k / 3^k / ln(2)

Nå, så har vi skovlen under, hvad der foregår: Vi skal beregne log2(y). Det gør vi, lidt indirekte må man sige, for det er jo egentlig ln(y) vi vil beregne, ved først at omskrive log2(y) = ln(y)/ln(2). Dernæst indskrænker vi vores argumentinterval ved omskrivningen

ln(y) = ln((3/2)*(2/3)*y) = ln(3/2) + ln((2/3)*y)

I termer af t = 2*y - 3, så y = (t+3)/2, er

(2/3)*y = (2/3)*(t+3)/2 = t/3 + 1

så

ln(y) = ln(3/2) + ln(t/3 + 1)

Nu har vi jo tidligere set rækkeudviklingen

ln(x) = (x-1) - (x-1)^2/2 + (x-1)^3/3 - (x-1)^4/4 + ...

Hvis vi nu skriver u = x-1, så x = u+1, står der jo:

ln(u+1) = u - u^2/2 + u^3/3 - u^4/4 + ...

Og med u = t/3:

ln(t/3 + 1) = (t/3) - (t/3)^2/2 + (t/3)^3/3 - (t/3)^4/4 + ...
  = t/3 - t^2/2/3^2 + t^3/3/3^3 - t^4/3/3^4 + ...

Og kniber vi øjnene sammen, kan vi jo godt se koefficienten til t^k: Den er (-1)^(k-1)/k/3^k og dermed = d[k]*ln(2), som vi fandt tidligere.

Il grande finale:

log2(y) = ln(y)/ln(2) = [ln(3/2) + ln(t/3 + 1)]/ln(2)
  = [ln(3/2) + t/3 - t^2/2/3^2 + t^3/3/3^3 - t^4/3/3^4 + ...]/ln(2)
  = ln(3/2)/ln(2) + t/3/ln(2) - t^2/2/3^2/ln(2) + t^3/3/3^3/ln(2) - t^4/3/3^4/ln(2) + ...

Så:

d[0] = ln(3/2)/ln(2)
d[1] = 1/3/ln(2)
d[2] = -1/2/3^2/ln(2)
...
d[k] = (-1)^(k-1)/k/3^k/ln(2)
...

Dermed er vi nået frem til samme d[k], koefficienten til t^k i polynomiet p(t), ved at regne forfra, som vi gættede på ved numeriske betragtninger.

Haskell trickeri

For at forklare i detaljer hvordan den optrevlede udgave, p(t), af det approksimerende 25.-grads polynomium fremkommer, gengiver jeg her det Haskell modul, FuncEval1, der indeholder de anvendte beregningsfunktioner:

  1 -- FuncEval1.hs: Function evaluation experiments, inspired by the GAIIIDobbeltPraecision double ln() questions
  2 -- 2019-Jan-27 11.34 / TN
  3 
  4 {-
  5 
  6 Basis: original/dp1.asc, typed in and adjusted by MK:
  7 
  8 algol<
  9 _b_e_g_i_n
    ...
210 _e_n_d;
211 t<
212 
213 -}
214 
215 module Main (
216   module Main
217 
218   , module Data.Maybe
219   , module Numeric
220   , module Data.List
221   , module Data.Ratio
222 
223   , module Primtal
224   , module Poly
225   ) where
226 
227   import Data.Maybe
228   import Numeric
229   import Data.List
230   import Data.Ratio
231 
232   import qualified Primtal
233   import Poly
234 
235   progName :: String
236   progName = "FuncEval1"
237   progStamp :: String
238   progStamp = "2020-May-24 05.51"
239 
240   -- Present a rational number as a decimal fraction with a given number of decimals:
241 
242   feShow1 decs val
243     = (if val < 0 then "-" else "") ++ integerPart ++ "." ++ fractionPart
244       where
245       (q,r) = (abs $ numerator val) `divMod` (denominator val)
246       integerPart = show q
247       b = 10 ^ decs
248       fractionPart = tail ( show ( b + r*b `div` (denominator val) ) )
249 
250   -- With rounding:
251 
252   feShow2 decs val = feShow1 decs (val + (signum val) / (2*10^decs))
253 
254   -- Taylor series for arctan:
255 
256   arctanSeries :: Ratio Integer -> [Ratio Integer]
257   arctanSeries x = zipWith (/) ( iterate (*(negate x * x)) x ) [ 1, 3 .. ]
258 
259   -- ln(x), 0 < x <= 2:
260 
261   -- lnSeries1 :: Ratio Integer -> [Ratio Integer]
262   lnSeries1 x = zipWith (/) ( iterate (*(1-x)) (x-1) ) [ 1, 2 .. ]
263 
264   -- ln(x), x > 0:
265 
266   -- lnSeries2 :: Ratio Integer -> [Ratio Integer]
267   lnSeries2 x = zipWith (/) ( iterate (*((x-1)/(x+1))^2) (2*(x-1)/(x+1)) ) [ 1, 3 .. ]
268 
269   -- ln(x), 0 < x <= 2:
270 
271   lnSeries3 x = zipWith (*) ( iterate (*([1]-x)) (x-[1]) ) $ map ((:[]) . (1%)) [ 1, 2 .. ]
272 
273   -- Converting from rational to fractional:
274 
275   toFractional r = (fromIntegral $ numerator r) / (fromIntegral $ denominator r)
276 
277   -- ln(x) using truncated lnseries2:
278 
279   lnEps eps x = sum $ takeWhile ( (>= eps) . abs . toFractional ) $ lnSeries2 x
280 
281   -- Relative difference:
282 
283   relDiff a b = 2 * abs ( a - b ) / ( abs a + abs b )
284 
285   -- Adjusting numbers in a list (of lists (of lists)) of numbers.
286 
287   length1 = genericLength
288   length2 = (sum . (map length1))
289   length3 = (sum . (map length2))
290 
291   update1 f n [] = []
292   update1 f n (x:xs) = (if n == 0 then f x else x) : update1 f (n-1) xs
293   update2 f n [] = []
294   update2 f n (xs:xss) = update1 f n xs : update2 f (n - length1 xs) xss
295   update3 f n [] = []
296   update3 f n (xss:xsss) = update2 f n xss : update3 f (n - length2 xss) xsss
297 
298   -- Constants from the code, typein:
299 
300   ln2List1 = [ 177, 456, 383, 500, 231, 755, 350, 316 ]
301   ln2List2 = [ 177, 456, 383, 500 ] ++ map (*2) [ 231, 755, 350, 316 ]
302 
303   -- Extract from the Mogens doubleln() code, based on:
304   --
305   --   thorkilnaur@DESKTOP-MLR81P9 /cygdrive/e/tn/DDHF/GIER/GAIIIDobbeltPraecision
306   --   $ awk '{if(match($0,/qq *([0-9]+)\.9\+ *([0-9]+)\.19\+ *([0-9]+)\.29\+ *([0-9]+)\.39/,a)){if(match($0,/;.*(a[0-9]+|b2)/,b)){d=e "        " f "[ [";e="] ] -- " b[1]"\n";f=", "};printf d a[1]","a[2]","a[3]","a[4]};d="], ["};END{printf e}' original/dp1.asc
307   --
308 
309   dp1Constants
310     = [
311         [ [0,556,144,566], [53,381,566,463] ] -- a0
312         , [ [0,253,479,229], [279,1016,744,971] ] -- a1
313         , [ [1023,1002,262,25], [358,69,715,6] ] -- a2
314         , [ [0,2,498,844], [351,747,646,94] ] -- a3
315         , [ [1023,1023,696,282], [30,184,831,56] ] -- a4
316         , [ [0,0,44,1006], [258,507,538,494] ] -- a5
317         , [ [1023,1023,1017,582], [52,752,508,222] ] -- a6
318         , [ [0,0,0,968], [275,594,951,406] ] -- a7
319         , [ [1023,1023,1023,878], [305,664,256,8] ] -- a8
320         , [ [0,0,0,22], [89,773,283,1007] ] -- a9
321         , [ [1023,1023,1023,1020], [294,822,514,497] ] -- a10
322         , [ [273,466,405,750] ] -- a11
323         , [ [980,1016,70,409] ] -- a12
324         , [ [6,830,842,171] ] -- a13
325         , [ [1022,936,806,455] ] -- a14
326         , [ [0,177,966,310] ] -- a15
327         , [ [1023,995,386,929] ] -- a16
328         , [ [0,4,636,863] ] -- a17
329         , [ [1023,1023,257,87] ] -- a18
330         , [ [0,0,124,672] ] -- a19
331         , [ [1023,1023,1003,698] ] -- a20
332         , [ [0,0,3,328] ] -- a21
333         , [ [1023,1023,1023,467] ] -- a22
334         , [ [0,0,0,91] ] -- a23
335         , [ [1023,1023,1023,1009] ] -- a24
336         , [ [0,0,0,3] ] -- a25
337         , [ [177,456,383,500], [231,755,350,316] ] -- b2
338     ]
339 
340   as = init dp1Constants
341 
342   ln2 = last dp1Constants
343 
344   -- Adjusted series:
345 
346   as2 = update3 (\x -> x-50) 103 as
347   as3 = update3 (\x -> x+90) 103 as
348   as4 = update3 (\x -> x+60) 23 as
349 
350   -- Interpretations:
351 
352   base1 = 2^10
353 
354   toRational1 (ds@(d0:_))
355     = if 2*d0 < base1 then c else c - 4
356     where
357     c = sum $ zipWith (*) (map fromIntegral ds) (iterate (/fromIntegral base1) (4 % fromIntegral base1))
358 
359   base2 = 2^39
360 
361   toRational2 [ds] = toRational1 ds / fromIntegral base2
362   toRational2 [ds0,ds1] = toRational1 ds0 + toRational1 ds1 / fromIntegral base2
363 
364   -- Polynomium type:
365 
366   instance Num a => Num [a]
367     where
368     (+) x y = polAdd (+) x y
369     (-) x y = polAdd (-) x y
370     (*) x y = polMult (+) (*) x y
371     fromInteger x = if x /= 0 then [fromIntegral x] else []
372 
373   -- Calculations:
374 
375   feCalc fromScalar t (b,b1,_) ap = ( b * t * fromScalar 2 - b1 + fromScalar ap, b, b1 )
376 
377   feLogBase_0 fromScalar as t
378     = let
379         (b,b1,b2) = foldl (feCalc fromScalar t) (0,0,0) (map ((*512) . toRational2) $ reverse $ as)
380       in
381         (b-b2)*(fromScalar $ 1%2)
382 
383   feLogBase_1 = feLogBase_0 id
384 
385   feLogBase_2 = feLogBase_0 (:[])
386 
387   feLog2_1 as x = feLogBase_1 as ((x-1)*2-1)
388   feLog2_2 as x = feLogBase_2 as ((x-1)*2-1)
389 
390   feLn_1 as x = feLog2_1 as x * toRational2 ln2
391 
392   -- Chebyshev polynomials:
393 
394   t = [1] : [0,1] : zipWith (-) (tail $ map (*[0,2]) t) t
395 
396   -- Ralston, A first course in numerical analysis, chapter 7, problem {27}:
397 
398   -- ...
399 
400   -- List of ( test, arg, double ln( arg ) ) from gaIIIdp2.pdf:
401 
402   doubleLns
403     = [
404         ( 95%100, -512932943875505333542%10^22 ) -- Test 10.127 double ln, double exp, exlog
405         , ( 147%100, 385262400790644933614%10^21 ) -- Test 10.127 double ln, double exp, exlog
406         , ( 193%100, 657520002916794183895%10^21 ) -- Test 10.127 double ln, double exp, exlog
407         , ( 1%1, -338813178901720135627%10^41 ) -- Test 10.128 double ln
408       ]
409 
410   -- Main:
411 
412   main :: IO ()
413   main
414     = putStrLn $ progName ++ " (" ++ progStamp ++ ")"

For rationelle tal implementeres det tidligere omtalte grundskridt i beregningen af det approksimerende polynomium

 ( b, b1, b2 ) := ( 2*t*b - b1 + a[p], b, b1 )

af funktionen feCalc:

375   feCalc fromScalar t (b,b1,_) ap = ( b * t * fromScalar 2 - b1 + fromScalar ap, b, b1 )

Parameteren fromScalar er en funktion, der konverterer skalarer (rationelle tal) til den type værdier, der indgår i beregningen. Og det er i første ombæring jo også bare rationelle tal, så man kan til en start tænke på fromScalar som blot identitetsfunktionen, der returnerer sit argument uændret.

Parameteren t til feCalc er det omskalerede y*2-3 og funktionen (feCalc fromScalar t) bringer således (b,b1,ligegyldigt) over i (2*t*b-b1+ap,b,b1) ved anvendelse af en enkelt koefficient, ap, i listen af a'er.

I funktionen feLogBase_0

377   feLogBase_0 fromScalar as t
378     = let
379         (b,b1,b2) = foldl (feCalc fromScalar t) (0,0,0) (map ((*512) . toRational2) $ reverse $ as)
380       in
381         (b-b2)*(fromScalar $ 1%2)

der tager fromScalar, a'erne og t som parametre, itereres (feCalc fromScalar t) så med standard funktionen foldl over (b,b1,b2), startende med (0,0,0), og med det ene ap efter det andet fra listen af a'er. Denne liste, parameteren as, er som udgangspunkt blot den globale

340   as = init dp1Constants

der indeholder alle, undtagen den sidste, af konstanterne fra den ligeledes globale

309   dp1Constants
310     = [
311         [ [0,556,144,566], [53,381,566,463] ] -- a0
            ...
336         , [ [0,0,0,3] ] -- a25
337         , [ [177,456,383,500], [231,755,350,316] ] -- b2
338     ]

der jo ultimativt stammer fra Mogens' indtastede double ln() kode. Der er også forskellige, justerede, udgaver af a'erne:

346   as2 = update3 (\x -> x-50) 103 as
347   as3 = update3 (\x -> x+90) 103 as
348   as4 = update3 (\x -> x+60) 23 as

Ikke mindst as4, der repræsenterer den korrigerede udgave.

Med reverse $ as arrangeres, at a'erne fiskes ud af listen i den rigtige rækkefølge (a25, a24, ...). Hvert a udsættes med ((*512) . toRational2) for at blive konverteret til et rationelt tal og dernæst multipliceret med 512, som vi tidligere har fundet gav det rigtige resultat.

På den måde beregnes talsættene

 ( 0, 0, 0 ), ( a25, 0, 0 ), ( 2*t*a25+a24, a25, 0 ),
   ( 2*t*(2*t*a25+a24)-a25+a23, 2*t*a25+a24, a25 ), ...

Resultatet af feLogBase_0 er (b-b2)/2 for den sidste af disse talsæt.

Som eksempel lad os beregne en approksimation til log2( y=1.23 ). Vi får hjælp af

383   feLogBase_1 = feLogBase_0 id

der kalder feLogBase_0 med fromScalar sat til den meget indviklede standardfunktion id x = x, der blot returnerer sit argument uændret.

Til beregningen af log2( y=1.23 ), skal vi først finde t = y*2-3 = 1.23*2-3 = -0.54 = -54/100. Med det t kan vi beregne

 Main> feLogBase_1 as (-54%100)
 3914574273767213354065209997885358443043136709029164551
   % 13107200000000000000000000000000000000000000000000000000

("Main>" er Hugs' prompt under interaktive beregninger.) Og det ser jo rædderligt ud, men konverterer vi til decimal med 25 decimaler (feShow1 25), får vi

 Main> feShow1 25 $ feLogBase_1 as (-54%100)
 "0.2986583155645151789905708"

i god overensstemmelse med

 Main> let log2 x = log x / log 2 in log2 1.23
 0.298658315564515

("log" er Haskell's naturlige logaritme, loge() eller ln()).

Polynomier repræsenteres ved en bar liste af tal, koefficienterne, så f.eks. polynomiet 2*x^2 + 3*x + 4 skrives [4,3,2]. Første tal i listen, 4, er således konstantleddet, koefficienten til x^0; det andet tal, 3, er koefficienten til x (altså x^1); og tredie tal, 2, er koefficienten til x^2.

Til beregninger med polynomier har vi polAdd (module Poly, ikke medtaget her) til addition, f.eks.

Main> polAdd (+) [16,8,4] [2,1]
[18,9,4]

til beregningen (4*x^2 + 8*x + 16) + (x + 2) = 4*x^2 + 9*x + 18. polAdd tager additionsoperationen til koefficienterne (+) som parameter. Det betyder bl.a., at vi også kan subtrahere. F.eks. kan vi med

Main> polAdd (-) [16,8,4] [2,1]
[14,7,4]

beregne (4*x^2 + 8*x + 16) - (x + 2) = 4*x^2 + 7*x + 14.

Multiplikation af polynomier sker med polMult, f.eks.

Main> polMult (+) (*) [-1,1] [-2,1]
[2,-3,1]

som beregner (x - 1)*(x - 2) = x^2 - 3*x + 2 og videre

Main> polMult (+) (*) [2,-3,1] [-3,1]
[-6,11,-6,1]

så vi får (x - 1)*(x - 2)*(x - 3) = x^3 - 6*x^2 + 11*x - 6, et trediegradspolynomie med de tre rødder (nulpunkter) 1, 2 og 3. polMult bruger parametrene (+) og (*) til at udføre de nødvendige operationer på koefficienterne.

Haskell er udstyret med en mekanisme, type klasser, til at definere hvad der kaldes ad-hoc polymorfi, altså dette, at samme funktion eller operator kan bruges på forskellige typer. Et eksempel er +, addition, der kan bruges på både heltal og flydende tal. (Haskell klasser er noget andet end, og skal ikke forveksles med, klassebegrebet fra objekt orienteret programmering.)

Som eksempel kan vi tage Haskell klassen Num, der omfatter tal i forskellige varianter. Eller mere præcist, ting som man kan addere, subtrahere eller gange sammen og i alle tilfælde få noget af samme slags som resultat. Division er specifikt sat udenfor døren her og er ikke med til festen.

Num defineres omtrent på denne måde:

class Num a where
  (+), (-), (*) :: a -> a -> a
  fromInteger   :: Integer -> a

"Num a" skal læses som udsagnet "typen a er i klassen Num". Og betingelserne, der skal sikre det, remses så op: Der skal være 3 funktioner, operatorer er det, +, - og *, der hver tager to værdier af typen a som parametre og producerer en værdi af typen a som resultat.

Endelig skal der være en funktion, fromInteger, der tager et heltal og konverterer det til en værdi af typen a.

Så kommer magien: Vi vil gerne passe vores polynomier ind i klassen Num og det sker med erklæringen:

366   instance Num a => Num [a]
367     where
368     (+) x y = polAdd (+) x y
369     (-) x y = polAdd (-) x y
370     (*) x y = polMult (+) (*) x y
371     fromInteger x = if x /= 0 then [fromIntegral x] else []

Ved "Num [a]" erklæres [a], typen af lister med elementer af type a, at tilhøre Num klassen. Den foranstillede "Num a =>" angiver betingelsen for, at Num [a] gælder, nemlig at typen a selv tilhører Num.

Og derefter angives så de konkrete funktioner (polAdd osv.), der skal understøtte typen [a] som medlem af Num klassen. Og en passende fromInteger, der håndterer konvertering af skalarer.

Nu kan vi skrive

Main> [-1,1]*[-2,1]*[-3,1]
[-6,11,-6,1]

og få udført den tidligere beregning af trediegradspolynomiet med rødderne 1, 2 og 3 på en lidt mere elegant måde. Og ikke nok med det, hele double ln() beregningen i feLogBase_0 kan gennemføres med polynomier i stedet for skalarer: Vi kan simpelthen sætte parameteren t til feLogBase_0 til at være et polynomium, altså en liste af koefficienter, i stedet for blot et tal. Det vil sige, helt af sig selv går det ikke, og det er så grunden til, at vi må hjælpe Haskell lidt med parameteren fromScalar på et par steder for at fortælle, at her skal en skalar opfattes på en bestemt måde.

Til brug for beregninger med polynomier i stedet for skalarer får vi hjælp af

385   feLogBase_2 = feLogBase_0 (:[])

der bruger fromScalar = (:[]) til konvertering af skalarer til 0'te-grads polynomier. (:[]) er en lidt skummel Haskell-specialitet: I Haskell er x:xs udtryk for en liste, der starter med elementet x, som påhæftes forrest på listen xs. For eksempel 1:[2,3] = [1,2,3]. Og (:[]) er så en såkaldt "section", hvor den ene operand af en to-operand operator, her :, er udeladt. Så (:[]) er en funktion, der tager et element og klistrer det foran en tom liste, eksempel (:[]) 1 = (1:[]) = [1]; med andre ord, (:[]) x skaber en liste [x] med blot det ene element, x. Så konverterer en skalar til et 0'te-grads polynomium med skalaren som konstantled, koefficienten til x^0.

Hjælpefunktionen fromScalar kan jo ses ses anvendt et par steder i feCalc og feLogBase_0. Og det er jo steder, hvor der skal ske konvertering fra et rationelt tal til et 0'te-grads polynomium. Andre steder, som f.eks. 0'erne i konstanten (0,0,0) i feLogBase_0, sker konverteringen i stilhed, først via fromInteger i min Num [a] instance, og dernæst via fromInteger i den instance-erklæring, der understøtter, at rationelle tal er med i Num klassen.

Men lad os se magien udfolde sig. I første omgang kan vi genberegne log(1.23), t = -0.54, nu som 0'te-grads polynomium:

Main> feLogBase_2 as [-54%100]
[3914574273767213354065209997885358443043136709029164551
  % 13107200000000000000000000000000000000000000000000000000]

Og for at få hoved og hale i det, med 25 decimaler:

Main> map (feShow1 25) $ feLogBase_2 as [-54%100]
["0.2986583155645151789905708"]

Og hvad er det så, vi gerne vil? Jo, vi vil gerne beregne double ln() af polynomiet x, så vi kan se, hvad er det egentlig for nogle koefficienter, der bruges på de enkelte potenser af t. Og polynomiet x, det hedder jo [0,1]:

Main> map (feShow1 25) $ feLogBase_2 as [0,1]
["0.5849625007211561817129218"
,"0.4808983469629878026568034"
,"-0.0801497244938313002101488"
,"0.0178110498875180551500490"
,"-0.0044527624718795391341261"
,"0.0011874033258348311368861"
,"-0.0003298342571758053776848"
,"0.0000942383591894804199112"
,"-0.0000274861881033373368909"
,"0.0000081440557606494295983"
,"-0.0000024432166824371603919"
,"0.0000007403685703850915178"
,"-0.0000002262239376452090944"
,"0.0000000696077350481605350"
,"-0.0000000215446253060136427"
,"0.0000000067020082261670665"
,"-0.0000000020956579804831676"
,"0.0000000006585383349744233"
,"-0.0000000002055307035675468"
,"0.0000000000638848973721906"
,"-0.0000000000219060325434838"
,"0.0000000000075743855632026"
,"-0.0000000000013997691894473"
,"0.0000000000002273736754432"
,"-0.0000000000004263256414560"
,"0.0000000000001705302565824"]

Og det er netop det p(t), som vi har set tidligere.

Som kontrol af magien, lad os bruge evaluering af polynomiet p til at beregne endnu en approksimation, log2( y=1.89 ). Udover polAdd og polMult har vi også polEval, der kan evaluere værdien af et polynomium for en given værdi af den indgående variabel. For eksempel, for det tidligere omtalte q(x) = x^3 - 6*x^2 + 11*x - 6 med nulpunkter 1, 2 og 3, kan vi beregne q(4) = 6:

Main> let q = [-1,1]*[-2,1]*[-3,1] in polEval (+) (*) q 4
6

Til log2( y=1.89 ), skal vi bruge t = 2*y-3 = 2*1.89-3 = 0.78 = 78/100. Og vi kan beregne:

Main> let p = feLogBase_2 as [0,1] in feShow1 25 $ polEval (+) (*) p (78%100)
"0.9183862344463479561130693"

Og det ser jo fint ud:

Main> let log2 x = log x / log 2 in log2 1.89
0.918386234446348

Hokuspokus filihankat.

Konklusion: Hvad med Chebyshev?

Beregningerne tager jo tilsyneladende noget af en omvej: Først evaluering af et polynomium, på kringlet vis, for log2(), der er fremkommet via en kendt rækkeudvikling for ln(), hver af koefficienterne så divideret med ln(2). Og endelig resultatet, log2(y), ganget med ln(2) for at få det endelige resultat, ln(y) for 1 <= y < 2. Turen omkring log2() er utvivlsomt begrundet i den lettere håndtering af den indledende skalering.

Men, som sagt, det hele er jo baseret på en kendt rækkeudvikling, Taylor-rækken ville man nok kalde den, så Chebyshev blev der ikke behov for i denne omgang. Det må blive en anden god gang.