// Copyright (c) 2023-2026 Chris (boreddevnl)
// This software is released under the GNU General Public License v3.0. See LICENSE file for details.
// This header needs to maintain in any file it is present in, as per the GPL license terms.

// BOREDOS_APP_DESC: 3/2D Graphing and plotting utility.
// BOREDOS_APP_ICONS: /Library/images/icons/colloid/se.sjoerd.Graphs.png;/Library/images/icons/colloid/app-icon-preview.png
#include "syscall.h"
#include "libui.h"
#include "utf-8.h"
#include "../../wm/libwidget.h"
#include <stdbool.h>
#include "stdlib.h"
#include "math.h"

// External syscalls from libc
extern void sys_parallel_run(void (*fn)(void*), void **args, int count);

// =========
// Constants
// =========
// Adjust the below variables to suit your system specification and preferences.
// --- User Configuration ---
#define ROTATE 1  // Set to 0 to disable auto-rotation in 3D mode.
#define GRID_3D 41 // Grid resolution. Adjust on how much you want your PC to die (lmao)
// --------------------------
#define TOOLBAR_H 30
#define STATUSBAR_H 30
#define GRAPH_Y TOOLBAR_H
#define CLIENT_H (win_h - 20)

static int win_w = 750;
static int win_h = 600;
static int graph_w = 750;
static int graph_h = 520;
static int fb_capacity = 0;
static uint32_t *graph_fb = NULL;
static int32_t *graph_zb = NULL; 
static uint64_t *graph_czb = NULL;  

#define COLOR_BG        0xFF1A1A2E
#define COLOR_GRID      0xFF2A2A4A
#define COLOR_AXIS      0xFF6A6A8A
#define COLOR_CURVE     0xFF00FF88
#define COLOR_TEXT       0xFFE0E0E0
#define COLOR_TOOLBAR_BG 0xFF2D2D2D
#define COLOR_STATUS_BG  0xFF252535
#define COLOR_DARK_BG    0xFF1E1E1E
#define COLOR_DARK_PANEL 0xFF2D2D2D
#define COLOR_DARK_TEXT  0xFFF0F0F0
#define MAX_TOKENS 256
#define MAX_NODES  128

#define MODE_2D 0
#define MODE_3D 1


// =================
// Expression Parser
// =================
enum {
    TOK_NUM, TOK_VAR_X, TOK_VAR_Y, TOK_VAR_Z,
    TOK_PLUS, TOK_MINUS, TOK_STAR, TOK_SLASH, TOK_CARET,
    TOK_LPAREN, TOK_RPAREN, TOK_EQUALS,
    TOK_FN_SIN, TOK_FN_COS, TOK_FN_TAN, TOK_FN_SQRT, TOK_FN_ABS, TOK_FN_LOG,
    TOK_END
};

enum {
    NODE_NUM, NODE_VAR,
    NODE_ADD, NODE_SUB, NODE_MUL, NODE_DIV, NODE_POW,
    NODE_NEG,
    NODE_SIN, NODE_COS, NODE_TAN, NODE_SQRT, NODE_ABS, NODE_LOG
};

typedef struct { int type; double value; } Token;
typedef struct { int type; double value; int var_idx; int left, right; } ASTNode;

enum {
    OP_PUSH_NUM, OP_PUSH_VAR,
    OP_ADD, OP_SUB, OP_MUL, OP_DIV, OP_POW,
    OP_NEG,
    OP_SIN, OP_COS, OP_TAN, OP_SQRT, OP_ABS, OP_LOG
};

typedef struct {
    int op;
    double val;
    int var_idx;
} Instruction;

#define MAX_BC_SIZE 256
static Instruction lhs_bc[MAX_BC_SIZE], rhs_bc[MAX_BC_SIZE];
static int lhs_bc_len = 0, rhs_bc_len = 0;



static bool is_alpha(char c) { return (c >= 'a' && c <= 'z') || (c >= 'A' && c <= 'Z'); }
static bool is_digit(char c) { return c >= '0' && c <= '9'; }
static bool is_value_tok(int t) {
    return t == TOK_NUM || t == TOK_VAR_X || t == TOK_VAR_Y || t == TOK_VAR_Z || t == TOK_RPAREN;
}

static int tokenize(const char *input, Token *tokens) {
    int count = 0;
    const char *p = input;

    while (*p && count < MAX_TOKENS - 1) {
        while (*p == ' ' || *p == '\t') p++;
        if (!*p) break;

        if (count > 0 && is_value_tok(tokens[count-1].type)) {
            if (is_digit(*p) || *p == '.' || is_alpha(*p) || *p == '(') {
                tokens[count].type = TOK_STAR; tokens[count].value = 0; count++;
            }
        }

        if (is_digit(*p) || (*p == '.' && is_digit(p[1]))) {
            double val = 0;
            while (is_digit(*p)) { val = val * 10 + (*p - '0'); p++; }
            if (*p == '.') {
                p++; double frac = 0.1;
                while (is_digit(*p)) { val += (*p - '0') * frac; frac *= 0.1; p++; }
            }
            tokens[count].type = TOK_NUM; tokens[count].value = val; count++;
        } else if (is_alpha(*p)) {
            char buf[32]; int len = 0;
            while (is_alpha(*p) && len < 31) buf[len++] = *p++;
            buf[len] = 0;

            if (len == 1 && (buf[0] == 'x' || buf[0] == 'X')) tokens[count].type = TOK_VAR_X;
            else if (len == 1 && (buf[0] == 'y' || buf[0] == 'Y')) tokens[count].type = TOK_VAR_Y;
            else if (len == 1 && (buf[0] == 'z' || buf[0] == 'Z')) tokens[count].type = TOK_VAR_Z;
            else if (strcmp(buf, "sin") == 0) tokens[count].type = TOK_FN_SIN;
            else if (strcmp(buf, "cos") == 0) tokens[count].type = TOK_FN_COS;
            else if (strcmp(buf, "tan") == 0) tokens[count].type = TOK_FN_TAN;
            else if (strcmp(buf, "sqrt") == 0) tokens[count].type = TOK_FN_SQRT;
            else if (strcmp(buf, "abs") == 0) tokens[count].type = TOK_FN_ABS;
            else if (strcmp(buf, "log") == 0) tokens[count].type = TOK_FN_LOG;
            else if (strcmp(buf, "pi") == 0 || strcmp(buf, "PI") == 0) {
                tokens[count].type = TOK_NUM; tokens[count].value = M_PI;
            } else { tokens[count].type = TOK_NUM; tokens[count].value = 0; }
            count++;
        } else {
            switch (*p) {
                case '+': tokens[count].type = TOK_PLUS; break;
                case '-': tokens[count].type = TOK_MINUS; break;
                case '*': tokens[count].type = TOK_STAR; break;
                case '/': tokens[count].type = TOK_SLASH; break;
                case '^': tokens[count].type = TOK_CARET; break;
                case '(': {
                    if (count > 0 && is_value_tok(tokens[count-1].type)) {
                        tokens[count].type = TOK_STAR; count++;
                    }
                    tokens[count].type = TOK_LPAREN; break;
                }
                case ')': tokens[count].type = TOK_RPAREN; break;
                case '=': tokens[count].type = TOK_EQUALS; break;
                default: p++; continue;
            }
            p++; count++;
        }
    }
    tokens[count].type = TOK_END;
    return count;
}

static int parse_expr(Token *t, int *p, ASTNode *n, int *nc);

static int parse_atom(Token *t, int *p, ASTNode *n, int *nc) {
    if (t[*p].type == TOK_NUM) {
        int i = (*nc)++;
        n[i].type = NODE_NUM; n[i].value = t[*p].value; n[i].left = n[i].right = -1;
        (*p)++; return i;
    }
    if (t[*p].type >= TOK_VAR_X && t[*p].type <= TOK_VAR_Z) {
        int i = (*nc)++;
        n[i].type = NODE_VAR; n[i].var_idx = t[*p].type - TOK_VAR_X;
        n[i].left = n[i].right = -1; (*p)++; return i;
    }
    if (t[*p].type == TOK_LPAREN) {
        (*p)++;
        int inner = parse_expr(t, p, n, nc);
        if (t[*p].type == TOK_RPAREN) (*p)++;
        return inner;
    }
    int i = (*nc)++;
    n[i].type = NODE_NUM; n[i].value = 0; n[i].left = n[i].right = -1;
    return i;
}

static int parse_unary(Token *t, int *p, ASTNode *n, int *nc) {
    if (t[*p].type == TOK_MINUS) {
        (*p)++;
        int op = parse_unary(t, p, n, nc);
        int i = (*nc)++; n[i].type = NODE_NEG; n[i].left = op; n[i].right = -1; return i;
    }
    if (t[*p].type >= TOK_FN_SIN && t[*p].type <= TOK_FN_LOG) {
        int fn_type = t[*p].type; (*p)++;
        int arg;
        if (t[*p].type == TOK_LPAREN) {
            (*p)++;
            arg = parse_expr(t, p, n, nc);
            if (t[*p].type == TOK_RPAREN) (*p)++;
        } else {
            arg = parse_unary(t, p, n, nc);
        }
        int node_type = NODE_SIN;
        if (fn_type == TOK_FN_COS) node_type = NODE_COS;
        else if (fn_type == TOK_FN_TAN) node_type = NODE_TAN;
        else if (fn_type == TOK_FN_SQRT) node_type = NODE_SQRT;
        else if (fn_type == TOK_FN_ABS) node_type = NODE_ABS;
        else if (fn_type == TOK_FN_LOG) node_type = NODE_LOG;
        int i = (*nc)++; n[i].type = node_type; n[i].left = arg; n[i].right = -1; return i;
    }
    return parse_atom(t, p, n, nc);
}

