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.