/* Draw a pythagoras triangle with no trig calls */
#include <math.h>
#include <stdio.h>
#include <stdint.h>
#include <assert.h>

#if 1
// Implementation for 16.16 fixed point
#include <inttypes.h> // PRId32
typedef int32_t F; // or long if stdint not available
#define M(name) (name ## fix)
#define ZERO (0)
#define ONE (INT32_C(1) << 16)
#define PRIF(a, b) #a PRId32
#define MUL(a,b) mulfix((a), (b))
#define RECIP(x) ((0x80000000 / (x)) << 1)
#define PROLOGUE "1 65536 div dup scale"

// From "Hacker's Delight" figure 11-1 adapted for 16.16 fixed point
int32_t mulfix(int32_t u, int32_t v) {
    uint32_t u0, v0, w0;
    int32_t u1, v1, w1, w2, t;
    u0 = u & 0xFFFF; u1 = u >> 16;
    v0 = v & 0xFFFF; v1 = v >> 16;
    w0 = u0*v0;
    t = u1*v0 + (w0 >> 16);
    w1 = t & 0xFFFF;
    w2 = t >> 16;
    w1 = u0*v1 + w1;
    return ((u1*v1 + w2) << 16) + w1;
}


// From https://github.com/chmike/fpsqrt/blob/master/fpsqrt.c MIT license
F sqrtfix(F v) {
    uint32_t t, q, b, r;
    r = (int32_t)v; 
    q = 0;          
    b = 0x40000000UL;
    if( r < 0x40000200 )
    {
        while( b != 0x40 )
        {
            t = q + b;
            if( r >= t )
            {
                r -= t;
                q = t + b; // equivalent to q += 2*b
            }
            r <<= 1;
            b >>= 1;
        }
        q >>= 8;
        return q;
    }
    while( b > 0x40 )
    {
        t = q + b;
        if( r >= t )
        {
            r -= t;
            q = t + b; // equivalent to q += 2*b
        }
        if( (r & 0x80000000) != 0 )
        {
            q >>= 1;
            b >>= 1;
            r >>= 1;
            while( b > 0x20 )
            {
                t = q + b;
                if( r >= t )
                {
                    r -= t;
                    q = t + b;
                }
                r <<= 1;
                b >>= 1;
            }
            q >>= 7;
            return q;
        }
        r <<= 1;
        b >>= 1;
    }
    q >>= 8;
    return q;
}
#else
typedef float F;
#define M(name) (name ## fix)
#define ZERO (0.f)
#define ONE (1.f)
#define PRIF(a, b) #a "." #b "f"
#define MUL(a,b) (a*b)
#define RECIP(x) (1.f/x)
#define PROLOGUE ""
#endif

#define MOVETO "moveto"
#define LINETO "lineto"
#define LINETO_CLOSEPATH_STROKE "lineto closepath stroke"

typedef struct { F a, b, c, d, e, f; } PSMatrix;

void matmul(PSMatrix *result, PSMatrix *m1, PSMatrix *m2) {
    PSMatrix r;
    r.a = MUL(m2->a, m1->a) + MUL(m2->b, m1->c);
    r.b = MUL(m2->a, m1->b) + MUL(m2->b, m1->d);
    r.c = MUL(m2->c, m1->a) + MUL(m2->d, m1->c);
    r.d = MUL(m2->c, m1->b) + MUL(m2->d, m1->d);
    r.e = MUL(m2->e, m1->a) + MUL(m2->f, m1->c) + m1->e;
    r.f = MUL(m2->e, m1->b) + MUL(m2->f, m1->d) + m1->f;
    *result = r;
}

F A = 4 * ONE, B = 3 * ONE, C, C_recip;
F cos_theta, sin_theta;

// Each matrix moves from the old baseline [bottom of square] to the baseline of the square on the B or C side
// of the triangle.
PSMatrix toB, toC, identity = {ONE, ZERO, ZERO, ONE, ZERO, ZERO};

void make_matrix(PSMatrix *m, F x1, F y1, F x2, F y2) {
    F dx = x2-x1, dy = y2-y1;
    F ux = MUL(dx, C_recip);
    F uy = MUL(dy, C_recip);

    PSMatrix mm = { ux, uy, -uy, ux, x1, y1 };

    *m = mm;
}

void op_transformed(const char *op, PSMatrix *m, F x, F y) {
    F tx = MUL(m->a, x) + MUL(m->c, y) + m->e,
      ty = MUL(m->b, x) + MUL(m->d, y) + m->f;
    printf("%" PRIF(,1) " %" PRIF(,1) " %s\n", tx, ty, op);
}

void recurse(PSMatrix *m, int n) {
    op_transformed(MOVETO, m, 0, 0);
    op_transformed(LINETO, m, 0, C);
    op_transformed(LINETO, m, C, C);
    op_transformed(LINETO_CLOSEPATH_STROKE, m, C, 0);

    if (n > 0) {
        PSMatrix m1;
        matmul(&m1, m, &toB);
        recurse(&m1, n-1);
        matmul(&m1, m, &toC);
        recurse(&m1, n-1);
    }

}

int main() {
    C = M(sqrt)(MUL(A, A) + MUL(B, B));
    C_recip = RECIP(C);
    F dy = MUL(MUL(A, B), C_recip);
    F dx = M(sqrt)(MUL(A, A) - MUL(dy, dy));
    make_matrix(&toB, 0, C, dx, C+dy);
    make_matrix(&toC, dx, C+dy, C, C);

    fprintf(stderr, "% " PRIF(6,2) " % " PRIF(6,2) " % " PRIF(6,2) "\n% " PRIF(6,2) " % " PRIF(6,2) " % " PRIF(6,2) "\n\n",
            toB.a, toB.b, toB.c,
            toB.d, toB.e, toB.f);

    PSMatrix initial = {ONE * 72 / 8, ZERO, ZERO, ONE*72 / 8, ONE*72*3, ONE*8};

    printf("%s\n", PROLOGUE);
    recurse(&initial, 6);
    printf("showpage\n");
}