static int parse_power(Token *t, int *p, ASTNode *n, int *nc) {
    int left = parse_unary(t, p, n, nc);
    if (t[*p].type == TOK_CARET) {
        (*p)++;
        int right = parse_power(t, p, n, nc); 
        int i = (*nc)++; n[i].type = NODE_POW; n[i].left = left; n[i].right = right; return i;
    }
    return left;
}

static int parse_term(Token *t, int *p, ASTNode *n, int *nc) {
    int left = parse_power(t, p, n, nc);
    while (t[*p].type == TOK_STAR || t[*p].type == TOK_SLASH) {
        int op = (t[*p].type == TOK_STAR) ? NODE_MUL : NODE_DIV; (*p)++;
        int right = parse_power(t, p, n, nc);
        int i = (*nc)++; n[i].type = op; n[i].left = left; n[i].right = right; left = i;
    }
    return left;
}

static int parse_expr(Token *t, int *p, ASTNode *n, int *nc) {
    int left = parse_term(t, p, n, nc);
    while (t[*p].type == TOK_PLUS || t[*p].type == TOK_MINUS) {
        int op = (t[*p].type == TOK_PLUS) ? NODE_ADD : NODE_SUB; (*p)++;
        int right = parse_term(t, p, n, nc);
        int i = (*nc)++; n[i].type = op; n[i].left = left; n[i].right = right; left = i;
    }
    return left;
}

static double eval_ast(ASTNode *n, int idx, double x, double y, double z) {
    if (idx < 0) return 0;
    ASTNode *nd = &n[idx];
    switch (nd->type) {
        case NODE_NUM: return nd->value;
        case NODE_VAR: return nd->var_idx == 0 ? x : nd->var_idx == 1 ? y : z;
        case NODE_ADD: return eval_ast(n, nd->left, x, y, z) + eval_ast(n, nd->right, x, y, z);
        case NODE_SUB: return eval_ast(n, nd->left, x, y, z) - eval_ast(n, nd->right, x, y, z);
        case NODE_MUL: return eval_ast(n, nd->left, x, y, z) * eval_ast(n, nd->right, x, y, z);
        case NODE_DIV: {
            double d = eval_ast(n, nd->right, x, y, z);
            return fabs(d) < 1e-15 ? 1e15 : eval_ast(n, nd->left, x, y, z) / d;
        }
        case NODE_POW: {
            double b = eval_ast(n, nd->left, x, y, z);
            double e = eval_ast(n, nd->right, x, y, z);
            if (e == 2.0) return b * b;
            if (e == 3.0) return b * b * b;
            return pow(b, e);
        }
        case NODE_NEG: return -eval_ast(n, nd->left, x, y, z);
        case NODE_SIN: return sin(eval_ast(n, nd->left, x, y, z));
        case NODE_COS: return cos(eval_ast(n, nd->left, x, y, z));
        case NODE_TAN: return tan(eval_ast(n, nd->left, x, y, z));
        case NODE_SQRT: return sqrt(eval_ast(n, nd->left, x, y, z));
        case NODE_ABS: return fabs(eval_ast(n, nd->left, x, y, z));
        case NODE_LOG: return log(eval_ast(n, nd->left, x, y, z));
    }
    return 0;
}

// =================
// Bytecode Compiler
// =================
static int compile_ast(ASTNode *nodes, int idx, Instruction *bc, int *len) {
    if (idx < 0 || *len >= MAX_BC_SIZE) return 0;
    ASTNode *n = &nodes[idx];
    
    compile_ast(nodes, n->left, bc, len);
    compile_ast(nodes, n->right, bc, len);
    
    Instruction *inst = &bc[(*len)++];
    switch (n->type) {
        case NODE_NUM: inst->op = OP_PUSH_NUM; inst->val = n->value; break;
        case NODE_VAR: inst->op = OP_PUSH_VAR; inst->var_idx = n->var_idx; break;
        case NODE_ADD: inst->op = OP_ADD; break;
        case NODE_SUB: inst->op = OP_SUB; break;
        case NODE_MUL: inst->op = OP_MUL; break;
        case NODE_DIV: inst->op = OP_DIV; break;
        case NODE_POW: inst->op = OP_POW; break;
        case NODE_NEG: inst->op = OP_NEG; break;
        case NODE_SIN: inst->op = OP_SIN; break;
        case NODE_COS: inst->op = OP_COS; break;
        case NODE_TAN: inst->op = OP_TAN; break;
        case NODE_SQRT: inst->op = OP_SQRT; break;
        case NODE_ABS: inst->op = OP_ABS; break;
        case NODE_LOG: inst->op = OP_LOG; break;
    }
    return 1;
}

static double run_bc(Instruction *bc, int len, double x, double y, double z) {
    if (len == 0) return 0;
    double stack[32];
    int sp = 0;
    
    for (int i = 0; i < len; i++) {
        Instruction *inst = &bc[i];
        switch (inst->op) {
            case OP_PUSH_NUM: stack[sp++] = inst->val; break;
            case OP_PUSH_VAR: stack[sp++] = (inst->var_idx == 0 ? x : inst->var_idx == 1 ? y : z); break;
            case OP_ADD: { double b = stack[--sp]; double a = stack[--sp]; stack[sp++] = a + b; break; }
            case OP_SUB: { double b = stack[--sp]; double a = stack[--sp]; stack[sp++] = a - b; break; }
            case OP_MUL: { double b = stack[--sp]; double a = stack[--sp]; stack[sp++] = a * b; break; }
            case OP_DIV: { 
                double b = stack[--sp]; double a = stack[--sp]; 
                stack[sp++] = (fabs(b) < 1e-15) ? 1e15 : a / b; break; 
            }
            case OP_POW: {
                double b = stack[--sp]; double a = stack[--sp];
                if (b == 2.0) stack[sp++] = a * a;
                else if (b == 3.0) stack[sp++] = a * a * a;
                else stack[sp++] = pow(a, b);
                break;
            }
            case OP_NEG: stack[sp-1] = -stack[sp-1]; break;
            case OP_SIN: stack[sp-1] = sin(stack[sp-1]); break;
            case OP_COS: stack[sp-1] = cos(stack[sp-1]); break;
            case OP_TAN: stack[sp-1] = tan(stack[sp-1]); break;
            case OP_SQRT: stack[sp-1] = sqrt(stack[sp-1]); break;
            case OP_ABS:  stack[sp-1] = fabs(stack[sp-1]); break;
            case OP_LOG:  stack[sp-1] = log(stack[sp-1]); break;
        }
    }
    return sp > 0 ? stack[sp-1] : 0;
}

// Check which variables an AST subtree uses
static void ast_find_vars(ASTNode *n, int idx, bool *has_x, bool *has_y, bool *has_z) {
    if (idx < 0) return;
    if (n[idx].type == NODE_VAR) {
        if (n[idx].var_idx == 0) *has_x = true;
        if (n[idx].var_idx == 1) *has_y = true;
        if (n[idx].var_idx == 2) *has_z = true;
    }
    ast_find_vars(n, n[idx].left, has_x, has_y, has_z);
    ast_find_vars(n, n[idx].right, has_x, has_y, has_z);
}

// =================
// Application State
// =================
static ui_window_t win_graph;

static char eq_buffer[256];
static int eq_len = 0;

// Parsed equation
static ASTNode lhs_nodes[MAX_NODES], rhs_nodes[MAX_NODES];
static int lhs_root = -1, rhs_root = -1;
static int lhs_nc = 0, rhs_nc = 0;
static bool eq_valid = false;
static int graph_mode = MODE_2D;

// 2D view
static double view_x_min = -10, view_x_max = 10;
static double view_y_min = -6.4, view_y_max = 6.4;

// 3D view
static double rot_x = 0.5, rot_y = 0.6; 
static double range_3d = 5.0;
static double zoom_3d = 1.0; 
static bool filled_mode = false;
static bool is_explicit_3d = false;
static bool is_explicit_2d = false;

static bool right_dragging = false;
static int drag_last_x = 0, drag_last_y = 0;

#define MAX_Z_PER_POINT 4
typedef struct {
    double z[MAX_Z_PER_POINT];
    float nx[MAX_Z_PER_POINT], ny[MAX_Z_PER_POINT], nz[MAX_Z_PER_POINT];
    int sx[MAX_Z_PER_POINT], sy[MAX_Z_PER_POINT], dz[MAX_Z_PER_POINT];
    int count;  
} surf_point_t;

static double surf_x[GRID_3D][GRID_3D], surf_y_3d[GRID_3D][GRID_3D];
static surf_point_t surf[GRID_3D][GRID_3D];

static bool surface_needs_eval = true;
static double rot_cx, rot_sx, rot_cy, rot_sy;

static widget_textbox_t tb_equation;
static widget_button_t btn_plot;
static widget_button_t btn_mode;
static widget_button_t btn_range_minus;
static widget_button_t btn_range_plus;

static const char *preset_labels[] = {
    "y = sin(x)", "y = x^2", "y = cos(x)*x",
    "z = sin(x)*cos(y)", "z = x^2 - y^2", "x^2+y^2+z^2=25", "x^2+y^2=16"
};
#define NUM_PRESETS 7

static bool presets_open = false;

// Widget context
static void gfx_draw_rect(void *ud, int x, int y, int w, int h, uint32_t c) {
    ui_draw_rect((ui_window_t)ud, x, y, w, h, c);
}
static void gfx_draw_rr(void *ud, int x, int y, int w, int h, int r, uint32_t c) {
    ui_draw_rounded_rect_filled((ui_window_t)ud, x, y, w, h, r, c);
}
static void gfx_draw_str(void *ud, int x, int y, const char *s, uint32_t c) {
    ui_draw_string((ui_window_t)ud, x, y, s, c);
}
static widget_context_t wctx = { 0, gfx_draw_rect, gfx_draw_rr, gfx_draw_str, NULL, NULL, false };

