Init Project

pkg/mathkit/calculus.go

@@ -0,0 +1,115 @@
+package mathkit
+
+import "fmt"
+
+// defaultStep is the finite-difference step size used when the
+// caller doesn't specify one. 1e-5 is the standard sweet spot for
+// float64 central differences: small enough for truncation error to
+// be negligible, large enough that it doesn't get swallowed by
+// float64 rounding error (the classic bias-variance trade-off of
+// numerical differentiation).
+const defaultStep = 1e-5
+
+// Derivative estimates d/d(variable) of expr at the point `at`
+// (variable -> value, including any other free variables the
+// expression needs) using the central-difference formula
+// (f(x+h)-f(x-h))/(2h), which is second-order accurate — far better
+// than the naive forward difference (f(x+h)-f(x))/h for the same h.
+// h<=0 uses defaultStep.
+func Derivative(expr *Expr, variable string, at map[string]float64, h float64) (float64, error) {
+	if h <= 0 {
+		h = defaultStep
+	}
+	x0, ok := at[variable]
+	if !ok {
+		return 0, fmt.Errorf("derivative: %q is not bound in \"at\"", variable)
+	}
+	plus := cloneVars(at)
+	plus[variable] = x0 + h
+	minus := cloneVars(at)
+	minus[variable] = x0 - h
+	fPlus, err := expr.Eval(plus)
+	if err != nil {
+		return 0, fmt.Errorf("derivative: evaluating at %s+h: %w", variable, err)
+	}
+	fMinus, err := expr.Eval(minus)
+	if err != nil {
+		return 0, fmt.Errorf("derivative: evaluating at %s-h: %w", variable, err)
+	}
+	return (fPlus - fMinus) / (2 * h), nil
+}
+
+// Gradient estimates ∇f at the point `at` with respect to each name
+// in variables, in order — the central operation every ML training
+// loop runs millions of times (just normally via backprop's chain
+// rule rather than finite differences; this is the numerical
+// stand-in that's useful for checking a hand-derived gradient, or for
+// expressions where deriving the symbolic gradient isn't worth it).
+func Gradient(expr *Expr, variables []string, at map[string]float64, h float64) ([]float64, error) {
+	if len(variables) == 0 {
+		return nil, fmt.Errorf("gradient: \"variables\" must list at least one variable")
+	}
+	grad := make([]float64, len(variables))
+	for i, v := range variables {
+		d, err := Derivative(expr, v, at, h)
+		if err != nil {
+			return nil, err
+		}
+		grad[i] = d
+	}
+	return grad, nil
+}
+
+func cloneVars(m map[string]float64) map[string]float64 {
+	out := make(map[string]float64, len(m))
+	for k, v := range m {
+		out[k] = v
+	}
+	return out
+}
+
+// Integrate estimates the definite integral of expr with respect to
+// variable over [from, to] using Simpson's rule, which is exact for
+// up to cubic polynomials and converges much faster than the
+// trapezoid rule for smooth functions. steps (the number of
+// subintervals) is rounded up to the nearest even number, since
+// Simpson's rule needs an even number of subintervals; steps<=0
+// defaults to 1000.
+func Integrate(expr *Expr, variable string, from, to float64, steps int, fixed map[string]float64) (float64, error) {
+	if steps <= 0 {
+		steps = 1000
+	}
+	if steps%2 != 0 {
+		steps++
+	}
+	h := (to - from) / float64(steps)
+	vars := cloneVars(fixed)
+
+	eval := func(x float64) (float64, error) {
+		vars[variable] = x
+		return expr.Eval(vars)
+	}
+
+	f0, err := eval(from)
+	if err != nil {
+		return 0, fmt.Errorf("integral: %w", err)
+	}
+	fn, err := eval(to)
+	if err != nil {
+		return 0, fmt.Errorf("integral: %w", err)
+	}
+	sum := f0 + fn
+	for i := 1; i < steps; i++ {
+		x := from + float64(i)*h
+		fx, err := eval(x)
+		if err != nil {
+			return 0, fmt.Errorf("integral: %w", err)
+		}
+		if i%2 == 0 {
+			sum += 2 * fx
+		} else {
+			sum += 4 * fx
+		}
+	}
+	return sum * h / 3, nil
+}

pkg/mathkit/expr.go

@@ -0,0 +1,421 @@
+// Package mathkit is a dependency-free math engine for the problems
+// that actually come up when building and debugging ML code: linear
+// algebra (the shape of every tensor op), calculus (gradients are how
+// everything trains), probability/statistics (data, metrics, losses),
+// numerical optimization (gradient descent, root finding), and the
+// handful of ML-specific primitives (activations, loss functions,
+// least-squares linear regression) that sit on top of those four
+// pillars. It has no third-party dependency — see the comment on
+// pkg/envkit/dotenv.go for why this module keeps go.mod minimal — so
+// the whole engine is just the Go standard library's math/sort/fmt.
+//
+// Every sub-area (linalg.go, calculus.go, statistics.go, optimize.go,
+// ml.go) is usable directly from Go code. solve.go adds one more
+// layer on top: Solve(operation, params) dispatches a string
+// operation name and a loosely-typed params map (exactly the shape
+// a JSON tool-call argument object decodes to) to the right function
+// and renders a human-readable report — that's the layer the agent
+// service's solve_ml_math tool calls (see
+// services/agent/internal/infrastructure/toolexec/mlmath.go), because
+// an LLM tool call arrives as JSON, not as typed Go arguments.
+package mathkit
+
+import (
+	"fmt"
+	"math"
+	"strings"
+	"unicode"
+)
+
+// Expr is a parsed arithmetic expression that can be evaluated
+// repeatedly against different variable bindings without re-parsing —
+// exactly what Gradient and GradientDescent need, since they evaluate
+// the same expression at many points.
+type Expr struct {
+	root node
+	src  string
+}
+
+// String returns the original source text the expression was parsed
+// from, for error messages and report rendering.
+func (e *Expr) String() string { return e.src }
+
+// node is the parsed AST. A tiny tagged union kept as an interface
+// rather than a struct-with-a-kind-field: each node type carries only
+// the fields it needs, and eval is a single recursive method.
+type node interface {
+	eval(vars map[string]float64) (float64, error)
+}
+
+type numNode float64
+
+func (n numNode) eval(map[string]float64) (float64, error) { return float64(n), nil }
+
+type varNode string
+
+func (n varNode) eval(vars map[string]float64) (float64, error) {
+	v, ok := vars[string(n)]
+	if !ok {
+		return 0, fmt.Errorf("unbound variable %q (pass it in \"variables\" or \"at\")", string(n))
+	}
+	return v, nil
+}
+
+type binNode struct {
+	op       byte // '+', '-', '*', '/', '^'
+	lhs, rhs node
+}
+
+func (n binNode) eval(vars map[string]float64) (float64, error) {
+	l, err := n.lhs.eval(vars)
+	if err != nil {
+		return 0, err
+	}
+	r, err := n.rhs.eval(vars)
+	if err != nil {
+		return 0, err
+	}
+	switch n.op {
+	case '+':
+		return l + r, nil
+	case '-':
+		return l - r, nil
+	case '*':
+		return l * r, nil
+	case '/':
+		if r == 0 {
+			return 0, fmt.Errorf("division by zero")
+		}
+		return l / r, nil
+	case '^':
+		return math.Pow(l, r), nil
+	}
+	return 0, fmt.Errorf("internal: unknown operator %q", n.op)
+}
+
+type negNode struct{ x node }
+
+func (n negNode) eval(vars map[string]float64) (float64, error) {
+	v, err := n.x.eval(vars)
+	return -v, err
+}
+
+type callNode struct {
+	fn   string
+	args []node
+}
+
+// unaryFuncs are functions of exactly one argument — the overwhelming
+// majority of what shows up in ML math (activations, logs, norms).
+var unaryFuncs = map[string]func(float64) float64{
+	"sin": math.Sin, "cos": math.Cos, "tan": math.Tan,
+	"asin": math.Asin, "acos": math.Acos, "atan": math.Atan,
+	"sinh": math.Sinh, "cosh": math.Cosh, "tanh": math.Tanh,
+	"exp": math.Exp, "sqrt": math.Sqrt, "abs": math.Abs,
+	"ln": math.Log, "log": math.Log, "log2": math.Log2, "log10": math.Log10,
+	"floor": math.Floor, "ceil": math.Ceil, "round": math.Round,
+	"sigmoid": Sigmoid, "relu": ReLU, "sign": func(x float64) float64 { return math.Copysign(1, x) * boolF(x != 0) },
+}
+
+func boolF(b bool) float64 {
+	if b {
+		return 1
+	}
+	return 0
+}
+
+// binaryFuncs are functions of exactly two arguments.
+var binaryFuncs = map[string]func(float64, float64) float64{
+	"max": math.Max, "min": math.Min, "pow": math.Pow,
+	"mod": math.Mod, "atan2": math.Atan2,
+}
+
+func (n callNode) eval(vars map[string]float64) (float64, error) {
+	switch len(n.args) {
+	case 1:
+		if f, ok := unaryFuncs[n.fn]; ok {
+			v, err := n.args[0].eval(vars)
+			if err != nil {
+				return 0, err
+			}
+			return f(v), nil
+		}
+	case 2:
+		if f, ok := binaryFuncs[n.fn]; ok {
+			a, err := n.args[0].eval(vars)
+			if err != nil {
+				return 0, err
+			}
+			b, err := n.args[1].eval(vars)
+			if err != nil {
+				return 0, err
+			}
+			return f(a, b), nil
+		}
+	}
+	return 0, fmt.Errorf("unknown function %q with %d argument(s)", n.fn, len(n.args))
+}
+
+// Eval evaluates the expression against a variable binding. Missing
+// variables are a reported error (never silently treated as zero),
+// since a silent zero would turn a typo into a wrong-but-plausible
+// answer — exactly the failure mode numerical tools must not have.
+func (e *Expr) Eval(vars map[string]float64) (float64, error) {
+	return e.root.eval(vars)
+}
+
+// ParseExpr parses a math expression written in ordinary infix
+// notation: + - * / ^ (right-associative power), unary minus,
+// parentheses, named variables, the constants pi/e, and function
+// calls (sin(x), sigmoid(x), max(a,b), ...). See unaryFuncs and
+// binaryFuncs for the full function list.
+func ParseExpr(src string) (*Expr, error) {
+	p := &parser{toks: tokenize(src)}
+	n, err := p.parseExpr()
+	if err != nil {
+		return nil, fmt.Errorf("parsing %q: %w", src, err)
+	}
+	if p.pos != len(p.toks) {
+		return nil, fmt.Errorf("parsing %q: unexpected trailing input near %q", src, p.toks[p.pos].text)
+	}
+	return &Expr{root: n, src: src}, nil
+}
+
+// --- tokenizer ---
+
+type tokKind int
+
+const (
+	tokNum tokKind = iota
+	tokIdent
+	tokOp
+	tokLParen
+	tokRParen
+	tokComma
+)
+
+type token struct {
+	kind tokKind
+	text string
+	num  float64
+}
+
+func tokenize(src string) []token {
+	var toks []token
+	r := []rune(src)
+	i := 0
+	for i < len(r) {
+		c := r[i]
+		switch {
+		case unicode.IsSpace(c):
+			i++
+		case c == '(':
+			toks = append(toks, token{kind: tokLParen, text: "("})
+			i++
+		case c == ')':
+			toks = append(toks, token{kind: tokRParen, text: ")"})
+			i++
+		case c == ',':
+			toks = append(toks, token{kind: tokComma, text: ","})
+			i++
+		case strings.ContainsRune("+-*/^", c):
+			toks = append(toks, token{kind: tokOp, text: string(c)})
+			i++
+		case unicode.IsDigit(c) || c == '.':
+			j := i
+			for j < len(r) && (unicode.IsDigit(r[j]) || r[j] == '.') {
+				j++
+			}
+			// scientific notation, e.g. 1e-3
+			if j < len(r) && (r[j] == 'e' || r[j] == 'E') {
+				k := j + 1
+				if k < len(r) && (r[k] == '+' || r[k] == '-') {
+					k++
+				}
+				if k < len(r) && unicode.IsDigit(r[k]) {
+					j = k
+					for j < len(r) && unicode.IsDigit(r[j]) {
+						j++
+					}
+				}
+			}
+			text := string(r[i:j])
+			var f float64
+			fmt.Sscanf(text, "%g", &f)
+			toks = append(toks, token{kind: tokNum, text: text, num: f})
+			i = j
+		case unicode.IsLetter(c) || c == '_':
+			j := i
+			for j < len(r) && (unicode.IsLetter(r[j]) || unicode.IsDigit(r[j]) || r[j] == '_') {
+				j++
+			}
+			toks = append(toks, token{kind: tokIdent, text: string(r[i:j])})
+			i = j
+		default:
+			// Unknown character: emit as its own operator token so the
+			// parser reports a clear "unexpected token" error instead
+			// of the tokenizer silently dropping it.
+			toks = append(toks, token{kind: tokOp, text: string(c)})
+			i++
+		}
+	}
+	return toks
+}
+
+// --- recursive-descent parser ---
+//
+// expr   := term (('+' | '-') term)*
+// term   := unary (('*' | '/') unary)*
+// unary  := ('-' | '+') unary | power
+// power  := primary ('^' unary)?         -- right-associative; the
+//           exponent recurses through unary (not power) so the
+//           precedence matches ordinary math/Python convention:
+//           unary minus binds LOOSER than '^' (-2^2 == -(2^2) == -4)
+//           but the exponent itself may carry its own sign
+//           (2^-2 == 0.25).
+// primary:= number | const | ident '(' args ')' | ident | '(' expr ')'
+// args   := expr (',' expr)*
+
+type parser struct {
+	toks []token
+	pos  int
+}
+
+func (p *parser) peek() (token, bool) {
+	if p.pos >= len(p.toks) {
+		return token{}, false
+	}
+	return p.toks[p.pos], true
+}
+
+func (p *parser) next() (token, bool) {
+	t, ok := p.peek()
+	if ok {
+		p.pos++
+	}
+	return t, ok
+}
+
+func (p *parser) parseExpr() (node, error) {
+	lhs, err := p.parseTerm()
+	if err != nil {
+		return nil, err
+	}
+	for {
+		t, ok := p.peek()
+		if !ok || t.kind != tokOp || (t.text != "+" && t.text != "-") {
+			return lhs, nil
+		}
+		p.next()
+		rhs, err := p.parseTerm()
+		if err != nil {
+			return nil, err
+		}
+		lhs = binNode{op: t.text[0], lhs: lhs, rhs: rhs}
+	}
+}
+
+func (p *parser) parseTerm() (node, error) {
+	lhs, err := p.parseUnary()
+	if err != nil {
+		return nil, err
+	}
+	for {
+		t, ok := p.peek()
+		if !ok || t.kind != tokOp || (t.text != "*" && t.text != "/") {
+			return lhs, nil
+		}
+		p.next()
+		rhs, err := p.parseUnary()
+		if err != nil {
+			return nil, err
+		}
+		lhs = binNode{op: t.text[0], lhs: lhs, rhs: rhs}
+	}
+}
+
+func (p *parser) parseUnary() (node, error) {
+	if t, ok := p.peek(); ok && t.kind == tokOp && (t.text == "-" || t.text == "+") {
+		p.next()
+		x, err := p.parseUnary()
+		if err != nil {
+			return nil, err
+		}
+		if t.text == "-" {
+			return negNode{x: x}, nil
+		}
+		return x, nil
+	}
+	return p.parsePower()
+}
+
+func (p *parser) parsePower() (node, error) {
+	lhs, err := p.parsePrimary()
+	if err != nil {
+		return nil, err
+	}
+	if t, ok := p.peek(); ok && t.kind == tokOp && t.text == "^" {
+		p.next()
+		rhs, err := p.parseUnary() // right-associative; lets the exponent carry its own sign
+		if err != nil {
+			return nil, err
+		}
+		return binNode{op: '^', lhs: lhs, rhs: rhs}, nil
+	}
+	return lhs, nil
+}
+
+func (p *parser) parsePrimary() (node, error) {
+	t, ok := p.next()
+	if !ok {
+		return nil, fmt.Errorf("unexpected end of expression")
+	}
+	switch t.kind {
+	case tokNum:
+		return numNode(t.num), nil
+	case tokLParen:
+		n, err := p.parseExpr()
+		if err != nil {
+			return nil, err
+		}
+		if c, ok := p.next(); !ok || c.kind != tokRParen {
+			return nil, fmt.Errorf("missing closing ')'")
+		}
+		return n, nil
+	case tokIdent:
+		switch strings.ToLower(t.text) {
+		case "pi":
+			return numNode(math.Pi), nil
+		case "e":
+			return numNode(math.E), nil
+		}
+		if nx, ok := p.peek(); ok && nx.kind == tokLParen {
+			p.next()
+			var args []node
+			if c, ok := p.peek(); !ok || c.kind != tokRParen {
+				for {
+					a, err := p.parseExpr()
+					if err != nil {
+						return nil, err
+					}
+					args = append(args, a)
+					c, ok := p.peek()
+					if !ok {
+						return nil, fmt.Errorf("missing closing ')' in call to %s", t.text)
+					}
+					if c.kind == tokComma {
+						p.next()
+						continue
+					}
+					break
+				}
+			}
+			if c, ok := p.next(); !ok || c.kind != tokRParen {
+				return nil, fmt.Errorf("missing closing ')' in call to %s", t.text)
+			}
+			return callNode{fn: strings.ToLower(t.text), args: args}, nil
+		}
+		return varNode(t.text), nil
+	}
+	return nil, fmt.Errorf("unexpected token %q", t.text)
+}

