math.kai

// TODO: @overloading
#foreign llvm #callconv "c" {

    pow  :: fn(x, power: f32) -> f32                                                         #linkname "llvm.pow.f32"
//  pow  :: fn(x, power: f64) -> f64                                                         #linkname "llvm.pow.f64"

//  pow  :: fn(x: f32, power: i32) -> f32                                                    #linkname "llvm.powi.f32"
//  pow  :: fn(x: f64, power: i32) -> f32                                                    #linkname "llvm.powi.f64"

    sqrt :: fn(x: f32) -> f32                                                                #linkname "llvm.sqrt.f32"
//  sqrt :: fn(x: f64) -> f64                                                                #linkname "llvm.sqrt.f64"

    exp  :: fn(x: f32) -> f32                                                                #linkname "llvm.exp.f32"
//  exp  :: fn(x: f64) -> f64                                                                #linkname "llvm.exp.f64"

    exp2  :: fn(x: f32) -> f32                                                               #linkname "llvm.exp2.f32"
//  exp2  :: fn(x: f64) -> f64                                                               #linkname "llvm.exp2.f64"

    log  :: fn(x: f32) -> f32                                                                #linkname "llvm.log.f32"
//  log  :: fn(x: f64) -> f64                                                                #linkname "llvm.log.f64"

    log10  :: fn(x: f32) -> f32                                                              #linkname "llvm.log10.f32"
//  log10  :: fn(x: f64) -> f64                                                              #linkname "llvm.log10.f64"

    log2  :: fn(x: f32) -> f32                                                               #linkname "llvm.log2.f32"
//  log2  :: fn(x: f64) -> f64                                                               #linkname "llvm.log2.f64"

    sin  :: fn(x: f32) -> f32                                                                #linkname "llvm.sin.f32"
//  sin  :: fn(x: f64) -> f64                                                                #linkname "llvm.sin.f64"

    cos  :: fn(x: f32) -> f32                                                                #linkname "llvm.cos.f32"
//  cos  :: fn(x: f64) -> f64                                                                #linkname "llvm.cos.f64"

    fma :: fn(a, b, c: f32) -> f32                                                           #linkname "llvm.fma.f32"
//  fma :: fn(a, b, c: f64) -> f64                                                           #linkname "llvm.fma.f64"

    fabs :: fn(x: f32) -> f32                                                                #linkname "llvm.fabs.f32"
//  fabs :: fn(x: f64) -> f64                                                                #linkname "llvm.fabs.f64"

    min :: fn(x: f32) -> f32                                                                 #linkname "llvm.minnum.f32"
//  min :: fn(x: f64) -> f64                                                                 #linkname "llvm.minnum.f64"

    max :: fn(x: f32) -> f32                                                                 #linkname "llvm.maxnum.f32"
//  max :: fn(x: f64) -> f64                                                                 #linkname "llvm.maxnum.f64"

    sign :: fn(x: f32) -> f32                                                                #linkname "llvm.copysign.f32"
//  sign :: fn(x: f64) -> f64                                                                #linkname "llvm.copysign.f64"

    floor :: fn(x: f32) -> f32                                                               #linkname "llvm.floor.f32"
//  floor :: fn(x: f64) -> f64                                                               #linkname "llvm.floor.f64"

    ceil :: fn(x: f32) -> f32                                                                #linkname "llvm.ceil.f32"
//  ceil :: fn(x: f64) -> f64                                                                #linkname "llvm.ceil.f64"

    trunc :: fn(x: f32) -> f32                                                               #linkname "llvm.trunc.f32"
//  trunc :: fn(x: f64) -> f64                                                               #linkname "llvm.trunc.f64"

    rint :: fn(x: f32) -> f32                                                                #linkname "llvm.rint.f32"
//  rint :: fn(x: f64) -> f64                                                                #linkname "llvm.rint.f64"

    nearbyint :: fn(x: f32) -> f32                                                           #linkname "llvm.nearbyint.f32"
//  nearbyint :: fn(x: f64) -> f64                                                           #linkname "llvm.nearbyint.f64"

    round :: fn(x: f32) -> f32                                                               #linkname "llvm.round.f32"
//  round :: fn(x: f64) -> f64                                                               #linkname "llvm.round.f64"

    // MARK: Specialized Arithmetic Intrinsics

    canonicalize :: fn(f32) -> f32                                                           #linkname "llvm.canonicalize.f32"
//  canonicalize :: fn(f64) -> f64                                                           #linkname "llvm.canonicalize.f64"

    fmuladd32  :: fn(a, b, c: f32) -> f32                                                    #linkname "llvm.fmuladd.f32"