// ================
// Graphics helpers
// ================
static void gfb_clear(uint32_t c) {
    int total = graph_w * graph_h;
    uint64_t clear_val = ((uint64_t)1000000LL << 32) | c;
    for (int i = 0; i < total; i++) {
        graph_fb[i] = c;
        graph_zb[i] = 1000000;
        if (graph_czb) graph_czb[i] = clear_val;
    }
}

static void gfb_pixel(int x, int y, uint32_t c) {
    if (x >= 0 && x < graph_w && y >= 0 && y < graph_h)
        graph_fb[y * graph_w + x] = c;
}

static void gfb_pixel_z(int x, int y, int z, uint32_t c) {
    if (x < 0 || x >= graph_w || y < 0 || y >= graph_h) return;
    int idx = y * graph_w + x;
    
    if (graph_czb) {
        uint64_t new_val = ((uint64_t)z << 32) | c;
        uint64_t old_val;
        while (z < (int32_t)((old_val = __atomic_load_n(&graph_czb[idx], __ATOMIC_RELAXED)) >> 32)) {
            if (__atomic_compare_exchange_n(&graph_czb[idx], &old_val, new_val, false, __ATOMIC_SEQ_CST, __ATOMIC_RELAXED))
                break;
        }
    } else {
        int32_t old_z;
        while (z < (old_z = __atomic_load_n(&graph_zb[idx], __ATOMIC_RELAXED))) {
            if (__atomic_compare_exchange_n(&graph_zb[idx], &old_z, z, false, __ATOMIC_SEQ_CST, __ATOMIC_RELAXED)) {
                graph_fb[idx] = c;
                break;
            }
        }
    }
}



static void gfb_line(int x0, int y0, int x1, int y1, uint32_t c) {
    int dx = fabs((double)(x1 - x0)), dy = fabs((double)(y1 - y0));
    int sx = x0 < x1 ? 1 : -1, sy = y0 < y1 ? 1 : -1;
    int err = dx - dy;
    for (int i = 0; i < 2000; i++) {
        gfb_pixel(x0, y0, c);
        if (x0 == x1 && y0 == y1) break;
        int e2 = 2 * err;
        if (e2 > -dy) { err -= dy; x0 += sx; }
        if (e2 < dx) { err += dx; y0 += sy; }
    }
}

static void gfb_line_z(int x0, int y0, int z0, int x1, int y1, int z1, uint32_t c) {
    int dx = fabs((double)(x1 - x0)), dy = fabs((double)(y1 - y0));
    int sx = x0 < x1 ? 1 : -1, sy = y0 < y1 ? 1 : -1;
    int err = dx - dy;
    int steps = (dx > dy ? dx : dy);
    if (steps == 0) { gfb_pixel_z(x0, y0, z0, c); return; }
    
    for (int i = 0; i <= steps && i < 2000; i++) {
        int cz = z0 + (int)((long)(z1 - z0) * i / steps);
        gfb_pixel_z(x0, y0, cz, c);
        if (x0 == x1 && y0 == y1) break;
        int e2 = 2 * err;
        if (e2 > -dy) { err -= dy; x0 += sx; }
        if (e2 < dx) { err += dx; y0 += sy; }
    }
}

static void gfb_triangle_z(int x0, int y0, int z0, int x1, int y1, int z1, int x2, int y2, int z2, uint32_t c) {
    if (y1 < y0) { int t; t=x0; x0=x1; x1=t; t=y0; y0=y1; y1=t; t=z0; z0=z1; z1=t; }
    if (y2 < y0) { int t; t=x0; x0=x2; x2=t; t=y0; y0=y2; y2=t; t=z0; z0=z2; z2=t; }
    if (y2 < y1) { int t; t=x1; x1=x2; x2=t; t=y1; y1=y2; y2=t; t=z1; z1=z2; z2=t; }

    if (y0 == y2) return;

    for (int y = y0; y <= y2; y++) {
        if (y < 0 || y >= graph_h) continue;
        bool second_half = y > y1 || y1 == y0;
        int h1 = second_half ? (y2 - y1) : (y1 - y0);
        int h2 = y2 - y0;
        if (h1 == 0) h1 = 1;
        if (h2 == 0) h2 = 1;

        float alpha = (float)(y - y0) / h2;
        float beta  = (float)(y - (second_half ? y1 : y0)) / h1;

        int ax = x0 + (int)((x2 - x0) * alpha);
        int bx = second_half ? x1 + (int)((x2 - x1) * beta) : x0 + (int)((x1 - x0) * beta);
        int az = z0 + (int)((z2 - z0) * alpha);
        int bz = second_half ? z1 + (int)((z2 - z1) * beta) : z0 + (int)((z1 - z0) * beta);

        if (ax > bx) { int t; t=ax; ax=bx; bx=t; t=az; az=bz; bz=t; }
        int span = bx - ax;
        int x_start = ax < 0 ? 0 : ax;
        int x_end   = bx >= graph_w ? graph_w - 1 : bx;
        for (int x = x_start; x <= x_end; x++) {
            float phi = (span == 0) ? 0.5f : (float)(x - ax) / span;
            int cz = az + (int)((bz - az) * phi);
            gfb_pixel_z(x, y, cz, c);
        }
    }
}

static uint32_t color_by_height(double z, double zmin, double zmax) {
    if (zmax <= zmin) return COLOR_CURVE;
    double t = (z - zmin) / (zmax - zmin);
    if (t < 0) t = 0;
    if (t > 1) t = 1;
    int r, g, b;
    if (t < 0.25) { double s = t/0.25; r=0; g=(int)(s*255); b=255; }
    else if (t < 0.5) { double s=(t-0.25)/0.25; r=0; g=255; b=(int)((1-s)*255); }
    else if (t < 0.75) { double s=(t-0.5)/0.25; r=(int)(s*255); g=255; b=0; }
    else { double s=(t-0.75)/0.25; r=255; g=(int)((1-s)*255); b=0; }
    return 0xFF000000 | (r<<16) | (g<<8) | b;
}

static uint32_t apply_shading(uint32_t color, double nx, double ny, double nz) {
    double len = sqrt(nx*nx + ny*ny + nz*nz);
    if (len > 1e-9) { nx /= len; ny /= len; nz /= len; }
    else { nx = 0; ny = 1; nz = 0; }

    double lx = 0.577, ly = 0.707, lz = 0.408;
    double dot = nx * lx + ny * ly + nz * lz;
    if (dot < 0) dot = -dot * 0.2; 
    
    double intensity = 0.3 + 0.7 * dot;
    if (intensity > 1.0) intensity = 1.0;
    
    int r = (color >> 16) & 0xFF;
    int g = (color >> 8) & 0xFF;
    int b = color & 0xFF;
    
    r = (int)(r * intensity);
    g = (int)(g * intensity);
    b = (int)(b * intensity);
    if (r > 255) r = 255; if (g > 255) g = 255; if (b > 255) b = 255;
    
    return 0xFF000000 | (r << 16) | (g << 8) | b;
}

// ========================
// 2D coordinate transforms
// ========================
static int screen_x_2d(double wx) {
    return (int)((wx - view_x_min) / (view_x_max - view_x_min) * graph_w);
}
static int screen_y_2d(double wy) {
    return (int)((view_y_max - wy) / (view_y_max - view_y_min) * graph_h);
}
static double world_x_2d(int px) {
    return view_x_min + (px / (double)graph_w) * (view_x_max - view_x_min);
}

// =============
// 3D projection
// =============
static void project_3d(double px, double py, double pz, int *sx, int *sy, int *sz) {
    double nx = px * rot_cy + pz * rot_sy, nz = -px * rot_sy + pz * rot_cy;
    px = nx; pz = nz;
    double ny = py * rot_cx - pz * rot_sx;
    nz = py * rot_sx + pz * rot_cx;
    py = ny; pz = nz;
    
    double view_dim = graph_h < graph_w ? graph_h : graph_w;
    double base_scale = view_dim * 0.35 / range_3d * zoom_3d;
    
    double d = range_3d * 5;
    double zd = pz + d;
    if (zd < d * 0.1) zd = d * 0.1;
    double persp = d / zd;
    *sx = (int)(px * base_scale * persp) + graph_w / 2;
    *sy = (int)(-py * base_scale * persp) + graph_h / 2;
    *sz = (int)(pz * 10000); // Fixed-point depth
}
// istg math is scary
// ====================================================
// Evaluate the implicit function: f(x,y,z) = LHS - RHS
// ====================================================
static double eval_implicit(double x, double y, double z) {
    return run_bc(lhs_bc, lhs_bc_len, x, y, z) - run_bc(rhs_bc, rhs_bc_len, x, y, z);
}

static double eval_rhs_only(double x, double y, double z) {
    return run_bc(rhs_bc, rhs_bc_len, x, y, z);
}