pkg/mathkit/linalg.go

@@ -0,0 +1,409 @@
+package mathkit
+
+import (
+	"fmt"
+	"math"
+)
+
+// Matrix is a dense matrix stored row-major: Matrix[i][j] is row i,
+// column j. Every row must have the same length — functions that take
+// a Matrix validate that with Matrix.Dims before doing any math, so a
+// ragged input is reported as an error rather than panicking.
+type Matrix [][]float64
+
+// Dims returns the matrix's row and column counts and an error if the
+// rows aren't all the same length (a ragged "matrix" isn't one).
+func (m Matrix) Dims() (rows, cols int, err error) {
+	rows = len(m)
+	if rows == 0 {
+		return 0, 0, fmt.Errorf("empty matrix")
+	}
+	cols = len(m[0])
+	for i, row := range m {
+		if len(row) != cols {
+			return 0, 0, fmt.Errorf("ragged matrix: row 0 has %d columns, row %d has %d", cols, i, len(row))
+		}
+	}
+	return rows, cols, nil
+}
+
+// Identity returns the n×n identity matrix.
+func Identity(n int) Matrix {
+	m := make(Matrix, n)
+	for i := range m {
+		m[i] = make([]float64, n)
+		m[i][i] = 1
+	}
+	return m
+}
+
+func sameDims(a, b Matrix) (rows, cols int, err error) {
+	ra, ca, err := a.Dims()
+	if err != nil {
+		return 0, 0, err
+	}
+	rb, cb, err := b.Dims()
+	if err != nil {
+		return 0, 0, err
+	}
+	if ra != rb || ca != cb {
+		return 0, 0, fmt.Errorf("shape mismatch: %dx%d vs %dx%d", ra, ca, rb, cb)
+	}
+	return ra, ca, nil
+}
+
+// MatAdd returns a+b, element-wise.
+func MatAdd(a, b Matrix) (Matrix, error) {
+	rows, cols, err := sameDims(a, b)
+	if err != nil {
+		return nil, fmt.Errorf("matrix_add: %w", err)
+	}
+	out := make(Matrix, rows)
+	for i := 0; i < rows; i++ {
+		out[i] = make([]float64, cols)
+		for j := 0; j < cols; j++ {
+			out[i][j] = a[i][j] + b[i][j]
+		}
+	}
+	return out, nil
+}
+
+// MatSub returns a-b, element-wise.
+func MatSub(a, b Matrix) (Matrix, error) {
+	rows, cols, err := sameDims(a, b)
+	if err != nil {
+		return nil, fmt.Errorf("matrix_sub: %w", err)
+	}
+	out := make(Matrix, rows)
+	for i := 0; i < rows; i++ {
+		out[i] = make([]float64, cols)
+		for j := 0; j < cols; j++ {
+			out[i][j] = a[i][j] - b[i][j]
+		}
+	}
+	return out, nil
+}
+
+// MatMul returns the matrix product a×b (standard inner-dimension
+// matching: a is m×n, b must be n×p).
+func MatMul(a, b Matrix) (Matrix, error) {
+	ra, ca, err := a.Dims()
+	if err != nil {
+		return nil, fmt.Errorf("matrix_mul: left operand: %w", err)
+	}
+	rb, cb, err := b.Dims()
+	if err != nil {
+		return nil, fmt.Errorf("matrix_mul: right operand: %w", err)
+	}
+	if ca != rb {
+		return nil, fmt.Errorf("matrix_mul: inner dimensions must match: %dx%d times %dx%d", ra, ca, rb, cb)
+	}
+	out := make(Matrix, ra)
+	for i := 0; i < ra; i++ {
+		out[i] = make([]float64, cb)
+		for j := 0; j < cb; j++ {
+			var sum float64
+			for k := 0; k < ca; k++ {
+				sum += a[i][k] * b[k][j]
+			}
+			out[i][j] = sum
+		}
+	}
+	return out, nil
+}
+
+// Scale multiplies every element of m by k.
+func Scale(m Matrix, k float64) Matrix {
+	out := make(Matrix, len(m))
+	for i, row := range m {
+		out[i] = make([]float64, len(row))
+		for j, v := range row {
+			out[i][j] = v * k
+		}
+	}
+	return out
+}
+
+// Transpose returns m^T.
+func Transpose(m Matrix) (Matrix, error) {
+	rows, cols, err := m.Dims()
+	if err != nil {
+		return nil, fmt.Errorf("matrix_transpose: %w", err)
+	}
+	out := make(Matrix, cols)
+	for j := 0; j < cols; j++ {
+		out[j] = make([]float64, rows)
+		for i := 0; i < rows; i++ {
+			out[j][i] = m[i][j]
+		}
+	}
+	return out, nil
+}
+
+// Trace returns the sum of the diagonal of a square matrix.
+func Trace(m Matrix) (float64, error) {
+	rows, cols, err := m.Dims()
+	if err != nil {
+		return 0, fmt.Errorf("matrix_trace: %w", err)
+	}
+	if rows != cols {
+		return 0, fmt.Errorf("matrix_trace: matrix must be square, got %dx%d", rows, cols)
+	}
+	var sum float64
+	for i := 0; i < rows; i++ {
+		sum += m[i][i]
+	}
+	return sum, nil
+}
+
+// clone makes a deep copy so Gaussian-elimination routines never
+// mutate the caller's matrix.
+func (m Matrix) clone() Matrix {
+	out := make(Matrix, len(m))
+	for i, row := range m {
+		out[i] = append([]float64(nil), row...)
+	}
+	return out
+}
+
+// luDecompose runs Gaussian elimination with partial pivoting on a
+// clone of m, returning the resulting upper-triangular-ish matrix and
+// the sign of the permutation (needed to get the determinant's sign
+// right) and the pivots used. Both Determinant and Inverse build on
+// this one routine rather than duplicating elimination logic.
+func luDecompose(m Matrix) (a Matrix, sign float64, err error) {
+	rows, cols, err := m.Dims()
+	if err != nil {
+		return nil, 0, err
+	}
+	if rows != cols {
+		return nil, 0, fmt.Errorf("matrix must be square, got %dx%d", rows, cols)
+	}
+	n := rows
+	a = m.clone()
+	sign = 1
+	for col := 0; col < n; col++ {
+		// Partial pivot: swap in the row with the largest absolute
+		// value in this column, for numerical stability.
+		pivotRow := col
+		best := math.Abs(a[col][col])
+		for r := col + 1; r < n; r++ {
+			if v := math.Abs(a[r][col]); v > best {
+				best, pivotRow = v, r
+			}
+		}
+		if best < 1e-12 {
+			return a, 0, nil // singular; determinant is 0, caller decides
+		}
+		if pivotRow != col {
+			a[col], a[pivotRow] = a[pivotRow], a[col]
+			sign = -sign
+		}
+		for r := col + 1; r < n; r++ {
+			factor := a[r][col] / a[col][col]
+			for c := col; c < n; c++ {
+				a[r][c] -= factor * a[col][c]
+			}
+		}
+	}
+	return a, sign, nil
+}
+
+// Determinant computes det(m) via LU (Gaussian elimination with
+// partial pivoting), O(n^3) — the standard numerically-stable
+// approach, versus the O(n!) cofactor expansion a direct reading of
+// the definition would suggest.
+func Determinant(m Matrix) (float64, error) {
+	a, sign, err := luDecompose(m)
+	if err != nil {
+		return 0, fmt.Errorf("matrix_determinant: %w", err)
+	}
+	n := len(a)
+	det := sign
+	for i := 0; i < n; i++ {
+		det *= a[i][i]
+	}
+	return det, nil
+}
+
+// Inverse computes m^-1 via Gauss-Jordan elimination on [m | I].
+// Returns an error for a non-square or singular (determinant ≈ 0)
+// matrix — ML code hitting a singular matrix (e.g. a rank-deficient
+// design matrix in least squares) usually wants to know that, not a
+// silently garbage answer.
+func Inverse(m Matrix) (Matrix, error) {
+	rows, cols, err := m.Dims()
+	if err != nil {
+		return nil, fmt.Errorf("matrix_inverse: %w", err)
+	}
+	if rows != cols {
+		return nil, fmt.Errorf("matrix_inverse: matrix must be square, got %dx%d", rows, cols)
+	}
+	n := rows
+	aug := make(Matrix, n)
+	for i := 0; i < n; i++ {
+		aug[i] = make([]float64, 2*n)
+		copy(aug[i], m[i])
+		aug[i][n+i] = 1
+	}
+	for col := 0; col < n; col++ {
+		pivotRow := col
+		best := math.Abs(aug[col][col])
+		for r := col + 1; r < n; r++ {
+			if v := math.Abs(aug[r][col]); v > best {
+				best, pivotRow = v, r
+			}
+		}
+		if best < 1e-12 {
+			return nil, fmt.Errorf("matrix_inverse: matrix is singular (not invertible)")
+		}
+		aug[col], aug[pivotRow] = aug[pivotRow], aug[col]
+		pivot := aug[col][col]
+		for c := 0; c < 2*n; c++ {
+			aug[col][c] /= pivot
+		}
+		for r := 0; r < n; r++ {
+			if r == col {
+				continue
+			}
+			factor := aug[r][col]
+			for c := 0; c < 2*n; c++ {
+				aug[r][c] -= factor * aug[col][c]
+			}
+		}
+	}
+	out := make(Matrix, n)
+	for i := 0; i < n; i++ {
+		out[i] = append([]float64(nil), aug[i][n:]...)
+	}
+	return out, nil
+}
+
+// Rank computes the matrix rank as the number of non-negligible
+// pivots found by Gaussian elimination (works for non-square
+// matrices too, unlike Determinant/Inverse).
+func Rank(m Matrix) (int, error) {
+	rows, cols, err := m.Dims()
+	if err != nil {
+		return 0, fmt.Errorf("matrix_rank: %w", err)
+	}
+	a := m.clone()
+	rank := 0
+	for col := 0; col < cols && rank < rows; col++ {
+		pivotRow := -1
+		best := 1e-9
+		for r := rank; r < rows; r++ {
+			if v := math.Abs(a[r][col]); v > best {
+				best, pivotRow = v, r
+			}
+		}
+		if pivotRow == -1 {
+			continue
+		}
+		a[rank], a[pivotRow] = a[pivotRow], a[rank]
+		for r := rank + 1; r < rows; r++ {
+			factor := a[r][col] / a[rank][col]
+			for c := col; c < cols; c++ {
+				a[r][c] -= factor * a[rank][c]
+			}
+		}
+		rank++
+	}
+	return rank, nil
+}
+
+// MatVecMul returns m×v for a column vector v.
+func MatVecMul(m Matrix, v []float64) ([]float64, error) {
+	rows, cols, err := m.Dims()
+	if err != nil {
+		return nil, fmt.Errorf("matrix-vector multiply: %w", err)
+	}
+	if cols != len(v) {
+		return nil, fmt.Errorf("matrix-vector multiply: matrix has %d columns but vector has %d entries", cols, len(v))
+	}
+	out := make([]float64, rows)
+	for i := 0; i < rows; i++ {
+		var sum float64
+		for j := 0; j < cols; j++ {
+			sum += m[i][j] * v[j]
+		}
+		out[i] = sum
+	}
+	return out, nil
+}
+
+// Dot returns the dot (inner) product of two equal-length vectors.
+func Dot(a, b []float64) (float64, error) {
+	if len(a) != len(b) {
+		return 0, fmt.Errorf("dot_product: vectors have different lengths (%d vs %d)", len(a), len(b))
+	}
+	var sum float64
+	for i := range a {
+		sum += a[i] * b[i]
+	}
+	return sum, nil
+}
+
+// Norm computes the Lp norm of a vector. p=2 is Euclidean (the
+// default any ML codebase means when it just says "norm"); p=1 is
+// Manhattan; p=math.Inf(1) is the max-absolute-value norm.
+func Norm(v []float64, p float64) float64 {
+	if math.IsInf(p, 1) {
+		var m float64
+		for _, x := range v {
+			if a := math.Abs(x); a > m {
+				m = a
+			}
+		}
+		return m
+	}
+	var sum float64
+	for _, x := range v {
+		sum += math.Pow(math.Abs(x), p)
+	}
+	return math.Pow(sum, 1/p)
+}
+
+// DominantEigen estimates the eigenvalue of largest absolute value
+// and its eigenvector via power iteration — the workhorse behind
+// things like PageRank and quick PCA-component estimates, and much
+// simpler/more numerically robust to implement from scratch than a
+// full QR eigendecomposition, which is the right trade-off for a
+// hand-rolled numerical toolkit. iterations<=0 defaults to 200, which
+// converges for any matrix with a well-separated dominant eigenvalue
+// to within float64 precision.
+func DominantEigen(m Matrix, iterations int) (eigenvalue float64, eigenvector []float64, err error) {
+	rows, cols, err := m.Dims()
+	if err != nil {
+		return 0, nil, fmt.Errorf("dominant_eigen: %w", err)
+	}
+	if rows != cols {
+		return 0, nil, fmt.Errorf("dominant_eigen: matrix must be square, got %dx%d", rows, cols)
+	}
+	if iterations <= 0 {
+		iterations = 200
+	}
+	n := rows
+	v := make([]float64, n)
+	for i := range v {
+		v[i] = 1 / math.Sqrt(float64(n)) // deterministic, reproducible start
+	}
+	var lambda float64
+	for iter := 0; iter < iterations; iter++ {
+		next, err := MatVecMul(m, v)
+		if err != nil {
+			return 0, nil, err
+		}
+		norm := Norm(next, 2)
+		if norm < 1e-15 {
+			return 0, nil, fmt.Errorf("dominant_eigen: did not converge (matrix may be nilpotent or zero)")
+		}
+		for i := range next {
+			next[i] /= norm
+		}
+		v = next
+		mv, _ := MatVecMul(m, v)
+		lambda, _ = Dot(v, mv)
+	}
+	return lambda, v, nil
+}