//  fmuladd64  :: fn(a, b, c: f64) -> f64                                                    #linkname "llvm.fmuladd.f64"

    // TODO: Arithmetic with Overflow Intrinsics
}

tau          :: 6.28318530717958647692528676655900576
pi           :: 3.14159265358979323846264338327950288

e            :: 2.71828182845904523536

epsilon      :: 1.19209290e-7

Vec3 :: #vector(3)f32
Vec4 :: #vector(4)f32
Mat4 :: [4][4]f32

tan  :: fn(x: f32) -> f32 { return sin(x) / cos(x) }

dot :: fn(a, b: Vec3) -> f32 {
    c := a * b
    return c.x + c.y + c.z
}

dot_4 :: fn(a, b: Vec4) -> f32 {
    c := a * b
    return c.x + c.y + c.z + c.w
}

mag :: fn(v: Vec3) -> f32 {
    return sqrt(dot(v, v))
}

mag_4 :: fn(v: Vec4) -> f32 {
    return sqrt(dot_4(v, v))
}

norm :: fn(v: Vec3) -> Vec3 {
    m := mag(v)
    mv := Vec3 { m, m, m }
    return v / mv
}

norm_4 :: fn(v: Vec4) -> Vec4 {
    m := mag_4(v)
    mv := Vec4 { m, m, m, m }
    return v / mv
}

cross :: fn(x, y: Vec3) -> Vec3 {
    a := x.yzx * y.zxy
    b := x.zxy * y.yzx
    return a - b
}

identity_mat :: fn() -> Mat4 {
    return Mat4 {
        {1, 0, 0, 0},
        {0, 1, 0, 0},
        {0, 0, 1, 0},
        {0, 0, 0, 1}
    }
}

transpose :: fn(m: Mat4) -> Mat4 {
    for j := 0; j < 4; j += 1 {
        for i := 0; i < 4; i += 1 {
            m[i][j], m[j][i] = m[j][i], m[i][j]
        }
    }
    return m
}

mul :: fn(a,b: Mat4) -> Mat4 {
    c: Mat4
    for j := 0; j < 4; j += 1 {
        for i := 0; i < 4; i += 1 {
            c[j][i] = a[0][i] * b[j][0] +
                      a[1][i] * b[j][1] +
                      a[2][i] * b[j][2] +
                      a[3][i] * b[j][3]
        }
    }
    return c
}

inverse :: fn(m: Mat4) -> Mat4 {
    o: Mat4

    sf00 := m[2][2] * m[3][3] - m[3][2] * m[2][3]
    sf01 := m[2][1] * m[3][3] - m[3][1] * m[2][3]
    sf02 := m[2][1] * m[3][2] - m[3][1] * m[2][2]
    sf03 := m[2][0] * m[3][3] - m[3][0] * m[2][3]
    sf04 := m[2][0] * m[3][2] - m[3][0] * m[2][2]
    sf05 := m[2][0] * m[3][1] - m[3][0] * m[2][1]
    sf06 := m[1][2] * m[3][3] - m[3][2] * m[1][3]
    sf07 := m[1][1] * m[3][3] - m[3][1] * m[1][3]
    sf08 := m[1][1] * m[3][2] - m[3][1] * m[1][2]
    sf09 := m[1][0] * m[3][3] - m[3][0] * m[1][3]
    sf10 := m[1][0] * m[3][2] - m[3][0] * m[1][2]
    sf11 := m[1][1] * m[3][3] - m[3][1] * m[1][3]
    sf12 := m[1][0] * m[3][1] - m[3][0] * m[1][1]
    sf13 := m[1][2] * m[2][3] - m[2][2] * m[1][3]
    sf14 := m[1][1] * m[2][3] - m[2][1] * m[1][3]
    sf15 := m[1][1] * m[2][2] - m[2][1] * m[1][2]
    sf16 := m[1][0] * m[2][3] - m[2][0] * m[1][3]
    sf17 := m[1][0] * m[2][2] - m[2][0] * m[1][2]
    sf18 := m[1][0] * m[2][1] - m[2][0] * m[1][1]

    o[0][0] = +(m[1][1] * sf00 - m[1][2] * sf01 + m[1][3] * sf02)
    o[0][1] = -(m[1][0] * sf00 - m[1][2] * sf03 + m[1][3] * sf04)
    o[0][2] = +(m[1][0] * sf01 - m[1][1] * sf03 + m[1][3] * sf05)
    o[0][3] = -(m[1][0] * sf02 - m[1][1] * sf04 + m[1][2] * sf05)
    o[1][0] = -(m[0][1] * sf00 - m[0][2] * sf01 + m[0][3] * sf02)
    o[1][1] = +(m[0][0] * sf00 - m[0][2] * sf03 + m[0][3] * sf04)
    o[1][2] = -(m[0][0] * sf01 - m[0][1] * sf03 + m[0][3] * sf05)
    o[1][3] = +(m[0][0] * sf02 - m[0][1] * sf04 + m[0][2] * sf05)
    o[2][0] = +(m[0][1] * sf06 - m[0][2] * sf07 + m[0][3] * sf08)
    o[2][1] = -(m[0][0] * sf06 - m[0][2] * sf09 + m[0][3] * sf10)
    o[2][2] = +(m[0][0] * sf11 - m[0][1] * sf09 + m[0][3] * sf12)
    o[2][3] = -(m[0][0] * sf08 - m[0][1] * sf10 + m[0][2] * sf12)
    o[3][0] = -(m[0][1] * sf13 - m[0][2] * sf14 + m[0][3] * sf15)
    o[3][1] = +(m[0][0] * sf13 - m[0][2] * sf16 + m[0][3] * sf17)
    o[3][2] = -(m[0][0] * sf14 - m[0][1] * sf16 + m[0][3] * sf18)
    o[3][3] = +(m[0][0] * sf15 - m[0][1] * sf17 + m[0][2] * sf18)

    ood := 1.0 / (m[0][0] * o[0][0] +
                  m[0][1] * o[0][1] +
                  m[0][2] * o[0][2] +
                  m[0][3] * o[0][3])

    o[0][0] *= ood
    o[0][1] *= ood
    o[0][2] *= ood
    o[0][3] *= ood
    o[1][0] *= ood
    o[1][1] *= ood
    o[1][2] *= ood
    o[1][3] *= ood
    o[2][0] *= ood
    o[2][1] *= ood
    o[2][2] *= ood
    o[2][3] *= ood
    o[3][0] *= ood
    o[3][1] *= ood
    o[3][2] *= ood
    o[3][3] *= ood

    return o
}