// ===========================
// Parse and classify equation
// ===========================
static void parse_equation(void) {
    surface_needs_eval = true;
    eq_valid = false;
    lhs_nc = 0; rhs_nc = 0;
    lhs_root = -1; rhs_root = -1;

    // Find '='
    int eq_pos = -1;
    for (int i = 0; eq_buffer[i]; i++) {
        if (eq_buffer[i] == '=') { eq_pos = i; break; }
    }

    if (eq_pos >= 0) {
        // Split into LHS and RHS
        char lhs_str[256], rhs_str[256];
        memcpy(lhs_str, eq_buffer, eq_pos); lhs_str[eq_pos] = 0;
        strcpy(rhs_str, eq_buffer + eq_pos + 1);

        Token lt[MAX_TOKENS], rt[MAX_TOKENS];
        tokenize(lhs_str, lt); tokenize(rhs_str, rt);
        int lp = 0, rp = 0;
        lhs_root = parse_expr(lt, &lp, lhs_nodes, &lhs_nc);
        rhs_root = parse_expr(rt, &rp, rhs_nodes, &rhs_nc);
    } else {
        Token tt[MAX_TOKENS];
        tokenize(eq_buffer, tt);
        int tp = 0;
        rhs_root = parse_expr(tt, &tp, rhs_nodes, &rhs_nc);
        // LHS = y (for display as y=f(x) if no y/z in expr)
        bool hx=false, hy=false, hz=false;
        ast_find_vars(rhs_nodes, rhs_root, &hx, &hy, &hz);
        if (!hy && !hz) {
            lhs_nc = 1;
            lhs_nodes[0].type = NODE_VAR; lhs_nodes[0].var_idx = 1;
            lhs_nodes[0].left = lhs_nodes[0].right = -1;
            lhs_root = 0;
        } else {
            lhs_root = rhs_root;
            memcpy(lhs_nodes, rhs_nodes, sizeof(ASTNode) * rhs_nc);
            lhs_nc = rhs_nc;
            rhs_nc = 1;
            rhs_nodes[0].type = NODE_NUM; rhs_nodes[0].value = 0;
            rhs_nodes[0].left = rhs_nodes[0].right = -1;
            rhs_root = 0;
        }
    }

    // Determine mode
    bool lhx=false, lhy=false, lhz=false, rhx=false, rhy=false, rhz=false;
    ast_find_vars(lhs_nodes, lhs_root, &lhx, &lhy, &lhz);
    ast_find_vars(rhs_nodes, rhs_root, &rhx, &rhy, &rhz);

    bool has_z = lhz || rhz;
    is_explicit_2d = false;
    is_explicit_3d = false;

    if (has_z) {
        graph_mode = MODE_3D;
        if (lhs_nc >= 1 && lhs_nodes[lhs_root].type == NODE_VAR &&
            lhs_nodes[lhs_root].var_idx == 2 && !rhz) {
            is_explicit_3d = true;
        }
    } else {
        graph_mode = MODE_2D;
        if (lhs_nc >= 1 && lhs_nodes[lhs_root].type == NODE_VAR &&
            lhs_nodes[lhs_root].var_idx == 1 && !rhy) {
            is_explicit_2d = true;
        }
    }

    eq_valid = true;

    // Compile to bytecode
    lhs_bc_len = 0; rhs_bc_len = 0;
    if (lhs_root >= 0) compile_ast(lhs_nodes, lhs_root, lhs_bc, &lhs_bc_len);
    if (rhs_root >= 0) compile_ast(rhs_nodes, rhs_root, rhs_bc, &rhs_bc_len);
}

// =========
// Rendering
// =========
static double get_nice_step(double range, int target_divisions) {
    if (range <= 0) return 1.0;
    double approx = range / target_divisions;
    
    // Find magnitude (10^n)
    double mag = 1.0;
    if (approx >= 1.0) {
        while (approx >= 10.0) { approx /= 10.0; mag *= 10.0; }
    } else {
        while (approx < 1.0) { approx *= 10.0; mag /= 10.0; }
    }
    
    // Pick nice residual
    double res;
    if (approx < 1.5) res = 1.0;
    else if (approx < 3.0) res = 2.0;
    else if (approx < 7.0) res = 5.0;
    else res = 10.0;
    
    return res * mag;
}

static void apply_aspect_ratio(void) {
    if (graph_w <= 0 || graph_h <= 0) return;
    double cy = (view_y_min + view_y_max) / 2.0;
    double x_range = view_x_max - view_x_min;
    double target_y_range = x_range * (double)graph_h / (double)graph_w;
    view_y_min = cy - target_y_range / 2.0;
    view_y_max = cy + target_y_range / 2.0;
}

static void autofit_2d_view(void) {
    if (!is_explicit_2d) {
        apply_aspect_ratio();
        return;
    }
    double y_min_data = 1e30, y_max_data = -1e30;
    bool found = false;
    for (int px = 0; px < graph_w; px += 2) {
        double wx = world_x_2d(px);
        double wy = eval_rhs_only(wx, 0, 0);
        if (wy == wy && fabs(wy) < 1e10) {
            if (wy < y_min_data) y_min_data = wy;
            if (wy > y_max_data) y_max_data = wy;
            found = true;
        }
    }
    if (found) {
        if (y_min_data * y_max_data <= 0 || fabs(y_min_data) < (y_max_data - y_min_data) * 0.5) {
            double max_abs = fabs(y_min_data);
            if (fabs(y_max_data) > max_abs) max_abs = fabs(y_max_data);
            double pad = max_abs * 0.15;
            if (pad < 0.5) pad = 0.5;
            view_y_min = -(max_abs + pad);
            view_y_max = (max_abs + pad);
        } else {
            double pad = (y_max_data - y_min_data) * 0.15;
            if (pad < 0.5) pad = 0.5;
            view_y_min = y_min_data - pad;
            view_y_max = y_max_data + pad;
        }

        double x_range = view_x_max - view_x_min;
        double target_y_range = x_range * (double)graph_h / (double)graph_w;
        double current_y_range = view_y_max - view_y_min;

        if (current_y_range < target_y_range) {
            double cy = (view_y_min + view_y_max) / 2.0;
            if (fabs(cy) < current_y_range * 0.1) cy = 0;
            view_y_min = cy - target_y_range / 2.0;
            view_y_max = cy + target_y_range / 2.0;
        } else {
            double target_x_range = current_y_range * (double)graph_w / (double)graph_h;
            double cx = (view_x_min + view_x_max) / 2.0;
            if (fabs(cx) < x_range * 0.1) cx = 0;
            view_x_min = cx - target_x_range / 2.0;
            view_x_max = cx + target_x_range / 2.0;
        }
    } else {
        apply_aspect_ratio();
    }
}

static void render_2d_grid(void) {
    // Grid intervals
    double x_step = get_nice_step(view_x_max - view_x_min, 8);
    double y_step = get_nice_step(view_y_max - view_y_min, 6);

    // X grid lines
    double x_start = (int)(view_x_min / x_step) * x_step;
    int safety = 0;
    for (double wx = x_start - x_step; wx <= view_x_max + x_step && safety < 200; wx += x_step, safety++) {
        int sx = screen_x_2d(wx);
        if (sx >= 0 && sx < graph_w) {
            for (int y = 0; y < graph_h; y++) gfb_pixel(sx, y, COLOR_GRID);
        }
    }

    // Y grid lines
    double y_start = (int)(view_y_min / y_step) * y_step;
    safety = 0;
    for (double wy = y_start - y_step; wy <= view_y_max + y_step && safety < 200; wy += y_step, safety++) {
        int sy = screen_y_2d(wy);
        if (sy >= 0 && sy < graph_h) {
            for (int x = 0; x < graph_w; x++) gfb_pixel(x, sy, COLOR_GRID);
        }
    }

    // Axes
    int ax = screen_x_2d(0), ay = screen_y_2d(0);
    if (ax >= 0 && ax < graph_w)
        for (int y = 0; y < graph_h; y++) gfb_pixel(ax, y, COLOR_AXIS);
    if (ay >= 0 && ay < graph_h)
        for (int x = 0; x < graph_w; x++) gfb_pixel(x, ay, COLOR_AXIS);
}

static void render_2d_explicit(void) {
    int prev_sx = -1, prev_sy = -1;
    bool prev_valid = false;

    for (int px = 0; px < graph_w; px++) {
        double wx = world_x_2d(px);
        double wy = eval_rhs_only(wx, 0, 0);
        if (wy != wy || fabs(wy) > 1e10) { prev_valid = false; continue; }
        int sy = screen_y_2d(wy);
        if (prev_valid && fabs((double)(sy - prev_sy)) < graph_h) {
            gfb_line(prev_sx, prev_sy, px, sy, COLOR_CURVE);
        }
        prev_sx = px; prev_sy = sy; prev_valid = true;
    }
}

static void render_2d_implicit(void) {
    // Marching squares for f(x,y) = 0
    int grid_x = 200, grid_y = 130;
    double dx = (view_x_max - view_x_min) / grid_x;
    double dy = (view_y_max - view_y_min) / grid_y;

    for (int gy = 0; gy < grid_y; gy++) {
        for (int gx = 0; gx < grid_x; gx++) {
            double x0 = view_x_min + gx * dx;
            double y0 = view_y_max - gy * dy;
            double x1 = x0 + dx, y1 = y0 - dy;

            double f00 = eval_implicit(x0, y0, 0);
            double f10 = eval_implicit(x1, y0, 0);
            double f01 = eval_implicit(x0, y1, 0);
            double f11 = eval_implicit(x1, y1, 0);

            // Check edges for sign changes
            int sx0 = screen_x_2d(x0), sx1 = screen_x_2d(x1);
            int sy0 = screen_y_2d(y0), sy1 = screen_y_2d(y1);

            if ((f00 > 0) != (f10 > 0)) {
                double t = f00 / (f00 - f10);
                int mx = sx0 + (int)(t * (sx1 - sx0));
                gfb_pixel(mx, sy0, COLOR_CURVE);
                gfb_pixel(mx+1, sy0, COLOR_CURVE);
            }
            if ((f00 > 0) != (f01 > 0)) {
                double t = f00 / (f00 - f01);
                int my = sy0 + (int)(t * (sy1 - sy0));
                gfb_pixel(sx0, my, COLOR_CURVE);
                gfb_pixel(sx0, my+1, COLOR_CURVE);
            }
            if ((f10 > 0) != (f11 > 0)) {
                double t = f10 / (f10 - f11);
                int my = sy0 + (int)(t * (sy1 - sy0));
                gfb_pixel(sx1, my, COLOR_CURVE);
            }
            if ((f01 > 0) != (f11 > 0)) {
                double t = f01 / (f01 - f11);
                int mx = sx0 + (int)(t * (sx1 - sx0));
                gfb_pixel(mx, sy1, COLOR_CURVE);
            }
        }
    }
}

