Wiskundigen maken er vaak een sport van om zoveel mogelijk decimalen van wiskundige constanten te berekenen, en in de Online Encyclopedia of Integer Sequences (de OEIS) vinden we daar veel voorbeelden van. De eerste 20000 cijfers van $\sqrt2$ vinden we in A002193 van de OEIS. Natuurlijk is het zo dat we het getal $\sqrt2$ nooit helemaal precies kunnen opschrijven in een decimale representatie, maar met een beetje inspanning kunnen we wel zoveel decimale cijfers produceren als we zelf willen. De programmeertaal Python blijkt daar uitermate geschikt voor te zijn.

De constante van Pythagoras is nauw verbonden met de stelling van Pythagoras. Het probleem destijds was: hoe lang is de diagonaal van een vierkant met zijde 1. Of: hoe lang is de zijde van een vierkant met oppervlakte $2$. Op beide vragen is het antwoord $\sqrt2$.

De Mesopotamiërs hadden daar al meer dan drieduizend jaar geleden een oplossing voor gevonden. Zij stelden het volgende. Laten we beginnen met een rechthoek met lengte $2$ en breedte $1$, zodat de oppervlakte van die rechthoek $2$ is. Neem vervolgens een nieuwe lengte voor een nieuwe rechthoek, door het rekenkundig gemiddelde te nemen van de lengte en de breedte van de vorige rechthoek. In de eerste stap is dat dus de lengte $\frac{2 + 1}{2} = \frac32$. Om de breedte te bepalen moeten we er voor zorgen dat de oppervlakte van de rechthoek $2$ blijft; in de informatica noemen ze zoiets een invariant, maar daar komen we zo meteen op terug. De nieuwe breedte wordt dus $\frac{2}{\frac32} = \frac43$. Met deze nieuwe rechthoek kunnen we weer hetzelfde procedé herhalen. In de tweede stap vinden we dus de lengte $\frac{\frac32 + \frac43}{2} = \frac{17}{12}$, en de bijbehorende breedte wordt $\frac{2}{\frac{17}{12}} = \frac{24}{17}$. Dit proces kunnen we eindeloos herhalen waarbij we telkens rationale getallen vinden die steeds dichter bij de werkelijke waarde van $\sqrt2$ komen te liggen. De Mesopotamiërs moesten al die berekeningen met de hand uitvoeren op een kleitablet, maar desondanks waren ze in staat om het getal $\sqrt2$ tot zes cijfers achter de komma te bepalen, wat voor die tijd een hele prestatie was. Tegenwoordig hoeven we dat niet meer met de hand te doen maar hebben we daar computers voor, waardoor we al snel heel wat stappen met die rechthoeken kunnen maken. Natuurlijk moeten we daarvoor het rekenprocedé opschrijven in een programmeertaal, in dit geval Python. Daarvoor schrijven we de lengte en de breedte van de rechthoek als verhoudingsgetallen. Voor de teller van de lengte gebruiken we de variabele tl (van teller lengte). Evenzo gebruiken we nl voor noemer lengte, tb voor teller breedte en nb voor noemer breedte. De eerste rechthoek kunnen we dus in Python beschrijven met de volgende toekenning aan de variabelen:

tl = 2
nl = 1
tb = 1
nb = 1

Voor de programmeertekst gebruiken we de volgende layout. Het =-teken is hierbij geen “is gelijk aan”, maar een “wordt gelijk aan”; dat wil zeggen dat de variabele links van het =-teken de waarde krijgt van de uitdrukking die rechts van het =-teken staat. We spreken het =-teken in dit geval uit als “wordt”. In Python kunnen we het bovenstaande ook korter schrijven, namelijk:

tl, nl, tb, nb = 2, 1, 1, 1 