translate_mat :: fn(v: Vec3) -> Mat4 {
    m := identity_mat()
    m[3][0] = v.x
    m[3][1] = v.y
    m[3][2] = v.z
    m[3][3] = 1
    return m
}

//rotate_mat :: fn(v: Vec3, rads: f32) -> Mat4 {
//    c := cos(rads)
//    s := sin(rads)
//
//    a := norm(v)
//    b := 1.0 - c
//    vc := Vec3 { b, b, b }
//    t := a * vc
//
//    rot := identity_mat()
//
//    rot[0][0] = c + t.x * a.x
//    rot[0][1] = 0.0 + t.x * a.y + s * a.z
//    rot[0][2] = 0.0 + t.x * a.z - s * a.y
//    rot[0][3] = 0.0
//
//    rot[1][0] = 0.0 + t.y * a.x - s * a.z
//    rot[1][1] = c + t.y * a.y
//    rot[1][2] = 0.0 + t.y * a.z + s * a.x
//    rot[1][3] = 0.0
//
//    rot[2][0] = 0.0 + t.z * a.x + s * a.y
//    rot[2][1] = 0.0 + t.z * a.y - s * a.x
//    rot[2][2] = c + t.z * a.z
//    rot[2][3] = 0.0
//
//    return rot
//}

scale :: fn(m: Mat4, s: f32) -> Mat4 {
    m[0][0] *= s
    m[1][1] *= s
    m[2][2] *= s
    return m
}

look_at :: fn(eye, centre, up: Vec3) -> Mat4 {
    f := norm(centre - eye)
    s := norm(cross(f, up))
    u := cross(s, f)

    return Mat4 {
        { s.x, u.x, -f.x, 0 },
        { s.y, u.y, -f.y, 0 },
        { s.z, u.z, -f.z, 0 },
        { -dot(s, eye), -dot(u, eye), dot(f, eye), 1}
    }
}

perspective :: fn(fovy, aspect, near, far: f32) -> Mat4 {
    m := Mat4 {
        {0, 0, 0, 0},
        {0, 0, 0, 0},
        {0, 0, 0, 1},
        {0, 0, 0, 0}
    }
    tan_half_fovy := tan(0.5 * fovy)

    m[0][0] = 1.0 / (aspect * tan_half_fovy)
    m[1][1] = 1.0 / (tan_half_fovy)
    m[2][2] = -(far + near) / (far - near)
    m[2][3] = -1.0
    m[3][2] = -2.0 * far * near / (far - near)
    return m
}