static void render_3d_axes(void) {
    int ox, oy, oz;
    project_3d(0, 0, 0, &ox, &oy, &oz);

    int ax, ay, az;
    project_3d(range_3d, 0, 0, &ax, &ay, &az); gfb_line_z(ox, oy, oz, ax, ay, az, 0xFFFF4444);
    project_3d(0, range_3d, 0, &ax, &ay, &az); gfb_line_z(ox, oy, oz, ax, ay, az, 0xFF44FF44);
    project_3d(0, 0, range_3d, &ax, &ay, &az); gfb_line_z(ox, oy, oz, ax, ay, az, 0xFF4444FF);
}

// =======================
// Parallel Evaluation Job
// =======================

typedef struct {
    float c[MAX_Z_PER_POINT];  
    float nx[MAX_Z_PER_POINT], ny[MAX_Z_PER_POINT], nz[MAX_Z_PER_POINT]; 
    int   count;
} eval_cache_entry_t;

static eval_cache_entry_t eval_cache[3][GRID_3D][GRID_3D];

typedef struct {
    int start_j, end_j;
    double range;
    double step;
    double z_scale;
    int march_axis;
} eval_job_t;

static void eval_3d_explicit_job(void *arg) {
    eval_job_t *job = (eval_job_t *)arg;
    for (int j = job->start_j; j < job->end_j; j++) {
        for (int i = 0; i < GRID_3D; i++) {
            double wx = -job->range + i * job->step;
            double wy = -job->range + j * job->step;
            double wz = eval_rhs_only(wx, wy, 0);
            surf_x[j][i] = wx;
            surf_y_3d[j][i] = wy;
            surf[j][i].count = 0;
            
            if (fabs(wz) > 1e10 || wz != wz) {
                continue;
            }
            
            surf[j][i].z[0] = wz;
            
            double eps = 0.001;
            double dfx = 0.5 * (eval_rhs_only(wx+eps, wy, 0) - eval_rhs_only(wx-eps, wy, 0)) / eps;
            double dfy = 0.5 * (eval_rhs_only(wx, wy+eps, 0) - eval_rhs_only(wx, wy-eps, 0)) / eps;
            
            double nx = -dfx;
            double ny = -dfy;
            double nz = 1.0;
            double len = sqrt(nx*nx + ny*ny + nz*nz);
            if (len > 1e-9) { nx /= len; ny /= len; nz /= len; }
            else { nx = 0; ny = 0; nz = 1; }
            
            surf[j][i].nx[0] = (float)nx;
            surf[j][i].ny[0] = (float)ny;
            surf[j][i].nz[0] = (float)nz;
            surf[j][i].count = 1;
        }
    }
}

static void eval_3d_implicit_job(void *arg) {
    eval_job_t *job = (eval_job_t *)arg;
    const int axis        = job->march_axis;
    const int march_steps = 170;  
    const double mstep    = job->range * 2.0 / march_steps;

    for (int j = job->start_j; j < job->end_j; j++) {
        for (int i = 0; i < GRID_3D; i++) {
            double a = -job->range + i * job->step; 
            double b = -job->range + j * job->step; 

            eval_cache_entry_t *ce = &eval_cache[axis][j][i];
            ce->count = 0;

            double prev_f;
            switch (axis) {
                case 1:  prev_f = eval_implicit(-job->range, a, b); break;
                case 2:  prev_f = eval_implicit(a, -job->range, b); break;
                default: prev_f = eval_implicit(a, b, -job->range); break;
            }

            double found[MAX_Z_PER_POINT]; int nf = 0;
            for (int k = 1; k <= march_steps && nf < MAX_Z_PER_POINT; k++) {
                double c = -job->range + k * mstep;
                double cur_f;
                switch (axis) {
                    case 1:  cur_f = eval_implicit(c, a, b); break;
                    case 2:  cur_f = eval_implicit(a, c, b); break;
                    default: cur_f = eval_implicit(a, b, c); break;
                }
                if ((prev_f > 0) != (cur_f > 0) && fabs(prev_f) < 1e10 && fabs(cur_f) < 1e10) {
                    double ca = c - mstep, cb = c, fa = prev_f, fb = cur_f; (void)fb;
                    for (int bi = 0; bi < 15; bi++) {
                        double cm = (ca+cb)*0.5, fm;
                        switch (axis) {
                            case 1:  fm = eval_implicit(cm, a, b); break;
                            case 2:  fm = eval_implicit(a, cm, b); break;
                            default: fm = eval_implicit(a, b, cm); break;
                        }
                        if ((fa>0)!=(fm>0)) { cb=cm; fb=fm; } else { ca=cm; fa=fm; }
                    }
                    found[nf++] = (ca+cb)*0.5;
                }
                prev_f = cur_f;
            }
            for (int r = 0; r < nf-1; r++)
                for (int s = r+1; s < nf; s++)
                    if (found[r] > found[s]) { double t=found[r]; found[r]=found[s]; found[s]=t; }

            double eps = 0.001;
            for (int r = 0; r < nf; r++) {
                double c = found[r];
                double nx, ny, nz;
                switch (axis) {
                    case 1:
                        nx = eval_implicit(c+eps,a,b)-eval_implicit(c-eps,a,b);
                        ny = eval_implicit(c,a+eps,b)-eval_implicit(c,a-eps,b);
                        nz = eval_implicit(c,a,b+eps)-eval_implicit(c,a,b-eps);
                        break;
                    case 2:
                        nx = eval_implicit(a+eps,c,b)-eval_implicit(a-eps,c,b);
                        ny = eval_implicit(a,c+eps,b)-eval_implicit(a,c-eps,b);
                        nz = eval_implicit(a,c,b+eps)-eval_implicit(a,c,b-eps);
                        break;
                    default:
                        nx = eval_implicit(a+eps,b,c)-eval_implicit(a-eps,b,c);
                        ny = eval_implicit(a,b+eps,c)-eval_implicit(a,b-eps,c);
                        nz = eval_implicit(a,b,c+eps)-eval_implicit(a,b,c-eps);
                        break;
                }
                double d = sqrt(nx*nx+ny*ny+nz*nz);
                if (d > 1e-12) { nx/=d; ny/=d; nz/=d; } else { nx=0; ny=1; nz=0; }
                ce->c[r]  = (float)c;
                ce->nx[r] = (float)nx;
                ce->ny[r] = (float)ny;
                ce->nz[r] = (float)nz;
            }
            ce->count = nf;
        }
    }
}

static void eval_3d_project_job(void *arg) {
    eval_job_t *job = (eval_job_t *)arg;
    const int axis = job->march_axis;
    for (int j = job->start_j; j < job->end_j; j++) {
        for (int i = 0; i < GRID_3D; i++) {
            eval_cache_entry_t *ce = &eval_cache[axis][j][i];
            double a = -job->range + i * job->step;
            double b = -job->range + j * job->step;

            surf_x[j][i]    = a;
            surf_y_3d[j][i] = b;
            surf[j][i].count = ce->count;
            for (int s = 0; s < ce->count; s++) {
                double c = ce->c[s];
                surf[j][i].z[s]  = (float)c;
                surf[j][i].nx[s] = ce->nx[s];
                surf[j][i].ny[s] = ce->ny[s];
                surf[j][i].nz[s] = ce->nz[s];
                switch (axis) {
                    case 1: project_3d(c, b, a, &surf[j][i].sx[s], &surf[j][i].sy[s], &surf[j][i].dz[s]); break;
                    case 2: project_3d(a, b, c, &surf[j][i].sx[s], &surf[j][i].sy[s], &surf[j][i].dz[s]); break;
                    default:project_3d(a, c, b, &surf[j][i].sx[s], &surf[j][i].sy[s], &surf[j][i].dz[s]); break;
                }
            }
        }
    }
}

static void eval_3d_explicit_project_job(void *arg) {
    eval_job_t *job = (eval_job_t *)arg;
    for (int j = job->start_j; j < job->end_j; j++) {
        for (int i = 0; i < GRID_3D; i++) {
            double a = -job->range + i * job->step;
            double b = -job->range + j * job->step;
            
            for (int s = 0; s < surf[j][i].count; s++) {
                double c = surf[j][i].z[s];
                project_3d(a, c * job->z_scale, b, &surf[j][i].sx[s], &surf[j][i].sy[s], &surf[j][i].dz[s]);
            }
        }
    }
}

typedef struct {
    int start_j, end_j;
    double zmin, zmax;
    int march_axis;   
    int normal_axis; 
} draw_job_t;

