Een raytracer schrijven in Brainfuck

"Vaak is het juiste antwoord om een tool te schrijven in een algemene programmeertaal die het probleem oplost, en CMake te leren hoe hij die tool moet aanroepen als onderdeel van het buildproces. Code-generatie, cryptografische handtekening-utiliteiten en zelfs raytracers zijn geschreven in de CMake-taal, maar dit is geen aanbevolen praktijk."

Omdat ik eerder al een raytracer had geschreven (en deze later had herschreven voor de GPU), wekte deze uitspraak mijn nieuwsgierigheid. Ik vroeg me af in welke taal het nóg uitdagender zou zijn om een raytracer te schrijven.

De vorige herschrijving was gebaseerd op code die weinig gebruikmaakte van basisprincipes en vooral afhankelijk was van een relatief complex geheel aan API's. Daarom koos ik voor de simpelste taal die ik kende: Brainfuck (BF). Een eenvoudige taal resulteert immers logischerwijs in een zeer eenvoudige codebase; BF-programma's zijn vaak slechts enkele regels lang. Bovendien zette een opmerking in de README van Muller me aan om een tegenvoorbeeld te laten zien.

De code is beschikbaar op mTvare6/rayfuck.

Inleiding tot Brainfuck

BF is een beslist eenvoudige taal die slechts uit acht operaties bestaat en één "datastructuur" gebruikt: een oneindig lint van cellen, waarbij elke cel een u8 (unsigned 8-bit integer) kan opslaan.

• Bij het teken > beweegt de datapointed (die naar een cel op het lint wijst) naar rechts, en bij < naar links.
• I/O wordt beheerd via , en .. De eerste slaat het input-byte op in de cel waar de datapointed naar wijst, en de tweede print deze waarde uit.
• De enige primitieven naast I/O waarmee een waarde kan worden gewijzigd, zijn het incrementen (+) en decrementen (-) van de waarde op de huidige positie van de datapointed.

De beperkingen zijn duidelijk: er zijn geen registers ($n > 1$) zoals andere machines hebben, geen instructies die op meer dan één cel tegelijk opereren, en geen instructies voor optelling of vermenigvuldiging.