pkg/mathkit/mathkit_test.go

@@ -0,0 +1,364 @@
+package mathkit
+
+import (
+	"math"
+	"strings"
+	"testing"
+)
+
+func almostEqual(a, b, tol float64) bool { return math.Abs(a-b) <= tol }
+
+func mustExpr(t *testing.T, src string) *Expr {
+	t.Helper()
+	e, err := ParseExpr(src)
+	if err != nil {
+		t.Fatalf("ParseExpr(%q): %v", src, err)
+	}
+	return e
+}
+
+func TestExprEval(t *testing.T) {
+	cases := []struct {
+		src  string
+		vars map[string]float64
+		want float64
+	}{
+		{"2 + 3 * 4", nil, 14},
+		{"(2 + 3) * 4", nil, 20},
+		{"2 ^ 3 ^ 2", nil, 512}, // right-associative: 2^(3^2) = 2^9
+		{"-2 ^ 2", nil, -4},     // unary minus binds looser than ^: -(2^2)
+		{"x * x + 2*x + 1", map[string]float64{"x": 3}, 16},
+		{"sigmoid(0)", nil, 0.5},
+		{"relu(-5)", nil, 0},
+		{"relu(5)", nil, 5},
+		{"max(3, 7)", nil, 7},
+		{"sqrt(16)", nil, 4},
+		{"pi", nil, math.Pi},
+	}
+	for _, c := range cases {
+		e := mustExpr(t, c.src)
+		got, err := e.Eval(c.vars)
+		if err != nil {
+			t.Fatalf("Eval(%q): %v", c.src, err)
+		}
+		if !almostEqual(got, c.want, 1e-9) {
+			t.Errorf("Eval(%q) = %v, want %v", c.src, got, c.want)
+		}
+	}
+}
+
+func TestExprUnboundVariable(t *testing.T) {
+	e := mustExpr(t, "x + 1")
+	if _, err := e.Eval(nil); err == nil {
+		t.Fatal("expected an error for an unbound variable, got nil")
+	}
+}
+
+func TestMatMul(t *testing.T) {
+	a := Matrix{{1, 2}, {3, 4}}
+	b := Matrix{{5, 6}, {7, 8}}
+	got, err := MatMul(a, b)
+	if err != nil {
+		t.Fatal(err)
+	}
+	want := Matrix{{19, 22}, {43, 50}}
+	for i := range want {
+		for j := range want[i] {
+			if got[i][j] != want[i][j] {
+				t.Errorf("MatMul[%d][%d] = %v, want %v", i, j, got[i][j], want[i][j])
+			}
+		}
+	}
+}
+
+func TestMatMulShapeMismatch(t *testing.T) {
+	a := Matrix{{1, 2, 3}}
+	b := Matrix{{1, 2}}
+	if _, err := MatMul(a, b); err == nil {
+		t.Fatal("expected a shape-mismatch error")
+	}
+}
+
+func TestDeterminantAndInverse(t *testing.T) {
+	a := Matrix{{4, 7}, {2, 6}}
+	det, err := Determinant(a)
+	if err != nil {
+		t.Fatal(err)
+	}
+	if !almostEqual(det, 10, 1e-9) {
+		t.Errorf("det = %v, want 10", det)
+	}
+	inv, err := Inverse(a)
+	if err != nil {
+		t.Fatal(err)
+	}
+	// A * A^-1 should be the identity.
+	prod, err := MatMul(a, inv)
+	if err != nil {
+		t.Fatal(err)
+	}
+	for i := 0; i < 2; i++ {
+		for j := 0; j < 2; j++ {
+			want := 0.0
+			if i == j {
+				want = 1
+			}
+			if !almostEqual(prod[i][j], want, 1e-9) {
+				t.Errorf("A*A^-1[%d][%d] = %v, want %v", i, j, prod[i][j], want)
+			}
+		}
+	}
+}
+
+func TestInverseSingular(t *testing.T) {
+	a := Matrix{{1, 2}, {2, 4}} // row 2 = 2 * row 1
+	if _, err := Inverse(a); err == nil {
+		t.Fatal("expected a singular-matrix error")
+	}
+}
+
+func TestDominantEigen(t *testing.T) {
+	// [[2,0],[0,1]] has eigenvalues 2 and 1 — dominant is 2, eigenvector (1,0).
+	a := Matrix{{2, 0}, {0, 1}}
+	lambda, v, err := DominantEigen(a, 0)
+	if err != nil {
+		t.Fatal(err)
+	}
+	if !almostEqual(lambda, 2, 1e-6) {
+		t.Errorf("eigenvalue = %v, want 2", lambda)
+	}
+	if !almostEqual(math.Abs(v[0]), 1, 1e-6) || !almostEqual(v[1], 0, 1e-6) {
+		t.Errorf("eigenvector = %v, want ~(1,0)", v)
+	}
+}
+
+func TestDerivative(t *testing.T) {
+	// d/dx[x^2] at x=3 is 6.
+	e := mustExpr(t, "x^2")
+	d, err := Derivative(e, "x", map[string]float64{"x": 3}, 0)
+	if err != nil {
+		t.Fatal(err)
+	}
+	if !almostEqual(d, 6, 1e-4) {
+		t.Errorf("derivative = %v, want 6", d)
+	}
+}
+
+func TestGradient(t *testing.T) {
+	// f(x,y) = x^2 + y^2 ; gradient at (3,4) is (6,8).
+	e := mustExpr(t, "x^2 + y^2")
+	g, err := Gradient(e, []string{"x", "y"}, map[string]float64{"x": 3, "y": 4}, 0)
+	if err != nil {
+		t.Fatal(err)
+	}
+	if !almostEqual(g[0], 6, 1e-4) || !almostEqual(g[1], 8, 1e-4) {
+		t.Errorf("gradient = %v, want [6 8]", g)
+	}
+}
+
+func TestIntegrate(t *testing.T) {
+	// ∫0..1 x^2 dx = 1/3.
+	e := mustExpr(t, "x^2")
+	v, err := Integrate(e, "x", 0, 1, 100, nil)
+	if err != nil {
+		t.Fatal(err)
+	}
+	if !almostEqual(v, 1.0/3.0, 1e-9) {
+		t.Errorf("integral = %v, want 1/3", v)
+	}
+}
+
+func TestGradientDescentConvergesToMinimum(t *testing.T) {
+	// f(x) = (x-3)^2 has its minimum at x=3.
+	e := mustExpr(t, "(x - 3)^2")
+	res, err := GradientDescent(e, []string{"x"}, []float64{0}, 0.1, 200, 1e-10)
+	if err != nil {
+		t.Fatal(err)
+	}
+	final := res.Path[len(res.Path)-1][0]
+	if !almostEqual(final, 3, 1e-3) {
+		t.Errorf("gradient descent converged to x=%v, want ~3", final)
+	}
+}
+
+func TestNewtonRootFindsSquareRoot(t *testing.T) {
+	// Root of x^2 - 2 = 0 is sqrt(2).
+	e := mustExpr(t, "x^2 - 2")
+	res, err := NewtonRoot(e, "x", 1, 0, 0, nil)
+	if err != nil {
+		t.Fatal(err)
+	}
+	root := res.Path[len(res.Path)-1]
+	if !almostEqual(root, math.Sqrt2, 1e-6) {
+		t.Errorf("newton root = %v, want %v", root, math.Sqrt2)
+	}
+}
+
+func TestStatistics(t *testing.T) {
+	data := []float64{2, 4, 4, 4, 5, 5, 7, 9}
+	mean, _ := Mean(data)
+	if !almostEqual(mean, 5, 1e-9) {
+		t.Errorf("mean = %v, want 5", mean)
+	}
+	popVar, _ := Variance(data, true)
+	if !almostEqual(popVar, 4, 1e-9) {
+		t.Errorf("population variance = %v, want 4", popVar)
+	}
+	std, _ := StdDev(data, true)
+	if !almostEqual(std, 2, 1e-9) {
+		t.Errorf("population std_dev = %v, want 2", std)
+	}
+}
+
+func TestCorrelationPerfect(t *testing.T) {
+	x := []float64{1, 2, 3, 4, 5}
+	y := []float64{2, 4, 6, 8, 10}
+	c, err := Correlation(x, y)
+	if err != nil {
+		t.Fatal(err)
+	}
+	if !almostEqual(c, 1, 1e-9) {
+		t.Errorf("correlation = %v, want 1", c)
+	}
+}
+
+func TestNormalDistribution(t *testing.T) {
+	// Standard normal at x=0: pdf = 1/sqrt(2*pi), cdf = 0.5.
+	pdf, _ := NormalPDF(0, 0, 1)
+	if !almostEqual(pdf, 1/math.Sqrt(2*math.Pi), 1e-9) {
+		t.Errorf("pdf(0) = %v", pdf)
+	}
+	cdf, _ := NormalCDF(0, 0, 1)
+	if !almostEqual(cdf, 0.5, 1e-9) {
+		t.Errorf("cdf(0) = %v, want 0.5", cdf)
+	}
+}
+
+func TestSigmoidNumericallyStable(t *testing.T) {
+	if v := Sigmoid(-1000); v != 0 {
+		t.Errorf("sigmoid(-1000) = %v, want 0", v)
+	}
+	if v := Sigmoid(1000); !almostEqual(v, 1, 1e-12) {
+		t.Errorf("sigmoid(1000) = %v, want ~1", v)
+	}
+	if v := Sigmoid(0); v != 0.5 {
+		t.Errorf("sigmoid(0) = %v, want 0.5", v)
+	}
+}
+
+func TestSoftmaxSumsToOne(t *testing.T) {
+	out := Softmax([]float64{1, 2, 3})
+	var sum float64
+	for _, v := range out {
+		sum += v
+	}
+	if !almostEqual(sum, 1, 1e-9) {
+		t.Errorf("softmax sums to %v, want 1", sum)
+	}
+	// Largest logit should get the largest probability.
+	if out[2] <= out[1] || out[1] <= out[0] {
+		t.Errorf("softmax ordering wrong: %v", out)
+	}
+}
+
+func TestMSELoss(t *testing.T) {
+	v, err := MSELoss([]float64{1, 2, 3}, []float64{1, 2, 4})
+	if err != nil {
+		t.Fatal(err)
+	}
+	// errors: 0, 0, -1 -> squared: 0,0,1 -> mean = 1/3
+	if !almostEqual(v, 1.0/3.0, 1e-9) {
+		t.Errorf("mse = %v, want 1/3", v)
+	}
+}
+
+func TestLinearRegressionPerfectFit(t *testing.T) {
+	// y = 2 + 3x, exactly.
+	X := Matrix{{1}, {2}, {3}, {4}, {5}}
+	y := []float64{5, 8, 11, 14, 17}
+	res, err := LinearRegression(X, y)
+	if err != nil {
+		t.Fatal(err)
+	}
+	if !almostEqual(res.Coefficients[0], 2, 1e-6) {
+		t.Errorf("intercept = %v, want 2", res.Coefficients[0])
+	}
+	if !almostEqual(res.Coefficients[1], 3, 1e-6) {
+		t.Errorf("slope = %v, want 3", res.Coefficients[1])
+	}
+	if !almostEqual(res.RSquared, 1, 1e-6) {
+		t.Errorf("R^2 = %v, want 1", res.RSquared)
+	}
+}
+
+// --- Solve dispatcher: the JSON-boundary layer the agent tool calls ---
+
+func TestSolveEvaluate(t *testing.T) {
+	out, err := Solve("evaluate", map[string]any{
+		"expression": "x^2 + y^2",
+		"variables":  map[string]any{"x": 3.0, "y": 4.0},
+	})
+	if err != nil {
+		t.Fatal(err)
+	}
+	if !strings.Contains(out, "25") {
+		t.Errorf("Solve(evaluate) = %q, want it to contain 25", out)
+	}
+}
+
+func TestSolveMatrixMul(t *testing.T) {
+	out, err := Solve("matrix_mul", map[string]any{
+		"a": []any{[]any{1.0, 2.0}, []any{3.0, 4.0}},
+		"b": []any{[]any{5.0, 6.0}, []any{7.0, 8.0}},
+	})
+	if err != nil {
+		t.Fatal(err)
+	}
+	if !strings.Contains(out, "19") || !strings.Contains(out, "50") {
+		t.Errorf("Solve(matrix_mul) = %q, missing expected entries", out)
+	}
+}
+
+func TestSolveUnknownOperation(t *testing.T) {
+	if _, err := Solve("not_a_real_operation", nil); err == nil {
+		t.Fatal("expected an error for an unknown operation")
+	}
+}
+
+func TestSolveGradientDescent(t *testing.T) {
+	out, err := Solve("gradient_descent", map[string]any{
+		"expression":    "(x - 3)^2",
+		"variables":     []any{"x"},
+		"start":         []any{0.0},
+		"learning_rate": 0.1,
+		"iterations":    200.0,
+	})
+	if err != nil {
+		t.Fatal(err)
+	}
+	if !strings.Contains(out, "converged=true") {
+		t.Errorf("Solve(gradient_descent) = %q, expected convergence", out)
+	}
+}
+
+func TestSolveMissingField(t *testing.T) {
+	_, err := Solve("mean", map[string]any{})
+	if err == nil {
+		t.Fatal("expected a missing-field error")
+	}
+	if !strings.Contains(err.Error(), "data") {
+		t.Errorf("error %q should mention the missing field", err.Error())
+	}
+}
+
+func TestOperationsListIsSorted(t *testing.T) {
+	ops := Operations()
+	if len(ops) == 0 {
+		t.Fatal("Operations() returned nothing")
+	}
+	for i := 1; i < len(ops); i++ {
+		if ops[i-1] > ops[i] {
+			t.Fatalf("Operations() not sorted: %q before %q", ops[i-1], ops[i])
+		}
+	}
+}