static void render_3d_draw_job(void *arg) {
    draw_job_t *job = (draw_job_t *)arg;
    for (int j = job->start_j; j < job->end_j; j++) {
        for (int i = 0; i < GRID_3D; i++) {
            if (surf[j][i].count == 0) continue;

            for (int s = 0; s < surf[j][i].count; s++) {
                int sx0 = surf[j][i].sx[s], sy0 = surf[j][i].sy[s], dz0 = surf[j][i].dz[s];

                double height_val = (job->march_axis == 0) ? surf[j][i].z[s] : surf_y_3d[j][i];
                uint32_t col = color_by_height(height_val, job->zmin, job->zmax);
                if (filled_mode && job->normal_axis >= 0) {
                    float anx = (float)fabs(surf[j][i].nx[s]);
                    float any = (float)fabs(surf[j][i].ny[s]);
                    float anz = (float)fabs(surf[j][i].nz[s]);
                    bool dominant;
                    switch (job->normal_axis) {
                        case 0: dominant = (anz >= anx - 0.05f) && (anz >= any - 0.05f); break;
                        case 1: dominant = (anx >= any - 0.05f) && (anx >= anz - 0.05f); break;
                        case 2: dominant = (any >= anx - 0.05f) && (any >= anz - 0.05f); break;
                        default: dominant = true; break;
                    }
                    if (!dominant) continue;
                    
                    // Depth bias to prevent Z-fighting between passes
                    int bias = (job->normal_axis == 0) ? 0 : (job->normal_axis == 1) ? 2 : 4;
                    dz0 += bias; 
                }

                // Refined neighbor selection: only connect if points are world-space neighbors
                float world_step = (float)(job->zmax - job->zmin) / (GRID_3D - 1); // rough scaling
                if (world_step < 0.1f) world_step = 0.5f;
                float thresh = world_step * 2.5f;

                int s_tr = -1;
                if (i+1 < GRID_3D) {
                    float mind = 1e30f;
                    for (int n=0; n < surf[j][i+1].count; n++) {
                        float d = (float)fabs(surf[j][i+1].z[n] - surf[j][i].z[s]);
                        if (d < mind) { mind = d; s_tr = n; }
                    }
                    if (mind > thresh) s_tr = -1;
                }
                int s_bl = -1;
                if (j+1 < GRID_3D) {
                    float mind = 1e30f;
                    for (int n=0; n < surf[j+1][i].count; n++) {
                        float d = (float)fabs(surf[j+1][i].z[n] - surf[j][i].z[s]);
                        if (d < mind) { mind = d; s_bl = n; }
                    }
                    if (mind > thresh) s_bl = -1;
                }
                int s_br = -1;
                if (i+1 < GRID_3D && j+1 < GRID_3D) {
                    float mind = 1e30f;
                    for (int n=0; n < surf[j+1][i+1].count; n++) {
                        float d = (float)fabs(surf[j+1][i+1].z[n] - surf[j][i].z[s]);
                        if (d < mind) { mind = d; s_br = n; }
                    }
                    if (mind > thresh) s_br = -1;
                }

                if (filled_mode) {
                    bool v_tr = (s_tr >= 0);
                    bool v_bl = (s_bl >= 0);
                    bool v_br = (s_br >= 0);

                    if (v_tr && v_bl && v_br) {
                        int bias = (job->normal_axis == 0) ? 0 : (job->normal_axis == 1) ? 2 : 4;
                        int sx_tr = surf[j][i+1].sx[s_tr], sy_tr = surf[j][i+1].sy[s_tr], dz_tr = surf[j][i+1].dz[s_tr] + bias;
                        int sx_bl = surf[j+1][i].sx[s_bl], sy_bl = surf[j+1][i].sy[s_bl], dz_bl = surf[j+1][i].dz[s_bl] + bias;
                        int sx_br = surf[j+1][i+1].sx[s_br], sy_br = surf[j+1][i+1].sy[s_br], dz_br = surf[j+1][i+1].dz[s_br] + bias;

                        float avg_nx = surf[j][i].nx[s] + surf[j][i+1].nx[s_tr] + surf[j+1][i].nx[s_bl] + surf[j+1][i+1].nx[s_br];
                        float avg_ny = surf[j][i].ny[s] + surf[j][i+1].ny[s_tr] + surf[j+1][i].ny[s_bl] + surf[j+1][i+1].ny[s_br];
                        float avg_nz = surf[j][i].nz[s] + surf[j][i+1].nz[s_tr] + surf[j+1][i].nz[s_bl] + surf[j+1][i+1].nz[s_br];
                        avg_nx *= 0.25f; avg_ny *= 0.25f; avg_nz *= 0.25f;

                        float nlen = (float)sqrt(avg_nx*avg_nx + avg_ny*avg_ny + avg_nz*avg_nz);
                        if (nlen > 1e-9) { avg_nx /= nlen; avg_ny /= nlen; avg_nz /= nlen; }
                        else { avg_nx = 0; avg_ny = 0; avg_nz = 1; }

                        uint32_t scol = apply_shading(col, avg_nx, avg_ny, avg_nz);
                        gfb_triangle_z(sx0, sy0, dz0, sx_tr, sy_tr, dz_tr, sx_bl, sy_bl, dz_bl, scol);
                        gfb_triangle_z(sx_tr, sy_tr, dz_tr, sx_br, sy_br, dz_br, sx_bl, sy_bl, dz_bl, scol);
                    } else if (v_tr && v_bl) {
                        uint32_t scol = apply_shading(col, surf[j][i].nx[s], surf[j][i].ny[s], surf[j][i].nz[s]);
                        gfb_triangle_z(sx0, sy0, dz0,
                                       surf[j][i+1].sx[s_tr], surf[j][i+1].sy[s_tr], surf[j][i+1].dz[s_tr],
                                       surf[j+1][i].sx[s_bl], surf[j+1][i].sy[s_bl], surf[j+1][i].dz[s_bl], scol);
                    } else if (v_tr && v_br) {
                        uint32_t scol = apply_shading(col, surf[j][i].nx[s], surf[j][i].ny[s], surf[j][i].nz[s]);
                        gfb_triangle_z(sx0, sy0, dz0,
                                       surf[j][i+1].sx[s_tr], surf[j][i+1].sy[s_tr], surf[j][i+1].dz[s_tr],
                                       surf[j+1][i+1].sx[s_br], surf[j+1][i+1].sy[s_br], surf[j+1][i+1].dz[s_br], scol);
                    } else if (v_bl && v_br) {
                        uint32_t scol = apply_shading(col, surf[j][i].nx[s], surf[j][i].ny[s], surf[j][i].nz[s]);
                        gfb_triangle_z(sx0, sy0, dz0,
                                       surf[j+1][i].sx[s_bl], surf[j+1][i].sy[s_bl], surf[j+1][i].dz[s_bl],
                                       surf[j+1][i+1].sx[s_br], surf[j+1][i+1].sy[s_br], surf[j+1][i+1].dz[s_br], scol);
                    }
                } else {
                    if (i + 1 < GRID_3D && s_tr >= 0) {
                        gfb_line_z(sx0, sy0, dz0, surf[j][i+1].sx[s_tr], surf[j][i+1].sy[s_tr], surf[j][i+1].dz[s_tr], col);
                    }
                    if (j + 1 < GRID_3D && s_bl >= 0) {
                        gfb_line_z(sx0, sy0, dz0, surf[j+1][i].sx[s_bl], surf[j+1][i].sy[s_bl], surf[j+1][i].dz[s_bl], col);
                    }
                }
            }


        }
    }
}

static void render_3d_explicit(void) {
    double step = range_3d * 2.0 * 1.05 / (GRID_3D - 1);
    double zmin = 1e30, zmax = -1e30;

    if (surface_needs_eval) {
        for (int j = 0; j < GRID_3D; j++) { 
            for (int i = 0; i < GRID_3D; i++) { 
                surf[j][i].count = 0;
            } 
        }

        int num_chunks = 4;
        eval_job_t jobs[4];
        void *job_args[4];
        int rows_per_chunk = GRID_3D / num_chunks;
        
        for (int c = 0; c < num_chunks; c++) {
            jobs[c].start_j = c * rows_per_chunk;
            jobs[c].end_j = (c == num_chunks - 1) ? GRID_3D : (c + 1) * rows_per_chunk;
            jobs[c].range = range_3d;
            jobs[c].step = step;
            jobs[c].march_axis = 0;
            job_args[c] = &jobs[c];
        }

        sys_parallel_run(eval_3d_explicit_job, job_args, num_chunks);
    }
    
    // Compute min/max for coloring
    for (int j = 0; j < GRID_3D; j++) {
        for (int i = 0; i < GRID_3D; i++) {
            for (int s = 0; s < surf[j][i].count; s++) {
                if (surf[j][i].z[s] < zmin) zmin = surf[j][i].z[s];
                if (surf[j][i].z[s] > zmax) zmax = surf[j][i].z[s];
            }
        }
    }

    double z_scale = 1.0;
    if (zmax > zmin && (zmax - zmin) > 0.001) {
        z_scale = (range_3d * 2.0) / (zmax - zmin);
    }

    {
        int num_chunks = 4;
        eval_job_t jobs[4];
        void *job_args[4];
        int rows_per_chunk = GRID_3D / num_chunks;
        for (int c = 0; c < num_chunks; c++) {
            jobs[c].start_j = c * rows_per_chunk;
            jobs[c].end_j = (c == num_chunks - 1) ? GRID_3D : (c + 1) * rows_per_chunk;
            jobs[c].range = range_3d;
            jobs[c].step = step;
            jobs[c].z_scale = z_scale;
            jobs[c].march_axis = 0;
            job_args[c] = &jobs[c];
        }
        sys_parallel_run(eval_3d_explicit_project_job, job_args, num_chunks);
    }
    
    {
        int num_chunks = 4;
        draw_job_t jobs[4];
        void *job_args[4];
        int rows_per_chunk = GRID_3D / num_chunks;
        for (int c = 0; c < num_chunks; c++) {
            jobs[c].start_j = c * rows_per_chunk;
            jobs[c].end_j = (c == num_chunks - 1) ? GRID_3D : (c + 1) * rows_per_chunk;
            jobs[c].zmin = zmin;
            jobs[c].zmax = zmax;
            jobs[c].march_axis = 0;
            jobs[c].normal_axis = -1;
            job_args[c] = &jobs[c];
        }
        sys_parallel_run(render_3d_draw_job, job_args, num_chunks);
    }
}