Dit wordt ook wel simultane toekenning genoemd. Nu we waarden hebben voor tl, nl, tb en nb kunnen we het procedé opschrijven om nieuwe waarden te berekenen voor deze variabelen. Daarbij weten we dat de variabelen tl, nl, tb en nb voldoen aan de volgende voorwaarde:

(tl/nl)*(tb/nb) = 2, 

oftewel

tl*tb = 2*nl*nb 

Verder weten we dat (vandaar ook de naamgeving lengte en breedte)

tl/nl > tb/nb 

oftewel

tl*nb > tb*nl 

De uitdrukking tl*tb = 2*nl*nb noemen we de eerste invariant, die we kortweg zullen aanduiden met I1. De uitdrukking tl*nb > tb*nl noemen we de tweede invariant en duiden we aan met I2. We zijn nu zover dat we de stap naar de volgende rechthoek in Python kunnen beschrijven. Daarvoor hebben we hierboven gezien dat de nieuwe lengte van de rechthoek het rekenkundig gemiddelde is van de lengte en de breedte van de vorige rechthoek, dus

tl/nl = (tl/nl + tb/nb)/2

Op deze manier kunnen we echter in Python geen waarden toekennen aan de variabelen tl en nl. Als een tekst nog niet in programmeerbare vorm staat noemen we dat pseudo-programmeertekst. De bovenstaande tekst is dus pseudo-programmeertekst. Overigens, als we met de Pythonbibliotheek “fractions” zouden werken, zouden we wel met een dergelijke toekenning kunnen werken. Python heeft wat dat betreft ongekende mogelijkheden.

We moeten er dus nu nog voor zorgen dat in het rechterdeel van de toekenning een uitdrukking staat met een teller en een noemer, want zoals hierboven weergegeven kunnen we dit niet gebruiken. Door herschrijven komen we tot

tl/nl = (tl*nb+tb*nl)/(2*nl*nb)

Ook dit is nog steeds pseudo-programmeertekst. De nieuwe waarden van tl en nl kunnen we echter toekennen volgens

tl, nl = tl*nb+tb*nl, 2*nl*nb

In deze vorm is het bruikbare programmeertekst geworden, en dus geen pseudo programmeertekst meer. In dit geval zouden we ook mogen schrijven

tl = tl*nb+tb*nl
nl = 2*nl*nb

Op de eerste plaats vinden we de simultane toekenning veel mooier, en op de tweede plaats kan dit hier alleen omdat de nieuwe waarde van tl niet voorkomt in de uitdrukking voor nl. Vervolgens willen we de eerste invariant, I1, behouden, zodat

tb/nb = 2/(tl/nl)

oftewel

tb = 2*nl
nb = tl

Ook deze twee kunnen we simultaan toekennen:

tb, nb = 2*nl, tl

Merk hierbij op dat we niet alle vier de variabelen tegelijk mogen toekennen, want voor de breedte (tb en nb) gebruiken we de nieuwe waarden van de lengte (tl en nl). De stap om tot de waarden van de volgende rechthoek te komen kunnen we nu in Python als volgt beschrijven:

Het #-teken heeft een bijzondere rol in Python. Het #-teken zelf en de tekst die op die regel volgt na het #-teken zijn commentaar, of annotaties: zij helpen ons bij herlezen de werking van het programma te begrijpen.

Zouden we nu de waarde van $\sqrt2$ precies willen bepalen, dan zouden we het volledige programma als volgt kunnen beschrijven:

Dit programma zal echter nooit stoppen. Het programma zou namelijk pas stoppen als de ontkenning van tl*nb > tb*nl waar zou zijn, dus tl*nb <= tb*nl, wat volgens de invariant I2 onmogelijk is. In principe zouden we de voorwaarde tl*nb <= tb*nl mogen versterken door te stellen dat tl*nb = tb*nl, maar dat zal nooit gebeuren. Zou dit wel mogelijk zijn, dan zou dit namelijk betekenen dat we een vierkant hebben gevonden met zijden die een rationaal getal zijn, oftewel dat $\sqrt2$ een rationaal getal is. We weten echter dat dit niet zo is. Daarom introduceren we maar een tellertje (t), zodat we een eindig aantal rechthoeken berekenen. We houden het maar even bij $10$ opeenvolgende rechthoeken, maar als je er meer wilt kan dat natuurlijk ook.