pkg/mathkit/ml.go

@@ -0,0 +1,191 @@
+package mathkit
+
+import (
+	"fmt"
+	"math"
+)
+
+// Sigmoid is the logistic function 1 / (1 + e^-x), computed in the
+// numerically stable branch-on-sign form so it doesn't overflow
+// math.Exp for very negative x (a plain 1/(1+exp(-x)) blows up
+// exp(-x) to +Inf for x below about -745 and silently returns 0
+// instead of propagating the overflow visibly).
+func Sigmoid(x float64) float64 {
+	if x >= 0 {
+		z := math.Exp(-x)
+		return 1 / (1 + z)
+	}
+	z := math.Exp(x)
+	return z / (1 + z)
+}
+
+// ReLU is max(0, x).
+func ReLU(x float64) float64 {
+	if x > 0 {
+		return x
+	}
+	return 0
+}
+
+// Softmax turns a vector of raw scores ("logits") into a probability
+// distribution: non-negative entries that sum to 1. Subtracts the max
+// before exponentiating (the standard numerically-stable formulation)
+// so large logits don't overflow math.Exp.
+func Softmax(x []float64) []float64 {
+	if len(x) == 0 {
+		return nil
+	}
+	max := x[0]
+	for _, v := range x[1:] {
+		if v > max {
+			max = v
+		}
+	}
+	out := make([]float64, len(x))
+	var sum float64
+	for i, v := range x {
+		e := math.Exp(v - max)
+		out[i] = e
+		sum += e
+	}
+	for i := range out {
+		out[i] /= sum
+	}
+	return out
+}
+
+// MSELoss returns the mean squared error between predictions and
+// targets: mean((pred - target)^2), the default regression loss.
+func MSELoss(predictions, targets []float64) (float64, error) {
+	if len(predictions) != len(targets) {
+		return 0, fmt.Errorf("mse_loss: predictions and targets have different lengths (%d vs %d)", len(predictions), len(targets))
+	}
+	if len(predictions) == 0 {
+		return 0, fmt.Errorf("mse_loss: empty input")
+	}
+	var sum float64
+	for i := range predictions {
+		d := predictions[i] - targets[i]
+		sum += d * d
+	}
+	return sum / float64(len(predictions)), nil
+}
+
+// crossEntropyEps clamps predicted probabilities away from exactly 0
+// or 1 before taking a log — ln(0) is -Inf, which would make a single
+// confidently-wrong prediction blow up the whole averaged loss to
+// +Inf instead of just contributing a large-but-finite penalty. Every
+// production cross-entropy implementation (e.g. scikit-learn's
+// log_loss) clips for exactly this reason.
+const crossEntropyEps = 1e-12
+
+// CrossEntropyLoss returns the mean binary cross-entropy between
+// predicted probabilities and 0/1 (or soft, in [0,1]) targets:
+// mean(-(y*ln(p) + (1-y)*ln(1-p))), the default classification loss.
+func CrossEntropyLoss(predictions, targets []float64) (float64, error) {
+	if len(predictions) != len(targets) {
+		return 0, fmt.Errorf("cross_entropy_loss: predictions and targets have different lengths (%d vs %d)", len(predictions), len(targets))
+	}
+	if len(predictions) == 0 {
+		return 0, fmt.Errorf("cross_entropy_loss: empty input")
+	}
+	var sum float64
+	for i, p := range predictions {
+		if p < 0 || p > 1 {
+			return 0, fmt.Errorf("cross_entropy_loss: predictions must be probabilities in [0, 1], got %g at index %d", p, i)
+		}
+		p = math.Min(math.Max(p, crossEntropyEps), 1-crossEntropyEps)
+		y := targets[i]
+		sum += -(y*math.Log(p) + (1-y)*math.Log(1-p))
+	}
+	return sum / float64(len(predictions)), nil
+}
+
+// LinearRegressionResult is a fitted ordinary-least-squares linear
+// model.
+type LinearRegressionResult struct {
+	// Coefficients[0] is the intercept, Coefficients[1:] line up with
+	// the feature columns of X in order.
+	Coefficients []float64
+	Predictions  []float64
+	Residuals    []float64
+	RSquared     float64
+}
+
+// LinearRegression fits y ≈ X·β + intercept by ordinary least squares
+// via the normal equations β = (XᵀX)^-1 Xᵀy, with an intercept column
+// of 1s prepended automatically (the caller supplies only the real
+// feature columns in X). The normal-equation approach is the direct,
+// closed-form solution to minimizing squared error — the textbook
+// ML-101 regression fit, as opposed to reaching it iteratively via
+// GradientDescent on the MSE loss.
+func LinearRegression(X Matrix, y []float64) (*LinearRegressionResult, error) {
+	rows, cols, err := X.Dims()
+	if err != nil {
+		return nil, fmt.Errorf("linear_regression: %w", err)
+	}
+	if rows != len(y) {
+		return nil, fmt.Errorf("linear_regression: X has %d rows but y has %d entries", rows, len(y))
+	}
+	if rows <= cols {
+		return nil, fmt.Errorf("linear_regression: need more observations (%d) than features+1 (%d) to fit uniquely", rows, cols+1)
+	}
+
+	// Prepend the intercept column of 1s.
+	design := make(Matrix, rows)
+	for i := 0; i < rows; i++ {
+		row := make([]float64, cols+1)
+		row[0] = 1
+		copy(row[1:], X[i])
+		design[i] = row
+	}
+
+	designT, err := Transpose(design)
+	if err != nil {
+		return nil, fmt.Errorf("linear_regression: %w", err)
+	}
+	xtx, err := MatMul(designT, design)
+	if err != nil {
+		return nil, fmt.Errorf("linear_regression: %w", err)
+	}
+	xtxInv, err := Inverse(xtx)
+	if err != nil {
+		return nil, fmt.Errorf("linear_regression: features are collinear (XᵀX is singular): %w", err)
+	}
+	xty, err := MatVecMul(designT, y)
+	if err != nil {
+		return nil, fmt.Errorf("linear_regression: %w", err)
+	}
+	beta, err := MatVecMul(xtxInv, xty)
+	if err != nil {
+		return nil, fmt.Errorf("linear_regression: %w", err)
+	}
+
+	preds, err := MatVecMul(design, beta)
+	if err != nil {
+		return nil, fmt.Errorf("linear_regression: %w", err)
+	}
+	residuals := make([]float64, rows)
+	var ssRes float64
+	for i := range preds {
+		residuals[i] = y[i] - preds[i]
+		ssRes += residuals[i] * residuals[i]
+	}
+	meanY, _ := Mean(y)
+	var ssTot float64
+	for _, v := range y {
+		d := v - meanY
+		ssTot += d * d
+	}
+	rSquared := 1.0
+	if ssTot > 0 {
+		rSquared = 1 - ssRes/ssTot
+	}
+
+	return &LinearRegressionResult{
+		Coefficients: beta,
+		Predictions:  preds,
+		Residuals:    residuals,
+		RSquared:     rSquared,
+	}, nil
+}