static void render_3d_implicit(void) {
    double step = range_3d * 2.0 * 1.05 / (GRID_3D - 1);

    if (surface_needs_eval) {
        int num_chunks = 4, rows_per_chunk = GRID_3D / num_chunks;
        eval_job_t jobs[4]; void *job_args[4];
        for (int axis = 0; axis < 3; axis++) {
            for (int c = 0; c < num_chunks; c++) {
                jobs[c].start_j    = c * rows_per_chunk;
                jobs[c].end_j      = (c == num_chunks-1) ? GRID_3D : (c+1)*rows_per_chunk;
                jobs[c].range      = range_3d;
                jobs[c].step       = step;
                jobs[c].z_scale    = 1.0;
                jobs[c].march_axis = axis;
                job_args[c] = &jobs[c];
            }
            sys_parallel_run(eval_3d_implicit_job, job_args, num_chunks);
        }
        surface_needs_eval = false;
    }

    double zmin = 1e30, zmax = -1e30;
    for (int j = 0; j < GRID_3D; j++)
        for (int i = 0; i < GRID_3D; i++) {
            int cnt = eval_cache[0][j][i].count;
            for (int s = 0; s < cnt; s++) {
                float z = eval_cache[0][j][i].c[s];
                if (z < zmin) zmin = z;
                if (z > zmax) zmax = z;
            }
        }
    if (zmin > zmax) { zmin = -(float)range_3d; zmax = (float)range_3d; }

    static const int normal_axes[3] = {0, 1, 2};
    int num_chunks = 4, rows_per_chunk = GRID_3D / num_chunks;

    for (int axis = 0; axis < 3; axis++) {
        {
            eval_job_t jobs[4]; void *job_args[4];
            for (int c = 0; c < num_chunks; c++) {
                jobs[c].start_j    = c * rows_per_chunk;
                jobs[c].end_j      = (c == num_chunks-1) ? GRID_3D : (c+1)*rows_per_chunk;
                jobs[c].range      = range_3d;
                jobs[c].step       = step;
                jobs[c].z_scale    = 1.0;
                jobs[c].march_axis = axis;
                job_args[c] = &jobs[c];
            }
            sys_parallel_run(eval_3d_project_job, job_args, num_chunks);
        }
        {
            draw_job_t jobs[4]; void *job_args[4];
            for (int c = 0; c < num_chunks; c++) {
                jobs[c].start_j    = c * rows_per_chunk;
                jobs[c].end_j      = (c == num_chunks-1) ? GRID_3D : (c+1)*rows_per_chunk;
                jobs[c].zmin       = zmin;
                jobs[c].zmax       = zmax;
                jobs[c].march_axis = axis;
                jobs[c].normal_axis = filled_mode ? normal_axes[axis] : -1;
                job_args[c] = &jobs[c];
            }
            sys_parallel_run(render_3d_draw_job, job_args, num_chunks);
        }
    }
}

static void double_to_str(double val, char *buf) {
    if (val != val) { strcpy(buf, "NaN"); return; }
    if (val > 1e15) { strcpy(buf, "inf"); return; }
    if (val < -1e15) { strcpy(buf, "-inf"); return; }
    
    if (val > 2e9) {
        strcpy(buf, ">2G");
        return;
    }
    if (val < -2e9) {
        strcpy(buf, "<-2G");
        return;
    }

    if (val < 0) { *buf++ = '-'; val = -val; }
    
    int ipart = (int)val;
    itoa(ipart, buf);
    while (*buf) buf++;
    double frac = val - (double)ipart;
    if (frac > 0.005) {
        *buf++ = '.';
        int d1 = (int)(frac * 10) % 10;
        int d2 = (int)(frac * 100) % 10;
        *buf++ = '0' + d1;
        if (d2) *buf++ = '0' + d2;
    }
    *buf = 0;
}

static void render_graph(void) {
    gfb_clear(COLOR_BG);

    if (!eq_valid) {
        if (graph_mode == MODE_2D) render_2d_grid();
        ui_draw_image(win_graph, 0, GRAPH_Y, graph_w, graph_h, graph_fb);
        return;
    }

    if (graph_mode == MODE_2D) {
        if (is_explicit_2d) render_2d_explicit(); 
        render_2d_grid();
        if (is_explicit_2d) {
            int prev_sx = -1, prev_sy = -1;
            bool prev_valid = false;
            for (int px = 0; px < graph_w; px++) {
                double wx = world_x_2d(px);
                double wy = eval_rhs_only(wx, 0, 0);
                if (wy != wy || fabs(wy) > 1e10) { prev_valid = false; continue; }
                int sy = screen_y_2d(wy);
                if (prev_valid && fabs((double)(sy - prev_sy)) < graph_h)
                    gfb_line(prev_sx, prev_sy, px, sy, COLOR_CURVE);
                prev_sx = px; prev_sy = sy; prev_valid = true;
            }
        } else render_2d_implicit();
    } else {
        render_3d_axes();
        if (is_explicit_3d) render_3d_explicit();
        else render_3d_implicit();
    }

    if (graph_mode == MODE_3D && graph_czb) {
        int total = graph_w * graph_h;
        for (int i = 0; i < total; i++) {
            graph_fb[i] = (uint32_t)(graph_czb[i] & 0xFFFFFFFF);
        }
    }

    ui_draw_image(win_graph, 0, GRAPH_Y, graph_w, graph_h, graph_fb);

    if (graph_mode == MODE_2D) {
        double x_step = get_nice_step(view_x_max - view_x_min, 8);
        double y_step = get_nice_step(view_y_max - view_y_min, 6);

        int axis_y = screen_y_2d(0);
        if (axis_y < 10) axis_y = 10;
        if (axis_y > graph_h - 20) axis_y = graph_h - 20;
        axis_y += GRAPH_Y;

        double x_start = (int)(view_x_min / x_step) * x_step;
        int safety = 0;
        for (double wx = x_start - x_step; wx <= view_x_max + x_step && safety < 100; wx += x_step, safety++) {
            if (fabs(wx) < x_step * 0.1) continue; 
            int sx = screen_x_2d(wx);
            if (sx > 20 && sx < graph_w - 40) {
                char buf[32]; double_to_str(wx, buf);
                ui_draw_string(win_graph, sx - 10, axis_y + 4, buf, COLOR_TEXT);
            }
        }

        int axis_x = screen_x_2d(0);
        if (axis_x < 5) axis_x = 5;
        if (axis_x > graph_w - 40) axis_x = graph_w - 40;

        double y_start = (int)(view_y_min / y_step) * y_step;
        safety = 0;
        for (double wy = y_start - y_step; wy <= view_y_max + y_step && safety < 100; wy += y_step, safety++) {
            if (fabs(wy) < y_step * 0.1) continue; 
            int sy = screen_y_2d(wy);
            if (sy > 20 && sy < graph_h - 20) {
                char buf[32]; double_to_str(wy, buf);
                ui_draw_string(win_graph, axis_x + 6, sy + GRAPH_Y - 5, buf, COLOR_TEXT);
            }
        }
        
        int zx = screen_x_2d(0), zy = screen_y_2d(0);
        if (zx >= 0 && zx < graph_w - 10 && zy >= 0 && zy < graph_h - 10) {
            ui_draw_string(win_graph, zx + 4, zy + GRAPH_Y + 4, "0", COLOR_TEXT);
        }
    }
}

// =====
// Paint
// =====
static void paint_all(void) {
    rot_cx = cos(rot_x); rot_sx = sin(rot_x);
    rot_cy = cos(rot_y); rot_sy = sin(rot_y);

    ui_draw_rect(win_graph, 0, 0, win_w, TOOLBAR_H, COLOR_TOOLBAR_BG);
    widget_textbox_draw(&wctx, &tb_equation);
    widget_button_draw(&wctx, &btn_plot);
    widget_button_draw(&wctx, &btn_mode);
    ui_draw_rounded_rect_filled(win_graph, win_w - 80, 4, 70, 22, 4, 0xFF3A3A5A);
    ui_draw_string(win_graph, win_w - 72, 8, "Presets", COLOR_DARK_TEXT);
    render_graph();

    // Status bar
    int sty = GRAPH_Y + graph_h;
    ui_draw_rect(win_graph, 0, sty, win_w, STATUSBAR_H, COLOR_STATUS_BG);
    ui_draw_string(win_graph, 10, sty + 8,
                   graph_mode == MODE_3D ? "3D | CPU INTENSIVE!!" :
                                           "2D | Scroll=Zoom", COLOR_TEXT);

    char range_buf[64];
    if (graph_mode == MODE_2D) {
        strcpy(range_buf, "x:[");
        char tmp[16]; double_to_str(view_x_min, tmp); strcat(range_buf, tmp);
        strcat(range_buf, ",");
        double_to_str(view_x_max, tmp); strcat(range_buf, tmp);
        strcat(range_buf, "]");
        ui_draw_string(win_graph, win_w - 150, sty + 8, range_buf, COLOR_TEXT);
    } else {
        strcpy(range_buf, "Range:");
        char tmp[16]; double_to_str(range_3d, tmp); strcat(range_buf, tmp);
        ui_draw_string(win_graph, win_w - 200, sty + 8, range_buf, COLOR_TEXT);
        widget_button_draw(&wctx, &btn_range_minus);
        widget_button_draw(&wctx, &btn_range_plus);
    }

    if (presets_open) {
        int px = win_w - 150, py = TOOLBAR_H;
        ui_draw_rounded_rect_filled(win_graph, px, py, 140, NUM_PRESETS * 20 + 4, 4, COLOR_DARK_PANEL);
        for (int i = 0; i < NUM_PRESETS; i++) {
            ui_draw_string(win_graph, px + 8, py + 4 + i * 20, preset_labels[i], COLOR_DARK_TEXT);
        }
    }

    ui_mark_dirty(win_graph, 0, 0, win_w, CLIENT_H);
    surface_needs_eval = false;
}