Als je dit programma laat uitvoeren zul je zien dat de getallen al gauw heel groot worden. De kracht van Python is onder meer dat het kan rekenen met willekeurig grote gehele getallen. Als we kijken naar de gevonden waarden voor de teller en de noemer van de lengte van de opeenvolgende rechthoeken, dan vinden we dezelfde getallen terug in de OEIS als de reeksen A001601 en A051009. Uit de naamgeving van A001601 blijkt dat we deze waarden ook hadden kunnen vinden door de teller van de eerste nieuw bepaalde rechthoek gelijk te stellen aan $3$, en vervolgens de volgende waarde te bepalen door tl = 2*tl**2-1, maar dit geheel terzijde. We geven dit alleen maar aan om aan te tonen dat je in de OEIS ook vaker alternatieve methoden kunt vinden om iets te berekenen. In sommige gevallen kan dat heel handig zijn. 

Uiteindelijk wilden we niet een breuk als benadering voor $\sqrt2$, maar een decimale representatie. Het getal voor de decimale komma kunnen we heel eenvoudig berekenen door te bepalen hoe vaak de noemer in de teller gaat. Dit gaat heel eenvoudig door geheeltallig te delen en er alvast de decimale komma achter te plaatsen:

print(tl//nl, “,”)

Om vervolgens te bepalen hoeveel decimale cijfers correct zijn, bepalen we de decimale representatie van zowel de lengte als van de breedte, en zolang die beide niet verschillen, gebruiken we de berekende cijfers als resultaat van de berekening.

We gaan vervolgens via een while-loop telkens een volgend cijfer bepalen. Dat doen we zo lang als de cijfers uit de lengte hetzelfde zijn als de cijfers uit de breedte van de bepaalde rechthoek. Maar ook nu moeten we een invariant voor die berekening bepalen. Laten we daarom eens kijken hoe we dat kunnen doen, door eerst de eerste stap te bekijken.

In de eerste stap bepalen we de cijfers voor de decimale komma door tl//nl te berekenen. Dat is iets anders dan de “gewone” deling $tl/nl$. In de wiskunde wordt de uitkomst van tl//nl veelal beschreven door $\lfloor tl/nl\rfloor$, wat we uitspreken als “vloer van tl/nl”. De vloer van $tl/nl$ is het grootste gehele getal $x$ waarvoor geldt dat $x = tl/nl$. De notatie $\lfloor tl/nl \rfloor$ is echter niet zo geschikt om dat als commentaar aan een programma toe te voegen. De bekende Nederlandse informaticus Edsger Dijkstra, uit de beginjaren van de hedendaagse informatica, introduceerde daarom de operator div, afgeleid van het Engelse woord “divide”, dus tl div nl is hetzelfde als $\lfloor tl/nl\rfloor$, en programmeren we in Python als tl//nl. Nu geldt echter dat (tl//nl)*nl <= tl, en wel zodanig dat 0 <= tl-(tl//nl)*nl < nl. Maar laten we dit wat algemener beschrijven met behulp van de getallen $x$ en $y$ en de operator $\mbox{ div }$. Er geldt dan dat $0 = x - (x \mbox{ div } y) \times y < y$.

Het verschil tussen $x$ en $(x \mbox{ div } y) \times y$ schrijft Dijkstra als “$x \mbox{ mod } y$”, waarbij $x \mbox{ mod } y$ de restwaarde is van de geheeltallige deling. Dat verklaart meteen waarom $x - (x \mbox{ div } y) \times y < y$ is. Als we de waarde van de geheeltallige deling $(x \mbox{ div } y)$ weer vermenigvuldigen met de noemer $(y)$, en er vervolgens de restwaarde $(x \mbox{ mod } y)$ bij optellen, krijgen we weer de oorspronkelijke waarde van de teller $(x)$ terug. Dus geldt er 

$$(x \mbox{ div } y) \times y + x \mbox{ mod } y = x$$

en voor de restwaarde $x \mbox{ mod } y$

$$0 = x \mbox{ mod } y < y.$$

In Python hebben we zowel operaties voor $x \mbox{ div } y$ als voor $x \mbox{ mod } y$. We zagen al dat $x \mbox{ div } y$ in Python geschreven wordt als x//y; $x \mbox{ mod } y$ wordt in Python geschreven als x%y. Misschien vraag je je af waarom we zowel “$x \mbox{ mod } y$” als “x%y” gebruiken. Het antwoord is heel simpel, de taal van de wiskunde is een andere taal dan de taal Python. De tweede invariant in het bovenstaande programma hadden we daarom wat sterker kunnen formuleren door te stellen dat

tl/nl > sqrt(2) > tb/nb

De waarheid van deze invariant volgt direct uit het feit dat tl/nl > tb/nb en de eerste invariant, (tl/nl)*(tb/nb) = 2. Dus in feite volgt deze sterkere formulering direct uit de twee eerdergenoemde invarianten. Dit schrijven we als volgt:

I1, I2 ==> tl/nl > sqrt(2) > tb/nb

te lezen als: uit invarianten I1 en I2 volgt dat tl/nl > sqrt(2) > tb/nb. Laten we terugkeren naar ons probleem om een groot aantal decimale cijfers te bepalen van $\sqrt2$. Het getal voor de komma (in dit geval slechts één cijfer) konden we bepalen door tl//nl te berekenen, zodat de cijfers achter de komma besloten liggen in de waarde van tl%nl. De restwaarde van $\sqrt2$ voor de cijfers achter de komma wordt dus bepaald door (tl mod nl) / nl. We beginnen daarom met de toekenning

tlr = tl%nl

waarbij tlr staat voor “teller lengte rest”. Merk verder op dat de noemer hetzelfde blijft.

Het eerstvolgende decimale cijfer kunnen we nu eenvoudig bepalen door tlr met $10$ te vermenigvuldigen en dan de geheeltallige deling door de noemer van de lengte (nl) te herhalen, en dan opnieuw de restwaarde voor de volgende decimale cijfers te bepalen:

cijfer = (10*tlr)//nl tlr = (10*tlr)%nl

Uiteindelijk gaan we dit niet alleen doen voor de lengte van de rechthoek in kwestie, maar ook voor de breedte. Zolang beide decimale representaties hetzelfde zijn, zullen de gevonden cijfers juist zijn. Voor beide cijfers zullen we de variabelen cl (voor cijfer lengte) en cb (voor cijfer breedte) gebruiken. 

We hebben hierbij alvast de eerste regel van de uitvoer bepaald met behulp van het string-type. De operatie str(cl) zet de waarde van cl om in karakters die het getal cl representeren in het decimale stelsel. Vervolgens kunnen we daar de decimale komma als tekst aanplakken. Dat “aanplakken” wordt concateneren genoemd, en de operatie om strings te concateneren is in Python de +. De cijfers achter de decimale komma hebben we in het bovenstaande programma wel al bepaald, maar die hebben we nog niet zichtbaar gemaakt in de uitvoer van het programma.

Om de resultaten leesbaar te houden, zullen we per regel $50$ cijfers achter de komma weergeven, met voor het leesgemak per regel het totaal aantal cijfers dat we tot daartoe hebben bepaald, achter een dubbele punt. Om te weten hoeveel cijfers we hebben zullen we wederom een tellertje moeten introduceren. We gebruiken hiervoor ct (voor cijfer teller), en maken meteen de grote stap die de beoogde while-loop weergeeft

De laatste stap die we zetten is dat we de berekening van de decimale representaties integreren in het oorspronkelijke programma. Als je het programma zelf uitvoert, zul je zien dat bij elke volgende rechthoek die we berekenen, het aantal correcte decimale cijfers achter de komma ongeveer verdubbelt. In het uiteindelijke programma hebben we aan de programmatekst ook commentaar meegegeven voor de while-loop die de decimale cijfers bepaalt. Deze “invarianten” zijn niet zo wiskundig opgeschreven als de eerdere invarianten, maar zijn wel voldoende om bij herlezing van de programmatekst te begrijpen wat er in die while-loop gebeurt.

Om af te sluiten geven we jullie nog wat opdrachten mee.

Veel plezier met programmeren, en bedenk dat je nu al aardig kunt rekenen als een Mesopotamiër. De Mesopotamiërs deden dat al ongeveer drie en een half duizend jaar geleden, maar met je computer kun je in een paar seconden net zoveel berekenen als de beste Mesopotamische rekenaars in een heel leven konden.

OPDRACHT 1

In A002193 van de OEIS kunnen we de eerste 20000 cijfers van $\sqrt2$ terugvinden. Controleer de uitkomsten die je gevonden hebt met die van de OEIS. Probeer eens zelf meer dan 20000 cijfers van $\sqrt2$ te bepalen. In plaats van $\sqrt2$ te berekenen, kun je het programma ook gebruiken om andere wortels te berekenen, bijvoorbeeld $\sqrt{\frac23}$, of een andere wortel; naar eigen keuze. Probeer zelf het programma zodanig aan te passen dat je de betreffende wortel kunt uitrekenen. Van sommige wortels zijn de cijfers ook te vinden in de OEIS, dus je kunt je eigen resultaten controleren.

OPDRACHT 2

We hebben het programma zo geschreven dat het de decimale representatie van $\sqrt2$ berekent. Kun je het programma zodanig veranderen dat het de representatie weergeeft in een grondtal van de getallenrepresentatie waarbij het grondtal hoogstens $10$ is, bijvoorbeeld binair (grondtal $2$) of octaal (grondtal $8$)?

Voor grondtallen die groter zijn dan $10$, zul je wellicht in eerste instantie een probleem tegenkomen. Allereerst, hoe geef je cijfers van 10 of meer weer? Voor het hexadecimaal getalstelsel (grondtal $16$) kun je de daarvoor gebruikelijke $10$ cijfers en de eerste zes letters van het alfabet gebruiken. Het is dan handig om eerst de string “cijfers” toe te kennen:
cijfers = "0123456789ABCDEF" Om het cijfer te bepalen gebruik je nu niet str(cl), maar cijfers[cl], waarmee je het juiste symbool voor het betreffende cijfer uit de string cijfers selecteert. De Mesopotamiërs gebruikten als grondtal zestig. Zulk een getallenstelsel noemen we een sexagesimaal getallenstelsel. Voor het sexagesimaal getallenstelsel lukt het niet meer om na het cijfer 9 letters te gaan gebruiken, omdat het alfabet daarvoor gewoonweg te weinig letters heeft. De oplossing is nu om de cijfers weer te geven met twee cijfers uit ons decimale getallenstelsel. Dit komt overigens sterk overeen met de manier waarop de Mesopotamiërs dat ook deden. Om de sexagesimaal cijfers goed van elkaar te onderscheiden is het nu verstandig er een spatie tussen te zetten. Je zou op die manier soortgelijk te werk kunnen gaan als bij het voorbeeld van de hexadecimale getallenrepresentatie. Een andere manier is om eerst het sexagesimaal cijfer cl te bepalen en vervolgens in plaats van s = s+str(cl) te schrijven s = s+str(cl//10)+str(cl%10)+str(" "). Probeer zelf beide methodes en kijk welke je het makkelijkst vindt.