pkg/mathkit/optimize.go

@@ -0,0 +1,135 @@
+package mathkit
+
+import (
+	"fmt"
+	"math"
+)
+
+// GradientDescentResult is the trace of one gradient-descent run.
+type GradientDescentResult struct {
+	// Path holds the point visited at the start of every iteration,
+	// Path[0] being the caller's starting point and the last entry
+	// being the final point — so len(Path) == Iterations+1. Kept in
+	// full (rather than just the final point) because seeing the
+	// trajectory is usually the point of running this: it's how you
+	// tell "converging nicely" apart from "diverging" or
+	// "oscillating because the learning rate is too high".
+	Path [][]float64
+	// Values holds f(Path[i]) for each point in Path, same length.
+	Values []float64
+	// Converged is true if the gradient norm dropped below the
+	// requested tolerance before Iterations ran out.
+	Converged bool
+}
+
+// GradientDescent runs plain (batch) gradient descent on expr, over
+// the variables listed (in order — that order defines the vector
+// space the descent happens in), starting at `start`, for at most
+// `iterations` steps of size `learningRate`. tol<=0 disables the
+// early-stop check (it always runs the full iteration count).
+//
+// This is deliberately vanilla gradient descent — no momentum, no
+// Adam-style adaptive rates — because the point of this tool is to
+// let a user see and reason about the bare update rule
+// x := x - lr * ∇f(x) that every fancier optimizer builds on, not to
+// re-implement a training framework.
+func GradientDescent(expr *Expr, variables []string, start []float64, learningRate float64, iterations int, tol float64) (*GradientDescentResult, error) {
+	if len(variables) != len(start) {
+		return nil, fmt.Errorf("gradient_descent: %d variable name(s) but start point has %d coordinate(s)", len(variables), len(start))
+	}
+	if learningRate <= 0 {
+		return nil, fmt.Errorf("gradient_descent: learning_rate must be positive")
+	}
+	if iterations <= 0 {
+		iterations = 100
+	}
+	point := append([]float64(nil), start...)
+	res := &GradientDescentResult{}
+	addPoint := func(p []float64) error {
+		vars := make(map[string]float64, len(variables))
+		for i, name := range variables {
+			vars[name] = p[i]
+		}
+		v, err := expr.Eval(vars)
+		if err != nil {
+			return err
+		}
+		res.Path = append(res.Path, append([]float64(nil), p...))
+		res.Values = append(res.Values, v)
+		return nil
+	}
+	if err := addPoint(point); err != nil {
+		return nil, fmt.Errorf("gradient_descent: evaluating start point: %w", err)
+	}
+	vars := make(map[string]float64, len(variables))
+	for iter := 0; iter < iterations; iter++ {
+		for i, name := range variables {
+			vars[name] = point[i]
+		}
+		grad, err := Gradient(expr, variables, vars, 0)
+		if err != nil {
+			return nil, fmt.Errorf("gradient_descent: iteration %d: %w", iter, err)
+		}
+		if tol > 0 && Norm(grad, 2) < tol {
+			res.Converged = true
+			break
+		}
+		for i := range point {
+			point[i] -= learningRate * grad[i]
+		}
+		if err := addPoint(point); err != nil {
+			return nil, fmt.Errorf("gradient_descent: iteration %d: %w", iter, err)
+		}
+		if math.IsNaN(point[0]) || math.IsInf(point[0], 0) {
+			return nil, fmt.Errorf("gradient_descent: diverged (values went to NaN/Inf) — try a smaller learning_rate")
+		}
+	}
+	return res, nil
+}
+
+// NewtonRootResult is the trace of one Newton's-method root-finding
+// run on a single-variable expression.
+type NewtonRootResult struct {
+	Path      []float64 // x_0, x_1, ..., x_final
+	Converged bool
+}
+
+// NewtonRoot finds a root of expr(variable) = 0 via Newton's method:
+// x_{n+1} = x_n - f(x_n)/f'(x_n), with f' estimated by central
+// difference (see Derivative). Any other free variables in expr must
+// be supplied in `fixed`. tol<=0 defaults to 1e-10.
+func NewtonRoot(expr *Expr, variable string, start float64, iterations int, tol float64, fixed map[string]float64) (*NewtonRootResult, error) {
+	if iterations <= 0 {
+		iterations = 100
+	}
+	if tol <= 0 {
+		tol = 1e-10
+	}
+	x := start
+	res := &NewtonRootResult{Path: []float64{x}}
+	vars := cloneVars(fixed)
+	for iter := 0; iter < iterations; iter++ {
+		vars[variable] = x
+		fx, err := expr.Eval(vars)
+		if err != nil {
+			return nil, fmt.Errorf("newton_root: iteration %d: %w", iter, err)
+		}
+		if math.Abs(fx) < tol {
+			res.Converged = true
+			break
+		}
+		dfx, err := Derivative(expr, variable, vars, 0)
+		if err != nil {
+			return nil, fmt.Errorf("newton_root: iteration %d: %w", iter, err)
+		}
+		if math.Abs(dfx) < 1e-14 {
+			return nil, fmt.Errorf("newton_root: derivative vanished near x=%g (flat spot) — try a different start", x)
+		}
+		x = x - fx/dfx
+		res.Path = append(res.Path, x)
+		if math.IsNaN(x) || math.IsInf(x, 0) {
+			return nil, fmt.Errorf("newton_root: diverged (x went to NaN/Inf) — try a different start")
+		}
+	}
+	return res, nil
+}

pkg/mathkit/solve.go

