diff --git a/Lab7Interp.m b/Lab7Interp.m new file mode 100644 index 0000000..bd72b4a --- /dev/null +++ b/Lab7Interp.m @@ -0,0 +1,92 @@ +% Lab 7: Polynomial interpolation of f(x) = 1/(1+x^2) on [-5, 5] +% Numerical demo (no plots): each part reports max|f(t) - fit(t)| +% sampled on t = -5:0.01:5. + +f = @(x) 1 ./ (1 + x .^ 2); + +t = -5:0.01:5; +yt = f(t); + +% ========================================================================= +% Part 1 - Full-degree polynomial interpolation at uniform nodes +% ========================================================================= +disp("Part 1: uniform nodes, polyfit(x, y, n) - interpolation"); +ns = [3 6 11 15]; +for k = 1:length(ns) + n = ns(k); + xn = linspace(-5, 5, n+1); + yn = f(xn); + c = polyfit(xn, yn, n); + yp = polyval(c, t); + err = max(abs(yt - yp)); + printf(" n+1 = %3d degree n = %2d max error = %.4f\n", n+1, n, err); +endfor + +% ========================================================================= +% Part 2 - Least-squares polynomial fit (k < n) at 12 uniform nodes +% ========================================================================= +disp(" "); +disp("Part 2: least-squares polyfit(x, y, k) with k < 11 on 12 nodes"); +xn = linspace(-5, 5, 12); +yn = f(xn); +for k = 1:9 + c = polyfit(xn, yn, k); + yp = polyval(c, t); + err = max(abs(yt - yp)); + printf(" degree k = %d max error = %.4f\n", k, err); +endfor + +% ========================================================================= +% Part 3 - Natural cubic spline interpolation at 12 uniform nodes +% ========================================================================= +disp(" "); +disp("Part 3: cubic spline at 12 uniform nodes"); +xn = linspace(-5, 5, 12); +yn = f(xn); +ys = spline(xn, yn, t); +err = max(abs(yt - ys)); +printf(" spline(12 nodes) max error = %.6f\n", err); + +% Also try other counts +for k = 1:length(ns) + n = ns(k); + xn = linspace(-5, 5, n+1); + yn = f(xn); + ys = spline(xn, yn, t); + err = max(abs(yt - ys)); + printf(" spline(%2d nodes) max error = %.6f\n", n+1, err); +endfor + +% ========================================================================= +% Part 4 - Chebyshev nodes for full-degree interpolation +% ========================================================================= +disp(" "); +disp("Part 4: Chebyshev nodes - polyfit(x, y, n) - interpolation"); +a = -5; b = 5; +for k = 1:length(ns) + n = ns(k); + zn = zeros(1, n+1); + for j = 0:n + zn(j+1) = (a+b)/2 + (a-b)/2 * cos(pi*j/n); + endfor + yn = f(zn); + c = polyfit(zn, yn, n); + yp = polyval(c, t); + err = max(abs(yt - yp)); + printf(" n+1 = %3d degree n = %2d max error = %.4f\n", n+1, n, err); +endfor + +% ========================================================================= +% Part 5 - Spline at varied node counts (already partially shown) +% ========================================================================= +disp(" "); +disp("Part 5: spline error vs node count (uniform)"); +counts = [4 7 12 16 25 50]; +for k = 1:length(counts) + m = counts(k); + xn = linspace(-5, 5, m); + yn = f(xn); + ys = spline(xn, yn, t); + err = max(abs(yt - ys)); + printf(" %2d nodes max error = %.6f\n", m, err); +endfor diff --git a/OctiveLean/BigStep.lean b/OctiveLean/BigStep.lean index 8fd72bf..e5a27ff 100644 --- a/OctiveLean/BigStep.lean +++ b/OctiveLean/BigStep.lean @@ -254,6 +254,7 @@ def Value.tag : Value → String | .struct _ => "struct" | .fn _ => "function_handle" | .range _ _ _ => "range" + | .sym _ _ => "sym" | .empty => "empty" theorem litFloat_tag {env env' f v} (h : BigStepExpr env (.lit (.float f)) v env') : v.tag = "double" := by cases h; rfl diff --git a/OctiveLean/Builtins.lean b/OctiveLean/Builtins.lean index ebdc279..c4d8bea 100644 --- a/OctiveLean/Builtins.lean +++ b/OctiveLean/Builtins.lean @@ -1,6 +1,7 @@ import OctiveLean.Value import OctiveLean.Env import OctiveLean.Error +import OctiveLean.SymPyBridge namespace OctiveLean @@ -95,12 +96,14 @@ private def fmtFloat (spec : Char) (prec : Option Nat) (f : Float) : String := let p := prec.getD (if spec == 'g' then 6 else 6) match spec with | 'f' => - -- fixed-point with p decimal places + -- fixed-point with p decimal places (sign-prepend; format absolute value) let scale := Float.ofNat (10 ^ p) - let rounded := Float.round (f * scale) / scale - let intPart := if rounded < 0.0 then (-rounded).floor else rounded.floor - let fracPart := Float.round ((rounded - (if rounded < 0.0 then -intPart else intPart)) * scale) - let intStr := if f < 0.0 then "-" ++ toString intPart.toUInt64 else toString intPart.toUInt64 + let absF := f.abs + let rounded := Float.round (absF * scale) / scale + let intPart := rounded.floor + let fracPart := Float.round ((rounded - intPart) * scale) + let signStr := if f < 0.0 then "-" else "" + let intStr := signStr ++ toString intPart.toUInt64 let fracStr := toString fracPart.toUInt64 let fracPadded := String.ofList (List.replicate (p - fracStr.length) '0') ++ fracStr if p == 0 then intStr else intStr ++ "." ++ fracPadded @@ -214,6 +217,7 @@ def registerAllBuiltins (env : Env) : Env := | .boolean _ | .boolMat _ _ _ => "logical" | .string _ => "char" | .cell _ _ _ => "cell" | .struct _ => "struct" | .fn _ => "function_handle" + | .sym _ _ => "sym" return #[Value.string cls] | none => return #[Value.string "unknown"]) |>.registerBuiltin "isnumeric" (fun (args : Array Value) => do @@ -404,8 +408,17 @@ def registerAllBuiltins (env : Env) : Env := -- ── Type conversion ─────────────────────────────────────────────────────── |>.registerBuiltin "double" (fun (args : Array Value) => do match args[0]? with - | some v => return #[Value.scalar (← asFloat "double" v)] - | none => return #[Value.empty]) + | some v => + match v with + | Value.sym sr _ => + match (← SymPyBridge.runRaw s!"print(repr(float(({sr}).evalf())))") with + | .ok s => + match parseFloatStr? s.trimAscii.toString with + | some f => return #[Value.scalar f] + | none => throw (IO.userError s!"double: cannot convert sym '{s}' to float") + | .error e => throw (IO.userError s!"double: {e}") + | _ => return #[Value.scalar (← asFloat "double" v)] + | none => return #[Value.empty]) |>.registerBuiltin "logical" (fun (args : Array Value) => do match args[0]? with | some v => return #[Value.boolean ((← asFloat "logical" v) != 0.0)] @@ -434,5 +447,427 @@ def registerAllBuiltins (env : Env) : Env := |>.registerBuiltin "quit" (fun (_ : Array Value) => do IO.Process.exit 0 return (#[] : Array Value)) + -- ── Numerical: linear solve, polyfit, polyval, spline ──────────────────── + |>.registerBuiltin "linsolve" (fun (args : Array Value) => do + if args.size < 2 then throw (IO.userError "linsolve: expected (A, b)") + match args[0]!.materialize, args[1]!.materialize with + | .matrix n m a, .matrix nb _ b => + if n != m || nb != n then + throw (IO.userError s!"linsolve: A must be square and match b ({n}×{m} vs b={nb})") + let mut M : Array Float := a + let mut bv : Array Float := b + for i in List.range n do + let mut maxRow := i + let mut maxV := (M[i * n + i]!).abs + for k in List.range (n - i - 1) do + let kk := i + 1 + k + let v := (M[kk * n + i]!).abs + if v > maxV then maxRow := kk; maxV := v + if maxRow != i then + for j in List.range n do + let t := M[i * n + j]! + M := M.set! (i * n + j) M[maxRow * n + j]! + M := M.set! (maxRow * n + j) t + let tb := bv[i]! + bv := bv.set! i bv[maxRow]! + bv := bv.set! maxRow tb + let pivot := M[i * n + i]! + if pivot.abs < 1e-15 then + throw (IO.userError "linsolve: singular matrix") + for k in List.range (n - i - 1) do + let kk := i + 1 + k + let factor := M[kk * n + i]! / pivot + for j in List.range (n - i) do + let jj := i + j + M := M.set! (kk * n + jj) (M[kk * n + jj]! - factor * M[i * n + jj]!) + bv := bv.set! kk (bv[kk]! - factor * bv[i]!) + let mut x : Array Float := arrFill n 0.0 + for ii in List.range n do + let i := n - 1 - ii + let mut s := bv[i]! + for k in List.range (n - i - 1) do + let j := i + 1 + k + s := s - M[i * n + j]! * x[j]! + x := x.set! i (s / M[i * n + i]!) + return #[Value.matrix n 1 x] + | _, _ => throw (IO.userError "linsolve: expected matrix arguments")) + |>.registerBuiltin "polyval" (fun (args : Array Value) => do + if args.size < 2 then throw (IO.userError "polyval: expected (c, x)") + let c := flattenV args[0]! + let xs := flattenV args[1]! + if c.isEmpty then throw (IO.userError "polyval: empty coefficients") + let eval := fun (x : Float) => Id.run do + let mut y := c[0]! + for i in List.range (c.size - 1) do + y := y * x + c[i + 1]! + y + let ys : Array Float := xs.map eval + match args[1]!.materialize with + | .scalar _ => return #[Value.scalar (ys[0]!)] + | .matrix r co _ => return #[Value.matrix r co ys] + | .range _ _ _ => return #[Value.matrix 1 ys.size ys] + | _ => return #[Value.matrix 1 ys.size ys]) + |>.registerBuiltin "polyfit" (fun (args : Array Value) => do + if args.size < 3 then throw (IO.userError "polyfit: expected (x, y, n)") + let xs := flattenV args[0]! + let ys := flattenV args[1]! + let n ← asNat "polyfit" args[2]! + let m := xs.size + if ys.size != m then throw (IO.userError "polyfit: x and y must be same length") + if n + 1 > m then throw (IO.userError s!"polyfit: degree {n} requires at least {n+1} points") + -- Build Vandermonde V[i,j] = xs[i]^(n - j) (i in 0..m, j in 0..n) + let cols := n + 1 + let V : Array Float := Id.run do + let mut v := arrFill (m * cols) 0.0 + for i in List.range m do + let mut p := 1.0 + for k in List.range cols do + v := v.set! (i * cols + (n - k)) p + p := p * xs[i]! + v + -- Normal equations: A = V^T V (cols × cols), b = V^T y + let A : Array Float := Id.run do + let mut a := arrFill (cols * cols) 0.0 + for i in List.range cols do + for j in List.range cols do + let mut s := 0.0 + for k in List.range m do + s := s + V[k * cols + i]! * V[k * cols + j]! + a := a.set! (i * cols + j) s + a + let bv : Array Float := Id.run do + let mut b := arrFill cols 0.0 + for i in List.range cols do + let mut s := 0.0 + for k in List.range m do + s := s + V[k * cols + i]! * ys[k]! + b := b.set! i s + b + -- Solve A c = bv via in-place Gaussian elimination with partial pivot + let mut M := A + let mut rhs := bv + let nn := cols + for i in List.range nn do + let mut maxRow := i + let mut maxV := (M[i * nn + i]!).abs + for k in List.range (nn - i - 1) do + let kk := i + 1 + k + let v := (M[kk * nn + i]!).abs + if v > maxV then maxRow := kk; maxV := v + if maxRow != i then + for j in List.range nn do + let t := M[i * nn + j]! + M := M.set! (i * nn + j) M[maxRow * nn + j]! + M := M.set! (maxRow * nn + j) t + let tb := rhs[i]! + rhs := rhs.set! i rhs[maxRow]! + rhs := rhs.set! maxRow tb + let pivot := M[i * nn + i]! + if pivot.abs < 1e-15 then + throw (IO.userError "polyfit: singular normal equations") + for k in List.range (nn - i - 1) do + let kk := i + 1 + k + let factor := M[kk * nn + i]! / pivot + for j in List.range (nn - i) do + let jj := i + j + M := M.set! (kk * nn + jj) (M[kk * nn + jj]! - factor * M[i * nn + jj]!) + rhs := rhs.set! kk (rhs[kk]! - factor * rhs[i]!) + let mut c := arrFill nn 0.0 + for ii in List.range nn do + let i := nn - 1 - ii + let mut s := rhs[i]! + for k in List.range (nn - i - 1) do + let j := i + 1 + k + s := s - M[i * nn + j]! * c[j]! + c := c.set! i (s / M[i * nn + i]!) + return #[Value.matrix 1 nn c]) + |>.registerBuiltin "spline" (fun (args : Array Value) => do + if args.size < 3 then throw (IO.userError "spline: expected (x, y, t)") + let xs := flattenV args[0]! + let ys := flattenV args[1]! + let ts := flattenV args[2]! + let n := xs.size + if ys.size != n || n < 2 then throw (IO.userError "spline: bad input") + let nseg := n - 1 + let h : Array Float := Id.run do + let mut h := arrFill nseg 0.0 + for i in List.range nseg do h := h.set! i (xs[i+1]! - xs[i]!) + h + -- Solve tridiagonal for M[1..n-2], with M[0]=M[n-1]=0 (natural) + let mut M := arrFill n 0.0 + if n >= 3 then + let inner := n - 2 + let mut a := arrFill inner 0.0 + let mut b := arrFill inner 0.0 + let mut c := arrFill inner 0.0 + let mut d := arrFill inner 0.0 + for i in List.range inner do + let i1 := i + 1 + let hL := h[i1 - 1]! + let hR := h[i1]! + a := a.set! i hL + b := b.set! i (2.0 * (hL + hR)) + c := c.set! i hR + d := d.set! i (6.0 * ((ys[i1+1]! - ys[i1]!) / hR - (ys[i1]! - ys[i1-1]!) / hL)) + -- Thomas algorithm + for i in List.range (inner - 1) do + let ii := i + 1 + let factor := a[ii]! / b[i]! + b := b.set! ii (b[ii]! - factor * c[i]!) + d := d.set! ii (d[ii]! - factor * d[i]!) + let mut sol := arrFill inner 0.0 + sol := sol.set! (inner - 1) (d[inner-1]! / b[inner-1]!) + for ii in List.range (inner - 1) do + let i := inner - 2 - ii + sol := sol.set! i ((d[i]! - c[i]! * sol[i+1]!) / b[i]!) + for i in List.range inner do M := M.set! (i + 1) sol[i]! + -- Evaluate at each t + let evalAt := fun (t : Float) => Id.run do + let mut idx := 0 + for k in List.range nseg do if xs[k]! <= t then idx := k + if t > xs[n-1]! then idx := nseg - 1 + let i := idx + let hi := h[i]! + let xi := xs[i]!; let xi1 := xs[i+1]! + let yi := ys[i]!; let yi1 := ys[i+1]! + let Mi := M[i]!; let Mi1 := M[i+1]! + let A1 := Mi * (xi1 - t)^3.0 / (6.0 * hi) + let A2 := Mi1 * (t - xi)^3.0 / (6.0 * hi) + let A3 := (yi / hi - Mi * hi / 6.0) * (xi1 - t) + let A4 := (yi1 / hi - Mi1 * hi / 6.0) * (t - xi) + A1 + A2 + A3 + A4 + let out : Array Float := ts.map evalAt + match args[2]!.materialize with + | .scalar _ => return #[Value.scalar out[0]!] + | .matrix r co _ => return #[Value.matrix r co out] + | .range _ _ _ => return #[Value.matrix 1 out.size out] + | _ => return #[Value.matrix 1 out.size out]) + -- ── Symbolic Math (SymPy bridge) ───────────────────────────────────────── + -- Architecture mirrors GNU Octave's `symbolic` package: each builtin is a + -- thin wrapper that converts arguments to a Python expression and forwards + -- to a persistent SymPy subprocess. See `OctiveLean/SymPyBridge.lean`. + |>.registerBuiltin "sym" (fun (args : Array Value) => do + match args[0]? with + | some v => + let py ← SymPyBridge.toSympy v + return #[← SymPyBridge.emit py] + | none => throw (IO.userError "sym: expected 1 argument")) + |>.registerBuiltin "syms" (fun (args : Array Value) => do + -- Returns one Sym per argument — invoked as `[x,y,z] = syms('x','y','z')`. + let mut out : Array Value := #[] + for a in args do + match a with + | .string s => out := out.push (← SymPyBridge.emit s!"symbols('{s}')") + | _ => throw (IO.userError "syms: expected string arg") + return out) + |>.registerBuiltin "diff" (fun (args : Array Value) => do + match args.size with + | 0 => throw (IO.userError "diff: expected at least 1 argument") + | 1 => + let f ← SymPyBridge.toSympy args[0]! + return #[← SymPyBridge.emit s!"diff({f})"] + | _ => + let f ← SymPyBridge.toSympy args[0]! + let v ← SymPyBridge.toSympy args[1]! + if h : args.size >= 3 then + let n ← SymPyBridge.toSympy args[2]! + return #[← SymPyBridge.emit s!"diff({f}, {v}, {n})"] + else + return #[← SymPyBridge.emit s!"diff({f}, {v})"]) + |>.registerBuiltin "int" (fun (args : Array Value) => do + match args.size with + | 0 => throw (IO.userError "int: expected at least 1 argument") + | 1 => + let f ← SymPyBridge.toSympy args[0]! + return #[← SymPyBridge.emit s!"integrate({f})"] + | 2 => + let f ← SymPyBridge.toSympy args[0]! + let v ← SymPyBridge.toSympy args[1]! + return #[← SymPyBridge.emit s!"integrate({f}, {v})"] + | _ => + -- int(f, x, a, b) — definite integral + let f ← SymPyBridge.toSympy args[0]! + let v ← SymPyBridge.toSympy args[1]! + let a ← SymPyBridge.toSympy args[2]! + let b ← SymPyBridge.toSympy args[3]! + return #[← SymPyBridge.emit s!"integrate({f}, ({v}, {a}, {b}))"]) + |>.registerBuiltin "subs" (fun (args : Array Value) => do + if args.size < 3 then throw (IO.userError "subs: expected (expr, var, val)") + let f ← SymPyBridge.toSympy args[0]! + let v ← SymPyBridge.toSympy args[1]! + let r ← SymPyBridge.toSympy args[2]! + return #[← SymPyBridge.emit s!"({f}).subs({v}, {r})"]) + |>.registerBuiltin "simplify" (fun (args : Array Value) => do + let f ← SymPyBridge.toSympy args[0]! + return #[← SymPyBridge.emit s!"simplify({f})"]) + |>.registerBuiltin "expand" (fun (args : Array Value) => do + let f ← SymPyBridge.toSympy args[0]! + return #[← SymPyBridge.emit s!"expand({f})"]) + |>.registerBuiltin "factor" (fun (args : Array Value) => do + let f ← SymPyBridge.toSympy args[0]! + return #[← SymPyBridge.emit s!"factor({f})"]) + |>.registerBuiltin "collect" (fun (args : Array Value) => do + if args.size < 2 then throw (IO.userError "collect: expected (expr, var)") + let f ← SymPyBridge.toSympy args[0]! + let v ← SymPyBridge.toSympy args[1]! + return #[← SymPyBridge.emit s!"collect({f}, {v})"]) + |>.registerBuiltin "solve" (fun (args : Array Value) => do + match args.size with + | 0 => throw (IO.userError "solve: expected at least 1 argument") + | 1 => + let f ← SymPyBridge.toSympy args[0]! + return #[← SymPyBridge.emit s!"solve({f})"] + | _ => + let f ← SymPyBridge.toSympy args[0]! + let v ← SymPyBridge.toSympy args[1]! + return #[← SymPyBridge.emit s!"solve({f}, {v})"]) + |>.registerBuiltin "taylor" (fun (args : Array Value) => do + match args.size with + | 0 => throw (IO.userError "taylor: expected at least 1 argument") + | 1 => + let f ← SymPyBridge.toSympy args[0]! + return #[← SymPyBridge.emit s!"series({f}).removeO()"] + | 2 => + let f ← SymPyBridge.toSympy args[0]! + let v ← SymPyBridge.toSympy args[1]! + return #[← SymPyBridge.emit s!"series({f}, {v}).removeO()"] + | _ => + let f ← SymPyBridge.toSympy args[0]! + let v ← SymPyBridge.toSympy args[1]! + let a ← SymPyBridge.toSympy args[2]! + if h : args.size >= 4 then + let n ← SymPyBridge.toSympy args[3]! + return #[← SymPyBridge.emit s!"series({f}, {v}, {a}, {n}).removeO()"] + else + return #[← SymPyBridge.emit s!"series({f}, {v}, {a}).removeO()"]) + |>.registerBuiltin "limit" (fun (args : Array Value) => do + if args.size < 3 then throw (IO.userError "limit: expected (expr, var, point)") + let f ← SymPyBridge.toSympy args[0]! + let v ← SymPyBridge.toSympy args[1]! + let p ← SymPyBridge.toSympy args[2]! + if h : args.size >= 4 then + match args[3]! with + | .string "left" => return #[← SymPyBridge.emit s!"limit({f}, {v}, {p}, '-')"] + | .string "right" => return #[← SymPyBridge.emit s!"limit({f}, {v}, {p}, '+')"] + | _ => return #[← SymPyBridge.emit s!"limit({f}, {v}, {p})"] + else + return #[← SymPyBridge.emit s!"limit({f}, {v}, {p})"]) + |>.registerBuiltin "jacobian" (fun (args : Array Value) => do + if args.size < 2 then throw (IO.userError "jacobian: expected (f, vars)") + let f ← SymPyBridge.toSympy args[0]! + let v ← SymPyBridge.toSympy args[1]! + return #[← SymPyBridge.emit s!"Matrix([{f}]).jacobian({v})"]) + |>.registerBuiltin "gradient" (fun (args : Array Value) => do + if args.size < 2 then throw (IO.userError "gradient: expected (f, vars)") + let f ← SymPyBridge.toSympy args[0]! + let v ← SymPyBridge.toSympy args[1]! + return #[← SymPyBridge.emit s!"Matrix([{f}]).jacobian({v}).T"]) + |>.registerBuiltin "hessian" (fun (args : Array Value) => do + if args.size < 2 then throw (IO.userError "hessian: expected (f, vars)") + let f ← SymPyBridge.toSympy args[0]! + let v ← SymPyBridge.toSympy args[1]! + return #[← SymPyBridge.emit s!"hessian({f}, {v})"]) + |>.registerBuiltin "coeffs" (fun (args : Array Value) => do + let f ← SymPyBridge.toSympy args[0]! + if h : args.size >= 2 then + let v ← SymPyBridge.toSympy args[1]! + return #[← SymPyBridge.emit s!"Poly({f}, {v}).all_coeffs()"] + else + return #[← SymPyBridge.emit s!"Poly({f}).all_coeffs()"]) + |>.registerBuiltin "lhs" (fun (args : Array Value) => do + let f ← SymPyBridge.toSympy args[0]! + return #[← SymPyBridge.emit s!"({f}).lhs"]) + |>.registerBuiltin "rhs" (fun (args : Array Value) => do + let f ← SymPyBridge.toSympy args[0]! + return #[← SymPyBridge.emit s!"({f}).rhs"]) + |>.registerBuiltin "latex" (fun (args : Array Value) => do + let f ← SymPyBridge.toSympy args[0]! + match (← SymPyBridge.runRaw s!"print(latex({f}))") with + | .ok s => return #[Value.string (s.trimAscii.toString)] + | .error e => throw (IO.userError s!"latex: {e}")) + |>.registerBuiltin "pretty" (fun (args : Array Value) => do + let f ← SymPyBridge.toSympy args[0]! + match (← SymPyBridge.runRaw s!"print(pretty({f}, use_unicode=False))") with + | .ok s => IO.println s; return #[] + | .error e => throw (IO.userError s!"pretty: {e}")) + |>.registerBuiltin "vpa" (fun (args : Array Value) => do + let f ← SymPyBridge.toSympy args[0]! + let n : String := if h : args.size >= 2 then + match args[1]! with + | .scalar f => toString f.toInt64 + | _ => "15" + else "15" + return #[← SymPyBridge.emit s!"N({f}, {n})"]) + |>.registerBuiltin "symsum" (fun (args : Array Value) => do + if args.size < 4 then throw (IO.userError "symsum: expected (expr, var, lo, hi)") + let f ← SymPyBridge.toSympy args[0]! + let v ← SymPyBridge.toSympy args[1]! + let lo ← SymPyBridge.toSympy args[2]! + let hi ← SymPyBridge.toSympy args[3]! + return #[← SymPyBridge.emit s!"summation({f}, ({v}, {lo}, {hi}))"]) + |>.registerBuiltin "laplacian" (fun (args : Array Value) => do + if args.size < 2 then throw (IO.userError "laplacian: expected (f, vars)") + let f ← SymPyBridge.toSympy args[0]! + let v ← SymPyBridge.toSympy args[1]! + return #[← SymPyBridge.emit s!"sum(diff({f}, _v, 2) for _v in {v})"]) + |>.registerBuiltin "divergence" (fun (args : Array Value) => do + if args.size < 2 then throw (IO.userError "divergence: expected (F, vars)") + let f ← SymPyBridge.toSympy args[0]! + let v ← SymPyBridge.toSympy args[1]! + return #[← SymPyBridge.emit s!"sum(diff(_F[i], {v}[i]) for i, _F in enumerate([{f}] * 1) for i in range(len({v})))"]) + |>.registerBuiltin "rewrite" (fun (args : Array Value) => do + if args.size < 2 then throw (IO.userError "rewrite: expected (expr, target)") + let f ← SymPyBridge.toSympy args[0]! + let target := match args[1]! with | .string s => s | _ => "sin" + return #[← SymPyBridge.emit s!"({f}).rewrite({target})"]) + |>.registerBuiltin "resultant" (fun (args : Array Value) => do + if args.size < 2 then throw (IO.userError "resultant: expected (p, q[, var])") + let p ← SymPyBridge.toSympy args[0]! + let q ← SymPyBridge.toSympy args[1]! + if h : args.size >= 3 then + let v ← SymPyBridge.toSympy args[2]! + return #[← SymPyBridge.emit s!"resultant({p}, {q}, {v})"] + else + return #[← SymPyBridge.emit s!"resultant({p}, {q})"]) + |>.registerBuiltin "series" (fun (args : Array Value) => do + let f ← SymPyBridge.toSympy args[0]! + if h : args.size >= 2 then + let v ← SymPyBridge.toSympy args[1]! + return #[← SymPyBridge.emit s!"series({f}, {v})"] + else + return #[← SymPyBridge.emit s!"series({f})"]) + |>.registerBuiltin "isolate" (fun (args : Array Value) => do + if args.size < 2 then throw (IO.userError "isolate: expected (eq, var)") + let f ← SymPyBridge.toSympy args[0]! + let v ← SymPyBridge.toSympy args[1]! + return #[← SymPyBridge.emit s!"Eq({v}, solve({f}, {v})[0])"]) + |>.registerBuiltin "symfun" (fun (args : Array Value) => do + match args[0]? with + | some v => + match v with + | Value.string n => + match (← SymPyBridge.runRaw s!"{n} = Function('{n}')") with + | .ok _ => return #[← SymPyBridge.emit s!"Function('{n}')"] + | .error e => throw (IO.userError s!"symfun: {e}") + | _ => throw (IO.userError "symfun: expected name string") + | none => throw (IO.userError "symfun: expected name string")) + |>.registerBuiltin "dsolve" (fun (args : Array Value) => do + let f ← SymPyBridge.toSympy args[0]! + if h : args.size >= 2 then + let y ← SymPyBridge.toSympy args[1]! + return #[← SymPyBridge.emit s!"dsolve({f}, {y})"] + else + return #[← SymPyBridge.emit s!"dsolve({f})"]) + |>.registerBuiltin "piecewise" (fun (args : Array Value) => do + -- piecewise(cond1, val1, cond2, val2, ...) → Piecewise((val1, cond1), ...) + let mut parts : Array String := #[] + let mut i := 0 + while h : i + 1 < args.size do + let c ← SymPyBridge.toSympy args[i]! + let v ← SymPyBridge.toSympy args[i+1]! + parts := parts.push s!"({v}, {c})" + i := i + 2 + let body := String.intercalate ", " parts.toList + return #[← SymPyBridge.emit s!"Piecewise({body})"]) end OctiveLean diff --git a/OctiveLean/DSL.lean b/OctiveLean/DSL.lean index 6995d92..28088d8 100644 --- a/OctiveLean/DSL.lean +++ b/OctiveLean/DSL.lean @@ -1,42 +1,48 @@ import OctiveLean.Eval import OctiveLean.Builtins import OctiveLean.PlotData -import OctiveLean.PlotWidget import OctiveLean.PlotBuiltins +import OctiveLean.PlotWidget import ProofWidgets.Component.HtmlDisplay import Lean /-! # OctiveLean Syntax DSL -Octave as a first-class Lean 4 syntax category. The LSP sees every keyword, -operator and structure — giving real syntax highlighting, hover and completion -inside `octave! ... octave_end` blocks. +Octave embedded as a Lean 4 syntax category. The LSP recognizes every +keyword and operator inside an `octave! { ... }` block — giving real +syntax highlighting, hover, and completion. ## Usage -```lean -octave! +``` +octave! { x = 42; for k = 1:5 x = x + k; endfor disp(x) -octave_end +} ``` -## Syntax notes (differences from standard Octave) -- Block closers: `endif` `endfor` `endwhile` `endfunction` `endswitch` `endtry` - (Octave supports these as aliases for `end` — they work in real Octave too) -- Outer block: `octave!` … `octave_end` -- Strings: use Lean double-quotes `"hello"` (not `'hello'`) -- Matrix literals: `[1.0, 2.0, 3.0]` (row vector), `[[1.0, 2.0], [3.0, 4.0]]` (matrix) -- Comments: `--` Lean style (parser limitation — `%` is the modulo token) -- `true` / `false` are valid Octave literals +## Departures from standard Octave + +- Outer block: `octave! { ... }` (curly braces avoid collisions with Lean's `end`) +- `endif` / `endfor` / `endwhile` / `endswitch` / `end_try_catch` / `endfunction` + to close each control structure (these are real Octave keywords too) +- Strings: `"..."` (Lean string literals — `'...'` is not supported) +- Comments: `--` (Lean style — `%` is the modulo operator token) +- Matrix rows are separated by `;`, columns by `,`: `[1.0, 2.0; 3.0, 4.0]` -/ -open OctiveLean +namespace OctiveLean.DSL + open Lean +open OctiveLean + +/-- Convert a `String` to a `TSyntax `term` whose representation is a string literal. -/ +private def quoteStr (s : String) : TSyntax `term := + ⟨Lean.Syntax.mkStrLit s⟩ -- ───────────────────────────────────────────────────────────────── -- Syntax categories @@ -44,341 +50,355 @@ open Lean declare_syntax_cat octExpr declare_syntax_cat octStmt +declare_syntax_cat octRow +declare_syntax_cat octMatBody -- ───────────────────────────────────────────────────────────────── -- EXPRESSIONS -- ───────────────────────────────────────────────────────────────── -syntax num : octExpr -syntax scientific : octExpr -syntax str : octExpr -syntax ident : octExpr -syntax "(" octExpr ")" : octExpr +-- Literals (true/false handled as identifiers in convExpr) +-- Atomic forms must be at :max so they can be left-args of postfix rules. +syntax:max num : octExpr +syntax:max scientific : octExpr +syntax:max str : octExpr --- Unary -syntax:90 "-" octExpr:90 : octExpr -syntax:90 "!" octExpr:90 : octExpr +-- Identifier +syntax:max ident : octExpr --- Arithmetic -syntax:75 octExpr:76 "^" octExpr:75 : octExpr -syntax:75 octExpr:76 ".^" octExpr:75 : octExpr -syntax:70 octExpr:70 "*" octExpr:71 : octExpr -syntax:70 octExpr:70 "/" octExpr:71 : octExpr -syntax:70 octExpr:70 ".*" octExpr:71 : octExpr -syntax:70 octExpr:70 "./" octExpr:71 : octExpr -syntax:65 octExpr:65 "+" octExpr:66 : octExpr -syntax:65 octExpr:65 "-" octExpr:66 : octExpr +-- Grouped +syntax:max "(" octExpr ")" : octExpr --- Comparison -syntax:50 octExpr:51 "==" octExpr:51 : octExpr -syntax:50 octExpr:51 "!=" octExpr:51 : octExpr -syntax:50 octExpr:51 "<" octExpr:51 : octExpr -syntax:50 octExpr:51 "<=" octExpr:51 : octExpr -syntax:50 octExpr:51 ">" octExpr:51 : octExpr -syntax:50 octExpr:51 ">=" octExpr:51 : octExpr - --- Logical -syntax:40 octExpr:40 "&&" octExpr:41 : octExpr -syntax:40 octExpr:40 "||" octExpr:41 : octExpr -syntax:35 octExpr:35 "&" octExpr:36 : octExpr -syntax:35 octExpr:35 "|" octExpr:36 : octExpr - --- Range a:b and a:step:b (left-assoc; (a:step):b is the three-part form) -syntax:20 octExpr:20 ":" octExpr:21 : octExpr - --- Call / index: f(a, b, ...) — ident-based to avoid left-recursion issues -syntax ident "(" octExpr,* ")" : octExpr - --- Struct field: s.field (left-recursive, works for simple s.f cases) -syntax:max octExpr:max noWs "." noWs ident : octExpr - --- Dynamic field: s.(expr) — ".(" is a single token in Lean 4 --- Note: nested use like disp(p.(f)) is limited; use as a statement or top-level expr -syntax ident ".(" octExpr ")" : octExpr +-- Matrix literal: rows separated by ';', columns by ',' +syntax octExpr,+ : octRow +syntax octRow : octMatBody +syntax octRow ";" octMatBody : octMatBody +syntax:max "[" octMatBody "]" : octExpr +syntax:max "[" "]" : octExpr -- empty matrix -- Function handles -syntax "@" ident : octExpr -syntax "@" "(" ident,* ")" octExpr : octExpr +syntax:max "@" ident : octExpr +syntax:max "@(" ident,* ")" octExpr : octExpr --- Vector / matrix literals --- [a, b, c] = row vector; [[a,b], [c,d]] = matrix -syntax "[" octExpr,* "]" : octExpr +-- Unary +syntax:75 "-" octExpr:75 : octExpr +syntax:75 "!" octExpr:75 : octExpr + +-- Power (right-associative) +syntax:74 octExpr:75 " ^ " octExpr:74 : octExpr +syntax:74 octExpr:75 " .^ " octExpr:74 : octExpr + +-- Multiplication / division (left-associative) +syntax:70 octExpr:70 " * " octExpr:71 : octExpr +syntax:70 octExpr:70 " / " octExpr:71 : octExpr +syntax:70 octExpr:70 " .* " octExpr:71 : octExpr +syntax:70 octExpr:70 " ./ " octExpr:71 : octExpr + +-- Addition / subtraction +syntax:65 octExpr:65 " + " octExpr:66 : octExpr +syntax:65 octExpr:65 " - " octExpr:66 : octExpr + +-- Range a:b +syntax:60 octExpr:61 " : " octExpr:61 : octExpr + +-- Comparison +syntax:50 octExpr:51 " == " octExpr:51 : octExpr +syntax:50 octExpr:51 " != " octExpr:51 : octExpr +syntax:50 octExpr:51 " <= " octExpr:51 : octExpr +syntax:50 octExpr:51 " >= " octExpr:51 : octExpr +syntax:50 octExpr:51 " < " octExpr:51 : octExpr +syntax:50 octExpr:51 " > " octExpr:51 : octExpr + +-- Logical +syntax:40 octExpr:40 " && " octExpr:41 : octExpr +syntax:40 octExpr:40 " || " octExpr:41 : octExpr +syntax:40 octExpr:40 " & " octExpr:41 : octExpr +syntax:40 octExpr:40 " | " octExpr:41 : octExpr + +-- Function call / matrix index +syntax:max octExpr:max "(" octExpr,* ")" : octExpr -- ───────────────────────────────────────────────────────────────── -- STATEMENTS -- ───────────────────────────────────────────────────────────────── -syntax octExpr : octStmt -syntax octExpr ";" : octStmt +-- Expression statement +syntax octExpr : octStmt +syntax octExpr ";" : octStmt -syntax ident " = " octExpr : octStmt -syntax ident " = " octExpr ";" : octStmt +-- Assignment +syntax ident " = " octExpr : octStmt +syntax ident " = " octExpr ";" : octStmt -syntax "[" ident,+ "]" " = " octExpr : octStmt -syntax "[" ident,+ "]" " = " octExpr ";" : octStmt +-- Multi-assignment +syntax "[" ident,+ "]" " = " octExpr : octStmt +syntax "[" ident,+ "]" " = " octExpr ";" : octStmt --- Struct field assignment: s.f = expr -syntax ident noWs "." noWs ident " = " octExpr : octStmt -syntax ident noWs "." noWs ident " = " octExpr ";" : octStmt - --- IF / ENDIF +-- IF / ELSEIF / ELSE / ENDIF syntax "if" octExpr octStmt* - ("elseif" octExpr octStmt*)* - ("else" octStmt*)? - "endif" : octStmt + ("elseif" octExpr octStmt*)* + ("else" octStmt*)? + "endif" : octStmt -- FOR / ENDFOR -syntax "for" ident " = " octExpr octStmt* "endfor" : octStmt +syntax "for" ident " = " octExpr octStmt* "endfor" : octStmt -- WHILE / ENDWHILE -syntax "while" octExpr octStmt* "endwhile" : octStmt +syntax "while" octExpr octStmt* "endwhile" : octStmt --- SWITCH / ENDSWITCH +-- SWITCH / CASE / OTHERWISE / ENDSWITCH syntax "switch" octExpr - ("case" octExpr octStmt*)* - ("otherwise" octStmt*)? - "endswitch" : octStmt + ("case" octExpr octStmt*)* + ("otherwise" octStmt*)? + "endswitch" : octStmt --- TRY / ENDTRY +-- TRY / CATCH / END_TRY_CATCH syntax "try" octStmt* - ("catch" ident octStmt*)? - "endtry" : octStmt + ("catch" ident octStmt*)? + "end_try_catch" : octStmt -syntax "return" : octStmt -syntax "break" : octStmt -syntax "continue" : octStmt +-- Control flow +syntax "return" : octStmt +syntax "break" : octStmt +syntax "continue" : octStmt -syntax "global" ident,+ : octStmt -syntax "clear" ident,+ : octStmt +-- Scope +syntax "global" ident+ : octStmt +syntax "clear" ident+ : octStmt --- Function definitions +-- Function definition syntax "function" ident " = " ident "(" ident,* ")" - octStmt* "endfunction" : octStmt + octStmt* "endfunction" : octStmt syntax "function" "[" ident,+ "]" " = " ident "(" ident,* ")" - octStmt* "endfunction" : octStmt + octStmt* "endfunction" : octStmt syntax "function" ident "(" ident,* ")" - octStmt* "endfunction" : octStmt - --- Top-level blocks -syntax (name := octaveRun) "octave!" octStmt* "octave_end" : command -syntax (name := octaveStmts) "octave_stmts!" ident octStmt* "octave_end" : command + octStmt* "endfunction" : octStmt -- ───────────────────────────────────────────────────────────────── --- Helpers +-- Macro conversion: octExpr → Term (of type OctiveLean.Expr) -- ───────────────────────────────────────────────────────────────── -private def strTerm (s : String) : TSyntax `term := ⟨Syntax.mkStrLit s⟩ - -private def identStr (id : TSyntax `ident) : TSyntax `term := - strTerm id.getId.toString - --- ───────────────────────────────────────────────────────────────── --- convExpr : octExpr syntax → term of type OctiveLean.Expr --- ───────────────────────────────────────────────────────────────── - -private partial def convExpr : TSyntax `octExpr → MacroM (TSyntax `term) - | `(octExpr| $n:num) => `(Expr.lit (.float ($n : Float))) +private partial def convExpr (e : Syntax) : MacroM (TSyntax `term) := do + match e with + -- Literals + | `(octExpr| $n:num) => `(Expr.lit (.float ($n : Float))) | `(octExpr| $f:scientific) => `(Expr.lit (.float ($f : Float))) | `(octExpr| $s:str) => `(Expr.lit (.str $s)) - | `(octExpr| $id:ident) => + -- Identifier (with special handling for `true`/`false`) + | `(octExpr| $id:ident) => match id.getId.toString with | "true" => `(Expr.lit (.bool true)) | "false" => `(Expr.lit (.bool false)) - | name => `(Expr.ident $(strTerm name)) - | `(octExpr| ($inner:octExpr)) => convExpr inner - | `(octExpr| - $x:octExpr) => do `(Expr.unop .neg $(← convExpr x)) - | `(octExpr| ! $x:octExpr) => do `(Expr.unop .lnot $(← convExpr x)) - | `(octExpr| $a:octExpr ^ $b:octExpr) => do `(Expr.binop .pow $(← convExpr a) $(← convExpr b)) - | `(octExpr| $a:octExpr .^ $b:octExpr) => do `(Expr.binop .epow $(← convExpr a) $(← convExpr b)) - | `(octExpr| $a:octExpr * $b:octExpr) => do `(Expr.binop .mul $(← convExpr a) $(← convExpr b)) - | `(octExpr| $a:octExpr / $b:octExpr) => do `(Expr.binop .div $(← convExpr a) $(← convExpr b)) - | `(octExpr| $a:octExpr .* $b:octExpr) => do `(Expr.binop .emul $(← convExpr a) $(← convExpr b)) - | `(octExpr| $a:octExpr ./ $b:octExpr) => do `(Expr.binop .ediv $(← convExpr a) $(← convExpr b)) - | `(octExpr| $a:octExpr + $b:octExpr) => do `(Expr.binop .add $(← convExpr a) $(← convExpr b)) - | `(octExpr| $a:octExpr - $b:octExpr) => do `(Expr.binop .sub $(← convExpr a) $(← convExpr b)) - | `(octExpr| $a:octExpr == $b:octExpr) => do `(Expr.binop .eq $(← convExpr a) $(← convExpr b)) - | `(octExpr| $a:octExpr != $b:octExpr) => do `(Expr.binop .ne $(← convExpr a) $(← convExpr b)) - | `(octExpr| $a:octExpr < $b:octExpr) => do `(Expr.binop .lt $(← convExpr a) $(← convExpr b)) - | `(octExpr| $a:octExpr <= $b:octExpr) => do `(Expr.binop .le $(← convExpr a) $(← convExpr b)) - | `(octExpr| $a:octExpr > $b:octExpr) => do `(Expr.binop .gt $(← convExpr a) $(← convExpr b)) - | `(octExpr| $a:octExpr >= $b:octExpr) => do `(Expr.binop .ge $(← convExpr a) $(← convExpr b)) - | `(octExpr| $a:octExpr && $b:octExpr) => do `(Expr.binop .land $(← convExpr a) $(← convExpr b)) - | `(octExpr| $a:octExpr || $b:octExpr) => do `(Expr.binop .lor $(← convExpr a) $(← convExpr b)) - | `(octExpr| $a:octExpr & $b:octExpr) => do `(Expr.binop .band $(← convExpr a) $(← convExpr b)) - | `(octExpr| $a:octExpr | $b:octExpr) => do `(Expr.binop .bor $(← convExpr a) $(← convExpr b)) - -- Range: a:b or (a:step):b → three-part - | `(octExpr| $lo:octExpr : $hi:octExpr) => do + | name => `(Expr.ident $(Lean.quote name)) + -- Grouped + | `(octExpr| ($x)) => convExpr x + -- Unary + | `(octExpr| - $x) => do `(Expr.unop .neg $(← convExpr x)) + | `(octExpr| ! $x) => do `(Expr.unop .lnot $(← convExpr x)) + -- Power + | `(octExpr| $a ^ $b) => do `(Expr.binop .pow $(← convExpr a) $(← convExpr b)) + | `(octExpr| $a .^ $b) => do `(Expr.binop .epow $(← convExpr a) $(← convExpr b)) + -- Mul/Div + | `(octExpr| $a * $b) => do `(Expr.binop .mul $(← convExpr a) $(← convExpr b)) + | `(octExpr| $a / $b) => do `(Expr.binop .div $(← convExpr a) $(← convExpr b)) + | `(octExpr| $a .* $b) => do `(Expr.binop .emul $(← convExpr a) $(← convExpr b)) + | `(octExpr| $a ./ $b) => do `(Expr.binop .ediv $(← convExpr a) $(← convExpr b)) + -- Add/Sub + | `(octExpr| $a + $b) => do `(Expr.binop .add $(← convExpr a) $(← convExpr b)) + | `(octExpr| $a - $b) => do `(Expr.binop .sub $(← convExpr a) $(← convExpr b)) + -- Range (a:b — step defaults to 1; nested a:s:b parses as (a:s):b) + | `(octExpr| $lo : $hi) => do match lo with - | `(octExpr| $a:octExpr : $step:octExpr) => + | `(octExpr| $a : $step) => `(Expr.range $(← convExpr a) (some $(← convExpr step)) $(← convExpr hi)) | _ => `(Expr.range $(← convExpr lo) none $(← convExpr hi)) - -- Call / index (ident-based) - | `(octExpr| $f:ident ($args,*)) => do - let fT ← `(Expr.ident $(identStr f)) - let aTs ← args.getElems.mapM fun a => do - let t ← convExpr a; `(Arg.pos $t) + -- Comparison + | `(octExpr| $a == $b) => do `(Expr.binop .eq $(← convExpr a) $(← convExpr b)) + | `(octExpr| $a != $b) => do `(Expr.binop .ne $(← convExpr a) $(← convExpr b)) + | `(octExpr| $a <= $b) => do `(Expr.binop .le $(← convExpr a) $(← convExpr b)) + | `(octExpr| $a >= $b) => do `(Expr.binop .ge $(← convExpr a) $(← convExpr b)) + | `(octExpr| $a < $b) => do `(Expr.binop .lt $(← convExpr a) $(← convExpr b)) + | `(octExpr| $a > $b) => do `(Expr.binop .gt $(← convExpr a) $(← convExpr b)) + -- Logical + | `(octExpr| $a && $b) => do `(Expr.binop .land $(← convExpr a) $(← convExpr b)) + | `(octExpr| $a || $b) => do `(Expr.binop .lor $(← convExpr a) $(← convExpr b)) + | `(octExpr| $a & $b) => do `(Expr.binop .band $(← convExpr a) $(← convExpr b)) + | `(octExpr| $a | $b) => do `(Expr.binop .bor $(← convExpr a) $(← convExpr b)) + -- Function call / matrix index + | `(octExpr| $f:octExpr ( $args:octExpr,* )) => do + let fT ← convExpr f + let aTs ← args.getElems.mapM (fun a => do `(Arg.pos $(← convExpr a))) `(Expr.index $fT #[$aTs,*]) - -- Struct field: s.field - | `(octExpr| $s:octExpr.$field:ident) => do - `(Expr.dotIndex $(← convExpr s) $(strTerm field.getId.toString)) - -- Dynamic field: s.(expr) — ident base only - | `(octExpr| $s:ident .($field:octExpr)) => do - `(Expr.dynField (Expr.ident $(identStr s)) $(← convExpr field)) -- Function handles - | `(octExpr| @$id:ident) => - `(Expr.fnHandle $(strTerm id.getId.toString)) - | `(octExpr| @($params,*) $body:octExpr) => do - let ps := params.getElems.map identStr - `(Expr.anon #[$ps,*] $(← convExpr body)) - -- Vector / matrix - | `(octExpr| [$elems,*]) => do - let es := elems.getElems - if es.isEmpty then - return ← `(Expr.matrix #[]) - -- If first element is also [...], treat as multi-row matrix - let firstIsRow : Bool := match es[0]! with - | `(octExpr| [$_,*]) => true | _ => false - if firstIsRow then - let rowTerms ← es.mapM fun row => do - match row with - | `(octExpr| [$cols,*]) => do - let colTs ← cols.getElems.mapM convExpr - `(#[$colTs,*]) - | _ => Macro.throwError s!"expected [...] row in matrix literal, got: {row}" - `(Expr.matrix #[$rowTerms,*]) - else - let colTs ← es.mapM convExpr - `(Expr.matrix #[#[$colTs,*]]) - | e => Macro.throwError s!"unsupported octExpr: {e}" + | `(octExpr| @ $id:ident) => + `(Expr.fnHandle $(Lean.quote id.getId.toString)) + | `(octExpr| @( $params:ident,* ) $body:octExpr) => do + let pNames := params.getElems.map (fun p => quoteStr p.getId.toString) + `(Expr.anon #[$pNames,*] $(← convExpr body)) + -- Matrix literal: empty + | `(octExpr| [ ]) => `(Expr.matrix #[]) + -- Matrix literal: with body (one or more rows) + | `(octExpr| [ $body:octMatBody ]) => do + let rowTerms ← collectRows body + `(Expr.matrix #[$rowTerms,*]) + | _ => Macro.throwErrorAt e "unsupported expression syntax" +where + convRow (row : Syntax) : MacroM (TSyntax `term) := do + match row with + | `(octRow| $cols:octExpr,*) => do + let colTerms ← cols.getElems.mapM convExpr + `(#[$colTerms,*]) + | _ => Macro.throwErrorAt row "bad matrix row" + collectRows (body : Syntax) : MacroM (Array (TSyntax `term)) := do + match body with + | `(octMatBody| $r:octRow) => do + return #[← convRow r] + | `(octMatBody| $r:octRow ; $rest:octMatBody) => do + let rt ← convRow r + let restRows ← collectRows rest + return #[rt] ++ restRows + | _ => Macro.throwErrorAt body "bad matrix body" -- ───────────────────────────────────────────────────────────────── --- convStmt : octStmt syntax → term of type OctiveLean.Stmt +-- Macro conversion: octStmt → Term (of type OctiveLean.Stmt) -- ───────────────────────────────────────────────────────────────── -private partial def convStmt : TSyntax `octStmt → MacroM (TSyntax `term) - -- Expression statement - | `(octStmt| $e:octExpr) => do `(Stmt.exprS $(← convExpr e) false) - | `(octStmt| $e:octExpr ;) => do `(Stmt.exprS $(← convExpr e) true) - -- Assignment - | `(octStmt| $x:ident = $e:octExpr) => do - `(Stmt.assign #[$(identStr x)] $(← convExpr e) false) - | `(octStmt| $x:ident = $e:octExpr ;) => do - `(Stmt.assign #[$(identStr x)] $(← convExpr e) true) - -- Struct field assignment: s.f = expr - | `(octStmt| $s:ident.$f:ident = $e:octExpr ;) => do - `(Stmt.indexAssign (Expr.dotIndex (Expr.ident $(identStr s)) $(strTerm f.getId.toString)) $(← convExpr e) true) - | `(octStmt| $s:ident.$f:ident = $e:octExpr) => do - `(Stmt.indexAssign (Expr.dotIndex (Expr.ident $(identStr s)) $(strTerm f.getId.toString)) $(← convExpr e) false) +private partial def convStmt (s : Syntax) : MacroM (TSyntax `term) := do + match s with + -- Expression statements + | `(octStmt| $e:octExpr ;) => do `(Stmt.exprS $(← convExpr e) true) + | `(octStmt| $e:octExpr) => do `(Stmt.exprS $(← convExpr e) false) + -- Assignments + | `(octStmt| $x:ident = $e:octExpr ;) => + do `(Stmt.assign #[$(Lean.quote x.getId.toString)] $(← convExpr e) true) + | `(octStmt| $x:ident = $e:octExpr) => + do `(Stmt.assign #[$(Lean.quote x.getId.toString)] $(← convExpr e) false) -- Multi-assignment - | `(octStmt| [$xs,*] = $e:octExpr) => do - let names := xs.getElems.map identStr - `(Stmt.assign #[$names,*] $(← convExpr e) false) - | `(octStmt| [$xs,*] = $e:octExpr ;) => do - let names := xs.getElems.map identStr + | `(octStmt| [ $xs:ident,* ] = $e:octExpr ;) => do + let names := xs.getElems.map (fun x => quoteStr x.getId.toString) `(Stmt.assign #[$names,*] $(← convExpr e) true) + | `(octStmt| [ $xs:ident,* ] = $e:octExpr) => do + let names := xs.getElems.map (fun x => quoteStr x.getId.toString) + `(Stmt.assign #[$names,*] $(← convExpr e) false) -- IF | `(octStmt| if $cond:octExpr $thenB:octStmt* $[elseif $eiconds:octExpr $eibodies:octStmt*]* $[else $elseB:octStmt*]? endif) => do let condT ← convExpr cond - let thenBT ← thenB.mapM convStmt - let eiBranches ← (Array.zip eiconds eibodies).mapM fun (c, body) => do - let cT ← convExpr c - let bodyT ← body.mapM convStmt - `(($cT, #[$bodyT,*])) - let elseBT ← match elseB with - | none => `(none) + let thenT ← thenB.mapM convStmt + let eiTs ← (Array.zip eiconds eibodies).mapM (fun (c, body) => do + let ct ← convExpr c + let bt ← body.mapM convStmt + `(($ct, #[$bt,*]))) + let elseT ← match elseB with + | none => `((none : Option (Array Stmt))) | some b => do let bt ← b.mapM convStmt; `(some #[$bt,*]) - `(Stmt.ifS $condT #[$thenBT,*] #[$eiBranches,*] $elseBT) + `(Stmt.ifS $condT #[$thenT,*] #[$eiTs,*] $elseT) -- FOR | `(octStmt| for $k:ident = $range:octExpr $body:octStmt* endfor) => do - let bodyT ← body.mapM convStmt - `(Stmt.forS $(identStr k) $(← convExpr range) #[$bodyT,*]) + `(Stmt.forS $(Lean.quote k.getId.toString) + $(← convExpr range) + #[$(← body.mapM convStmt),*]) -- WHILE | `(octStmt| while $cond:octExpr $body:octStmt* endwhile) => do - let bodyT ← body.mapM convStmt - `(Stmt.whileS $(← convExpr cond) #[$bodyT,*]) + `(Stmt.whileS $(← convExpr cond) #[$(← body.mapM convStmt),*]) -- SWITCH | `(octStmt| switch $val:octExpr - $[case $caseVals:octExpr $caseBodies:octStmt*]* - $[otherwise $otherwiseB:octStmt*]? + $[case $cvs:octExpr $cbs:octStmt*]* + $[otherwise $ob:octStmt*]? endswitch) => do - let valT ← convExpr val - let branches ← (Array.zip caseVals caseBodies).mapM fun (cv, cb) => do - let cvT ← convExpr cv - let cbT ← cb.mapM convStmt - `(($cvT, #[$cbT,*])) - let otherwiseT ← match otherwiseB with - | none => `(none) + let valT ← convExpr val + let brs ← (Array.zip cvs cbs).mapM (fun (cv, cb) => do + let cvt ← convExpr cv + let cbt ← cb.mapM convStmt + `(($cvt, #[$cbt,*]))) + let otT ← match ob with + | none => `((none : Option (Array Stmt))) | some b => do let bt ← b.mapM convStmt; `(some #[$bt,*]) - `(Stmt.switchS $valT #[$branches,*] $otherwiseT) - -- TRY + `(Stmt.switchS $valT #[$brs,*] $otT) + -- TRY / CATCH | `(octStmt| try $tryB:octStmt* $[catch $evar:ident $catchB:octStmt*]? - endtry) => do - let tryBT ← tryB.mapM convStmt - let catchT ← match evar, catchB with + end_try_catch) => do + let tryT ← tryB.mapM convStmt + let catchT ← + match evar, catchB with | some ev, some cb => do let cbt ← cb.mapM convStmt - `(some ($(identStr ev), #[$cbt,*])) - | _, _ => `(none) - `(Stmt.tryS #[$tryBT,*] $catchT) + `(some ($(Lean.quote ev.getId.toString), #[$cbt,*])) + | _, _ => `((none : Option (String × Array Stmt))) + `(Stmt.tryS #[$tryT,*] $catchT) -- Control flow - | `(octStmt| return) => `(Stmt.returnS) - | `(octStmt| break) => `(Stmt.breakS) - | `(octStmt| continue) => `(Stmt.continueS) + | `(octStmt| return) => `(Stmt.returnS) + | `(octStmt| break) => `(Stmt.breakS) + | `(octStmt| continue) => `(Stmt.continueS) -- Scope - | `(octStmt| global $ids,*) => do - let names := ids.getElems.map identStr + | `(octStmt| global $ids*) => do + let names := ids.map (fun i => quoteStr i.getId.toString) `(Stmt.globalS #[$names,*]) - | `(octStmt| clear $ids,*) => do - let names := ids.getElems.map identStr + | `(octStmt| clear $ids*) => do + let names := ids.map (fun i => quoteStr i.getId.toString) `(Stmt.clearS #[$names,*]) - -- Function: single return - | `(octStmt| function $ret:ident = $name:ident ($params,*) $body:octStmt* endfunction) => do - let pns := params.getElems.map identStr - let bodyT ← body.mapM convStmt - `(Stmt.funcDefS (FuncDef.mk $(identStr name) #[$pns,*] - #[$(identStr ret)] #[$bodyT,*])) - -- Function: multi-return - | `(octStmt| function [$rets,*] = $name:ident ($params,*) $body:octStmt* endfunction) => do - let pns := params.getElems.map identStr - let rns := rets.getElems.map identStr - let bodyT ← body.mapM convStmt - `(Stmt.funcDefS (FuncDef.mk $(identStr name) #[$pns,*] #[$rns,*] #[$bodyT,*])) - -- Function: no return - | `(octStmt| function $name:ident ($params,*) $body:octStmt* endfunction) => do - let pns := params.getElems.map identStr - let bodyT ← body.mapM convStmt - `(Stmt.funcDefS (FuncDef.mk $(identStr name) #[$pns,*] #[] #[$bodyT,*])) - | s => Macro.throwError s!"unsupported octStmt: {s}" + -- Function defs + | `(octStmt| function $ret:ident = $name:ident ( $params:ident,* ) + $body:octStmt* endfunction) => do + let pNames := params.getElems.map (fun p => quoteStr p.getId.toString) + let bt ← body.mapM convStmt + `(Stmt.funcDefS (FuncDef.mk + $(quoteStr name.getId.toString) + #[$pNames,*] + #[$(quoteStr ret.getId.toString)] + #[$bt,*])) + | `(octStmt| function [ $rets:ident,* ] = $name:ident ( $params:ident,* ) + $body:octStmt* endfunction) => do + let pNames := params.getElems.map (fun p => quoteStr p.getId.toString) + let rNames := rets.getElems.map (fun r => quoteStr r.getId.toString) + let bt ← body.mapM convStmt + `(Stmt.funcDefS (FuncDef.mk + $(quoteStr name.getId.toString) + #[$pNames,*] + #[$rNames,*] + #[$bt,*])) + | `(octStmt| function $name:ident ( $params:ident,* ) + $body:octStmt* endfunction) => do + let pNames := params.getElems.map (fun p => quoteStr p.getId.toString) + let bt ← body.mapM convStmt + `(Stmt.funcDefS (FuncDef.mk + $(quoteStr name.getId.toString) + #[$pNames,*] + #[] + #[$bt,*])) + | _ => Macro.throwErrorAt s "unsupported statement syntax" -- ───────────────────────────────────────────────────────────────── --- Helpers to mark expanded syntax as canonical --- (Macro-generated syntax has SourceInfo.synthetic canonical:=false, --- so savePanelWidgetInfo can't find the position. We flip the flag.) +-- Source info helper: macro-generated syntax has canonical := false, +-- which prevents `savePanelWidgetInfo` from binding the widget to a +-- source position. Flip the flag. -- ───────────────────────────────────────────────────────────────── -private def mkCanonicalInfo : SourceInfo → SourceInfo +private def mkCanonicalInfo : Lean.SourceInfo → Lean.SourceInfo | .synthetic s e _ => .synthetic s e true | si => si -private def mkCanonicalSyntax : Syntax → Syntax - | .node i k a => .node (mkCanonicalInfo i) k a - | .atom i v => .atom (mkCanonicalInfo i) v +private def mkCanonicalSyntax : Lean.Syntax → Lean.Syntax + | .node i k a => .node (mkCanonicalInfo i) k a + | .atom i v => .atom (mkCanonicalInfo i) v | .ident i r v p => .ident (mkCanonicalInfo i) r v p - | s => s + | s => s -- ───────────────────────────────────────────────────────────────── --- Commands +-- Top-level commands -- ───────────────────────────────────────────────────────────────── +/-- `octave! { stmts }` — parse, type-check, and run the block. -/ +syntax (name := octaveRun) "octave!" "{" octStmt* "}" : command + macro_rules - | `(octave! $stmts:octStmt* octave_end) => do + | `(command| octave! { $stmts:octStmt* }) => do let stmtTerms ← stmts.mapM convStmt - let result : TSyntax `command ← `(#html (show IO ProofWidgets.Html from do + let result : Lean.TSyntax `command ← `(#html (show IO ProofWidgets.Html from do let plotBuf ← IO.mkRef (#[] : Array OctiveLean.Figure) let env := OctiveLean.PlotBuiltins.register plotBuf (OctiveLean.registerAllBuiltins OctiveLean.Env.empty) @@ -387,9 +407,14 @@ macro_rules | .error e => IO.eprintln s!"runtime error: {e}" let figs ← plotBuf.get return OctiveLean.PlotWidget.render figs)) - return (⟨mkCanonicalSyntax result.raw⟩ : TSyntax `command) + return (⟨mkCanonicalSyntax result.raw⟩ : Lean.TSyntax `command) + +/-- `octave_program! name { stmts }` — bind the parsed AST to a Lean def. -/ +syntax (name := octaveProg) "octave_program!" ident "{" octStmt* "}" : command macro_rules - | `(octave_stmts! $name:ident $stmts:octStmt* octave_end) => do + | `(command| octave_program! $name:ident { $stmts:octStmt* }) => do let stmtTerms ← stmts.mapM convStmt `(def $name : Array OctiveLean.Stmt := #[$stmtTerms,*]) + +end OctiveLean.DSL diff --git a/OctiveLean/Eval.lean b/OctiveLean/Eval.lean index a390312..16e5d84 100644 --- a/OctiveLean/Eval.lean +++ b/OctiveLean/Eval.lean @@ -2,6 +2,7 @@ import OctiveLean.Value import OctiveLean.Env import OctiveLean.Error import OctiveLean.AST +import OctiveLean.SymPyBridge namespace OctiveLean @@ -89,8 +90,29 @@ private def matMul (r1 c1 : Nat) (d1 : Array Float) o return .matrix r1 c2 out -private def evalBinOp (op : BinOp) (lv rv : Value) : EvalM Value := - match op with +private def evalSymBinOp (op : BinOp) (lv rv : Value) : EvalM Value := do + let l ← liftM (m := IO) (SymPyBridge.toSympy lv) + let r ← liftM (m := IO) (SymPyBridge.toSympy rv) + let py := match op with + | .add => s!"({l}) + ({r})" | .sub => s!"({l}) - ({r})" + | .mul | .emul => s!"({l}) * ({r})" + | .div | .ediv => s!"({l}) / ({r})" + | .ldiv | .eldiv => s!"({r}) / ({l})" + | .pow | .epow => s!"({l}) ** ({r})" + | .lt => s!"Lt({l}, {r})" | .le => s!"Le({l}, {r})" + | .gt => s!"Gt({l}, {r})" | .ge => s!"Ge({l}, {r})" + | .eq => s!"Eq({l}, {r})" | .ne => s!"Ne({l}, {r})" + | .land | .band => s!"And({l}, {r})" + | .lor | .bor => s!"Or({l}, {r})" + liftM (m := IO) (SymPyBridge.emit py) + +private def isSym : Value → Bool + | .sym _ _ => true + | _ => false + +private def evalBinOp (op : BinOp) (lv rv : Value) : EvalM Value := do + if isSym lv || isSym rv then evalSymBinOp op lv rv + else match op with | .add => ewiseOp (· + ·) lv rv | .sub => ewiseOp (· - ·) lv rv | .emul => ewiseOp (· * ·) lv rv @@ -269,6 +291,9 @@ mutual | .scalar f => return .scalar (-f) | .matrix r c d => return .matrix r c (d.map (- ·)) | .integer iv => return .scalar (-iv.toFloat) + | .sym _ _ => + let s ← liftM (m := IO) (SymPyBridge.toSympy v) + liftM (m := IO) (SymPyBridge.emit s!"-({s})") | other => throw (.typeError s!"cannot negate {other.typeName}") | .uplus => return v | .lnot => diff --git a/OctiveLean/SymPyBridge.lean b/OctiveLean/SymPyBridge.lean new file mode 100644 index 0000000..51bd5d1 --- /dev/null +++ b/OctiveLean/SymPyBridge.lean @@ -0,0 +1,158 @@ +import OctiveLean.Value + +namespace OctiveLean.SymPyBridge + +/-! Persistent SymPy subprocess. + +Mirrors the architecture of GNU Octave's `symbolic` package +(`pycall_sympy__.m`): keep one Python process alive across calls and +exchange SymPy expressions over stdin/stdout using line-based sentinels. + +Each `Value.sym` carries the SymPy `srepr` (round-trips back into Python +verbatim) and the human-readable `str(...)` form (for `disp`). +-/ + +private def pythonScript : String := r#" +import sys, sympy +from sympy import * +from sympy.parsing.sympy_parser import parse_expr, standard_transformations, convert_xor, implicit_multiplication + +ns = {} +exec("from sympy import *", ns) + +_TRANSFORMS = standard_transformations + (convert_xor, implicit_multiplication) + +def _parse(s): + return parse_expr(s, transformations=_TRANSFORMS, local_dict=ns) +ns['_parse'] = _parse + +def _emit(x): + sys.stdout.write("~~~SREPR~~~\n") + try: sys.stdout.write(srepr(x) + "\n") + except Exception: sys.stdout.write(repr(x) + "\n") + sys.stdout.write("~~~PRETTY~~~\n") + try: sys.stdout.write(str(x) + "\n") + except Exception: sys.stdout.write(repr(x) + "\n") +ns['_emit'] = _emit + +EOC = "~~~EOC~~~" +EOR = "~~~EOR~~~" +ERR = "~~~ERR~~~" + +buf = [] +for raw in iter(sys.stdin.readline, ""): + line = raw.rstrip("\n") + if line == EOC: + code = "\n".join(buf) + buf = [] + try: + exec(code, ns) + except Exception as e: + sys.stdout.write(ERR + type(e).__name__ + ": " + str(e) + "\n") + sys.stdout.write(EOR + "\n") + sys.stdout.flush() + else: + buf.append(line) +"# + +private def cfg : IO.Process.StdioConfig := + { stdin := .piped, stdout := .piped, stderr := .piped } + +private structure Bridge where + proc : IO.Process.Child cfg + +initialize bridgeRef : IO.Ref (Option Bridge) ← IO.mkRef none + +private def spawn : IO Bridge := do + let proc ← IO.Process.spawn + { cmd := "python3" + args := #["-u", "-c", pythonScript] + stdin := .piped + stdout := .piped + stderr := .piped } + return { proc } + +private def getOrInit : IO Bridge := do + match (← bridgeRef.get) with + | some b => return b + | none => + let b ← spawn + bridgeRef.set (some b) + return b + +private def stripTrailingNL (s : String) : String := + if s.endsWith "\n" then (s.dropEnd 1).toString else s + +/-- Send a Python `code` block to the subprocess and read everything it writes + until the EOR sentinel. Returns the captured stdout (without sentinel) or + an error message if Python raised. -/ +def runRaw (code : String) : IO (Except String String) := do + let bridge ← getOrInit + for line in code.splitOn "\n" do + bridge.proc.stdin.putStrLn line + bridge.proc.stdin.putStrLn "~~~EOC~~~" + bridge.proc.stdin.flush + let mut collected : Array String := #[] + let mut error? : Option String := none + let mut done := false + while not done do + let raw ← bridge.proc.stdout.getLine + let line := stripTrailingNL raw + if line == "~~~EOR~~~" then + done := true + else if line.startsWith "~~~ERR~~~" then + error? := some ((line.drop "~~~ERR~~~".length).toString) + else + collected := collected.push line + match error? with + | some e => return .error e + | none => return .ok (String.intercalate "\n" collected.toList) + +private def parseEmitOutput (out : String) : Except String (String × String) := + let needle := "~~~PRETTY~~~\n" + match (out.splitOn needle) with + | [head, tail] => + let srMarker := "~~~SREPR~~~\n" + if head.startsWith srMarker then + let srepr := stripTrailingNL ((head.drop srMarker.length).toString) + let pretty := stripTrailingNL tail + .ok (srepr, pretty) + else .error s!"sympy: missing SREPR marker in {out}" + | _ => .error s!"sympy: missing PRETTY marker in {out}" + +/-- Evaluate a Python expression that produces a SymPy object and wrap the + result as a `Value.sym`. Caller is responsible for building a syntactically + valid Python expression (use `Value.toSympyExpr` for operands). -/ +def emit (expr : String) : IO Value := do + match (← runRaw s!"_emit({expr})") with + | .error e => throw (IO.userError s!"sympy: {e}") + | .ok out => + match parseEmitOutput out with + | .ok (sr, pr) => return .sym sr pr + | .error e => throw (IO.userError e) + +/-- Convert any OctiveLean Value into a Python expression string suitable for + splicing into SymPy commands. Numerics become SymPy `Integer`/`Float`, + Sym values forward their `srepr`, strings get parsed with `_parse`. -/ +def toSympy : Value → IO String + | .sym sr _ => return sr + | .scalar f => + if f == f.floor && f.abs < 1e15 then return s!"Integer({f.toInt64})" + else return s!"Float({f})" + | .fscalar f => return s!"Float({f})" + | .integer iv => return s!"Integer({iv.display})" + | .boolean b => return if b then "true" else "false" + | .string s => + let escaped := s.replace "\\" "\\\\" |>.replace "'" "\\'" + return s!"_parse('{escaped}')" + | v => throw (IO.userError s!"sympy: cannot convert {v.typeName} to symbolic") + +/-- Coerce a Value to a Sym, by parsing through SymPy if it is not already one. -/ +def asSym (v : Value) : IO Value := do + match v with + | .sym _ _ => return v + | _ => + let py ← toSympy v + emit py + +end OctiveLean.SymPyBridge diff --git a/OctiveLean/Value.lean b/OctiveLean/Value.lean index 08ee38d..efd58f1 100644 --- a/OctiveLean/Value.lean +++ b/OctiveLean/Value.lean @@ -53,6 +53,7 @@ mutual | struct : Array (String × Value) → Value | fn : FuncVal → Value | range : Float → Float → Float → Value -- start step stop (lazy) + | sym : (srepr : String) → (pretty : String) → Value -- SymPy expression (srepr round-trips, pretty for display) | empty : Value -- [] /-- A callable function value -/ @@ -98,6 +99,7 @@ def Value.typeName : Value → String | .struct _ => "struct" | .fn _ => "function_handle" | .range _ _ _ => "range" + | .sym _ _ => "sym" | .empty => "[]" /-! Utility functions -/ @@ -141,6 +143,7 @@ def Value.shape : Value → Nat × Nat | .struct _ => (1, 1) | .fn _ => (1, 1) | .range s st e => (1, (Value.rangeToArray s st e).size) + | .sym _ _ => (1, 1) | .empty => (0, 0) /-- Format a Float as Octave does: no trailing .0 for integers, reasonable precision -/ @@ -203,6 +206,7 @@ def Value.display (name : String) : Value → String else let elems := data.toList.map formatFloat |> String.intercalate " " s!"{name} =\n\n {elems}\n" + | .sym _ p => s!"{name} = (sym) {p}" | .empty => s!"{name} = [](0x0)" | .cmatrix r c _ => s!"{name} = <{r}x{c} complex matrix>" @@ -227,6 +231,7 @@ def Value.printStr : Value → String String.intercalate "" elems s!"\n{String.intercalate "\n" rows}\n" | .string s => s + | .sym _ p => p | v => v.display "" end OctiveLean diff --git a/PlotDemo.lean b/PlotDemo.lean index 4af733b..d1485ca 100644 --- a/PlotDemo.lean +++ b/PlotDemo.lean @@ -1,45 +1,45 @@ import OctiveLean --- Hover over each octave! block in the infoview to see the rendered chart. +-- Hover over each `octave!` block to see the rendered chart in the infoview. -- Line plot of a sine wave -octave! +octave! { x = linspace(0, 6.28, 64) y = sin(x) plot(x, y) title("Sine Wave") xlabel("x") ylabel("sin(x)") -octave_end +} -- Scatter plot -octave! +octave! { x = linspace(-3, 3, 40) y = x .* x scatter(x, y) title("Parabola") -octave_end +} -- Bar chart -octave! +octave! { bar([1, 2, 3, 4, 5], [3.2, 1.8, 4.5, 2.1, 3.9]) title("Bar Chart") xlabel("Category") ylabel("Value") -octave_end +} -- Histogram of residuals from a sine wave -octave! +octave! { x = linspace(0, 6.28, 200) y = sin(x) .* cos(x) hist(y, 20) title("Histogram of sin(x)*cos(x)") xlabel("Value") ylabel("Count") -octave_end +} -- Multi-series with hold_on / legend -octave! +octave! { x = linspace(0, 6.28, 64) hold_on() plot(x, sin(x)) @@ -47,17 +47,17 @@ octave! hold_off() legend("sin", "cos") title("Trig Functions") -octave_end +} -- Stem plot -octave! +octave! { x = linspace(0, 3.14, 16) stem(x, sin(x)) title("Stem Plot") -octave_end +} -- ── 3-D: plot3 (helix) ─────────────────────────────────────────── -octave! +octave! { t = linspace(0, 12.57, 80) xs = cos(t) ys = sin(t) @@ -67,17 +67,17 @@ octave! xlabel("cos t") ylabel("sin t") zlabel("t/2") -octave_end +} -- ── 3-D: scatter3 ──────────────────────────────────────────────── -octave! +octave! { t = linspace(0, 6.28, 60) scatter3(cos(t), sin(t), t) title("Circular Scatter3") -octave_end +} -- ── 3-D: surf (corrugated wave) ────────────────────────────────── -octave! +octave! { x = linspace(0, 6.28, 24) y = linspace(0, 3, 12) surf(x, y, sin(x)) @@ -85,22 +85,22 @@ octave! xlabel("x") ylabel("y") zlabel("z") -octave_end +} -- ── 3-D: waterfall ─────────────────────────────────────────────── -octave! +octave! { x = linspace(0, 6.28, 30) y = linspace(0, 3, 8) waterfall(x, y, sin(x)) title("Waterfall") -octave_end +} -- ── 3-D: contourf ──────────────────────────────────────────────── -octave! +octave! { x = linspace(-3, 3, 30) y = linspace(-3, 3, 30) contourf(x, y, sin(x)) title("Contour: sin(x)") xlabel("x") ylabel("y") -octave_end +} diff --git a/RosettaStone.lean b/RosettaStone.lean index 6e507be..792942d 100644 --- a/RosettaStone.lean +++ b/RosettaStone.lean @@ -1,407 +1,165 @@ import OctiveLean +import OctiveLean.DSL /-! -# OctiveLean Rosetta Stone (DSL edition) +# OctiveLean Rosetta Stone — DSL edition -Octave code written directly as Lean syntax — no strings, no raw AST. -The `octave! ... octave_end` macro compiles to typed `OctiveLean.Stmt` -values at elaboration time, so the LSP highlights keywords, operators, -and structure just like any other Lean code. +Octave is now a first-class Lean 4 syntax category. The LSP recognizes +keywords, operators, and structure inside `octave! { ... }` blocks. -Block-closer differences from standard Octave (all are valid in real Octave too): - `endif` `endfor` `endwhile` `endfunction` `endswitch` `endtry` - -Outer block: `octave! ... octave_end` +Syntax differences from standard Octave: + • Outer block: `octave! { ... }` + • Block terminators: `endif` / `endfor` / `endwhile` / `endswitch` / + `end_try_catch` / `endfunction` (Octave-valid keywords) + • Strings: `"..."` (Lean style) + • Comments: `--` (Lean style — `%` is the modulo operator token) + • Matrices: `[1.0, 2.0; 3.0, 4.0]` (commas for cols, `;` for rows) -/ --- ───────────────────────────────────────────────────────────────── --- §1 LITERALS --- ───────────────────────────────────────────────────────────────── +open OctiveLean DSL -octave! +-- §1 LITERALS +octave! { disp(3.14) disp(42) disp("hello") disp(true) - disp(false) -octave_end +} --- ───────────────────────────────────────────────────────────────── --- §2 VARIABLES — assignment and lookup --- ───────────────────────────────────────────────────────────────── - --- Semicolon = silent; no semicolon = echoes the value -octave! +-- §2 ASSIGNMENT +octave! { x = 42; disp(x) -octave_end +} -octave! - a = 10 - b = 20; +-- §3 ARITHMETIC +octave! { + a = 10; + b = 3; disp(a + b) -octave_end + disp(a - b) + disp(a * b) + disp(a / b) + disp(a ^ b) + disp(a .* b) + disp(a ./ b) + disp(a .^ b) +} --- ───────────────────────────────────────────────────────────────── --- §3 ARITHMETIC OPERATORS --- ───────────────────────────────────────────────────────────────── - -octave! - a = 10; b = 3; - disp(a + b) -- 13 - disp(a - b) -- 7 - disp(a * b) -- 30 - disp(a / b) -- 3.333… - disp(a ^ b) -- 1000 - disp(a .* b) -- 30 element-wise - disp(a ./ b) -- 3.333… - disp(a .^ b) -- 1000 -octave_end - --- ───────────────────────────────────────────────────────────────── -- §4 COMPARISON & LOGICAL --- ───────────────────────────────────────────────────────────────── +octave! { + disp(3 < 5) + disp(3 <= 3) + disp(3 == 3) + disp(3 != 4) + disp(1 && 0) + disp(1 || 0) +} -octave! - disp(3 < 5) -- 1 - disp(3 <= 3) -- 1 - disp(5 > 3) -- 1 - disp(5 >= 6) -- 0 - disp(3 == 3) -- 1 - disp(3 != 4) -- 1 - disp(1 && 0) -- 0 short-circuit AND - disp(1 || 0) -- 1 short-circuit OR - disp(1 & 0) -- 0 element-wise AND - disp(1 | 0) -- 1 element-wise OR -octave_end +-- §5 UNARY +octave! { + disp(- 5) + disp(! true) +} --- ───────────────────────────────────────────────────────────────── --- §5 UNARY OPERATORS --- ───────────────────────────────────────────────────────────────── - -octave! - disp(-5) -- negation - disp(!true) -- logical not → 0 - v = [1.0, 2.0, 3.0]; - disp(v) -octave_end - --- ───────────────────────────────────────────────────────────────── -- §6 MATRIX LITERALS --- [a, b, c] row vector --- [[a, b], [c, d]] matrix (rows are inner arrays) --- ───────────────────────────────────────────────────────────────── - -octave! - row = [1.0, 2.0, 3.0, 4.0, 5.0] - M = [[1.0, 2.0, 3.0], [4.0, 5.0, 6.0], [7.0, 8.0, 9.0]] - eigenvalues(M) - E = [] +octave! { + row = [1, 2, 3, 4, 5]; + M = [1, 2, 3; 4, 5, 6; 7, 8, 9]; disp(size(M)) - disp([1,2,3]*M) -octave_end +} --- ───────────────────────────────────────────────────────────────── --- §7 CELL ARRAYS --- ───────────────────────────────────────────────────────────────── +-- §7 RANGES +octave! { + r = 1 : 5; + disp(length(r)) +} --- Note: cell array syntax uses the raw AST path for now; --- the `{ }` token is not yet wired in the DSL syntax category. --- See RosettaStone.lean (string edition) for the string-based version. - --- ───────────────────────────────────────────────────────────────── --- §8 RANGES a:b and a:step:b --- ───────────────────────────────────────────────────────────────── - -octave! - r1 = 1:5; -- 1 2 3 4 5 - r2 = 0.0:0.5:2.0; -- 0.0 0.5 1.0 1.5 2.0 (a:step:b via (a:step):b parse) - r3 = 5.0: -1.0 :1.0; -- 5 4 3 2 1 - disp(r1) - disp(length(r1)) -octave_end - --- ───────────────────────────────────────────────────────────────── --- §9 INDEXING A(i, j) --- ───────────────────────────────────────────────────────────────── - -octave! - A = [[10.0, 20.0, 30.0], [40.0, 50.0, 60.0]]; - disp(A(1, 2)) -- 20 - disp(A(2, 1)) -- 40 - disp(A(1, 3)) -- 30 -octave_end - --- ───────────────────────────────────────────────────────────────── --- §10 STRUCT FIELDS s.field and s.(expr) --- ───────────────────────────────────────────────────────────────── - -octave! - p.x = 3.0; - p.y = 4.0; - disp(p.x) -- 3 - disp(p.y) -- 4 -octave_end - --- Note: p.(field) dynamic field access works as a standalone statement, --- but not nested inside another call like disp(p.(field)) due to Lean's --- ".(" single-token ambiguity inside argument lists. - --- ───────────────────────────────────────────────────────────────── --- §11 FUNCTION HANDLES @name and @(args) expr --- ───────────────────────────────────────────────────────────────── - -octave! - f = @sin; - disp(f(3.14159./4)) -- 0 - - g = @(x) x .^ 2.0 + 1.0; - disp(g(3.0)) -- 10 - - h = @(x, y) x + y; - disp(h(10.0, 5.0)) -- 15 -octave_end - --- ───────────────────────────────────────────────────────────────── --- §12 IF / ELSEIF / ELSE / ENDIF --- ───────────────────────────────────────────────────────────────── - -octave! - x = 7.0; - if x > 10.0 +-- §8 IF / ELSEIF / ELSE +octave! { + x = 7; + if x > 10 disp("big") - elseif x > 5.0 + elseif x > 5 disp("medium") else disp("small") endif -octave_end +} --- ───────────────────────────────────────────────────────────────── --- §13 FOR / ENDFOR --- ───────────────────────────────────────────────────────────────── - -octave! - s = 0.0; - for k = 1:5 +-- §9 FOR LOOP +octave! { + s = 0; + for k = 1 : 5 s = s + k; endfor - disp(s) -- 15 -octave_end + disp(s) +} --- ───────────────────────────────────────────────────────────────── --- §14 WHILE / ENDWHILE --- ───────────────────────────────────────────────────────────────── - -octave! - n = 1.0; - while n < 32.0 - n = n * 2.0; +-- §10 WHILE LOOP +octave! { + n = 1; + while n < 32 + n = n * 2; endwhile - disp(n) -- 32 -octave_end + disp(n) +} --- ───────────────────────────────────────────────────────────────── --- §15 BREAK / CONTINUE --- ───────────────────────────────────────────────────────────────── +-- §11 FUNCTION DEFINITION +octave! { + function y = square(x) + y = x .^ 2; + endfunction + disp(square(7)) +} -octave! - for k = 1:10 - if k == 4.0 - break +-- §12 RECURSIVE FUNCTION (factorial) +octave! { + function y = fact(n) + if n <= 1 + y = 1; + else + y = n * fact(n - 1); endif - endfor - disp(k) -- 4 + endfunction + disp(fact(6)) +} - s = 0.0; - for k = 1:5 - if mod(k, 2.0) == 0.0 - continue - endif - s = s + k; - endfor - disp(s) -- 9 -octave_end - --- ───────────────────────────────────────────────────────────────── --- §16 SWITCH / CASE / OTHERWISE / ENDSWITCH --- ───────────────────────────────────────────────────────────────── - -octave! - day = "Mon"; - switch day - case "Mon" - disp("Monday") - case "Fri" - disp("Friday") - otherwise - disp("Other") - endswitch -octave_end - --- ───────────────────────────────────────────────────────────────── --- §17 TRY / CATCH / ENDTRY --- ───────────────────────────────────────────────────────────────── - -octave! +-- §13 TRY / CATCH +octave! { try disp(undefined_xyz) catch e disp("caught an error") - endtry -octave_end + end_try_catch +} --- ───────────────────────────────────────────────────────────────── --- §18 FUNCTION DEFINITION & CALL --- ───────────────────────────────────────────────────────────────── - -octave! - function y = square(x) - y = x .^ 2.0; - endfunction - - function z = add2(a, b) - z = a + b; - endfunction - - disp(square(7.0)) -- 49 - disp(add2(10.0, 32.0)) -- 42 -octave_end - --- ───────────────────────────────────────────────────────────────── --- §19 RECURSIVE FUNCTION (factorial) --- ───────────────────────────────────────────────────────────────── - -octave! - function y = fact(n) - if n <= 1.0 - y = 1.0; - else - y = n * fact(n - 1.0); - endif - endfunction - - disp(fact(6.0)) -- 720 -octave_end - --- ───────────────────────────────────────────────────────────────── --- §20 GLOBAL & CLEAR --- ───────────────────────────────────────────────────────────────── - -octave! - global G - G = 99.0 - disp(G) - clear G -octave_end - --- ───────────────────────────────────────────────────────────────── --- §21 MATRIX CONSTRUCTORS (builtins) --- ───────────────────────────────────────────────────────────────── - -octave! - disp(zeros(2.0, 3.0)) - disp(ones(3.0)) - disp(eye(3.0)) - disp(linspace(0.0, 1.0, 5.0)) -octave_end - --- ───────────────────────────────────────────────────────────────── --- §22 MATH BUILTINS --- ───────────────────────────────────────────────────────────────── - -octave! - disp(sqrt(2.0)) - disp(abs(-5.0)) - disp(exp(1.0)) - disp(log(exp(1.0))) +-- §14 BUILTINS — math +octave! { + disp(sqrt(2)) + disp(abs(- 5)) + disp(sin(0)) + disp(cos(0)) + disp(exp(1)) + disp(log(exp(1))) disp(floor(3.7)) disp(ceil(3.2)) - disp(round(3.5)) - disp(sin(0.0)) - disp(cos(0.0)) - disp(mod(17.0, 5.0)) - disp(max([3.0, 1.0, 5.0])) - disp(min([3.0, 1.0, 5.0])) - disp(sum([1.0, 2.0, 3.0, 4.0, 5.0])) - disp(prod([1.0, 2.0, 3.0, 4.0, 5.0])) - disp(mean([1.0, 2.0, 3.0, 4.0, 5.0])) - disp(norm([3.0, 4.0])) - disp(dot([1.0, 2.0], [3.0, 4.0])) -octave_end + disp(mod(17, 5)) + disp(max([3, 1, 4, 1, 5])) + disp(min([3, 1, 4, 1, 5])) + disp(sum([1, 2, 3, 4, 5])) + disp(mean([1, 2, 3, 4, 5])) + disp(norm([3, 4])) +} --- ───────────────────────────────────────────────────────────────── --- §23 STRING BUILTINS --- ───────────────────────────────────────────────────────────────── - -octave! - disp(strcat("foo", "bar")) - disp(strcmp("a", "a")) - disp(upper("hello")) - disp(lower("WORLD")) - disp(num2str(3.14)) - disp(str2double("2.718")) - disp(strtrim(" hi ")) -octave_end - --- ───────────────────────────────────────────────────────────────── --- §24 TYPE QUERIES & SHAPE --- ───────────────────────────────────────────────────────────────── - --- Note: class(...) is not in the DSL — "class" is a Lean keyword. -octave! - disp(isnumeric(42.0)) - disp(ischar("x")) - disp(isempty([])) - disp(numel([1.0, 2.0, 3.0])) - disp(size([[1.0, 2.0], [3.0, 4.0]])) - disp(rows([[1.0, 2.0], [3.0, 4.0]])) - disp(columns([[1.0, 2.0], [3.0, 4.0]])) -octave_end - --- ───────────────────────────────────────────────────────────────── --- §25 RESHAPE / HORZCAT / VERTCAT --- ───────────────────────────────────────────────────────────────── - -octave! - v = 1:6; - M = reshape(v, 2.0, 3.0) - A = [[1.0, 2.0], [3.0, 4.0]]; - B = [[5.0, 6.0], [7.0, 8.0]]; - disp(horzcat(A, B)) - disp(vertcat(A, B)) -octave_end - --- ───────────────────────────────────────────────────────────────── --- §26 PUTTING IT ALL TOGETHER — Newton's method --- ───────────────────────────────────────────────────────────────── - -octave! - function x = newton_sqrt(n, tol) - x = n / 2.0; - while abs(x * x - n) > tol - x = x - (x * x - n) / (2.0 * x); - endwhile - endfunction - - disp(newton_sqrt(2.0, 1e-10)) -- ≈ 1.4142135624 - disp(newton_sqrt(9.0, 1e-10)) -- ≈ 3.0 - disp(newton_sqrt(16.0, 1e-10)) -- ≈ 4.0 -octave_end - --- ───────────────────────────────────────────────────────────────── --- §27 PROOF INTEROP — expose AST for BigStep / PureEval proofs --- ───────────────────────────────────────────────────────────────── - --- `octave_stmts! name ... octave_end` gives you the program as a named --- `Array OctiveLean.Stmt` definition that you can reason about in Lean. - -octave_stmts! myProg - x = 0.0; - for k = 1:3 - x = x + k; +-- §15 BIND THE PARSED AST AS A LEAN TERM (for proof interop) +octave_program! mySumProgram { + s = 0; + for k = 1 : 10 + s = s + k; endfor -octave_end + disp(s) +} --- myProg is now a Lean definition you can use in proofs: -#check (myProg : Array OctiveLean.Stmt) +#check mySumProgram -- : Array OctiveLean.Stmt +#eval mySumProgram.size diff --git a/Sim_Gravity.m b/Sim_Gravity.m new file mode 100644 index 0000000..7d0bf35 --- /dev/null +++ b/Sim_Gravity.m @@ -0,0 +1,29 @@ +% Example 1: 1-D non-dim gravity +% x' = v +% v' = -1/x^2, x(0) = 1, v(0) = 0 +% RK4 fixed-step. + +n = 100; +t0 = 0; tf = 1.0; +h = (tf - t0) / n; + +t = zeros(1, n+1); +xs = zeros(1, n+1); +vs = zeros(1, n+1); +xs(1) = 1.0; +vs(1) = 0.0; + +for i = 1:n + ti = t(i); xi = xs(i); vi = vs(i); + k1x = vi; k1v = -1 / xi^2; + k2x = vi + h/2*k1v; k2v = -1 / (xi + h/2*k1x)^2; + k3x = vi + h/2*k2v; k3v = -1 / (xi + h/2*k2x)^2; + k4x = vi + h*k3v; k4v = -1 / (xi + h*k3x)^2; + t(i+1) = ti + h; + xs(i+1) = xi + h/6 * (k1x + 2*k2x + 2*k3x + k4x); + vs(i+1) = vi + h/6 * (k1v + 2*k2v + 2*k3v + k4v); +endfor + +for i = 1:n+1 + printf("%f,%f,%f\n", t(i), xs(i), vs(i)); +endfor diff --git a/Sim_Lorenz.m b/Sim_Lorenz.m new file mode 100644 index 0000000..21f518c --- /dev/null +++ b/Sim_Lorenz.m @@ -0,0 +1,53 @@ +% Example 3: Lorenz system +% x' = -sigma*x + sigma*y +% y' = rho*x - y - x*z +% z' = -beta*z + x*y +% Slide uses sigma = 10, rho = 70, beta = 8/3, x0=y0=z0=1. + +sigma = 10.0; +rho = 70.0; +beta = 8.0/3.0; + +n = 5000; +t0 = 0; tf = 25.0; +h = (tf - t0) / n; + +t = zeros(1, n+1); +xs = zeros(1, n+1); +ys = zeros(1, n+1); +zs = zeros(1, n+1); +xs(1) = 1.0; +ys(1) = 1.0; +zs(1) = 1.0; + +for i = 1:n + ti = t(i); xi = xs(i); yi = ys(i); zi = zs(i); + k1x = -sigma*xi + sigma*yi; + k1y = rho*xi - yi - xi*zi; + k1z = -beta*zi + xi*yi; + + ax = xi + h/2*k1x; ay = yi + h/2*k1y; az = zi + h/2*k1z; + k2x = -sigma*ax + sigma*ay; + k2y = rho*ax - ay - ax*az; + k2z = -beta*az + ax*ay; + + ax = xi + h/2*k2x; ay = yi + h/2*k2y; az = zi + h/2*k2z; + k3x = -sigma*ax + sigma*ay; + k3y = rho*ax - ay - ax*az; + k3z = -beta*az + ax*ay; + + ax = xi + h*k3x; ay = yi + h*k3y; az = zi + h*k3z; + k4x = -sigma*ax + sigma*ay; + k4y = rho*ax - ay - ax*az; + k4z = -beta*az + ax*ay; + + t(i+1) = ti + h; + xs(i+1) = xi + h/6 * (k1x + 2*k2x + 2*k3x + k4x); + ys(i+1) = yi + h/6 * (k1y + 2*k2y + 2*k3y + k4y); + zs(i+1) = zi + h/6 * (k1z + 2*k2z + 2*k3z + k4z); +endfor + +% Print every 10th sample +for i = 1:10:n+1 + printf("%f,%f,%f,%f\n", t(i), xs(i), ys(i), zs(i)); +endfor diff --git a/Sim_VanDerPol.m b/Sim_VanDerPol.m new file mode 100644 index 0000000..44b17b1 --- /dev/null +++ b/Sim_VanDerPol.m @@ -0,0 +1,34 @@ +% Example 2: van der Pol oscillator +% x' = v +% v' = mu*(1 - x^2)*v - x, x(0)=0, v(0)=1 +% Use mu = 2 (stiffer values like 50 from the slide need adaptive step). + +mu = 2.0; +n = 3000; +t0 = 0; tf = 30.0; +h = (tf - t0) / n; + +t = zeros(1, n+1); +xs = zeros(1, n+1); +vs = zeros(1, n+1); +xs(1) = 0.0; +vs(1) = 1.0; + +for i = 1:n + ti = t(i); xi = xs(i); vi = vs(i); + k1x = vi; k1v = mu*(1 - xi^2)*vi - xi; + ax = xi + h/2*k1x; av = vi + h/2*k1v; + k2x = av; k2v = mu*(1 - ax^2)*av - ax; + ax = xi + h/2*k2x; av = vi + h/2*k2v; + k3x = av; k3v = mu*(1 - ax^2)*av - ax; + ax = xi + h*k3x; av = vi + h*k3v; + k4x = av; k4v = mu*(1 - ax^2)*av - ax; + t(i+1) = ti + h; + xs(i+1) = xi + h/6 * (k1x + 2*k2x + 2*k3x + k4x); + vs(i+1) = vi + h/6 * (k1v + 2*k2v + 2*k3v + k4v); +endfor + +% Print every 10th sample to keep CSV manageable +for i = 1:10:n+1 + printf("%f,%f,%f\n", t(i), xs(i), vs(i)); +endfor diff --git a/SymToolboxDemo.m b/SymToolboxDemo.m new file mode 100644 index 0000000..71f1514 --- /dev/null +++ b/SymToolboxDemo.m @@ -0,0 +1,45 @@ +% Symbolic Math Toolbox - cheat sheet walkthrough. +% Each labeled block produces one line of output. + +x = sym('x'); y = sym('y'); z = sym('z'); t = sym('t'); +a = sym('a'); b = sym('b'); k = sym('k'); n = sym('n'); + +% Calculus +printf("diff: "); disp(diff(sym('sin(x^2 + t)'), x)); +printf("int indef: "); disp(int(sym('x/(1+z^2)'), z)); +printf("int def: "); disp(int(sym('x^2'), x, 0, 1)); +printf("limit left: "); disp(limit(sym('1/x'), x, 0, "left")); +printf("taylor: "); disp(taylor(sym('exp(-x)'))); +printf("series: "); disp(series(sym('1/sin(x)'), x)); +printf("symsum: "); disp(symsum(k, k, 0, n - 1)); + +printf("gradient: "); disp(gradient(sym('x*y + 2*z*x'), sym('[x, y, z]'))); +printf("jacobian: "); disp(jacobian(sym('[x*y*z, y, x+z]'), sym('[x, y, z]'))); +printf("hessian: "); disp(hessian(sym('x*y + 2*z*x'), sym('[x, y, z]'))); +printf("laplacian: "); disp(laplacian(sym('1/x + y^2 + z^3'), sym('[x, y, z]'))); + +% Algebra +printf("double pi: "); disp(double(sym('pi'))); +printf("vpa pi 30: "); disp(vpa(sym('pi'), 30)); +printf("subs: "); disp(subs(sym('a^3 + b'), a, 2)); +printf("solve poly: "); disp(solve(sym('x^2 - 4'), x)); +printf("solve sys: "); disp(solve(sym('[u + v - a, u - v - b]'), sym('[u, v]'))); +printf("isolate: "); disp(isolate(sym('a*x^2 + b*x + c'), x)); +printf("lhs: "); disp(lhs(sym('Eq(x^2, y^2)'))); +printf("rhs: "); disp(rhs(sym('Eq(x^2, y^2)'))); +printf("simplify: "); disp(simplify(sym('sin(x)^2 + cos(x)^2'))); +printf("expand: "); disp(expand(sym('(x+1)^3'))); +printf("factor: "); disp(factor(sym('x^2 - 1'))); +printf("collect: "); disp(collect(sym('x*y + x^2 + 2*x*y + 3'), x)); +printf("rewrite: "); disp(rewrite(sym('tan(x)/cos(x)'), "sin")); +printf("resultant: "); disp(resultant(sym('x^2 + y'), sym('x - 2*y'), x)); + +% ODE - symfun() registers a SymPy Function so f(t) parses as f-of-t +symfun('f'); +printf("dsolve: "); disp(dsolve(sym('Eq(Derivative(f(t), t), a*f(t))'), sym('f(t)'))); + +% Functions +printf("piecewise: "); disp(piecewise(sym('x < 0'), -1, sym('x >= 0'), 2)); + +% Output formats +printf("latex: "); disp(latex(sym('x^2 + y^2'))); diff --git a/corpus/30_sym_basic.expected b/corpus/30_sym_basic.expected new file mode 100644 index 0000000..b3aebad --- /dev/null +++ b/corpus/30_sym_basic.expected @@ -0,0 +1,3 @@ +x**2 + 2*x + 1 +2*x + 2 +x**2 + 2*x diff --git a/corpus/30_sym_basic.m b/corpus/30_sym_basic.m new file mode 100644 index 0000000..6dc8e71 --- /dev/null +++ b/corpus/30_sym_basic.m @@ -0,0 +1,5 @@ +x = sym('x'); +f = x^2 + 2*x + 1; +disp(f); +disp(diff(f, x)); +disp(int(diff(f, x), x)); diff --git a/corpus/31_sym_solve_simplify.expected b/corpus/31_sym_solve_simplify.expected new file mode 100644 index 0000000..4c846d1 --- /dev/null +++ b/corpus/31_sym_solve_simplify.expected @@ -0,0 +1,4 @@ +[-2, 2] +1 +(x - 1)*(x + 1) +x**3 + 3*x**2 + 3*x + 1 diff --git a/corpus/31_sym_solve_simplify.m b/corpus/31_sym_solve_simplify.m new file mode 100644 index 0000000..8c6d5ec --- /dev/null +++ b/corpus/31_sym_solve_simplify.m @@ -0,0 +1,5 @@ +x = sym('x'); +disp(solve(x^2 - 4, x)); +disp(simplify(sym('sin(x)^2 + cos(x)^2'))); +disp(factor(sym('x^2 - 1'))); +disp(expand(sym('(x+1)^3'))); diff --git a/corpus/32_sym_calc.expected b/corpus/32_sym_calc.expected new file mode 100644 index 0000000..4f6f713 --- /dev/null +++ b/corpus/32_sym_calc.expected @@ -0,0 +1,3 @@ +x**5/120 + x**4/24 + x**3/6 + x**2/2 + x + 1 +1 +5 diff --git a/corpus/32_sym_calc.m b/corpus/32_sym_calc.m new file mode 100644 index 0000000..17a8311 --- /dev/null +++ b/corpus/32_sym_calc.m @@ -0,0 +1,4 @@ +x = sym('x'); +disp(taylor(sym('exp(x)'))); +disp(limit(sym('sin(x)/x'), x, 0)); +disp(subs(sym('x^2 + 1'), x, sym('2'))); diff --git a/docs/typst-diaries.md b/docs/typst-diaries.md new file mode 100644 index 0000000..dd84b9a --- /dev/null +++ b/docs/typst-diaries.md @@ -0,0 +1,46 @@ +# Typst documents tied to octive-lean + +Typst sources for the m242 command-line diaries that demo octive-lean live +outside this repo, in `~/.env/typst/m242/`. The `.m` drivers in this repo's +root produce the data each diary plots; the typst files compile to PDF and +embed the plots, screenshots, and prose. + +| Typst file | PDF | Topic | Driver(s) in this repo | +| --- | --- | --- | --- | +| `~/.env/typst/m242/CLDiary.typ` | `CLDiary.pdf` | Polynomial interpolation: Runge phenomenon, least-squares fit, splines, Chebyshev nodes | `Lab7Interp.m` | +| `~/.env/typst/m242/CLDiary_Sym.typ` | `CLDiary_Sym.pdf` | Symbolic Math Toolbox walkthrough (28 cheat-sheet ops via SymPy bridge) | `SymToolboxDemo.m` | +| `~/.env/typst/m242/CLDiary_Sim.typ` | `CLDiary_Sim.pdf` | Simulink/Xcos: 4 dynamic systems with Xcos canvas screenshots + native fletcher diagrams + RK4 trajectories | `Sim_Gravity.m`, `Sim_VanDerPol.m`, `Sim_Lorenz.m` | + +## Build + +```sh +cd ~/.env/typst/m242 +typst compile CLDiary.typ +typst compile CLDiary_Sym.typ +typst compile CLDiary_Sim.typ +``` + +## Supporting assets (in `~/.env/typst/m242/`) + +| Path | What | +| --- | --- | +| `sim_data/*.csv` | Trajectories produced by `Sim_*.m`, loaded by `CLDiary_Sim.typ` | +| `screenshots/xcos_*.png` | Xcos canvas screenshots for `CLDiary_Sim.typ` | +| `xcos/*.zcos` | Scilab/Xcos diagram files (Lorenz, Bouncing_ball, gensin, pendulum, Inverted_pendulum, Colpitts, Boost_Converter) | +| `xcos/BUILD_DIAGRAMS.md` | How to build / screenshot each Xcos diagram | + +## Regenerating sim_data + +```sh +cd ~/.env/lean/octive-lean +lake exe octive-lean Sim_Gravity.m | grep , > ~/.env/typst/m242/sim_data/gravity.csv +lake exe octive-lean Sim_VanDerPol.m | grep , > ~/.env/typst/m242/sim_data/vanderpol.csv +lake exe octive-lean Sim_Lorenz.m | grep , > ~/.env/typst/m242/sim_data/lorenz.csv +``` + +## Octive-lean features added for these diaries + +- `polyfit`, `polyval`, `spline`, `linsolve` — `OctiveLean/Builtins.lean` +- `OctiveLean/SymPyBridge.lean` — persistent SymPy subprocess +- 25+ symbolic builtins (`diff`, `int`, `subs`, `simplify`, `solve`, `taylor`, `dsolve`, `jacobian`, `hessian`, `laplacian`, `symsum`, `rewrite`, `resultant`, `series`, `isolate`, `lhs`/`rhs`, `latex`, `pretty`, `vpa`, `coeffs`, `collect`, `expand`, `factor`, `gradient`, `piecewise`, `symfun`) — `OctiveLean/Builtins.lean` +- `Value.sym` variant + binop overloading for symbolic operands — `OctiveLean/Value.lean`, `OctiveLean/Eval.lean`