// ====
// Zoom
// ====
static void reset_view(void) {
    if (graph_mode == MODE_2D) {
        view_x_min = -10.0; view_x_max = 10.0;
        view_y_min = -6.4; view_y_max = 6.4;
        apply_aspect_ratio();
    } else {
        zoom_3d = 1.0;
        rot_x = -0.5;
        rot_y = 0.5;
        range_3d = 10.0;
    }
}

static void zoom_2d(double factor) {
    if (factor <= 0) return;
    
    double cx = (view_x_min + view_x_max) / 2.0;
    double cy = (view_y_min + view_y_max) / 2.0;
    
    double half_x = (view_x_max - view_x_min) / 2.0;
    double half_y = (view_y_max - view_y_min) / 2.0;
    
    half_x *= factor;
    half_y *= factor;
    
    // Safety caps for zoom range
    if (half_x < 1e-12) half_x = 1e-12;
    if (half_x > 1e12) half_x = 1e12;
    if (half_y < 1e-12) half_y = 1e-12;
    if (half_y > 1e12) half_y = 1e12;

    view_x_min = cx - half_x;
    view_x_max = cx + half_x;
    view_y_min = cy - half_y;
    view_y_max = cy + half_y;
    
    apply_aspect_ratio();
}

static void handle_scroll(int dz) {
    if (graph_mode == MODE_2D) {
        surface_needs_eval = true;
        if (dz > 0) zoom_2d(0.85);
        else zoom_2d(1.18);
    } else {
        if (dz > 0) zoom_3d *= 1.15;
        else zoom_3d *= 0.87;
    }
}

static void update_widget_layout(void) {
    tb_equation.w = win_w - 220;
    btn_plot.x = win_w - 200;
    widget_button_init(&btn_mode, win_w - 145, 4, 60, 22, filled_mode ? "Wire" : "Filled");
    
    int sty = win_h - STATUSBAR_H - 20;
    widget_button_init(&btn_range_minus, win_w - 95, sty + 4, 30, 22, "-");
    widget_button_init(&btn_range_plus, win_w - 55, sty + 4, 30, 22, "+");
}

// ====
// Main
// ====
int main(void) {
    win_graph = ui_window_create("Grapher", 80, 60, win_w, win_h);
    ui_window_set_resizable(win_graph, true);
    wctx.user_data = (void *)win_graph;

    fb_capacity = graph_w * graph_h;
    graph_fb = (uint32_t *)malloc(fb_capacity * sizeof(uint32_t));
    graph_zb = (int32_t *)malloc(fb_capacity * sizeof(int32_t));
    graph_czb = (uint64_t *)malloc(fb_capacity * sizeof(uint64_t));
    if (!graph_fb || !graph_zb || !graph_czb) return 1;

    memset(eq_buffer, 0, sizeof(eq_buffer));
    strcpy(eq_buffer, "y = sin(x)");
    eq_len = strlen(eq_buffer);

    widget_textbox_init(&tb_equation, 10, 4, 340, 22, eq_buffer, 255);
    tb_equation.cursor_pos = eq_len;
    widget_button_init(&btn_plot, 360, 4, 50, 22, "Plot");
    update_widget_layout();

    // Parse initial equation
    parse_equation();

    paint_all();

    gui_event_t ev;
    bool needs_repaint = false;
    while (1) {
        bool got_event = false;
        while (ui_get_event(win_graph, &ev)) {
            got_event = true;
            if (ev.type == GUI_EVENT_CLOSE) {
                sys_exit(0);
            } else if (ev.type == GUI_EVENT_PAINT) {
                needs_repaint = true;
            } else if (ev.type == GUI_EVENT_KEY) {
                uint32_t cp = (uint32_t)ev.arg4;
                char c = (char)ev.arg1;
                if (tb_equation.focused) {
                    if (cp == '\n') {
                        eq_len = (int)strlen(eq_buffer);
                        parse_equation();
                        if (graph_mode == MODE_2D) autofit_2d_view();
                        needs_repaint = true;
                    } else {
                        widget_textbox_handle_key(&tb_equation, cp, (int)ev.arg1, NULL);
                        eq_len = (int)strlen(eq_buffer);
                        needs_repaint = true;
                    }
                } else if ((c == 'r' || c == 'R') && ev.arg3 == 1) { 
                    reset_view();
                    surface_needs_eval = true;
                    needs_repaint = true;
                } else if ((c == 'f' || c == 'F')) {
                    filled_mode = !filled_mode;
                    update_widget_layout();
                    needs_repaint = true;
                }
            } else if (ev.type == GUI_EVENT_RESIZE) {
                win_w = ev.arg1;
                win_h = ev.arg2;
                graph_w = win_w;
                graph_h = (win_h - 20) - TOOLBAR_H - STATUSBAR_H;
                if (graph_h < 50) graph_h = 50;
                
                int req_cap = graph_w * graph_h;
                if (req_cap > fb_capacity) {
                    if (graph_fb) free(graph_fb);
                    if (graph_zb) free(graph_zb);
                    if (graph_czb) free(graph_czb);
                    fb_capacity = (int)(req_cap * 1.5);
                    graph_fb = (uint32_t *)malloc(fb_capacity * sizeof(uint32_t));
                    graph_zb = (int32_t *)malloc(fb_capacity * sizeof(int32_t));
                    graph_czb = (uint64_t *)malloc(fb_capacity * sizeof(uint64_t));
                }
                
                update_widget_layout();
                
                if (graph_mode == MODE_2D) {
                    apply_aspect_ratio();
                }
                surface_needs_eval = true;
                needs_repaint = true;
            } else if (ev.type == GUI_EVENT_CLICK) {
                int mx = ev.arg1, my = ev.arg2;

                if (presets_open) {
                    int px = win_w - 150, py = TOOLBAR_H;
                    if (mx >= px && mx < px + 140 && my >= py && my < py + NUM_PRESETS * 20 + 4) {
                        int idx = (my - py - 2) / 20;
                        if (idx >= 0 && idx < NUM_PRESETS) {
                            strcpy(eq_buffer, preset_labels[idx]);
                            eq_len = strlen(eq_buffer);
                            tb_equation.cursor_pos = eq_len;
                            parse_equation();
                            if (graph_mode == MODE_2D) autofit_2d_view();
                        }
                    }
                    presets_open = false;
                    needs_repaint = true;
                    continue;
                }

                if (mx >= win_w - 80 && mx < win_w - 10 && my >= 4 && my < 26) {
                    presets_open = !presets_open;
                    needs_repaint = true;
                    continue;
                }

                if (widget_button_handle_mouse(&btn_plot, mx, my, false, true, NULL)) {
                    parse_equation();
                    if (graph_mode == MODE_2D) autofit_2d_view();
                    needs_repaint = true;
                    continue;
                }
                
                if (widget_button_handle_mouse(&btn_mode, mx, my, false, true, NULL)) {
                    filled_mode = !filled_mode;
                    update_widget_layout();
                    needs_repaint = true;
                    continue;
                }
                
                if (graph_mode == MODE_3D) {
                    if (widget_button_handle_mouse(&btn_range_plus, mx, my, false, true, NULL)) {
                        range_3d *= 1.25;
                        surface_needs_eval = true;
                        needs_repaint = true;
                        continue;
                    }
                    if (widget_button_handle_mouse(&btn_range_minus, mx, my, false, true, NULL)) {
                        range_3d *= 0.8;
                        surface_needs_eval = true;
                        needs_repaint = true;
                        continue;
                    }
                }

                widget_textbox_handle_mouse(&wctx, &tb_equation, mx, my, true, NULL);
                needs_repaint = true;

            } else if (ev.type == GUI_EVENT_MOUSE_DOWN) {
                int mx = ev.arg1, my = ev.arg2;
                widget_button_handle_mouse(&btn_plot, mx, my, true, false, NULL);
                widget_button_handle_mouse(&btn_mode, mx, my, true, false, NULL);
                if (graph_mode == MODE_3D) {
                    widget_button_handle_mouse(&btn_range_plus, mx, my, true, false, NULL);
                    widget_button_handle_mouse(&btn_range_minus, mx, my, true, false, NULL);
                }
                needs_repaint = true;

            } else if (ev.type == GUI_EVENT_MOUSE_UP) {
                int mx = ev.arg1, my = ev.arg2;
                widget_button_handle_mouse(&btn_plot, mx, my, false, false, NULL);
                widget_button_handle_mouse(&btn_mode, mx, my, false, false, NULL);
                if (graph_mode == MODE_3D) {
                    widget_button_handle_mouse(&btn_range_plus, mx, my, false, false, NULL);
                    widget_button_handle_mouse(&btn_range_minus, mx, my, false, false, NULL);
                }
                needs_repaint = true;

            } else if (ev.type == GUI_EVENT_RIGHT_CLICK) {
                if (graph_mode == MODE_3D) {
                    right_dragging = true;
                    drag_last_x = ev.arg1;
                    drag_last_y = ev.arg2;
                }
            } else if (ev.type == GUI_EVENT_MOUSE_MOVE) {
                int mx = ev.arg1, my = ev.arg2;
                int buttons = ev.arg3;

                if (graph_mode == MODE_3D && (buttons & 2)) {
                    double dx = mx - drag_last_x;
                    double dy = my - drag_last_y;
                    rot_y += dx * 0.01;
                    rot_x += dy * 0.01;
                    drag_last_x = mx;
                    drag_last_y = my;
                    needs_repaint = true;
                } else {
                    right_dragging = false;
                }
            } else if (ev.type == 9) { 
                handle_scroll(ev.arg1);
                if (eq_valid) needs_repaint = true;
            }
        }
        
        if (graph_mode == MODE_3D && ROTATE == 1) {
            rot_y += 0.01;
            needs_repaint = true;
        }

        if (needs_repaint) {
            paint_all();
            needs_repaint = false;
        }

        if (!got_event) {
            sleep(16);
        }
    }

    free(graph_fb);
    sys_exit(0);
    return 0;
}