@@ -0,0 +1,835 @@
+package mathkit
+
+import (
+	"fmt"
+	"math"
+	"sort"
+	"strings"
+)
+
+// Solve dispatches a math problem by string operation name against a
+// loosely-typed params map — exactly the shape a JSON object (an LLM
+// tool-call argument blob, decoded with encoding/json into
+// map[string]any) takes once decoded. It exists as one layer above
+// the typed functions in linalg.go/calculus.go/statistics.go/
+// optimize.go/ml.go specifically so a caller at the JSON boundary
+// (services/agent/internal/infrastructure/toolexec/mlmath.go) has one
+// function to call and one error type to handle, instead of a type
+// switch per operation.
+//
+// Returns a human-readable report (not raw JSON or a bare float) on
+// success, because the caller here is ultimately an LLM deciding what
+// to do next — a labeled, explained answer is directly usable in its
+// response, where a bare number would have to be re-explained.
+func Solve(operation string, params map[string]any) (string, error) {
+	op := strings.ToLower(strings.TrimSpace(operation))
+	switch op {
+	// --- general ---
+	case "evaluate":
+		return solveEvaluate(params)
+
+	// --- linear algebra ---
+	case "matrix_add", "matrix_sub", "matrix_mul":
+		return solveMatrixBinOp(op, params)
+	case "matrix_transpose":
+		return solveMatrixTranspose(params)
+	case "matrix_determinant":
+		return solveMatrixDeterminant(params)
+	case "matrix_inverse":
+		return solveMatrixInverse(params)
+	case "matrix_rank":
+		return solveMatrixRank(params)
+	case "matrix_trace":
+		return solveMatrixTrace(params)
+	case "dot_product":
+		return solveDotProduct(params)
+	case "vector_norm":
+		return solveVectorNorm(params)
+	case "dominant_eigen":
+		return solveDominantEigen(params)
+
+	// --- calculus ---
+	case "derivative":
+		return solveDerivative(params)
+	case "gradient":
+		return solveGradient(params)
+	case "integral":
+		return solveIntegral(params)
+
+	// --- optimization ---
+	case "gradient_descent":
+		return solveGradientDescent(params)
+	case "newton_root":
+		return solveNewtonRoot(params)
+
+	// --- statistics & probability ---
+	case "mean":
+		return solveReduce1(params, "mean", func(d []float64) (float64, error) { return Mean(d) })
+	case "variance":
+		return solveReduce1Pop(params, "variance", Variance)
+	case "std_dev":
+		return solveReduce1Pop(params, "std_dev", StdDev)
+	case "median":
+		return solveReduce1(params, "median", Median)
+	case "covariance":
+		return solveReduce2Pop(params, "covariance", Covariance)
+	case "correlation":
+		return solveReduce2(params, "correlation", Correlation)
+	case "z_score_normalize":
+		return solveNormalize(params, "z_score_normalize", ZScoreNormalize)
+	case "min_max_normalize":
+		return solveNormalize(params, "min_max_normalize", MinMaxNormalize)
+	case "normal_pdf":
+		return solveNormalDist(params, "normal_pdf", NormalPDF)
+	case "normal_cdf":
+		return solveNormalDist(params, "normal_cdf", NormalCDF)
+	case "binomial_pmf":
+		return solveBinomialPMF(params)
+
+	// --- ML primitives ---
+	case "sigmoid":
+		return solveElementwise(params, "sigmoid", Sigmoid)
+	case "relu":
+		return solveElementwise(params, "relu", ReLU)
+	case "tanh":
+		return solveElementwise(params, "tanh", math.Tanh)
+	case "softmax":
+		return solveSoftmax(params)
+	case "mse_loss":
+		return solveLoss(params, "mse_loss", MSELoss)
+	case "cross_entropy_loss":
+		return solveLoss(params, "cross_entropy_loss", CrossEntropyLoss)
+	case "linear_regression":
+		return solveLinearRegression(params)
+
+	default:
+		return "", fmt.Errorf("unknown operation %q — see the tool description for the full list of supported operations", operation)
+	}
+}
+
+// ---------- param extraction helpers ----------
+//
+// Every helper below reports a specific, actionable error (which
+// field was wrong and why) rather than a generic type-assertion
+// panic recovery message, because these errors are fed straight back
+// to the calling model as the tool result — see toolexec's rpcFail /
+// IsError convention — and "field X must be a number, got a string"
+// is something the model can act on; "interface conversion failed"
+// is not.
+
+func reqFloat(params map[string]any, key string) (float64, error) {
+	v, ok := params[key]
+	if !ok {
+		return 0, fmt.Errorf("missing required field %q", key)
+	}
+	f, ok := toFloat(v)
+	if !ok {
+		return 0, fmt.Errorf("field %q must be a number", key)
+	}
+	return f, nil
+}
+
+func optFloat(params map[string]any, key string, fallback float64) float64 {
+	v, ok := params[key]
+	if !ok {
+		return fallback
+	}
+	if f, ok := toFloat(v); ok {
+		return f
+	}
+	return fallback
+}
+
+func optInt(params map[string]any, key string, fallback int) int {
+	v, ok := params[key]
+	if !ok {
+		return fallback
+	}
+	if f, ok := toFloat(v); ok {
+		return int(f)
+	}
+	return fallback
+}
+
+func optBool(params map[string]any, key string, fallback bool) bool {
+	v, ok := params[key]
+	if !ok {
+		return fallback
+	}
+	if b, ok := v.(bool); ok {
+		return b
+	}
+	return fallback
+}
+
+func toFloat(v any) (float64, bool) {
+	switch n := v.(type) {
+	case float64:
+		return n, true
+	case int:
+		return float64(n), true
+	case int64:
+		return float64(n), true
+	}
+	return 0, false
+}
+
+func reqString(params map[string]any, key string) (string, error) {
+	v, ok := params[key]
+	if !ok {
+		return "", fmt.Errorf("missing required field %q", key)
+	}
+	s, ok := v.(string)
+	if !ok || strings.TrimSpace(s) == "" {
+		return "", fmt.Errorf("field %q must be a non-empty string", key)
+	}
+	return s, nil
+}
+
+func optString(params map[string]any, key, fallback string) string {
+	v, ok := params[key]
+	if !ok {
+		return fallback
+	}
+	s, ok := v.(string)
+	if !ok {
+		return fallback
+	}
+	return s
+}
+
+func reqFloatSlice(params map[string]any, key string) ([]float64, error) {
+	v, ok := params[key]
+	if !ok {
+		return nil, fmt.Errorf("missing required field %q (expected an array of numbers)", key)
+	}
+	raw, ok := v.([]any)
+	if !ok {
+		return nil, fmt.Errorf("field %q must be an array of numbers", key)
+	}
+	out := make([]float64, len(raw))
+	for i, x := range raw {
+		f, ok := toFloat(x)
+		if !ok {
+			return nil, fmt.Errorf("field %q: element %d is not a number", key, i)
+		}
+		out[i] = f
+	}
+	return out, nil
+}
+
+func reqStringSlice(params map[string]any, key string) ([]string, error) {
+	v, ok := params[key]
+	if !ok {
+		return nil, fmt.Errorf("missing required field %q (expected an array of strings)", key)
+	}
+	raw, ok := v.([]any)
+	if !ok {
+		return nil, fmt.Errorf("field %q must be an array of strings", key)
+	}
+	out := make([]string, len(raw))
+	for i, x := range raw {
+		s, ok := x.(string)
+		if !ok {
+			return nil, fmt.Errorf("field %q: element %d is not a string", key, i)
+		}
+		out[i] = s
+	}
+	return out, nil
+}
+
+func reqMatrix(params map[string]any, key string) (Matrix, error) {
+	v, ok := params[key]
+	if !ok {
+		return nil, fmt.Errorf("missing required field %q (expected a 2D array of numbers)", key)
+	}
+	rows, ok := v.([]any)
+	if !ok {
+		return nil, fmt.Errorf("field %q must be a 2D array of numbers, e.g. [[1,2],[3,4]]", key)
+	}
+	m := make(Matrix, len(rows))
+	for i, r := range rows {
+		rowRaw, ok := r.([]any)
+		if !ok {
+			return nil, fmt.Errorf("field %q: row %d is not an array", key, i)
+		}
+		row := make([]float64, len(rowRaw))
+		for j, x := range rowRaw {
+			f, ok := toFloat(x)
+			if !ok {
+				return nil, fmt.Errorf("field %q: row %d, column %d is not a number", key, i, j)
+			}
+			row[j] = f
+		}
+		m[i] = row
+	}
+	return m, nil
+}
+
+// varMap builds the variable bindings an expression needs from an
+// "at" (or "variables") object field — {"x": 1, "y": 2} — falling
+// back to an empty map if the field is absent (a single-variable
+// expression evaluated via the "variable"/"at" pair still works).
+func varMap(params map[string]any, key string) (map[string]float64, error) {
+	v, ok := params[key]
+	if !ok {
+		return map[string]float64{}, nil
+	}
+	raw, ok := v.(map[string]any)
+	if !ok {
+		return nil, fmt.Errorf("field %q must be an object mapping variable names to numbers, e.g. {\"x\": 1.5}", key)
+	}
+	out := make(map[string]float64, len(raw))
+	for k, x := range raw {
+		f, ok := toFloat(x)
+		if !ok {
+			return nil, fmt.Errorf("field %q: value for %q is not a number", key, k)
+		}
+		out[k] = f
+	}
+	return out, nil
+}
+
+func parseExprField(params map[string]any) (*Expr, error) {
+	src, err := reqString(params, "expression")
+	if err != nil {
+		return nil, err
+	}
+	expr, err := ParseExpr(src)
+	if err != nil {
+		return nil, err
+	}
+	return expr, nil
+}
+
+func fmtFloat(f float64) string {
+	if math.IsNaN(f) {
+		return "NaN"
+	}
+	if math.IsInf(f, 1) {
+		return "+Inf"
+	}
+	if math.IsInf(f, -1) {
+		return "-Inf"
+	}
+	return strings.TrimRight(strings.TrimRight(fmt.Sprintf("%.10f", f), "0"), ".")
+}
+
+func fmtFloats(fs []float64) string {
+	parts := make([]string, len(fs))
+	for i, f := range fs {
+		parts[i] = fmtFloat(f)
+	}
+	return "[" + strings.Join(parts, ", ") + "]"
+}
+
+func fmtMatrix(m Matrix) string {
+	var sb strings.Builder
+	for _, row := range m {
+		sb.WriteString(fmtFloats(row))
+		sb.WriteByte('\n')
+	}
+	return strings.TrimRight(sb.String(), "\n")
+}
+
+// ---------- operation implementations ----------
+
+func solveEvaluate(params map[string]any) (string, error) {
+	expr, err := parseExprField(params)
+	if err != nil {
+		return "", err
+	}
+	vars, err := varMap(params, "variables")
+	if err != nil {
+		return "", err
+	}
+	result, err := expr.Eval(vars)
+	if err != nil {
+		return "", err
+	}
+	return fmt.Sprintf("%s = %s", expr.String(), fmtFloat(result)), nil
+}
+
+func solveMatrixBinOp(op string, params map[string]any) (string, error) {
+	a, err := reqMatrix(params, "a")
+	if err != nil {
+		return "", err
+	}
+	b, err := reqMatrix(params, "b")
+	if err != nil {
+		return "", err
+	}
+	var out Matrix
+	switch op {
+	case "matrix_add":
+		out, err = MatAdd(a, b)
+	case "matrix_sub":
+		out, err = MatSub(a, b)
+	case "matrix_mul":
+		out, err = MatMul(a, b)
+	}
+	if err != nil {
+		return "", err
+	}
+	return fmt.Sprintf("%s result:\n%s", op, fmtMatrix(out)), nil
+}
+
+func solveMatrixTranspose(params map[string]any) (string, error) {
+	a, err := reqMatrix(params, "a")
+	if err != nil {
+		return "", err
+	}
+	out, err := Transpose(a)
+	if err != nil {
+		return "", err
+	}
+	return fmt.Sprintf("transpose result:\n%s", fmtMatrix(out)), nil
+}
+
+func solveMatrixDeterminant(params map[string]any) (string, error) {
+	a, err := reqMatrix(params, "a")
+	if err != nil {
+		return "", err
+	}
+	det, err := Determinant(a)
+	if err != nil {
+		return "", err
+	}
+	return fmt.Sprintf("det(A) = %s", fmtFloat(det)), nil
+}
+
+func solveMatrixInverse(params map[string]any) (string, error) {
+	a, err := reqMatrix(params, "a")
+	if err != nil {
+		return "", err
+	}
+	inv, err := Inverse(a)
+	if err != nil {
+		return "", err
+	}
+	return fmt.Sprintf("A^-1:\n%s", fmtMatrix(inv)), nil
+}
+
+func solveMatrixRank(params map[string]any) (string, error) {
+	a, err := reqMatrix(params, "a")
+	if err != nil {
+		return "", err
+	}
+	r, err := Rank(a)
+	if err != nil {
+		return "", err
+	}
+	return fmt.Sprintf("rank(A) = %d", r), nil
+}
+
+func solveMatrixTrace(params map[string]any) (string, error) {
+	a, err := reqMatrix(params, "a")
+	if err != nil {
+		return "", err
+	}
+	tr, err := Trace(a)
+	if err != nil {
+		return "", err
+	}
+	return fmt.Sprintf("trace(A) = %s", fmtFloat(tr)), nil
+}
+
+func solveDotProduct(params map[string]any) (string, error) {
+	a, err := reqFloatSlice(params, "a")
+	if err != nil {
+		return "", err
+	}
+	b, err := reqFloatSlice(params, "b")
+	if err != nil {
+		return "", err
+	}
+	d, err := Dot(a, b)
+	if err != nil {
+		return "", err
+	}
+	return fmt.Sprintf("a · b = %s", fmtFloat(d)), nil
+}
+
+func solveVectorNorm(params map[string]any) (string, error) {
+	a, err := reqFloatSlice(params, "a")
+	if err != nil {
+		return "", err
+	}
+	p := optFloat(params, "p", 2)
+	n := Norm(a, p)
+	label := fmt.Sprintf("L%g", p)
+	if math.IsInf(p, 1) {
+		label = "L∞"
+	}
+	return fmt.Sprintf("%s norm = %s", label, fmtFloat(n)), nil
+}
+
+func solveDominantEigen(params map[string]any) (string, error) {
+	a, err := reqMatrix(params, "a")
+	if err != nil {
+		return "", err
+	}
+	iters := optInt(params, "iterations", 0)
+	lambda, v, err := DominantEigen(a, iters)
+	if err != nil {
+		return "", err
+	}
+	return fmt.Sprintf("dominant eigenvalue ≈ %s\neigenvector ≈ %s", fmtFloat(lambda), fmtFloats(v)), nil
+}
+
+func solveDerivative(params map[string]any) (string, error) {
+	expr, err := parseExprField(params)
+	if err != nil {
+		return "", err
+	}
+	variable, err := reqString(params, "variable")
+	if err != nil {
+		return "", err
+	}
+	at, err := varMap(params, "at")
+	if err != nil {
+		return "", err
+	}
+	if _, ok := at[variable]; !ok {
+		at[variable] = optFloat(params, "value", 0)
+	}
+	h := optFloat(params, "h", 0)
+	d, err := Derivative(expr, variable, at, h)
+	if err != nil {
+		return "", err
+	}
+	return fmt.Sprintf("d/d%s[%s] at %s=%s ≈ %s", variable, expr.String(), variable, fmtFloat(at[variable]), fmtFloat(d)), nil
+}
+
+func solveGradient(params map[string]any) (string, error) {
+	expr, err := parseExprField(params)
+	if err != nil {
+		return "", err
+	}
+	variables, err := reqStringSlice(params, "variables")
+	if err != nil {
+		return "", err
+	}
+	at, err := varMap(params, "at")
+	if err != nil {
+		return "", err
+	}
+	h := optFloat(params, "h", 0)
+	grad, err := Gradient(expr, variables, at, h)
+	if err != nil {
+		return "", err
+	}
+	return fmt.Sprintf("∇[%s] with respect to %s ≈ %s", expr.String(), strings.Join(variables, ","), fmtFloats(grad)), nil
+}
+
+func solveIntegral(params map[string]any) (string, error) {
+	expr, err := parseExprField(params)
+	if err != nil {
+		return "", err
+	}
+	variable, err := reqString(params, "variable")
+	if err != nil {
+		return "", err
+	}
+	from, err := reqFloat(params, "from")
+	if err != nil {
+		return "", err
+	}
+	to, err := reqFloat(params, "to")
+	if err != nil {
+		return "", err
+	}
+	steps := optInt(params, "steps", 0)
+	fixed, err := varMap(params, "variables")
+	if err != nil {
+		return "", err
+	}
+	result, err := Integrate(expr, variable, from, to, steps, fixed)
+	if err != nil {
+		return "", err
+	}
+	return fmt.Sprintf("∫[%s to %s] %s d%s ≈ %s", fmtFloat(from), fmtFloat(to), expr.String(), variable, fmtFloat(result)), nil
+}
+
+func solveGradientDescent(params map[string]any) (string, error) {
+	expr, err := parseExprField(params)
+	if err != nil {
+		return "", err
+	}
+	variables, err := reqStringSlice(params, "variables")
+	if err != nil {
+		return "", err
+	}
+	start, err := reqFloatSlice(params, "start")
+	if err != nil {
+		return "", err
+	}
+	lr, err := reqFloat(params, "learning_rate")
+	if err != nil {
+		return "", err
+	}
+	iterations := optInt(params, "iterations", 100)
+	tol := optFloat(params, "tolerance", 1e-8)
+	res, err := GradientDescent(expr, variables, start, lr, iterations, tol)
+	if err != nil {
+		return "", err
+	}
+	last := len(res.Path) - 1
+	var sb strings.Builder
+	fmt.Fprintf(&sb, "gradient descent on f(%s) = %s\n", strings.Join(variables, ","), expr.String())
+	fmt.Fprintf(&sb, "start: %s -> f = %s\n", fmtFloats(res.Path[0]), fmtFloat(res.Values[0]))
+	fmt.Fprintf(&sb, "ran %d step(s), converged=%v\n", last, res.Converged)
+	fmt.Fprintf(&sb, "final: %s -> f = %s\n", fmtFloats(res.Path[last]), fmtFloat(res.Values[last]))
+	// Show a short trajectory sample so the model can see the descent
+	// shape without a huge wall of numbers for large iteration counts.
+	sampleN := 5
+	if last+1 <= sampleN*2 {
+		for i := 0; i <= last; i++ {
+			fmt.Fprintf(&sb, "  step %d: x=%s f=%s\n", i, fmtFloats(res.Path[i]), fmtFloat(res.Values[i]))
+		}
+	} else {
+		for i := 0; i < sampleN; i++ {
+			fmt.Fprintf(&sb, "  step %d: x=%s f=%s\n", i, fmtFloats(res.Path[i]), fmtFloat(res.Values[i]))
+		}
+		sb.WriteString("  ...\n")
+		for i := last - sampleN + 1; i <= last; i++ {
+			fmt.Fprintf(&sb, "  step %d: x=%s f=%s\n", i, fmtFloats(res.Path[i]), fmtFloat(res.Values[i]))
+		}
+	}
+	return strings.TrimRight(sb.String(), "\n"), nil
+}
+
+func solveNewtonRoot(params map[string]any) (string, error) {
+	expr, err := parseExprField(params)
+	if err != nil {
+		return "", err
+	}
+	variable, err := reqString(params, "variable")
+	if err != nil {
+		return "", err
+	}
+	start, err := reqFloat(params, "start")
+	if err != nil {
+		return "", err
+	}
+	iterations := optInt(params, "iterations", 0)
+	tol := optFloat(params, "tolerance", 0)
+	fixed, err := varMap(params, "variables")
+	if err != nil {
+		return "", err
+	}
+	res, err := NewtonRoot(expr, variable, start, iterations, tol, fixed)
+	if err != nil {
+		return "", err
+	}
+	last := res.Path[len(res.Path)-1]
+	return fmt.Sprintf("newton's method on %s(%s)=0: converged=%v after %d step(s), root ≈ %s=%s",
+		expr.String(), variable, res.Converged, len(res.Path)-1, variable, fmtFloat(last)), nil
+}
+
+func solveReduce1(params map[string]any, name string, f func([]float64) (float64, error)) (string, error) {
+	data, err := reqFloatSlice(params, "data")
+	if err != nil {
+		return "", err
+	}
+	v, err := f(data)
+	if err != nil {
+		return "", err
+	}
+	return fmt.Sprintf("%s(%s) = %s", name, fmtFloats(data), fmtFloat(v)), nil
+}
+
+func solveReduce1Pop(params map[string]any, name string, f func([]float64, bool) (float64, error)) (string, error) {
+	data, err := reqFloatSlice(params, "data")
+	if err != nil {
+		return "", err
+	}
+	population := optBool(params, "population", false)
+	v, err := f(data, population)
+	if err != nil {
+		return "", err
+	}
+	kind := "sample"
+	if population {
+		kind = "population"
+	}
+	return fmt.Sprintf("%s %s(%s) = %s", kind, name, fmtFloats(data), fmtFloat(v)), nil
+}
+
+func solveReduce2(params map[string]any, name string, f func([]float64, []float64) (float64, error)) (string, error) {
+	x, err := reqFloatSlice(params, "x")
+	if err != nil {
+		return "", err
+	}
+	y, err := reqFloatSlice(params, "y")
+	if err != nil {
+		return "", err
+	}
+	v, err := f(x, y)
+	if err != nil {
+		return "", err
+	}
+	return fmt.Sprintf("%s(x, y) = %s", name, fmtFloat(v)), nil
+}
+
+func solveReduce2Pop(params map[string]any, name string, f func([]float64, []float64, bool) (float64, error)) (string, error) {
+	x, err := reqFloatSlice(params, "x")
+	if err != nil {
+		return "", err
+	}
+	y, err := reqFloatSlice(params, "y")
+	if err != nil {
+		return "", err
+	}
+	population := optBool(params, "population", false)
+	v, err := f(x, y, population)
+	if err != nil {
+		return "", err
+	}
+	kind := "sample"
+	if population {
+		kind = "population"
+	}
+	return fmt.Sprintf("%s %s(x, y) = %s", kind, name, fmtFloat(v)), nil
+}
+
+func solveNormalize(params map[string]any, name string, f func([]float64) ([]float64, error)) (string, error) {
+	data, err := reqFloatSlice(params, "data")
+	if err != nil {
+		return "", err
+	}
+	out, err := f(data)
+	if err != nil {
+		return "", err
+	}
+	return fmt.Sprintf("%s(%s) = %s", name, fmtFloats(data), fmtFloats(out)), nil
+}
+
+func solveNormalDist(params map[string]any, name string, f func(float64, float64, float64) (float64, error)) (string, error) {
+	x, err := reqFloat(params, "x")
+	if err != nil {
+		return "", err
+	}
+	mean := optFloat(params, "mean", 0)
+	std := optFloat(params, "std_dev", 1)
+	v, err := f(x, mean, std)
+	if err != nil {
+		return "", err
+	}
+	return fmt.Sprintf("%s(x=%s, mean=%s, std_dev=%s) = %s", name, fmtFloat(x), fmtFloat(mean), fmtFloat(std), fmtFloat(v)), nil
+}
+
+func solveBinomialPMF(params map[string]any) (string, error) {
+	k := optInt(params, "k", -1)
+	if k < 0 {
+		return "", fmt.Errorf("missing required field %q", "k")
+	}
+	n := optInt(params, "n", -1)
+	if n < 0 {
+		return "", fmt.Errorf("missing required field %q", "n")
+	}
+	p, err := reqFloat(params, "p")
+	if err != nil {
+		return "", err
+	}
+	v, err := BinomialPMF(k, n, p)
+	if err != nil {
+		return "", err
+	}
+	return fmt.Sprintf("P(X=%d) for X~Binomial(n=%d, p=%s) = %s", k, n, fmtFloat(p), fmtFloat(v)), nil
+}
+
+func solveElementwise(params map[string]any, name string, f func(float64) float64) (string, error) {
+	if raw, ok := params["x"].([]any); ok {
+		data := make([]float64, len(raw))
+		for i, v := range raw {
+			fv, ok := toFloat(v)
+			if !ok {
+				return "", fmt.Errorf("field \"x\": element %d is not a number", i)
+			}
+			data[i] = fv
+		}
+		out := make([]float64, len(data))
+		for i, v := range data {
+			out[i] = f(v)
+		}
+		return fmt.Sprintf("%s(%s) = %s", name, fmtFloats(data), fmtFloats(out)), nil
+	}
+	x, err := reqFloat(params, "x")
+	if err != nil {
+		return "", fmt.Errorf(`field "x" must be a number or an array of numbers: %w`, err)
+	}
+	return fmt.Sprintf("%s(%s) = %s", name, fmtFloat(x), fmtFloat(f(x))), nil
+}
+
+func solveSoftmax(params map[string]any) (string, error) {
+	x, err := reqFloatSlice(params, "x")
+	if err != nil {
+		return "", err
+	}
+	out := Softmax(x)
+	sum := 0.0
+	for _, v := range out {
+		sum += v
+	}
+	return fmt.Sprintf("softmax(%s) = %s (sums to %s)", fmtFloats(x), fmtFloats(out), fmtFloat(sum)), nil
+}
+
+func solveLoss(params map[string]any, name string, f func([]float64, []float64) (float64, error)) (string, error) {
+	pred, err := reqFloatSlice(params, "predictions")
+	if err != nil {
+		return "", err
+	}
+	target, err := reqFloatSlice(params, "targets")
+	if err != nil {
+		return "", err
+	}
+	v, err := f(pred, target)
+	if err != nil {
+		return "", err
+	}
+	return fmt.Sprintf("%s(predictions, targets) = %s", name, fmtFloat(v)), nil
+}
+
+func solveLinearRegression(params map[string]any) (string, error) {
+	X, err := reqMatrix(params, "x")
+	if err != nil {
+		return "", err
+	}
+	y, err := reqFloatSlice(params, "y")
+	if err != nil {
+		return "", err
+	}
+	res, err := LinearRegression(X, y)
+	if err != nil {
+		return "", err
+	}
+	var sb strings.Builder
+	fmt.Fprintf(&sb, "ordinary least squares fit: y ≈ %s", fmtFloat(res.Coefficients[0]))
+	for i, c := range res.Coefficients[1:] {
+		fmt.Fprintf(&sb, " + %s*x%d", fmtFloat(c), i+1)
+	}
+	fmt.Fprintf(&sb, "\nR² = %s\n", fmtFloat(res.RSquared))
+	fmt.Fprintf(&sb, "predictions: %s\n", fmtFloats(res.Predictions))
+	fmt.Fprintf(&sb, "residuals: %s", fmtFloats(res.Residuals))
+	return sb.String(), nil
+}
+
+// Operations returns the sorted list of every operation name Solve
+// accepts — used by the agent's tool-catalog description so the
+// supported operation list is generated from the real dispatch table
+// instead of hand-copied prose that can drift out of sync with it.
+func Operations() []string {
+	ops := []string{
+		"evaluate",
+		"matrix_add", "matrix_sub", "matrix_mul", "matrix_transpose", "matrix_determinant",
+		"matrix_inverse", "matrix_rank", "matrix_trace", "dot_product", "vector_norm", "dominant_eigen",
+		"derivative", "gradient", "integral",
+		"gradient_descent", "newton_root",
+		"mean", "variance", "std_dev", "median", "covariance", "correlation",
+		"z_score_normalize", "min_max_normalize", "normal_pdf", "normal_cdf", "binomial_pmf",
+		"sigmoid", "relu", "tanh", "softmax", "mse_loss", "cross_entropy_loss", "linear_regression",
+	}
+	sort.Strings(ops)
+	return ops
+}