Het laatste ingrediënt om BF Turing-compleet te maken is de loop, geschreven met [ en ]. Wanneer het token [ wordt aangetroffen, controleert de runtime de cel bij de datapointed: als deze nul is, springt de uitvoering naar het bijbehorende ]. Anders wordt de loop betreden. Bij ] keert de uitvoering terug naar het bijbehorende [ als de cel niet-nul is, en verlaat hij de loop als de cel nul is.

Een goede oefening is om een 'cat'-programma te schrijven; probeer er een te maken met slechts vijf tekens om een gevoel te krijgen voor de omgeving.

Voorbereiding

Voordat ik begon, besloot ik om geen resultaten of implementatiegegevens op te zoeken en zoveel mogelijk vanuit basisprincipes op te bouwen. Om de scope minimaal te houden en het programma een onmiskenbare raytracer te maken, besloot ik exact hetzelfde beeld te renderen als in de Metal-sectie van RIW.

De C-code was echter te complex, en het schrijven van een C-parser lag duidelijk buiten de scope. Ik besloot dat elk double (en andere datatypes zoals bool) zou worden gerepresenteerd door het combineren van cellen, waarbij de helft van de bits het fractionele deel vertegenwoordigt en de andere helft het gehele getal. Dit plaatst effectief een vast binair punt tussen beide delen, wat bekend staat als een Q-formaat.

Een signed Q8.8 zou een resolutie van $1/256$ en een bereik van ongeveer $[-128, 128)$ geven. Dat was onvoldoende, aangezien de bol die als grond in de scène dient een straal van $r=1000$ moest hebben om vlak te lijken. Daarom koos ik voor het duurdere signed Q16.16-formaat. Dit heeft een resolutie van $1/2^{16}$ en een bereik van $[-2^{15}, 2^{15})$, wat voldoende is.

Ik besloot de code om te zetten naar een SSA-achtig formaat (Static Single Assignment), waarbij recursieve code iteratief werd gemaakt en variabelen in functies werden voorzien van een Hongaarse notatie om naamconflicten tijdens het opzoeken van adressen te voorkomen.

Tevens was het noodzakelijk om het parsen en de codegeneratie te scheiden. De complexiteit werd verdeeld over twee code-regio's, met een tussenliggende "DSL" als Intermediate Representation (IR). Deze DSL bevatte eenvoudige operaties zoals: abs, add, and, call, copy, div, else, end, eq, func, ge, gt, if, int, le, lt, mul, neg, not, or, print2, print3, set, sqrt, sub, text, var en while.

Het volgende lastige punt waren enkele bibliotheekaanroepen: sqrt, rand en abs. Voor rand overwoog ik aanvankelijk een twee-status oplossing, maar de meeste varianten hebben een slechte periode. Ik koos uiteindelijk voor de eenvoudigste versie: A = (5 * A + 1) % 256 Dit garandeert dat de reeks pas herhaalt na een volledige sequentie van 256 waarden, wat acceptabel is voor dit gebruik (supersampling anti-aliasing).

Voor sqrt was de formule van Heron een kandidaat, maar die zou inefficiënt zijn vanwege de benodigde deling. Herhaalde aftrekking zou minder gegenereerde code opleveren (wat beter is, omdat de interpreter minder hoeft te bewegen), maar bleef relatief duur.

Andere kandidaten waren de Taylor-reeks en de "School Method" (staartdeling). Voor de Taylor-reeks geldt: $\sqrt{1 + x} = 1 + \frac{x}{2} - \frac{x^2}{8}$ $y = 2^{16} \cdot x$ $\frac{\sqrt{2^{16} + y}}{2^8} = (1 + \frac{y}{2^{17}} - \frac{y^2}{2^{35}})$ $\frac{\sqrt{y}}{2^8} = (1 + \frac{y - 2^{16}}{2^{17}} - \frac{(y - 2^{16})^2}{2^{35}})$

Een plot in Desmos liet zien dat de fit slecht was onder een gecodeerde waarde van $20\text{k}$ (ongeveer $0,305$), wat een belangrijke regio was. Dit liet me over aan de long-division methode. Als de echte waarde $x$ is, dan is de gerepresenteerde waarde: $N = x \cdot 2^{16}$ Om $\sqrt{x}$ te representeren, hebben we nodig: $N' = \sqrt{x} \cdot 2^{16}$ $\text{isqrt}(N) = \sqrt{x} \cdot 2^8$ $\text{isqrt}(2^{16} \cdot N) = \sqrt{x} \cdot 2^{16} = N'$

Hier is isqrt gerechtvaardigd, omdat een verschil van één in het gecodeerde resultaat de gedecodeerde vierkantswortel met minder dan $1/2^{16}$ (ongeveer $0,00001526$) verandert.

Ten slotte zou de volledige variabelenkaart en de bijbehorende BF-adressen worden beheerd met een dictionary.

Implementatie

Met de theoretische basis op orde, bleven alleen de implementatiedetails over. Twee belangrijke primitieven waren move en copy.

Move en Copy

Een move operatie werkt door een waarde continu te verlagen tot de cel bij de initiële datapointed nul wordt, terwijl de andere cel in gelijke mate wordt verhoogd. Als one-liner: [->+<]

Een copy operatie werkt als volgt: uitgaande van het array [a, 0, 0], gebruikt de code [->+>+<<] om dit om te zetten naar [0, a, a]. Indien nodig kan de laatste a vervolgens weer naar binnen worden verplaatst.

Omdat optellingen en later delingen herhaald gebruik maken van tijdelijke waarden, heeft elke waarde in de map nabijgelegen slots voor tijdelijke variabelen. Dit voorkomt dat de datapointed te veel heen en weer moet bewegen en biedt ruimte voor carry-cellen.

Vermenigvuldiging

Vermenigvuldiging is relatief eenvoudig en bestaat uit het vermenigvuldigen van elke cel, het opslaan van de resultaten en het later bij elkaar optellen. De vermenigvuldigingsstap wordt uitgevoerd via herhaalde optelling.

Voor het vermenigvuldigen van twee cellen wordt één cel gekopieerd naar een tijdelijke plek als outer loop, en de andere wordt bij elke iteratie gekopieerd om als inner loop te dienen. Elke iteratie van de inner loop verhoogt het resultaat met één.

Voor de waarden van vier cellen wordt elke cel van de ene set gepaard met elke cel van de andere. Een vermenigvuldiging van de cellen op positie $i$ en $j$ wordt toegevoegd aan $i+j$ in een resultaat van acht cellen. [a0, a1, a2, a3] * [b0, b1, b2, b3] $\text{result}[i+j] += ai \cdot bj$

Omdat beide inputs al $2^{16}$ in hun representatie hadden, worden de laagste twee cellen weggegooid bij het terugkopiëren van het resultaat: $N1 = x1 \cdot 2^{16}$ $N2 = x2 \cdot 2^{16}$ $(N1 \cdot N2) / 2^{16} = N = (x1 \cdot x2) \cdot 2^{16}$

Deling

Deling is minder direct, maar kan worden uitgevoerd via de methode van staartdeling. De dividend wordt gelezen vanaf de meest significante cel, en bij elke stap wordt het oude restant meegenomen naar de volgende cel. R = R * 256 + A[next]

De divisor wordt vervolgens herhaaldelijk van dit restant afgetrokken, waarna er één wordt toegevoegd aan de resultaatcel. Wanneer het restant negatief wordt, wordt de stap teruggedraaid en gaan we naar de volgende cel.

while R >= D:
    R -= D
    result += 1

Dit vereist maximaal 255 aftrekkingen per cel.

Net zoals de representatie bij vermenigvuldiging naar rechts wordt verschoven, verliest deling informatie door de verschuiving naar links. Daarom is er een verschuiving in de representatie nodig vóór het delen: $N1 = x1 \cdot 2^{16}$ $N2 = x2 \cdot 2^{16}$ $(N1 / N2) \cdot 2^{16} = (x1 / x2) \cdot 2^{16}$ (bits gaan verloren) $(N1 \cdot 2^{16}) / N2 = (x1 / x2) \cdot 2^{16}$

Bij al deze operaties worden eerst de tekens verwijderd, waarna het resultaat negatief wordt gemaakt als slechts één van de inputs negatief was.

Vergelijkingen en Booleaanse logica

Vergelijkingen gebruiken een kleinere operatie: twee cellen worden samen gedecimenteerd tot er minstens één nul wordt, totdat er een verschil is of de kopieën volledig nul zijn. Er is een kleine verwerking nodig waarbij $2^7$ wordt toegevoegd aan de meest significante cel; anders hebben negatieve getallen technisch gezien een hogere waarde wanneer ze puur als bytes worden bekeken:

• 00 ... 7f → positieve helft
• 80 ... ff → negatieve helft
• Na toevoeging van 128 en wrapping: 80 ... ff → positieve helft, 00 ... 7f → negatieve helft.

Booleaanse checks lezen de cellen en zetten de output op één als een van de cellen niet-nul is (om de "truthiness" van C te imiteren). In de representatie staat de laagste bit voor true en alle bits nul voor false. Operaties zoals and, or en not werken met deze bit-representatie.

Negatie maakt gebruik van de truc $-x = \sim x + 1$. Elke cel $x$ wordt gecomplementeerd (via $255-x$), waarna één wordt toegevoegd aan de laagste cel met carry-over. De abs-functie controleert alleen de hoogste bit van de laatste cel en voert deze negatie uit als deze is ingesteld.

Functies en Control Flow

Vanwege de beslissing om variabelen per functie te scopen en SSA-stijl code te gebruiken, waren functies zeer eenvoudig. De codegenerator noteert het functielichaam onder de naam en plaatst dit inline wanneer call func wordt aangetroffen.

if-statements zijn eenvoudig door het gebruik van loops: [- body ] Voor een else wordt een andere vlag op één gezet en gewist door het eerste body-blok. while-loops werken op soortgelijke wijze:

condition
[
  body
  move to condition and calculate
]

Het Resultaat

Het uiteindelijke programma was 23MB, wat groter is dan de afbeelding zelf (ongeveer 0,9MB). Dit maakt het een vrij slechte keuze als compressietechniek.

Uit ruwe berekeningen bleek dat het programma 100 ray-berekeningen per minuut deed, oftewel één pixel per minuut. Gezien de afbeelding $400 \times 225$ pixels is, zou de initiële schatting ongeveer 62,5 dagen op mijn laptop zijn zonder verdere optimalisaties. Ik realiseerde me echter dat ik alleen de lucht had gezien; het bouncen rond de bollen zou de ETA aanzienlijk verlengen. De afbeelding in dit artikel is daarom een benadering van wat er gerenderd zou worden, gemaakt met de C-code. Van de 1229 (van de 90.000) pixels die op het moment van schrijven waren gegenereerd, verschillen er slechts 10, meestal met een waarde van één.

Sommige optimalisaties zijn mogelijk. Het inleveren van precisie om het aantal geraakte cellen te verminderen is een optie als de grond minder nauwkeurig benaderd kan worden. Ook de normalisatiestap voor willekeurige vectoren kan worden overgeslagen, hoewel dat de scatter-distributie verandert en het programma niet langer exact de code zou uitvoeren die ik probeerde te reproduceren.

Update

Een reactie op mijn Reddit-thread vroeg hoe ik het zou kunnen verbeteren met fork/join primitieven. Omdat ik deze uitdaging interessant vond, ben ik mijn eigen JIT-interpreter gaan verbeteren, wat leidde tot een enorme prestatieverbetering. De uiteindelijke render ziet er nu een beetje uit als een Van Gogh-schilderij, waarschijnlijk door precisiefouten.