pkg/mathkit/statistics.go

@@ -0,0 +1,224 @@
+package mathkit
+
+import (
+	"fmt"
+	"math"
+)
+
+// Mean returns the arithmetic mean of data.
+func Mean(data []float64) (float64, error) {
+	if len(data) == 0 {
+		return 0, fmt.Errorf("mean: empty data set")
+	}
+	var sum float64
+	for _, v := range data {
+		sum += v
+	}
+	return sum / float64(len(data)), nil
+}
+
+// Variance returns the variance of data. population=true divides by
+// n (the population variance); population=false (the default for ML
+// work, where `data` is almost always a sample) divides by n-1
+// (Bessel's correction, the unbiased sample-variance estimator).
+func Variance(data []float64, population bool) (float64, error) {
+	n := len(data)
+	if n == 0 {
+		return 0, fmt.Errorf("variance: empty data set")
+	}
+	if !population && n < 2 {
+		return 0, fmt.Errorf("variance: sample variance needs at least 2 points (use population variance for a single point)")
+	}
+	mean, _ := Mean(data)
+	var ss float64
+	for _, v := range data {
+		d := v - mean
+		ss += d * d
+	}
+	denom := float64(n)
+	if !population {
+		denom = float64(n - 1)
+	}
+	return ss / denom, nil
+}
+
+// StdDev returns the standard deviation (sqrt of Variance).
+func StdDev(data []float64, population bool) (float64, error) {
+	v, err := Variance(data, population)
+	if err != nil {
+		return 0, err
+	}
+	return math.Sqrt(v), nil
+}
+
+// Median returns the median of data (average of the two middle
+// values for an even-length slice). Does not mutate the input.
+func Median(data []float64) (float64, error) {
+	if len(data) == 0 {
+		return 0, fmt.Errorf("median: empty data set")
+	}
+	sorted := append([]float64(nil), data...)
+	insertionSort(sorted)
+	n := len(sorted)
+	if n%2 == 1 {
+		return sorted[n/2], nil
+	}
+	return (sorted[n/2-1] + sorted[n/2]) / 2, nil
+}
+
+// insertionSort is O(n^2) but this package only ever sorts small
+// stats samples, where it beats sort.Float64s' overhead and keeps
+// the package allocation-free apart from the clone itself.
+func insertionSort(a []float64) {
+	for i := 1; i < len(a); i++ {
+		key := a[i]
+		j := i - 1
+		for j >= 0 && a[j] > key {
+			a[j+1] = a[j]
+			j--
+		}
+		a[j+1] = key
+	}
+}
+
+// Covariance returns the sample (or population) covariance of two
+// equal-length data sets — the building block of both Correlation and
+// a PCA covariance matrix.
+func Covariance(x, y []float64, population bool) (float64, error) {
+	if len(x) != len(y) {
+		return 0, fmt.Errorf("covariance: data sets have different lengths (%d vs %d)", len(x), len(y))
+	}
+	n := len(x)
+	if n == 0 {
+		return 0, fmt.Errorf("covariance: empty data set")
+	}
+	if !population && n < 2 {
+		return 0, fmt.Errorf("covariance: sample covariance needs at least 2 points")
+	}
+	mx, _ := Mean(x)
+	my, _ := Mean(y)
+	var sum float64
+	for i := 0; i < n; i++ {
+		sum += (x[i] - mx) * (y[i] - my)
+	}
+	denom := float64(n)
+	if !population {
+		denom = float64(n - 1)
+	}
+	return sum / denom, nil
+}
+
+// Correlation returns the Pearson correlation coefficient of x and y,
+// in [-1, 1].
+func Correlation(x, y []float64) (float64, error) {
+	cov, err := Covariance(x, y, true)
+	if err != nil {
+		return 0, fmt.Errorf("correlation: %w", err)
+	}
+	sx, err := StdDev(x, true)
+	if err != nil {
+		return 0, fmt.Errorf("correlation: %w", err)
+	}
+	sy, err := StdDev(y, true)
+	if err != nil {
+		return 0, fmt.Errorf("correlation: %w", err)
+	}
+	if sx == 0 || sy == 0 {
+		return 0, fmt.Errorf("correlation: a constant data set has undefined correlation")
+	}
+	return cov / (sx * sy), nil
+}
+
+// ZScoreNormalize returns (x - mean) / stddev for each point — the
+// standardization every ML feature-scaling pipeline applies before
+// training a distance- or gradient-sensitive model.
+func ZScoreNormalize(data []float64) ([]float64, error) {
+	mean, err := Mean(data)
+	if err != nil {
+		return nil, fmt.Errorf("z_score_normalize: %w", err)
+	}
+	std, err := StdDev(data, true)
+	if err != nil {
+		return nil, fmt.Errorf("z_score_normalize: %w", err)
+	}
+	if std == 0 {
+		return nil, fmt.Errorf("z_score_normalize: a constant data set has zero standard deviation")
+	}
+	out := make([]float64, len(data))
+	for i, v := range data {
+		out[i] = (v - mean) / std
+	}
+	return out, nil
+}
+
+// MinMaxNormalize rescales data to [0, 1].
+func MinMaxNormalize(data []float64) ([]float64, error) {
+	if len(data) == 0 {
+		return nil, fmt.Errorf("min_max_normalize: empty data set")
+	}
+	lo, hi := data[0], data[0]
+	for _, v := range data {
+		if v < lo {
+			lo = v
+		}
+		if v > hi {
+			hi = v
+		}
+	}
+	if hi == lo {
+		return nil, fmt.Errorf("min_max_normalize: a constant data set has zero range")
+	}
+	out := make([]float64, len(data))
+	for i, v := range data {
+		out[i] = (v - lo) / (hi - lo)
+	}
+	return out, nil
+}
+
+// NormalPDF is the probability density function of N(mean, stdDev^2)
+// at x.
+func NormalPDF(x, mean, stdDev float64) (float64, error) {
+	if stdDev <= 0 {
+		return 0, fmt.Errorf("normal_pdf: std_dev must be positive")
+	}
+	z := (x - mean) / stdDev
+	return math.Exp(-0.5*z*z) / (stdDev * math.Sqrt(2*math.Pi)), nil
+}
+
+// NormalCDF is the cumulative distribution function of N(mean,
+// stdDev^2) at x, computed via math.Erf (the standard closed form:
+// Φ(x) = (1 + erf((x-μ)/(σ√2))) / 2).
+func NormalCDF(x, mean, stdDev float64) (float64, error) {
+	if stdDev <= 0 {
+		return 0, fmt.Errorf("normal_cdf: std_dev must be positive")
+	}
+	return 0.5 * (1 + math.Erf((x-mean)/(stdDev*math.Sqrt2))), nil
+}
+
+// BinomialPMF returns P(X = k) for X ~ Binomial(n, p).
+func BinomialPMF(k, n int, p float64) (float64, error) {
+	if n < 0 || k < 0 || k > n {
+		return 0, fmt.Errorf("binomial_pmf: require 0 <= k <= n, got k=%d n=%d", k, n)
+	}
+	if p < 0 || p > 1 {
+		return 0, fmt.Errorf("binomial_pmf: p must be in [0, 1], got %g", p)
+	}
+	return binomialCoeff(n, k) * math.Pow(p, float64(k)) * math.Pow(1-p, float64(n-k)), nil
+}
+
+// binomialCoeff computes C(n, k) via the multiplicative formula,
+// which stays in float64 range far longer than computing n!/(k!(n-k)!)
+// from raw factorials would.
+func binomialCoeff(n, k int) float64 {
+	if k < 0 || k > n {
+		return 0
+	}
+	if k > n-k {
+		k = n - k
+	}
+	result := 1.0
+	for i := 0; i < k; i++ {
+		result *= float64(n-i) / float64(i+1)
+	}
+	return result
+}