diff --git a/benchmark/scripts/collect_results.py b/benchmark/scripts/collect_results.py
new file mode 100644
index 00000000000..b874876d871
--- /dev/null
+++ b/benchmark/scripts/collect_results.py
@@ -0,0 +1,62 @@
+#-------------------------------------------------------------
+#
+# Licensed to the Apache Software Foundation (ASF) under one
+# or more contributor license agreements.  See the NOTICE file
+# distributed with this work for additional information
+# regarding copyright ownership.  The ASF licenses this file
+# to you under the Apache License, Version 2.0 (the
+# "License"); you may not use this file except in compliance
+# with the License.  You may obtain a copy of the License at
+#
+#   http://www.apache.org/licenses/LICENSE-2.0
+#
+# Unless required by applicable law or agreed to in writing,
+# software distributed under the License is distributed on an
+# "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
+# KIND, either express or implied.  See the License for the
+# specific language governing permissions and limitations
+# under the License.
+#
+#-------------------------------------------------------------
+
+"""
+Parse per-run accuracy files into a single results.csv.
+
+Output columns: label, epsilon, private, accuracy
+"""
+import pathlib, csv, re
+
+RESULTS = pathlib.Path("benchmark/results")
+
+rows = []
+
+def parse_acc(path: pathlib.Path) -> float:
+    txt = path.read_text().strip()
+    # SystemDS writes a bare float.
+    return float(txt)
+
+# Non-private baseline.
+baseline_path = RESULTS / "acc_baseline.txt"
+if baseline_path.exists():
+    rows.append(dict(label="baseline", epsilon="inf",
+                     private=0, accuracy=parse_acc(baseline_path)))
+
+# DP runs.
+for eps in [0.5, 1, 4, 8]:
+    p = RESULTS / f"acc_eps_{eps}.txt"
+    if p.exists():
+        rows.append(dict(label=f"epsilon={eps}", epsilon=eps,
+                         private=1, accuracy=parse_acc(p)))
+    else:
+        print(f"Warning: {p} not found - skipping")
+
+out = RESULTS / "results.csv"
+with open(out, "w", newline="") as f:
+    w = csv.DictWriter(f, fieldnames=["label","epsilon","private","accuracy"])
+    w.writeheader()
+    w.writerows(rows)
+
+print(f"Wrote {out}")
+for r in rows:
+    print(f"  {r['label']:12s}  acc={r['accuracy']:.4f}")
+
diff --git a/benchmark/scripts/eval.dml b/benchmark/scripts/eval.dml
new file mode 100644
index 00000000000..a9ad9e179c0
--- /dev/null
+++ b/benchmark/scripts/eval.dml
@@ -0,0 +1,43 @@
+#-------------------------------------------------------------
+#
+# Licensed to the Apache Software Foundation (ASF) under one
+# or more contributor license agreements.  See the NOTICE file
+# distributed with this work for additional information
+# regarding copyright ownership.  The ASF licenses this file
+# to you under the Apache License, Version 2.0 (the
+# "License"); you may not use this file except in compliance
+# with the License.  You may obtain a copy of the License at
+#
+#   http://www.apache.org/licenses/LICENSE-2.0
+#
+# Unless required by applicable law or agreed to in writing,
+# software distributed under the License is distributed on an
+# "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
+# KIND, either express or implied.  See the License for the
+# specific language governing permissions and limitations
+# under the License.
+#
+#-------------------------------------------------------------
+
+# eval.dml — compute binary classification accuracy on held-out test set.
+# Arguments: data_dir, model_path, out_acc
+data_dir   = $data_dir;
+model_path = $model_path;
+out_acc    = $out_acc;
+
+X_test = read(data_dir + "/X_test.csv",
+              data_type="matrix", value_type="double", format="csv");
+y_test = read(data_dir + "/y_test.csv",
+              data_type="matrix", value_type="double", format="csv");
+w      = read(model_path,
+              data_type="matrix", value_type="double", format="csv");
+
+scores    = X_test %*% w;
+preds     = (scores > 0.0);          # threshold at 0 (log-odds)
+correct   = sum(preds == y_test);
+n_test    = nrow(y_test);
+accuracy  = correct / n_test;
+
+print("Accuracy: " + accuracy);
+write(accuracy, out_acc, format="csv");
+
diff --git a/benchmark/scripts/fedavg_dp.dml b/benchmark/scripts/fedavg_dp.dml
new file mode 100644
index 00000000000..1715f44758e
--- /dev/null
+++ b/benchmark/scripts/fedavg_dp.dml
@@ -0,0 +1,151 @@
+#-------------------------------------------------------------
+#
+# Licensed to the Apache Software Foundation (ASF) under one
+# or more contributor license agreements.  See the NOTICE file
+# distributed with this work for additional information
+# regarding copyright ownership.  The ASF licenses this file
+# to you under the Apache License, Version 2.0 (the
+# "License"); you may not use this file except in compliance
+# with the License.  You may obtain a copy of the License at
+#
+#   http://www.apache.org/licenses/LICENSE-2.0
+#
+# Unless required by applicable law or agreed to in writing,
+# software distributed under the License is distributed on an
+# "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
+# KIND, either express or implied.  See the License for the
+# specific language governing permissions and limitations
+# under the License.
+#
+#-------------------------------------------------------------
+
+# ── fedavg_dp.dml ────────────────────────────────────────────────────────────
+# Arguments (passed via -nvargs):
+#   data_dir   : path to benchmark/data/
+#   epsilon    : DP privacy budget ε  (use 9999 for non-private baseline)
+#   delta      : DP delta             (1e-5 can be used)
+#   clip_norm  : per-example gradient L2-norm clip bound (default 4.0) —
+#                sensitivity = clip_norm / n follows from this, since
+#                clipping each record's gradient contribution to clip_norm
+#                is what makes that bound actually hold.
+#   n_rounds   : number of FedAvg rounds (default 50)
+#   lr         : learning rate          (default 0.1)
+#   w1_rows    : rows in worker 1 shard
+#   w2_rows    : rows in worker 2 shard
+#   w3_rows    : rows in worker 3 shard
+#   w4_rows    : rows in worker 4 shard
+#   n_features : number of features
+#   out        : path to write final weights
+#   private    : 1 = apply DP noise (default), 0 = non-private baseline
+
+data_dir   = $data_dir;
+epsilon    = $epsilon;
+delta      = $delta;
+clip_norm  = ifdef($clip_norm,  4.0);
+n_rounds   = ifdef($n_rounds,   50);
+lr         = ifdef($lr,         0.1);
+private    = ifdef($private,    1);
+out_path   = $out;
+
+# Per-round epsilon: dp_gaussian is called once per round, and
+# DPBudgetAccountant composes the cost of every release against the single
+# total budget set by dp_set_budget(). Spending the full epsilon on every
+# round would exhaust the budget almost immediately, so split it evenly
+# across rounds instead.
+round_epsilon = epsilon / n_rounds;
+
+w1r = $w1_rows;
+w2r = $w2_rows;
+w3r = $w3_rows;
+w4r = $w4_rows;
+d   = $n_features;
+n   = w1r + w2r + w3r + w4r;
+
+# Sensitivity of the released (mean) gradient to a single record changing:
+# with each record's contribution clipped to L2-norm <= clip_norm below,
+# the sum can move by at most clip_norm, so the average moves by clip_norm/n.
+sensitivity = clip_norm / n;
+
+# ── Build federated matrix from 4 local workers ───────────────────────────
+# Row ranges are 0-based [begin, end) for each partition.
+r1s=0;        r1e=w1r;
+r2s=w1r;      r2e=w1r+w2r;
+r3s=w1r+w2r;  r3e=w1r+w2r+w3r;
+r4s=w1r+w2r+w3r; r4e=n;
+
+X = federated(
+    addresses=list(
+        "localhost:8301/" + data_dir + "/worker1/X_train.csv",
+        "localhost:8302/" + data_dir + "/worker2/X_train.csv",
+        "localhost:8303/" + data_dir + "/worker3/X_train.csv",
+        "localhost:8304/" + data_dir + "/worker4/X_train.csv"),
+    ranges=list(
+        list(r1s, 0), list(r1e, d),
+        list(r2s, 0), list(r2e, d),
+        list(r3s, 0), list(r3e, d),
+        list(r4s, 0), list(r4e, d)));
+
+y = federated(
+    addresses=list(
+        "localhost:8301/" + data_dir + "/worker1/y_train.csv",
+        "localhost:8302/" + data_dir + "/worker2/y_train.csv",
+        "localhost:8303/" + data_dir + "/worker3/y_train.csv",
+        "localhost:8304/" + data_dir + "/worker4/y_train.csv"),
+    ranges=list(
+        list(r1s, 0), list(r1e, 1),
+        list(r2s, 0), list(r2e, 1),
+        list(r3s, 0), list(r3e, 1),
+        list(r4s, 0), list(r4e, 1)));
+
+# ── Initialise weights ────────────────────────────────────────────────────
+w = matrix(0, rows=d, cols=1);
+
+if (private == 1) {
+        eps = dp_set_budget($epsilon, $delta);
+}
+
+# ── Training rounds ───────────────────────────────────────────────────────
+for (round in 1:n_rounds) {
+
+    # Forward pass — executes on federated workers.
+    scores = X %*% w;                           # (n × 1), federated
+    probs  = 1.0 / (1.0 + exp(-scores));        # (n × 1), federated
+
+    residuals = probs - y;                       # (n × 1), federated
+
+    # DP noise injection — only when private=1.
+    if (private == 1) {
+        # Per-example L2-norm clipping: each record's raw gradient
+        # contribution is X[i,:] * residual_i, with norm
+        # ||X[i,:]||_2 * |residual_i|. Scaling it down to clip_norm
+        # whenever it exceeds that bound is what makes
+        # sensitivity = clip_norm / n (set above) an actual, provable
+        # bound instead of an arbitrary constant.
+        row_norms         = sqrt(rowSums(X^2));               # (n × 1) ||X[i,:]||_2
+        contrib_norms     = abs(residuals) * row_norms;        # (n × 1) ||X[i,:]*residual_i||_2
+        clip_scale        = clip_norm / pmax(contrib_norms, clip_norm); # (n × 1), in (0,1]
+        clipped_residuals = residuals * clip_scale;            # (n × 1)
+
+        # Gradient aggregation — t(X) %*% clipped_residuals is a (d × 1)
+        # local sum that the coordinator collects in one federated
+        # aggregate instruction.
+        grad = t(X) %*% clipped_residuals / n;                # (d × 1), LOCAL after agg, clipped
+
+        noisy_grad = dp_gaussian(grad,
+                                 "identity",
+                                 sensitivity=sensitivity,
+                                 epsilon=round_epsilon,
+                                 delta=delta);
+        w = w - lr * noisy_grad;
+    } else {
+        # Gradient aggregation — t(X) %*% residuals is an (d × 1) local sum
+        # that the coordinator collects in one federated aggregate instruction.
+        grad = t(X) %*% residuals / n;                        # (d × 1), LOCAL after agg
+        w = w - lr * grad;
+    }
+}
+
+# ── Write model weights ───────────────────────────────────────────────────
+write(w, out_path, format="csv");
+print("FedAvg done. epsilon=" + epsilon + "  rounds=" + n_rounds);
+
diff --git a/benchmark/scripts/plot.py b/benchmark/scripts/plot.py
new file mode 100644
index 00000000000..7b160d065de
--- /dev/null
+++ b/benchmark/scripts/plot.py
@@ -0,0 +1,134 @@
+#-------------------------------------------------------------
+#
+# Licensed to the Apache Software Foundation (ASF) under one
+# or more contributor license agreements.  See the NOTICE file
+# distributed with this work for additional information
+# regarding copyright ownership.  The ASF licenses this file
+# to you under the Apache License, Version 2.0 (the
+# "License"); you may not use this file except in compliance
+# with the License.  You may obtain a copy of the License at
+#
+#   http://www.apache.org/licenses/LICENSE-2.0
+#
+# Unless required by applicable law or agreed to in writing,
+# software distributed under the License is distributed on an
+# "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
+# KIND, either express or implied.  See the License for the
+# specific language governing permissions and limitations
+# under the License.
+#
+#-------------------------------------------------------------
+"""
+Read results.csv and produce two figures:
+
+1. accuracy_vs_epsilon.png
+   Line plot: x = epsilon, y = accuracy.
+   Horizontal dashed line = non-private baseline.
+   Points labelled with accuracy values.
+
+2. privacy_cost.png
+   Bar chart showing accuracy loss relative to baseline (utility cost of DP).
+"""
+import pathlib
+import csv
+import matplotlib
+matplotlib.use("Agg")
+import matplotlib.pyplot as plt
+import matplotlib.ticker as ticker
+
+RESULTS = pathlib.Path("benchmark/results")
+
+# ── Load ──────────────────────────────────────────────────────────────────
+rows = []
+with open(RESULTS / "results.csv") as f:
+    for r in csv.DictReader(f):
+        rows.append({
+            "label":    r["label"],
+            "epsilon":  float(r["epsilon"]) if r["epsilon"] != "inf" else None,
+            "private":  int(r["private"]),
+            "accuracy": float(r["accuracy"]),
+        })
+
+baseline = next(r for r in rows if r["private"] == 0)
+dp_rows  = sorted([r for r in rows if r["private"] == 1],
+                  key=lambda r: r["epsilon"])
+
+eps_vals = [r["epsilon"] for r in dp_rows]
+acc_vals = [r["accuracy"] for r in dp_rows]
+baseline_acc = baseline["accuracy"]
+
+# ── Figure 1: Accuracy vs epsilon ───────────────────────────────────────────────
+fig, ax = plt.subplots(figsize=(7, 4.5))
+
+ax.plot(eps_vals, acc_vals, marker="o", linewidth=2,
+        color="#028090", label="DP-FedAvg (Gaussian)")
+ax.axhline(baseline_acc, linestyle="--", color="#1C3A5E",
+           linewidth=1.5, label=f"Non-private baseline ({baseline_acc:.3f})")
+
+# Annotate each DP point.
+for eps, acc in zip(eps_vals, acc_vals):
+    ax.annotate(f"{acc:.3f}", xy=(eps, acc),
+                xytext=(0, 8), textcoords="offset points",
+                ha="center", fontsize=9, color="#028090")
+
+ax.set_xscale("log")
+ax.set_xticks(eps_vals)
+ax.get_xaxis().set_major_formatter(ticker.ScalarFormatter())
+ax.set_xlabel("Privacy budget epsilon  (smaller = stronger privacy)", fontsize=11)
+ax.set_ylabel("Test accuracy", fontsize=11)
+ax.set_title("Accuracy vs. Privacy Budget - DP-FedAvg on Adult (4 workers)",
+             fontsize=12)
+ax.legend(fontsize=9)
+ax.set_ylim(max(0, min(acc_vals) - 0.05), min(1.0, baseline_acc + 0.05))
+ax.grid(True, which="both", linestyle=":", alpha=0.5)
+
+plt.tight_layout()
+out1 = RESULTS / "accuracy_vs_epsilon.png"
+fig.savefig(out1, dpi=150)
+print(f"Saved {out1}")
+plt.close()
+
+# ── Figure 2: Utility cost (accuracy drop) ────────────────────────────────
+fig, ax = plt.subplots(figsize=(6, 4))
+
+drops = [baseline_acc - acc for acc in acc_vals]
+colors = ["#B91C1C" if d > 0.02 else "#028090" for d in drops]
+# Position bars at their true epsilon value on a log-scaled x-axis (rather than
+# evenly-spaced categorical slots) so the visual spacing between 0.5 to 1 and
+# 4 to 8 reflects the same 2x ratio. Bar widths scale with x so they stay a
+# constant fraction of their slot in log space instead of shrinking/growing.
+widths = [e * 0.4 for e in eps_vals]
+bars = ax.bar(eps_vals, drops, color=colors, width=widths,
+              edgecolor="white")
+ax.set_xscale("log")
+ax.set_xticks(eps_vals)
+ax.get_xaxis().set_major_formatter(ticker.ScalarFormatter())
+
+for bar, drop in zip(bars, drops):
+    ax.text(bar.get_x() + bar.get_width() / 2,
+            bar.get_height() + 0.001,
+            f"{drop:.3f}", ha="center", va="bottom", fontsize=9)
+
+ax.legend(["drop > 0.02", "drop <= 0.02"])
+
+ax.axhline(0, color="black", linewidth=0.8)
+ax.set_xlabel("Privacy budget epsilon", fontsize=11)
+ax.set_ylabel("Accuracy drop vs. baseline", fontsize=11)
+ax.set_title("Utility Cost of Differential Privacy - DP-FedAvg on Adult",
+             fontsize=11)
+ax.grid(True, axis="y", linestyle=":", alpha=0.5)
+
+plt.tight_layout()
+out2 = RESULTS / "privacy_cost.png"
+fig.savefig(out2, dpi=150)
+print(f"Saved {out2}")
+plt.close()
+
+# ── Console summary table ─────────────────────────────────────────────────
+print()
+print(f"{'epsilon':>8}  {'accuracy':>10}  {'drop':>8}")
+print("-" * 34)
+print(f"{'baseline':>8}  {baseline_acc:10.4f}  {'-':>8}")
+for eps, acc, drop in zip(eps_vals, acc_vals, drops):
+    print(f"{eps:>8.1f}  {acc:10.4f}  {drop:8.4f}")
+
diff --git a/benchmark/scripts/prepare_data.py b/benchmark/scripts/prepare_data.py
new file mode 100644
index 00000000000..1a18e067bc7
--- /dev/null
+++ b/benchmark/scripts/prepare_data.py
@@ -0,0 +1,143 @@
+#-------------------------------------------------------------
+#
+# Licensed to the Apache Software Foundation (ASF) under one
+# or more contributor license agreements.  See the NOTICE file
+# distributed with this work for additional information
+# regarding copyright ownership.  The ASF licenses this file
+# to you under the Apache License, Version 2.0 (the
+# "License"); you may not use this file except in compliance
+# with the License.  You may obtain a copy of the License at
+#
+#   http://www.apache.org/licenses/LICENSE-2.0
+#
+# Unless required by applicable law or agreed to in writing,
+# software distributed under the License is distributed on an
+# "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
+# KIND, either express or implied.  See the License for the
+# specific language governing permissions and limitations
+# under the License.
+#
+#-------------------------------------------------------------
+"""
+Download the UCI Adult dataset, binarise labels, standardise features,
+split into 4 equal horizontal partitions for federated workers, and write
+SystemDS .mtd metadata files alongside each CSV shard.
+
+Outputs
+-------
+benchmark/data/worker{1..4}/X_train.csv  + X_train.csv.mtd
+benchmark/data/worker{1..4}/y_train.csv  + y_train.csv.mtd
+benchmark/data/X_test.csv                + X_test.csv.mtd
+benchmark/data/y_test.csv                + y_test.csv.mtd
+benchmark/data/meta.txt                  # n_train, n_test, n_features
+"""
+
+import json, os, pathlib
+import numpy as np
+import pandas as pd
+from sklearn.model_selection import train_test_split
+from sklearn.preprocessing import StandardScaler
+
+ADULT_TRAIN_URL = (
+    "https://archive.ics.uci.edu/ml/machine-learning-databases"
+    "/adult/adult.data"
+)
+ADULT_TEST_URL = (
+    "https://archive.ics.uci.edu/ml/machine-learning-databases"
+    "/adult/adult.test"
+)
+
+COLS = [
+    "age","workclass","fnlwgt","education","education_num","marital_status",
+    "occupation","relationship","race","sex","capital_gain","capital_loss",
+    "hours_per_week","native_country","label",
+]
+NUMERIC = ["age","fnlwgt","education_num","capital_gain",
+           "capital_loss","hours_per_week"]
+
+DATA_DIR = pathlib.Path("benchmark/data")
+N_WORKERS = 4
+
+def download(url, dest):
+    import urllib.request
+    if not dest.exists():
+        print(f"Downloading {url}")
+        urllib.request.urlretrieve(url, dest)
+
+def load_adult(path, skip_rows=0):
+    df = pd.read_csv(path, names=COLS, skipinitialspace=True,
+                     skiprows=skip_rows, na_values="?").dropna()
+    # binarise label: >50K => 1, else 0
+    df["label"] = (df["label"].str.strip().str.rstrip(".") == ">50K").astype(float)
+    # one-hot encode categoricals
+    cats = [c for c in COLS[:-1] if c not in NUMERIC]
+    df = pd.get_dummies(df, columns=cats, drop_first=True)
+    return df
+
+def write_csv_and_mtd(arr: np.ndarray, path: pathlib.Path, description: str):
+    """Write a CSV and a matching SystemDS .mtd metadata file."""
+    path.parent.mkdir(parents=True, exist_ok=True)
+    np.savetxt(path, arr, delimiter=",", fmt="%.8f")
+    rows, cols = arr.shape
+    mtd = {
+        "data_type": "matrix",
+        "value_type": "double",
+        "rows": rows,
+        "cols": cols,
+        "format": "csv",
+        "header": False,
+        "description": description,
+    }
+    with open(str(path) + ".mtd", "w") as f:
+        json.dump(mtd, f, indent=2)
+    print(f"  {path}  ({rows} × {cols})")
+
+# ── Download ─────────────────────────────────────────────────────────────────
+download(ADULT_TRAIN_URL, DATA_DIR / "raw" / "adult.data")
+download(ADULT_TEST_URL,  DATA_DIR / "raw" / "adult.test")
+
+train_df = load_adult(DATA_DIR / "raw" / "adult.data")
+test_df  = load_adult(DATA_DIR / "raw" / "adult.test", skip_rows=1)
+
+# Align columns (test may have different dummies after get_dummies).
+train_df, test_df = train_df.align(test_df, join="left", axis=1, fill_value=0)
+
+# ── Feature / label split ────────────────────────────────────────────────────
+feature_cols = [c for c in train_df.columns if c != "label"]
+X_train = train_df[feature_cols].values.astype(float)
+y_train = train_df["label"].values.reshape(-1, 1).astype(float)
+X_test  = test_df[feature_cols].values.astype(float)
+y_test  = test_df["label"].values.reshape(-1, 1).astype(float)
+
+# ── Standardise (fit on train only) ─────────────────────────────────────────
+scaler = StandardScaler()
+X_train = scaler.fit_transform(X_train)
+X_test  = scaler.transform(X_test)
+
+n_train, n_features = X_train.shape
+n_test = X_test.shape[0]
+print(f"Train: {n_train} × {n_features}  |  Test: {n_test} × {n_features}")
+
+# ── Partition across workers (equal horizontal splits) ──────────────────────
+indices = np.array_split(np.arange(n_train), N_WORKERS)
+for i, idx in enumerate(indices, start=1):
+    wdir = DATA_DIR / f"worker{i}"
+    write_csv_and_mtd(X_train[idx], wdir / "X_train.csv",
+                      f"Adult features, worker {i}")
+    write_csv_and_mtd(y_train[idx], wdir / "y_train.csv",
+                      f"Adult labels,   worker {i}")
+
+# ── Test set (coordinator-local) ────────────────────────────────────────────
+write_csv_and_mtd(X_test, DATA_DIR / "X_test.csv",  "Adult test features")
+write_csv_and_mtd(y_test, DATA_DIR / "y_test.csv",  "Adult test labels")
+
+# ── Metadata for DML scripts ─────────────────────────────────────────────────
+worker_rows = [len(idx) for idx in indices]
+with open(DATA_DIR / "meta.txt", "w") as f:
+    f.write(f"n_train={n_train}\n")
+    f.write(f"n_test={n_test}\n")
+    f.write(f"n_features={n_features}\n")
+    for i, r in enumerate(worker_rows, start=1):
+        f.write(f"worker{i}_rows={r}\n")
+print("Wrote benchmark/data/meta.txt")
+
diff --git a/benchmark/scripts/run_benchmark.sh b/benchmark/scripts/run_benchmark.sh
new file mode 100755
index 00000000000..52226c43f28
--- /dev/null
+++ b/benchmark/scripts/run_benchmark.sh
@@ -0,0 +1,36 @@
+#-------------------------------------------------------------
+#
+# Licensed to the Apache Software Foundation (ASF) under one
+# or more contributor license agreements.  See the NOTICE file
+# distributed with this work for additional information
+# regarding copyright ownership.  The ASF licenses this file
+# to you under the Apache License, Version 2.0 (the
+# "License"); you may not use this file except in compliance
+# with the License.  You may obtain a copy of the License at
+#
+#   http://www.apache.org/licenses/LICENSE-2.0
+#
+# Unless required by applicable law or agreed to in writing,
+# software distributed under the License is distributed on an
+# "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
+# KIND, either express or implied.  See the License for the
+# specific language governing permissions and limitations
+# under the License.
+#
+#-------------------------------------------------------------
+
+# 1. Prepare data (once).
+python benchmark/scripts/prepare_data.py
+
+# 2. Run the sweep (starts workers, trains, evaluates, stops workers).
+bash benchmark/scripts/run_sweep.sh
+
+# 3. Collect results into CSV.
+python benchmark/scripts/collect_results.py
+
+# 4. Generate plots.
+python benchmark/scripts/plot.py
+
+# 5. Confirm outputs exist.
+ls -lh benchmark/results/accuracy_vs_epsilon.png benchmark/results/privacy_cost.png
+
diff --git a/benchmark/scripts/run_sweep.sh b/benchmark/scripts/run_sweep.sh
new file mode 100755
index 00000000000..9c9ff3d9b60
--- /dev/null
+++ b/benchmark/scripts/run_sweep.sh
@@ -0,0 +1,90 @@
+#!/usr/bin/env bash
+#-------------------------------------------------------------
+#
+# Licensed to the Apache Software Foundation (ASF) under one
+# or more contributor license agreements.  See the NOTICE file
+# distributed with this work for additional information
+# regarding copyright ownership.  The ASF licenses this file
+# to you under the Apache License, Version 2.0 (the
+# "License"); you may not use this file except in compliance
+# with the License.  You may obtain a copy of the License at
+#
+#   http://www.apache.org/licenses/LICENSE-2.0
+#
+# Unless required by applicable law or agreed to in writing,
+# software distributed under the License is distributed on an
+# "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
+# KIND, either express or implied.  See the License for the
+# specific language governing permissions and limitations
+# under the License.
+#
+#-------------------------------------------------------------
+
+# Runs FedAvg for each epsilon value and the non-private baseline,
+# then evaluates accuracy. Results are appended to results/results.csv.
+set -euo pipefail
+
+REPO_ROOT="$(cd "$(dirname "$0")/../.." && pwd)"
+JAR="$REPO_ROOT/target/SystemDS.jar"
+SCRIPTS="$REPO_ROOT/benchmark/scripts"
+DATA="$REPO_ROOT/benchmark/data"
+RESULTS="$REPO_ROOT/benchmark/results"
+mkdir -p "$RESULTS"
+
+# Read dataset metadata written by prepare_data.py.
+source <(grep -E '^(n_train|n_test|n_features|worker[1-4]_rows)=' \
+              "$DATA/meta.txt" | sed 's/=/="/;s/$/"/')
+
+COMMON_ARGS="\
+  data_dir=$DATA \
+  n_features=$n_features \
+  w1_rows=$worker1_rows \
+  w2_rows=$worker2_rows \
+  w3_rows=$worker3_rows \
+  w4_rows=$worker4_rows \
+  n_rounds=300 \
+  lr=1.0 \
+  clip_norm=4.0 \
+  delta=1e-5"
+
+# ── Start workers ─────────────────────────────────────────────────────────
+echo "=== Starting federated workers ==="
+bash "$SCRIPTS/start_workers.sh"
+
+run_one() {
+  local label="$1"    # e.g. "eps_0.5" or "baseline"
+  local extra="$2"    # extra -nvargs for this run
+  local model="$RESULTS/model_${label}.csv"
+  local acc_file="$RESULTS/acc_${label}.txt"
+
+  echo ""
+  echo "--- Training: $label ---"
+  java --add-modules=jdk.incubator.vector -jar "$JAR" \
+      -f "$SCRIPTS/fedavg_dp.dml" \
+      -nvargs $COMMON_ARGS $extra out="$model" \
+      2>&1 | tee "$RESULTS/train_${label}.log" | grep -E "FedAvg|error|Error" || true
+
+  echo "  Evaluating …"
+  java --add-modules=jdk.incubator.vector -jar "$JAR" \
+      -f "$SCRIPTS/eval.dml" \
+      -nvargs data_dir="$DATA" model_path="$model" out_acc="$acc_file" \
+      2>&1 | grep "Accuracy:"
+}
+
+# ── Non-private baseline ──────────────────────────────────────────────────
+run_one "baseline" "private=0 epsilon=9999"
+
+# ── DP runs ───────────────────────────────────────────────────────────────
+for EPS in 0.5 1 4 8; do
+  LABEL="eps_${EPS}"
+  run_one "$LABEL" "private=1 epsilon=${EPS}"
+done
+
+# ── Stop workers ──────────────────────────────────────────────────────────
+echo ""
+echo "=== Stopping federated workers ==="
+bash "$SCRIPTS/stop_workers.sh"
+
+echo ""
+echo "All runs complete. Logs and model files in $RESULTS/"
+
diff --git a/benchmark/scripts/start_workers.sh b/benchmark/scripts/start_workers.sh
new file mode 100755
index 00000000000..d41e7a73944
--- /dev/null
+++ b/benchmark/scripts/start_workers.sh
@@ -0,0 +1,45 @@
+#!/usr/bin/env bash
+#-------------------------------------------------------------
+#
+# Licensed to the Apache Software Foundation (ASF) under one
+# or more contributor license agreements.  See the NOTICE file
+# distributed with this work for additional information
+# regarding copyright ownership.  The ASF licenses this file
+# to you under the Apache License, Version 2.0 (the
+# "License"); you may not use this file except in compliance
+# with the License.  You may obtain a copy of the License at
+#
+#   http://www.apache.org/licenses/LICENSE-2.0
+#
+# Unless required by applicable law or agreed to in writing,
+# software distributed under the License is distributed on an
+# "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
+# KIND, either express or implied.  See the License for the
+# specific language governing permissions and limitations
+# under the License.
+#
+#-------------------------------------------------------------
+
+# Start 4 local SystemDS federated workers on ports 8301-8304.
+# Each worker is given the absolute path to its data shard directory.
+set -e
+REPO_ROOT="$(cd "$(dirname "$0")/../.." && pwd)"
+JAR="$REPO_ROOT/target/SystemDS.jar"
+DATA_DIR="$REPO_ROOT/benchmark/data"
+LOG_DIR="$REPO_ROOT/benchmark/results"
+mkdir -p "$LOG_DIR"
+
+for i in 1 2 3 4; do
+  PORT=$((8300 + i))
+  echo "Starting worker $i on port $PORT …"
+  java --add-modules=jdk.incubator.vector -jar "$JAR" \
+      -w "$PORT" \
+      > "$LOG_DIR/worker${i}.log" 2>&1 &
+  echo $! > "$LOG_DIR/worker${i}.pid"
+done
+
+# Give workers time to bind their ports.
+sleep 3
+echo "Workers running. PIDs:"
+for i in 1 2 3 4; do cat "$LOG_DIR/worker${i}.pid"; done
+
diff --git a/benchmark/scripts/stop_workers.sh b/benchmark/scripts/stop_workers.sh
new file mode 100755
index 00000000000..bccd9a18486
--- /dev/null
+++ b/benchmark/scripts/stop_workers.sh
@@ -0,0 +1,32 @@
+#!/usr/bin/env bash
+#-------------------------------------------------------------
+#
+# Licensed to the Apache Software Foundation (ASF) under one
+# or more contributor license agreements.  See the NOTICE file
+# distributed with this work for additional information
+# regarding copyright ownership.  The ASF licenses this file
+# to you under the Apache License, Version 2.0 (the
+# "License"); you may not use this file except in compliance
+# with the License.  You may obtain a copy of the License at
+#
+#   http://www.apache.org/licenses/LICENSE-2.0
+#
+# Unless required by applicable law or agreed to in writing,
+# software distributed under the License is distributed on an
+# "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
+# KIND, either express or implied.  See the License for the
+# specific language governing permissions and limitations
+# under the License.
+#
+#-------------------------------------------------------------
+
+LOG_DIR="$(cd "$(dirname "$0")/../../benchmark/results" && pwd)"
+for i in 1 2 3 4; do
+  PID_FILE="$LOG_DIR/worker${i}.pid"
+  if [ -f "$PID_FILE" ]; then
+    PID=$(cat "$PID_FILE")
+    kill "$PID" 2>/dev/null && echo "Stopped worker $i (PID $PID)" || true
+    rm -f "$PID_FILE"
+  fi
+done
+
diff --git a/src/main/java/org/apache/sysds/common/Builtins.java b/src/main/java/org/apache/sysds/common/Builtins.java
index f5719641df7..6c3971f499f 100644
--- a/src/main/java/org/apache/sysds/common/Builtins.java
+++ b/src/main/java/org/apache/sysds/common/Builtins.java
@@ -127,6 +127,9 @@ public enum Builtins {
 	SETDIFF("setdiff", true),
 	DIST("dist", true),
 	DMV("dmv", true),
+	DP_GAUSSIAN("dp_gaussian", false, true),
+	DP_LAPLACE("dp_laplace", false, true),
+	DP_SET_BUDGET("dp_set_budget", false),
 	DROP_INVALID_TYPE("dropInvalidType", false),
 	DROP_INVALID_LENGTH("dropInvalidLength", false),
 	EIGEN("eigen", false, ReturnType.MULTI_RETURN),
diff --git a/src/main/java/org/apache/sysds/common/InstructionType.java b/src/main/java/org/apache/sysds/common/InstructionType.java
index e0e77c46c59..ed795d038e2 100644
--- a/src/main/java/org/apache/sysds/common/InstructionType.java
+++ b/src/main/java/org/apache/sysds/common/InstructionType.java
@@ -63,6 +63,7 @@ public enum InstructionType {
 	MMChain,
 	Union,
 	EINSUM,
+	DPBuiltin,
 
 	//SP Types
 	MAPMM,
diff --git a/src/main/java/org/apache/sysds/common/Opcodes.java b/src/main/java/org/apache/sysds/common/Opcodes.java
index 9a894dde13b..c614fb50460 100644
--- a/src/main/java/org/apache/sysds/common/Opcodes.java
+++ b/src/main/java/org/apache/sysds/common/Opcodes.java
@@ -194,6 +194,10 @@ public enum Opcodes {
 	LIST("list", InstructionType.BuiltinNary),
 	EINSUM("einsum", InstructionType.BuiltinNary),
 
+	//DP built-in functions
+	DP_GAUSSIAN("dp_gaussian", InstructionType.DPBuiltin),
+	DP_LAPLACE("dp_laplace", InstructionType.DPBuiltin),
+
 	//Parametrized builtin functions
 	AUTODIFF("autoDiff", InstructionType.ParameterizedBuiltin),
 	CONTAINS("contains", InstructionType.ParameterizedBuiltin),
diff --git a/src/main/java/org/apache/sysds/common/Types.java b/src/main/java/org/apache/sysds/common/Types.java
index 624c9eed3c6..fe6ac0f365c 100644
--- a/src/main/java/org/apache/sysds/common/Types.java
+++ b/src/main/java/org/apache/sysds/common/Types.java
@@ -806,10 +806,14 @@ public static ReOrgOp valueOfByOpcode(String opcode) {
 	
 	/** Parameterized operations that require named variable arguments */
 	public enum ParamBuiltinOp {
-		AUTODIFF, CDF, CONTAINS, INVALID, INVCDF, GROUPEDAGG, RMEMPTY, REPLACE, REXPAND,
-		LOWER_TRI, UPPER_TRI,
-		TRANSFORMAPPLY, TRANSFORMDECODE, TRANSFORMCOLMAP, TRANSFORMMETA,
-		TOKENIZE, TOSTRING, LIST, PARAMSERV
+		AUTODIFF, CDF, CONTAINS,
+		// DP_GAUSSIAN and DP_LAPLACE reuse this Hop/Lop family (ParameterizedBuiltinOp, ParameterizedBuiltin Lop)
+		// but route to a distinct CPType.DPBuiltin/DPBuiltinCPInstruction at the CP-instruction layer instead
+		// of ParameterizedBuiltinCPInstruction (see Opcodes.DP_GAUSSIAN, Opcodes.DP_LAPLACE).
+		DP_GAUSSIAN, DP_LAPLACE,
+		INVALID, INVCDF, GROUPEDAGG, RMEMPTY, REPLACE, REXPAND,
+		LOWER_TRI, UPPER_TRI, TRANSFORMAPPLY, TRANSFORMDECODE, TRANSFORMCOLMAP, TRANSFORMMETA, TOKENIZE, TOSTRING, LIST,
+		PARAMSERV
 	}
 	
 	/** Deep Neural Network specific operations */
diff --git a/src/main/java/org/apache/sysds/hops/ParameterizedBuiltinOp.java b/src/main/java/org/apache/sysds/hops/ParameterizedBuiltinOp.java
index 61a4b8b8f91..ae8742657c2 100644
--- a/src/main/java/org/apache/sysds/hops/ParameterizedBuiltinOp.java
+++ b/src/main/java/org/apache/sysds/hops/ParameterizedBuiltinOp.java
@@ -182,7 +182,7 @@ public Lop constructLops()
 			}
 			case CONTAINS:
 			case CDF:
-			case INVCDF: 
+			case INVCDF:
 			case REPLACE:
 			case LOWER_TRI:
 			case UPPER_TRI:
@@ -194,10 +194,12 @@ public Lop constructLops()
 			case TOSTRING:
 			case PARAMSERV:
 			case LIST:
-			case AUTODIFF:{
-				ParameterizedBuiltin pbilop = new ParameterizedBuiltin(
-					inputlops, _op, getDataType(), getValueType(), et);
-				if( isMultiThreadedOpType() )
+			case AUTODIFF:
+			case DP_GAUSSIAN:
+			case DP_LAPLACE: {
+				ParameterizedBuiltin pbilop = new ParameterizedBuiltin(inputlops, _op, getDataType(), getValueType(),
+					et);
+				if(isMultiThreadedOpType())
 					pbilop.setNumThreads(OptimizerUtils.getConstrainedNumThreads(_maxNumThreads));
 				setOutputDimensions(pbilop);
 				setLineNumbers(pbilop);
@@ -688,7 +690,17 @@ else if( _op == ParamBuiltinOp.TRANSFORMAPPLY ) {
 				return new MatrixCharacteristics(dc.getRows(), dc.getCols(), -1, dc.getLength());
 			}
 		}
-		
+		else if(_op == ParamBuiltinOp.DP_GAUSSIAN || _op == ParamBuiltinOp.DP_LAPLACE) {
+			if(dc.dimsKnown()) {
+				Hop query = getParameterHop("query");
+				String queryVal = (query instanceof LiteralOp) ? ((LiteralOp) query).getStringValue() : null;
+				if("colMeans".equals(queryVal) || "colSums".equals(queryVal))
+					ret = new MatrixCharacteristics(1, dc.getCols(), -1, dc.getCols());
+				else if("identity".equals(queryVal))
+					ret = new MatrixCharacteristics(dc.getRows(), dc.getCols(), -1, dc.getLength());
+			}
+		}
+
 		return ret;
 	}
 	@Override
@@ -755,10 +767,10 @@ && getTargetHop().areDimsBelowThreshold() ) {
 		// 2. For paramserv function, always be CP mode so that
 		// the parameter server could have a central instruction
 		// to determine the local or remote workers
-		if (_op == ParamBuiltinOp.TRANSFORMCOLMAP || _op == ParamBuiltinOp.TRANSFORMMETA
-				|| _op == ParamBuiltinOp.TOSTRING || _op == ParamBuiltinOp.LIST
-				|| _op == ParamBuiltinOp.CDF || _op == ParamBuiltinOp.INVCDF
-				|| _op == ParamBuiltinOp.PARAMSERV) {
+		if(_op == ParamBuiltinOp.TRANSFORMCOLMAP || _op == ParamBuiltinOp.TRANSFORMMETA ||
+			_op == ParamBuiltinOp.TOSTRING || _op == ParamBuiltinOp.LIST || _op == ParamBuiltinOp.CDF ||
+			_op == ParamBuiltinOp.INVCDF || _op == ParamBuiltinOp.PARAMSERV || _op == ParamBuiltinOp.DP_GAUSSIAN ||
+			_op == ParamBuiltinOp.DP_LAPLACE) {
 			_etype = ExecType.CP;
 		}
 
@@ -940,8 +952,14 @@ public boolean compare( Hop that )
 	{
 		if( !(that instanceof ParameterizedBuiltinOp) )
 			return false;
-		
-		ParameterizedBuiltinOp that2 = (ParameterizedBuiltinOp)that;	
+
+		// NOTE: dp_gaussian, dp_laplace draw fresh random noise on every call and record a
+		// privacy-budget charge as a side effect (see DPBuiltinCPInstruction), so two
+		// syntactically identical calls must never be merged into one execution.
+		if( _op == ParamBuiltinOp.DP_GAUSSIAN || _op == ParamBuiltinOp.DP_LAPLACE )
+			return false;
+
+		ParameterizedBuiltinOp that2 = (ParameterizedBuiltinOp)that;
 		boolean ret = (_op == that2._op
 					  && _paramIndexMap!=null && that2._paramIndexMap!=null
 					  && _paramIndexMap.size() == that2._paramIndexMap.size()
diff --git a/src/main/java/org/apache/sysds/lops/ParameterizedBuiltin.java b/src/main/java/org/apache/sysds/lops/ParameterizedBuiltin.java
index 3604121aac8..b71fdf3423b 100644
--- a/src/main/java/org/apache/sysds/lops/ParameterizedBuiltin.java
+++ b/src/main/java/org/apache/sysds/lops/ParameterizedBuiltin.java
@@ -204,7 +204,19 @@ public String getInstructions(String output)
 				compileGenericParamMap(sb, _inputParams);
 				break;
 			}
-			
+			case DP_GAUSSIAN: {
+				sb.append(Opcodes.DP_GAUSSIAN);
+				sb.append(OPERAND_DELIMITOR);
+				compileGenericParamMap(sb, _inputParams);
+				break;
+			}
+			case DP_LAPLACE: {
+				sb.append(Opcodes.DP_LAPLACE);
+				sb.append(OPERAND_DELIMITOR);
+				compileGenericParamMap(sb, _inputParams);
+				break;
+			}
+
 			default:
 				throw new LopsException(this.printErrorLocation() + "In ParameterizedBuiltin Lop, Unknown operation: " + _operation);
 		}
diff --git a/src/main/java/org/apache/sysds/parser/BuiltinFunctionExpression.java b/src/main/java/org/apache/sysds/parser/BuiltinFunctionExpression.java
index ab0c7993b4e..f14210d9b7c 100644
--- a/src/main/java/org/apache/sysds/parser/BuiltinFunctionExpression.java
+++ b/src/main/java/org/apache/sysds/parser/BuiltinFunctionExpression.java
@@ -2003,8 +2003,31 @@ else if(this.getOpCode() == Builtins.MAX_POOL || this.getOpCode() == Builtins.AV
 				output.setBlocksize (id.getBlocksize());
 				output.setValueType(id.getValueType());
 			}
-			else 
+			else
 				raiseValidateError("Local instruction not allowed in dml script");
+			break;
+		case DP_SET_BUDGET: {
+			checkNumParameters(2);
+			checkScalarParam(getFirstExpr());
+			checkScalarParam(getSecondExpr());
+			// resolved entirely at compile time (see DMLTranslator), so both
+			// arguments must be known before the HOP DAG is built.
+			if (!isConstant(getFirstExpr()) || !isConstant(getSecondExpr()))
+				raiseValidateError(getOpCode() + ": 'epsilon' and 'delta' must be compile-time numeric literals",
+					false, LanguageErrorCodes.INVALID_PARAMETERS);
+			double dpSetBudgetEpsilon = getDoubleValue(getFirstExpr());
+			double dpSetBudgetDelta = getDoubleValue(getSecondExpr());
+			if (dpSetBudgetEpsilon <= 0)
+				raiseValidateError(getOpCode() + ": epsilon must be > 0, got " + dpSetBudgetEpsilon,
+					false, LanguageErrorCodes.INVALID_PARAMETERS);
+			if ((dpSetBudgetDelta <= 0) || (dpSetBudgetDelta >= 1))
+				raiseValidateError(getOpCode() + ": delta must be in (0,1), got " + dpSetBudgetDelta,
+					false, LanguageErrorCodes.INVALID_PARAMETERS);
+			output.setDataType(DataType.SCALAR);
+			output.setValueType(ValueType.FP64);
+			output.setDimensions(0, 0);
+			break;
+		}
 		case COMPRESS:
 		case DECOMPRESS:
 			if(OptimizerUtils.ALLOW_SCRIPT_LEVEL_COMPRESS_COMMAND){
diff --git a/src/main/java/org/apache/sysds/parser/DMLProgram.java b/src/main/java/org/apache/sysds/parser/DMLProgram.java
index 2f69cb7ea0d..eed9f802aee 100644
--- a/src/main/java/org/apache/sysds/parser/DMLProgram.java
+++ b/src/main/java/org/apache/sysds/parser/DMLProgram.java
@@ -36,7 +36,18 @@ public class DMLProgram
 	private ArrayList<StatementBlock> _blocks;
 	private Map<String, FunctionDictionary<FunctionStatementBlock>> _namespaces;
 	private boolean _containsRemoteParfor;
-	
+
+	/**
+	 * Session-wide differential privacy budget resolved at compile time from a dp_set_budget(epsilon, delta)
+	 * call (its arguments must be numeric literals, see BuiltinFunctionExpression). Null until such a call is
+	 * encountered during HOP construction; consulted by ExecutionContext#getDPBudgetAccountant() in place of
+	 * its hardcoded default. This is deliberately a plain field (not a Hop/Lop/Instruction): since Program
+	 * holds a reference back to this DMLProgram (see Program#getDMLProg()), the value survives from
+	 * compile time through to runtime without needing a runtime instruction at all.
+	 */
+	private Double _dpBudgetEpsilon;
+	private Double _dpBudgetDelta;
+
 	public DMLProgram(){
 		_blocks = new ArrayList<>();
 		_namespaces = new HashMap<>();
@@ -67,6 +78,23 @@ public void setContainsRemoteParfor(boolean flag) {
 	public boolean containsRemoteParfor() {
 		return _containsRemoteParfor;
 	}
+
+	public void setDPBudget(double epsilon, double delta) {
+		_dpBudgetEpsilon = epsilon;
+		_dpBudgetDelta = delta;
+	}
+
+	public boolean hasDPBudget() {
+		return _dpBudgetEpsilon != null;
+	}
+
+	public double getDPBudgetEpsilon() {
+		return _dpBudgetEpsilon;
+	}
+
+	public double getDPBudgetDelta() {
+		return _dpBudgetDelta;
+	}
 	
 	public static boolean isInternalNamespace(String namespace) {
 		return DEFAULT_NAMESPACE.equals(namespace)
diff --git a/src/main/java/org/apache/sysds/parser/DMLTranslator.java b/src/main/java/org/apache/sysds/parser/DMLTranslator.java
index a8e1667d049..573e6697c13 100644
--- a/src/main/java/org/apache/sysds/parser/DMLTranslator.java
+++ b/src/main/java/org/apache/sysds/parser/DMLTranslator.java
@@ -2014,6 +2014,8 @@ private Hop processParameterizedBuiltinFunctionExpression(ParameterizedBuiltinFu
 			case TRANSFORMMETA:
 			case PARAMSERV:
 			case AUTODIFF:
+			case DP_GAUSSIAN:
+			case DP_LAPLACE:
 				currBuiltinOp = new ParameterizedBuiltinOp(target.getName(), target.getDataType(),
 					target.getValueType(), ParamBuiltinOp.valueOf(source.getOpCode().name()), paramHops);
 				break;
@@ -2589,6 +2591,25 @@ else if ( sop.equalsIgnoreCase(Opcodes.NOTEQUAL.toString()) )
 		case DECOMPRESS:
 			currBuiltinOp = new UnaryOp(target.getName(), target.getDataType(), ValueType.FP64, OpOp1.DECOMPRESS, expr);
 			break;
+		case DP_SET_BUDGET: {
+			// Resolved entirely at compile time: BuiltinFunctionExpression.validateExpression
+			// already enforced that both arguments are numeric literals, so 'expr'/'expr2' are
+			// guaranteed LiteralOps here. There is deliberately no runtime Hop/Lop/Instruction for
+			// this call - the budget is applied directly to the DMLProgram (reachable later from
+			// ExecutionContext via Program.getDMLProg(), see ExecutionContext.getDPBudgetAccountant())
+			// before any instruction executes, so there is nothing for the DAG linearizer to reorder
+			// or drop as dead code.
+			if (_dmlProg.hasDPBudget())
+				throw new LanguageException(source.getOpCode() + ": dp_set_budget may only be called once per "
+					+ "script (already set to epsilon=" + _dmlProg.getDPBudgetEpsilon()
+					+ ", delta=" + _dmlProg.getDPBudgetDelta() + ")");
+			if (!(expr instanceof LiteralOp) || !(expr2 instanceof LiteralOp))
+				throw new LanguageException(source.getOpCode()
+					+ ": epsilon and delta must be compile-time numeric literals");
+			_dmlProg.setDPBudget(((LiteralOp) expr).getDoubleValue(), ((LiteralOp) expr2).getDoubleValue());
+			currBuiltinOp = expr; // echo epsilon back as confirmation
+			break;
+		}
 		case QUANTIZE_COMPRESS:
 			currBuiltinOp = new BinaryOp(target.getName(), target.getDataType(), target.getValueType(), OpOp2.valueOf(source.getOpCode().name()), expr, expr2);
 			break;
diff --git a/src/main/java/org/apache/sysds/parser/ParameterizedBuiltinFunctionExpression.java b/src/main/java/org/apache/sysds/parser/ParameterizedBuiltinFunctionExpression.java
index 314440628e0..30622c85096 100644
--- a/src/main/java/org/apache/sysds/parser/ParameterizedBuiltinFunctionExpression.java
+++ b/src/main/java/org/apache/sysds/parser/ParameterizedBuiltinFunctionExpression.java
@@ -268,6 +268,11 @@ public void validateExpression(HashMap<String, DataIdentifier> ids, HashMap<Stri
 			validateUnique(output, conditional);
 			break;
 
+		case DP_GAUSSIAN:
+		case DP_LAPLACE:
+			validateDpMechanism(output, conditional);
+			break;
+
 		default: //always unconditional (because unsupported operation)
 			//handle common issue of transformencode
 			if( getOpCode()==Builtins.TRANSFORMENCODE )
@@ -770,6 +775,54 @@ private void validateReplace(DataIdentifier output, boolean conditional) {
 		output.setDimensions(target.getOutput().getDim1(), target.getOutput().getDim2());
 	}
 	
+	private void validateDpMechanism(DataIdentifier output, boolean conditional) {
+		// support X as the sole positional (unnamed) parameter, like lowerTri(X, ...)
+		HashMap<String, Expression> varParams = getVarParams();
+		if( varParams.containsKey(null) )
+			varParams.put("target", varParams.remove(null));
+
+		boolean gaussian = getOpCode() == Builtins.DP_GAUSSIAN;
+		Set<String> valid = gaussian ?
+			CollectionUtils.asSet("target", "query", "sensitivity", "epsilon", "delta") :
+			CollectionUtils.asSet("target", "query", "sensitivity", "epsilon");
+		checkInvalidParameters(getOpCode(), varParams, valid);
+
+		Expression target = getVarParam("target");
+		checkTargetParam(target, conditional);
+		for( String param : valid )
+			if( !param.equals("target") )
+				checkScalarParam(getOpCode().getName(), param, conditional);
+
+		long[] dims = getDPOutputDims(getVarParam("query"),
+			target.getOutput().getDim1(), target.getOutput().getDim2());
+		output.setDataType(DataType.MATRIX);
+		output.setValueType(ValueType.FP64);
+		output.setDimensions(dims[0], dims[1]);
+	}
+
+	/** Output dimensions of T %*% X for the given named query, X being n x d. */
+	private long[] getDPOutputDims(Expression queryExpr, long n, long d) {
+		// dp_gaussian/dp_laplace require the "query" parameter to be a compile-time string literal so that the output shape
+		// (and thus the transformation matrix T built at runtime) is known during validation.
+		if(!(queryExpr instanceof StringIdentifier))
+			raiseValidateError(getOpCode() + ": 'query' must be a string literal", false,
+				LanguageErrorCodes.INVALID_PARAMETERS);
+		String query = ((StringIdentifier) queryExpr).getValue();
+
+		switch(query) {
+			case "colMeans":
+			case "colSums":
+				return new long[] {1, d};
+			case "identity":
+				return new long[] {n, d};
+			default:
+				raiseValidateError(
+					getOpCode() + ": unknown query type '" + query + "' (expected colMeans, colSums, or identity)",
+					false, LanguageErrorCodes.INVALID_PARAMETERS);
+				return null; // unreachable
+		}
+	}
+
 	private void checkScalarParam(String group, String param, boolean conditional) {
 		Expression eparam = getVarParam(param);
 		if( eparam==null ) {
diff --git a/src/main/java/org/apache/sysds/runtime/controlprogram/context/ExecutionContext.java b/src/main/java/org/apache/sysds/runtime/controlprogram/context/ExecutionContext.java
index 67cda352a73..0055b9b521d 100644
--- a/src/main/java/org/apache/sysds/runtime/controlprogram/context/ExecutionContext.java
+++ b/src/main/java/org/apache/sysds/runtime/controlprogram/context/ExecutionContext.java
@@ -28,6 +28,7 @@
 import org.apache.sysds.conf.ConfigurationManager;
 import org.apache.sysds.hops.OptimizerUtils;
 import org.apache.sysds.hops.fedplanner.FTypes.FType;
+import org.apache.sysds.parser.DMLProgram;
 import org.apache.sysds.runtime.DMLRuntimeException;
 import org.apache.sysds.runtime.controlprogram.LocalVariableMap;
 import org.apache.sysds.runtime.controlprogram.Program;
@@ -62,6 +63,7 @@
 import org.apache.sysds.runtime.meta.MetaData;
 import org.apache.sysds.runtime.meta.MetaDataFormat;
 import org.apache.sysds.runtime.util.HDFSTool;
+import org.apache.sysds.runtime.privacy.dp.DPBudgetAccountant;
 import org.apache.sysds.utils.Statistics;
 
 import java.util.ArrayList;
@@ -90,6 +92,8 @@ public class ExecutionContext {
 
 	protected SEALClient _seal_client;
 
+	private DPBudgetAccountant _dpBudgetAccountant = null;
+
 	//parfor temporary functions (created by eval)
 	protected Set<String> _fnNames;
 
@@ -144,6 +148,21 @@ public void setLineage(Lineage lineage) {
 		_lineage = lineage;
 	}
 
+	/**
+	 * Returns the session-scoped {@link DPBudgetAccountant}, lazily initialised on first use. If the DML script called
+	 * dp_set_budget(epsilon, delta) with compile-time literal arguments, that value (resolved onto the
+	 * DMLProgram during HOP construction - see DMLTranslator's DP_SET_BUDGET case) is used
+	 * instead of the hardcoded defaults.
+	 */
+	public DPBudgetAccountant getDPBudgetAccountant() {
+		if(_dpBudgetAccountant == null) {
+			DMLProgram dmlProg = (_prog != null) ? _prog.getDMLProg() : null;
+			_dpBudgetAccountant = (dmlProg != null && dmlProg.hasDPBudget()) ? new DPBudgetAccountant(
+				dmlProg.getDPBudgetEpsilon(), dmlProg.getDPBudgetDelta()) : new DPBudgetAccountant();
+		}
+		return _dpBudgetAccountant;
+	}
+
 	public boolean isAutoCreateVars() {
 		return _autoCreateVars;
 	}
diff --git a/src/main/java/org/apache/sysds/runtime/instructions/CPInstructionParser.java b/src/main/java/org/apache/sysds/runtime/instructions/CPInstructionParser.java
index 92e11b425dd..ca6baf058b4 100644
--- a/src/main/java/org/apache/sysds/runtime/instructions/CPInstructionParser.java
+++ b/src/main/java/org/apache/sysds/runtime/instructions/CPInstructionParser.java
@@ -65,6 +65,7 @@
 import org.apache.sysds.runtime.instructions.cp.UnaryCPInstruction;
 import org.apache.sysds.runtime.instructions.cp.VariableCPInstruction;
 import org.apache.sysds.runtime.instructions.cp.UnionCPInstruction;
+import org.apache.sysds.runtime.instructions.cp.DPBuiltinCPInstruction;
 import org.apache.sysds.runtime.instructions.cp.EinsumCPInstruction;
 import org.apache.sysds.runtime.instructions.cpfile.MatrixIndexingCPFileInstruction;
 
@@ -226,7 +227,10 @@ public static CPInstruction parseSingleInstruction ( InstructionType cptype, Str
 			
 			case EINSUM:
 				return EinsumCPInstruction.parseInstruction(str);
-				
+
+			case DPBuiltin:
+				return DPBuiltinCPInstruction.parseInstruction(str);
+
 			default:
 				throw new DMLRuntimeException("Invalid CP Instruction Type: " + cptype );
 		}
diff --git a/src/main/java/org/apache/sysds/runtime/instructions/cp/CPInstruction.java b/src/main/java/org/apache/sysds/runtime/instructions/cp/CPInstruction.java
index b8d84ca3898..e45e17ac175 100644
--- a/src/main/java/org/apache/sysds/runtime/instructions/cp/CPInstruction.java
+++ b/src/main/java/org/apache/sysds/runtime/instructions/cp/CPInstruction.java
@@ -42,14 +42,11 @@ public enum CPType {
 		AggregateUnary, AggregateBinary, AggregateTernary,
 		Unary, Binary, Ternary, Quaternary, BuiltinNary, Ctable,
 		MultiReturnParameterizedBuiltin, ParameterizedBuiltin, MultiReturnBuiltin, MultiReturnComplexMatrixBuiltin,
-		Builtin, Reorg, Variable, FCall, Append, Rand, QSort, QPick, Local,
-		MatrixIndexing, MMTSJ, PMMJ, MMChain, Reshape, Partition, Compression, DeCompression, SpoofFused,
-		StringInit, CentralMoment, Covariance, UaggOuterChain, Dnn, Sql, Prefetch, Broadcast, TrigRemote,
-		EvictLineageCache, EINSUM,
-		NoOp,
-		Union,
-		QuantizeCompression
-	 }
+		Builtin, Reorg, Variable, FCall, Append, Rand, QSort, QPick, Local, MatrixIndexing, MMTSJ, PMMJ, MMChain,
+		Reshape, Partition, Compression, DeCompression, SpoofFused, StringInit, CentralMoment, Covariance,
+		UaggOuterChain, Dnn, Sql, Prefetch, Broadcast, TrigRemote, EvictLineageCache, EINSUM, NoOp, Union,
+		QuantizeCompression, DPBuiltin
+	}
 
 	protected final CPType _cptype;
 	protected final boolean _requiresLabelUpdate;
diff --git a/src/main/java/org/apache/sysds/runtime/instructions/cp/DPBuiltinCPInstruction.java b/src/main/java/org/apache/sysds/runtime/instructions/cp/DPBuiltinCPInstruction.java
new file mode 100755
index 00000000000..96763c11b92
--- /dev/null
+++ b/src/main/java/org/apache/sysds/runtime/instructions/cp/DPBuiltinCPInstruction.java
@@ -0,0 +1,426 @@
+/*
+ * Licensed to the Apache Software Foundation (ASF) under one
+ * or more contributor license agreements.  See the NOTICE file
+ * distributed with this work for additional information
+ * regarding copyright ownership.  The ASF licenses this file
+ * to you under the Apache License, Version 2.0 (the
+ * "License"); you may not use this file except in compliance
+ * with the License.  You may obtain a copy of the License at
+ *
+ *   http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing,
+ * software distributed under the License is distributed on an
+ * "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
+ * KIND, either express or implied.  See the License for the
+ * specific language governing permissions and limitations
+ * under the License.
+ */
+
+package org.apache.sysds.runtime.instructions.cp;
+
+import org.apache.commons.math3.distribution.LaplaceDistribution;
+import org.apache.commons.math3.distribution.NormalDistribution;
+import org.apache.commons.math3.random.Well1024a;
+
+import org.apache.sysds.runtime.DMLRuntimeException;
+import org.apache.sysds.runtime.controlprogram.context.ExecutionContext;
+import org.apache.sysds.runtime.functionobjects.Multiply;
+import org.apache.sysds.runtime.instructions.InstructionUtils;
+import org.apache.sysds.runtime.matrix.data.LibMatrixMult;
+import org.apache.sysds.runtime.matrix.data.LibMatrixReorg;
+import org.apache.sysds.runtime.matrix.data.MatrixBlock;
+import org.apache.sysds.runtime.matrix.operators.RightScalarOperator;
+import org.apache.sysds.runtime.privacy.dp.DPBudgetAccountant;
+
+import java.util.LinkedHashMap;
+
+/**
+ * CP instruction for differential-privacy release of a linear query over the original matrix.
+ *
+ * DML syntax (raw-matrix form):
+ * result = dp_laplace(X, query="colMeans", sensitivity=1.0, epsilon=0.5)
+ * result = dp_gaussian(X, query="colMeans", sensitivity=1.0, epsilon=0.5, delta=1e-5)
+ *
+ * The instruction receives the original n x d} matrix X, builds a transformation matrix T
+ * (k x n) from the named query (see {@link #buildTransform}), and returns a noisy release of
+ * T %*% X. The noise is not added as a separate elementwise pass over a materialised aggregate: it is injected
+ * by augmenting T with an identity block and X with the noise matrix, so that the noisy release is the
+ * result of a single {@link LibMatrixMult#matrixMult} call (see {@link #processInstruction} for the derivation).
+ *
+ * Sensitivity norm: sensitivity is not interchangeable between the two builtins. dp_laplace calibrates
+ * its noise scale to the L1 sensitivity of T %*% X to a single-record change; dp_gaussian calibrates
+ * its stdev to the L2 sensitivity. For a scalar release (e.g. query="colMeans" on single-column X) the two
+ * norms coincide, but for a vector- or matrix-valued release they generally differ - the caller is responsible for
+ * supplying the norm matching the builtin invoked (see {@link #sensitivityOf}).
+ *
+ * The {@link #sensitivityOf} method is deliberately separated from the noise-scale computation. It currently returns
+ * the caller-supplied constant. A future rewrite pass could replace the body of this single method with a static
+ * analysis that derives sensitivity from T's column norms and a declared per-record bound on X; every
+ * other line in this class would stay unchanged.
+ */
+public class DPBuiltinCPInstruction extends ComputationCPInstruction {
+
+	// -----------------------------------------------------------------------
+	// Constants
+	// -----------------------------------------------------------------------
+
+	/** Opcode registered in Builtins and CPInstructionParser. */
+	public static final String OPCODE_GAUSSIAN = "dp_gaussian";
+	public static final String OPCODE_LAPLACE = "dp_laplace";
+
+	// -----------------------------------------------------------------------
+	// Fields
+	// -----------------------------------------------------------------------
+
+	/**
+	 * Named parameters extracted from the serialised instruction string. Keys: "target", "query", "sensitivity",
+	 * "epsilon", "delta" (Gaussian only).
+	 *
+	 * Using the same LinkedHashMap<String,String> convention as ParameterizedBuiltinCPInstruction so that
+	 * CPInstructionParser can call the shared constructParameterMap() helper unchanged.
+	 */
+	private final LinkedHashMap<String, String> _params;
+
+	private static final NormalDistribution normal = new NormalDistribution();
+
+
+	// -----------------------------------------------------------------------
+	// Constructor (private – use parseInstruction)
+	// -----------------------------------------------------------------------
+
+	private DPBuiltinCPInstruction(CPOperand input, CPOperand output, String opcode, String istr,
+		LinkedHashMap<String, String> params) {
+		super(CPType.DPBuiltin, null, input, null, output, opcode, istr);
+		_params = params;
+	}
+
+	// -----------------------------------------------------------------------
+	// Static factory / parser
+	// -----------------------------------------------------------------------
+
+	/**
+	 * Reconstructs a DPBuiltinCPInstruction from its serialised instruction string produced by the LOP layer.
+	 *
+	 * Expected format (OPERAND_DELIM = '\u00b0'):
+	 * dp_gaussian°target=mVar1·MATRIX·FP64°query=colMeans·SCALAR·STRING·true
+	 * °sensitivity=1.0·SCALAR·FP64·true°epsilon=0.5·SCALAR·FP64·true °delta=1e-5·SCALAR·FP64·true°_mVar2·MATRIX·FP64
+	 *
+	 * The first token is always the opcode; the last token is always the output operand; the tokens in between are
+	 * key=value pairs. This matches the convention used by ParameterizedBuiltinCPInstruction exactly.
+	 */
+	public static DPBuiltinCPInstruction parseInstruction(String str) {
+		String[] parts = InstructionUtils.getInstructionPartsWithValueType(str);
+		InstructionUtils.checkNumFields(parts, 5, 6); // laplace=5, gaussian=6
+		String opcode = parts[0];
+
+		// Output operand is always the last token.
+		CPOperand output = new CPOperand(parts[parts.length - 1]);
+
+		// The "target" parameter holds the variable name of the input matrix.
+		// ParameterizedBuiltinCPInstruction.constructParameterMap strips the
+		// type suffixes and returns bare key=value pairs.
+		LinkedHashMap<String, String> params = ParameterizedBuiltinCPInstruction.constructParameterMap(parts);
+
+		// The target CPOperand is needed by ComputationCPInstruction's
+		// getInputs() / getLineageItem() machinery.
+		CPOperand input = new CPOperand(params.get("target"), org.apache.sysds.common.Types.ValueType.FP64,
+			org.apache.sysds.common.Types.DataType.MATRIX);
+
+		// Validate required keys.
+		if(!params.containsKey("query"))
+			throw new DMLRuntimeException(opcode + ": missing 'query'");
+		if(!params.containsKey("sensitivity"))
+			throw new DMLRuntimeException(opcode + ": missing 'sensitivity'");
+		if(!params.containsKey("epsilon"))
+			throw new DMLRuntimeException(opcode + ": missing 'epsilon'");
+		if(opcode.equals(OPCODE_GAUSSIAN) && !params.containsKey("delta"))
+			throw new DMLRuntimeException(opcode + ": missing 'delta'");
+
+		return new DPBuiltinCPInstruction(input, output, opcode, str, params);
+	}
+
+	// -----------------------------------------------------------------------
+	// Core execution
+	// -----------------------------------------------------------------------
+
+	/**
+	 * Executes the DP release.
+	 *
+	 * - Read the original {@link MatrixBlock} X from the variable table.
+	 * - Build the transformation matrix T (k x n) from query (see {@link #buildTransform}).
+	 * - Determine sensitivity via {@link #sensitivityOf}.
+	 * - Generate a noise {@link MatrixBlock} shaped k x d.
+	 * - Fuse T %*% X + noise into a single {@link LibMatrixMult#matrixMult} call (see below).
+	 * - Record the release with the session-scoped {@link DPBudgetAccountant}; throw if budget is exhausted.
+	 * - Write the noisy block back to the variable table and release the input pin.
+	 *
+	 * Fusion derivation: for T (k x n), X (n x d) and noise N (k x d),
+	 * let T' = [T | I_k] (k x (n+k)) and X' = [X ; N] ((n+k) x d). Then
+	 * T' %*% X' = T %*% X + I_k %*% N = T %*% X + N, computed as one matrix multiply instead of a multiply
+	 * followed by a separate elementwise add.
+	 */
+	@Override
+	public void processInstruction(ExecutionContext ec) {
+
+		// ── 1. Read original input matrix X ─────────────────────────────────
+		// getMatrixInput pins the block in memory and increments the
+		// reference count; we must call releaseMatrixInput afterwards.
+		MatrixBlock X = ec.getMatrixInput(_params.get("target"));
+
+		// ── 2. Parse DP parameters ──────────────────────────────────────────
+		double epsilon = parsePositiveDouble("epsilon");
+		double delta = instOpcode.equals(OPCODE_GAUSSIAN) ? parsePositiveDouble("delta") : 0.0;
+		String query = _params.get("query");
+
+		// ── 3. Build the transformation matrix T (k x n) ────────────────────
+		MatrixBlock T = buildTransform(query, X.getNumRows());
+
+		// ── 4. Determine sensitivity (caller-supplied constant) ─────────────
+		double sensitivity = sensitivityOf(T);
+
+		// ── 5. Generate noise shaped like the release T %*% X (k x d) ───────
+		MatrixBlock noiseBlock = generateNoise(T.getNumRows(), X.getNumColumns(), sensitivity, epsilon, delta);
+
+		// ── 6. Fuse T %*% X + noise into a single matrix multiply ───────────
+		MatrixBlock Ik = identity(T.getNumRows());
+		MatrixBlock Tp = T.append(Ik, null, true); // [T | I_k]
+		MatrixBlock Xp = X.append(noiseBlock, null, false); // [X ; noise]
+		MatrixBlock outBlock = LibMatrixMult.matrixMult(Tp, Xp);
+
+		// ── 7. Record release and enforce budget ────────────────────────────
+		// getDPBudgetAccountant() returns a lazy-initialised DPBudgetAccountant that is
+		// owned by this ExecutionContext (added in a companion EC patch).
+		DPBudgetAccountant accountant = ec.getDPBudgetAccountant();
+		accountant.compose(epsilon, delta, sensitivity); // throws on exhaustion
+
+		// ── 8. Write output and release input pin ───────────────────────────
+		ec.releaseMatrixInput(_params.get("target"));
+		ec.setMatrixOutput(output.getName(), outBlock);
+	}
+
+	// -----------------------------------------------------------------------
+	// Transformation matrix construction
+	// -----------------------------------------------------------------------
+
+	/**
+	 * Builds the k x n transformation matrix T for the given named query, to be left-multiplied against
+	 * the n x d input X as T %*% X.
+	 *
+	 * - "colMeans": T is 1 x n, filled with 1/n - T %*% X is the column-mean row vector.
+	 * - "colSums": T is 1 x n, filled with 1.0 - T %*% X is the column-sum row vector.
+	 * - "identity": T is the n x n identity (built sparsely via {@link #identity}) - T %*% X is X itself,
+	 * i.e. a noisy release of the raw matrix.
+	 *
+	 * Row-wise aggregates (rowMeans/rowSums) reduce across the feature axis of X, i.e. they are
+	 * naturally X %*% T' (right-multiply), not T %*% X, so they are intentionally not supported here.
+	 */
+	private static MatrixBlock buildTransform(String query, int n) {
+		switch(query) {
+			case "colMeans": {
+				MatrixBlock T = new MatrixBlock(1, n, false);
+				T.allocateDenseBlock();
+				double v = 1.0 / n;
+				for(int c = 0; c < n; c++)
+					T.set(0, c, v);
+				T.recomputeNonZeros();
+				return T;
+			}
+			case "colSums": {
+				MatrixBlock T = new MatrixBlock(1, n, false);
+				T.allocateDenseBlock();
+				for(int c = 0; c < n; c++)
+					T.set(0, c, 1.0);
+				T.recomputeNonZeros();
+				return T;
+			}
+			case "identity":
+				return identity(n);
+			default:
+				throw new DMLRuntimeException("dp_laplace/dp_gaussian: unknown query type '" + query
+					+ "' (expected colMeans, colSums, or identity)");
+		}
+	}
+
+	/**
+	 * Builds a k x k identity matrix, sparsely, by reusing the existing {@link LibMatrixReorg#diag} reorg
+	 * operator (the same runtime path DML's diag() builtin uses to expand a vector into a diagonal matrix).
+	 * Keeps memory O(k) rather than O(k^2), which matters for the query="identity" case where
+	 * k equals the number of rows of X.
+	 */
+	private static MatrixBlock identity(int k) {
+		MatrixBlock ones = new MatrixBlock(k, 1, false);
+		ones.allocateDenseBlock();
+		for(int i = 0; i < k; i++)
+			ones.set(i, 0, 1.0);
+		ones.recomputeNonZeros();
+		return LibMatrixReorg.diag(ones, new MatrixBlock(k, k, true));
+	}
+
+	// -----------------------------------------------------------------------
+	// Sensitivity seam
+	// -----------------------------------------------------------------------
+
+	/**
+	 * Returns the sensitivity of the release T %*% X to a single-record change, in the norm required by the
+	 * mechanism actually invoked: L1 for dp_laplace, L2 for dp_gaussian (see the class Javadoc). The
+	 * two only coincide when the release is scalar.
+	 *
+	 * Returns the caller-supplied literal from the DML script as-is, with no norm conversion or validation - the DML
+	 * author must compute the sensitivity in the correct norm for the builtin they call. A future rewrite pass could
+	 * replace this body with an analysis that derives sensitivity from T's column norms and a declared
+	 * per-record bound on X; no other line in this class would need to change.
+	 *
+	 * @param T the transformation matrix (unused for now; kept as the seam for a future sensitivity-derivation pass)
+	 * @return caller-supplied sensitivity constant, expected to already be in the L1 norm (Laplace) or L2 norm
+	 *         (Gaussian)
+	 */
+	private double sensitivityOf(MatrixBlock T) {
+		return parsePositiveDouble("sensitivity");
+	}
+
+	// -----------------------------------------------------------------------
+	// Noise generation
+	// -----------------------------------------------------------------------
+
+	/**
+	 * Generates a rows x cols noise {@link MatrixBlock} - matching the shape of the release T %*% X -
+	 * filled with samples from the mechanism-appropriate distribution calibrated to (sensitivity,
+	 * epsilon, delta).
+	 *
+	 * Both mechanisms produce a dense block. Sparsity exploitation is left for future work; for the releases targeted
+	 * here (e.g. column means, column sums) the noise is dense regardless.
+	 */
+	private MatrixBlock generateNoise(int rows, int cols, double sensitivity, double epsilon, double delta) {
+
+		MatrixBlock noise;
+
+		if(instOpcode.equals(OPCODE_LAPLACE)) {
+			// Laplace mechanism
+			// For a given epsilon, noise is drawn from the Laplace distribution at
+			// scale b = sensitivity / epsilon
+			noise = fillLaplaceNoise(rows, cols, sensitivity / epsilon);
+		}
+		else {
+			// Gaussian mechanism
+			// For a given epsilon and delta, noise is drawn from the Gaussian distribution
+			// N(0, sigma^2)
+			double sigma = computeGaussianSigma(sensitivity, epsilon, delta);
+			noise = fillGaussianNoise(rows, cols, sigma);
+		}
+
+		return noise;
+	}
+
+	// /**
+	//  * Gaussian mechanism: calibrate sigma for (epsilon, delta)-DP.
+	//  * Classical Gaussian mechanism calibration
+	//  *
+    //  * @param sensitivity L2 sensitivity
+    //  * @param epsilon     target epsilon
+    //  * @param delta       target delta
+    //  * @return optimal sigma
+	//  */
+	// public static double getGaussianSigma(double sensitivity, double epsilon, double delta) {
+	// 	double sigma = sensitivity * Math.sqrt(2.0 * Math.log(1.25 / delta)) / epsilon;
+	// 	return sigma;
+	// }
+
+    /**
+     * Compute the optimal sigma for the Analytic Gaussian Mechanism (Balle & Wang 2018).
+     * Returns the smallest sigma such that the Gaussian mechanism is (epsilon, delta)-DP.
+     *
+     * @param sensitivity L2 sensitivity
+     * @param epsilon     target epsilon
+     * @param delta       target delta
+     * @return optimal sigma
+     */
+	public static double computeGaussianSigma(double sensitivity, double epsilon, double delta) {
+
+		// Upper bound: classical Gaussian mechanism (loose but safe)
+		double sigmaHigh = (sensitivity * Math.sqrt(2 * Math.log(1.25 / delta))) / epsilon;
+		double sigmaLow = 1e-12;
+
+		for (int i = 0; i < 100; i++) {
+			double sigmaMid = 0.5 * (sigmaLow + sigmaHigh);
+
+			if (deltaUpperBound(sigmaMid, epsilon, delta, sensitivity) > delta) {
+				sigmaLow = sigmaMid;
+			} else {
+				sigmaHigh = sigmaMid;
+			}
+		}
+
+		return sigmaHigh;
+	}
+
+    /**
+     * Analytic Gaussian DP inequality from Balle & Wang (2018).
+     */
+	private static double deltaUpperBound(double sigma, double epsilon, double delta, double sensitivity) {
+		double c = sensitivity / (2 * sigma);
+
+		double term1 = normal.cumulativeProbability(c - epsilon * sigma / sensitivity);
+		double term2 = Math.exp(epsilon) *
+				normal.cumulativeProbability(-c - epsilon * sigma / sensitivity);
+
+		return term1 - term2;
+	}
+
+	/**
+	 * Fills block with i.i.d. Laplace(0, scale) samples.
+	 *
+	 * Draws from commons-math3's {@link LaplaceDistribution} (its inverse-CDF sampling, tested independently
+	 * of this class) seeded by {@link Well1024a}, the same long-period equidistributed generator
+	 * {@link org.apache.sysds.runtime.matrix.data.LibMatrixDatagen} uses for DML's rand() builtin.
+	 */
+	private static MatrixBlock fillLaplaceNoise(int rows, int cols, double scale) {
+		MatrixBlock noise = new MatrixBlock(rows, cols, false); // dense
+		noise.allocateDenseBlock();
+		LaplaceDistribution laplace = new LaplaceDistribution(new Well1024a(), 0, scale);
+		for(int r = 0; r < rows; r++) {
+			for(int c = 0; c < cols; c++) {
+				noise.set(r, c, laplace.sample());
+			}
+		}
+		return noise;
+	}
+
+	/**
+	 * Generates a rows x cols block of i.i.d. N(0, sigma^2) samples.
+	 *
+	 * Reuses the same Well1024a-seeded, Box-Muller normal generator that backs DML's rand(pdf="normal")
+	 * (see {@link MatrixBlock#randOperations}), so the noise gets the same long-period PRNG and block-parallel generation
+	 * as the rest of SystemDS's random matrix generation. randOperations produces standard N(0,1) samples
+	 * (pdf="normal" ignores min/max), so the sigma scaling is applied afterwards as a scalar multiply.
+	 */
+	private static MatrixBlock fillGaussianNoise(int rows, int cols, double sigma) {
+		MatrixBlock std = MatrixBlock.randOperations(rows, cols, 1.0, 0, 1, "normal", -1);
+		return std.scalarOperations(new RightScalarOperator(Multiply.getMultiplyFnObject(), sigma), null);
+	}
+
+	// -----------------------------------------------------------------------
+	// Helpers
+	// -----------------------------------------------------------------------
+
+	/**
+	 * Parses a parameter value as a positive double.
+	 *
+	 * @throws DMLRuntimeException if the key is absent, unparseable, or non-positive
+	 */
+	private double parsePositiveDouble(String key) {
+		String raw = _params.get(key);
+		if(raw == null)
+			throw new DMLRuntimeException(instOpcode + ": parameter '" + key + "' is missing");
+		double v;
+		try {
+			v = Double.parseDouble(raw);
+		}
+		catch(NumberFormatException e) {
+			throw new DMLRuntimeException(instOpcode + ": parameter '" + key + "' is not a valid number: " + raw);
+		}
+		if(!(v > 0.0))
+			throw new DMLRuntimeException(instOpcode + ": parameter '" + key + "' must be strictly positive, got " + v);
+		return v;
+	}
+}
diff --git a/src/main/java/org/apache/sysds/runtime/privacy/dp/DPBudgetAccountant.java b/src/main/java/org/apache/sysds/runtime/privacy/dp/DPBudgetAccountant.java
new file mode 100644
index 00000000000..c2fd4bf4d9f
--- /dev/null
+++ b/src/main/java/org/apache/sysds/runtime/privacy/dp/DPBudgetAccountant.java
@@ -0,0 +1,241 @@
+/*
+ * Licensed to the Apache Software Foundation (ASF) under one
+ * or more contributor license agreements.  See the NOTICE file
+ * distributed with this work for additional information
+ * regarding copyright ownership.  The ASF licenses this file
+ * to you under the Apache License, Version 2.0 (the
+ * "License"); you may not use this file except in compliance
+ * with the License.  You may obtain a copy of the License at
+ *
+ *   http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing,
+ * software distributed under the License is distributed on an
+ * "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
+ * KIND, either express or implied.  See the License for the
+ * specific language governing permissions and limitations
+ * under the License.
+ */
+
+package org.apache.sysds.runtime.privacy.dp;
+
+import org.apache.sysds.runtime.DMLRuntimeException;
+import org.apache.sysds.runtime.instructions.cp.DPBuiltinCPInstruction;
+
+/**
+ * Session-scoped differential privacy budget accountant.
+ *
+ * Tracks composition of DP releases across the lifetime of a DML script execution. Each call to {@link #compose}
+ * records one release and checks whether the cumulative privacy cost has exceeded the user-specified budget.
+ *
+ * The mechanism type (Laplace vs Gaussian) is inferred from the delta argument passed to {@link #compose}:
+ *
+ * - Laplace (delta == 0): pure epsilon-DP. The budget cost is tracked via basic composition: each release contributes
+ * exactly its epsilon to a running sum. This is the tightest possible bound for pure DP and avoids the looser estimate that
+ * results from routing Laplace through the RDP conversion path (which would introduce an unnecessary delta). Noise scale is
+ * calibrated to L1 sensitivity (see {@link #compose}).
+ * - Gaussian (delta > 0): (epsilon, delta)-DP via Renyi DP composition.
+ * Renyi divergences at a discrete set of orders alpha compose additively; the accumulated sum is converted to (epsilon, delta) at
+ * query time using the formula from Mironov 2017. This is substantially tighter than basic composition for repeated
+ * Gaussian releases, which is the common case in federated learning.
+ *
+ * When both mechanisms are used in the same script the total cost is:
+ * epsilon_total = epsilon_Laplace_sum + epsilon_Gaussian_RDP
+ * This follows from basic composition of a pure-DP mechanism with an approximate-DP mechanism, which is additive in epsilon.
+ *
+ * Renyi orders tracked (Gaussian path) alpha in {2, 4, 8, 16, 32, 64, 128, 256, 512, 1024}. At query time the minimum
+ * converted epsilon across all orders is taken as the tightest available bound.
+ *
+ * Gaussian RDP divergence For the Gaussian mechanism with noise scale sigma and L2 sensitivity delta_f:
+ * D_alpha = alpha * delta_f^2 / (2*sigma^2) sigma
+ * is back-derived from the caller's (epsilon, delta) via the standard calibration formula (see {@link #gaussianSigma}). Note that
+ * sensitivity cancels in the final expression, so the RDP cost depends only on the (epsilon, delta) parameters.
+ *
+ * RDP => (epsilon, delta) conversion (Mironov 2017, Proposition 3):
+ * epsilon(alpha) = R[alpha] + log(1/delta) / (alpha − 1)
+ *
+ * One instance is created per ExecutionContext (lazy init). It is garbage-collected with the context when the
+ * script finishes; no state leaks between script executions or between concurrent scripts.
+ *
+ * Not thread-safe. A single DML script executes instructions sequentially on one thread, so no synchronisation is
+ * needed.
+ *
+ * @see DPBuiltinCPInstruction
+ */
+public class DPBudgetAccountant {
+
+	// -----------------------------------------------------------------------
+	// Renyi orders used for Gaussian composition
+	// -----------------------------------------------------------------------
+
+	private static final double DEFAULT_EPSILON_BUDGET = 1.0;
+
+	private static final double DEFAULT_DELTA = 1e-5;
+
+	/**
+	 * Discrete set of Renyi orders alpha. All must be > 1. Finer grids give tighter bounds; this set covers the range
+	 * relevant for typical ML workloads.
+	 */
+	private static final double[] ORDERS = {2, 4, 8, 16, 32, 64, 128, 256, 512, 1024};
+
+	// -----------------------------------------------------------------------
+	// State
+	// -----------------------------------------------------------------------
+
+	/** Accumulated Renyi divergence at each order (Gaussian releases only). */
+	private final double[] _rdpSum = new double[ORDERS.length];
+
+	/**
+	 * Running sum of pure epsilon from Laplace releases.
+	 *
+	 * Laplace gives pure epsilon-DP (no delta). Basic composition is exact and tighter than the RDP conversion path for Laplace
+	 * (which would introduce an unnecessary delta and produce a looser bound). Each Laplace release adds its epsilon here; the
+	 * total is added directly in {@link #totalEpsilonSpent()}.
+	 */
+	private double _pureEpsilonSum = 0.0;
+
+	/** Total privacy budget (epsilon) for the script execution. */
+	private final double _epsilonBudget;
+
+	/** delta used for the Gaussian RDP-to-(epsilon,delta) conversion. */
+	private final double _delta;
+
+	/** Number of releases recorded so far (for error messages). */
+	private int _releaseCount = 0;
+
+	/** Whether at least one Gaussian release has been recorded. */
+	private boolean _hasGaussianReleases = false;
+
+	// -----------------------------------------------------------------------
+	// Constructors
+	// -----------------------------------------------------------------------
+
+	/**
+	 * Creates an accountant with the given global budget.
+	 *
+	 * Typical usage: the DML script sets the budget once at the top (future work: a
+	 * dp_set_budget(epsilon, delta) built-in), or the accountant is created with defaults and the budget is
+	 * checked after each release.
+	 *
+	 * @param epsilonBudget total epsilon budget for the script execution (must be > 0)
+	 * @param delta         delta used for the Gaussian RDP-to-(epsilon,delta) conversion (must be in (0,1))
+	 */
+	public DPBudgetAccountant(double epsilonBudget, double delta) {
+		if(!(epsilonBudget > 0))
+			throw new DMLRuntimeException("DPBudgetAccountant: epsilonBudget must be > 0, got " + epsilonBudget);
+		if(!(delta > 0 && delta < 1))
+			throw new DMLRuntimeException("DPBudgetAccountant: delta must be in (0,1), got " + delta);
+		_epsilonBudget = epsilonBudget;
+		_delta = delta;
+	}
+
+	/**
+	 * Convenience constructor using a liberal default delta = 1e-5. Suitable when the calling script does not specify delta
+	 * explicitly.
+	 */
+	public DPBudgetAccountant(double epsilonBudget) {
+		this(epsilonBudget, 1e-5);
+	}
+
+	/**
+	 * Default constructor using defaults. Suitable when the calling script does not specify epsilon, delta explicitly.
+	 */
+	public DPBudgetAccountant() {
+		this(DEFAULT_EPSILON_BUDGET, DEFAULT_DELTA);
+	}
+
+	// -----------------------------------------------------------------------
+	// Core API
+	// -----------------------------------------------------------------------
+
+	/**
+	 * Records one DP release and checks the budget.
+	 *
+	 * This method must be called before the result is written to the variable table. If the budget is exhausted it
+	 * throws and the caller's result is discarded, preventing an unaccounted release.
+	 *
+	 * Mechanism selection (see class-level Javadoc for details):
+	 * - delta == 0 => Laplace, pure epsilon-DP basic composition
+	 * - delta > 0 => Gaussian, Renyi DP composition
+	 *
+	 * @param epsilon     per-release epsilon parameter (must be >= 0)
+	 * @param delta       per-release delta parameter (0 for Laplace, >= 0 for Gaussian)
+	 * @param sensitivity sensitivity of the released quantity (must be > 0). The norm depends on the mechanism
+	 *                    selected by delta: callers must supply the L1 sensitivity when
+	 *                    delta == 0 (Laplace), and the L2 sensitivity when delta > 0
+	 *                    (Gaussian). The two coincide for scalar-valued releases but diverge for vector-valued ones, so
+	 *                    passing the wrong norm silently under- or over-calibrates the noise.
+	 * @throws DMLRuntimeException if the cumulative epsilon after this release would exceed the budget
+	 */
+	public void compose(double epsilon, double delta, double sensitivity) {
+		_releaseCount++;
+
+		if(delta == 0.0) {
+			// Laplace: pure epsilon-DP, basic composition - cost is exactly epsilon.
+			_pureEpsilonSum += epsilon;
+		}
+		else {
+			// Gaussian: accumulate Renyi divergence at each order, then convert.
+			_hasGaussianReleases = true;
+			for(int i = 0; i < ORDERS.length; i++) {
+				double sigma = DPBuiltinCPInstruction.computeGaussianSigma(sensitivity, epsilon, delta);
+				_rdpSum[i] += rdpGaussian(ORDERS[i], sensitivity, sigma);
+			}
+		}
+
+		double spentEpsilon = totalEpsilonSpent();
+		if(spentEpsilon > _epsilonBudget) {
+			throw new DMLRuntimeException(String.format(
+				"Privacy budget exhausted after %d release(s): " + "spent epsilon %.6f exceeds budget epsilon = %.6f (delta = %.2e). "
+					+ "Reduce the number of releases or widen the budget.",
+				_releaseCount, spentEpsilon, _epsilonBudget, _delta));
+		}
+	}
+
+	// -----------------------------------------------------------------------
+	// Inspection
+	// -----------------------------------------------------------------------
+
+	/**
+	 * Returns the current total privacy cost as an epsilon value.
+	 *
+	 * Total = Laplace pure-epsilon sum + Gaussian RDP-converted epsilon (clamped to zero when no Gaussian releases have been
+	 * recorded).
+	 */
+	public double totalEpsilonSpent() {
+		if(!_hasGaussianReleases)
+			return _pureEpsilonSum;
+
+		// Take min_alpha(epsilon_alpha) as the current total privacy cost
+		double gaussianEps = Double.MAX_VALUE;
+		for(int i = 0; i < ORDERS.length; i++) {
+			double alpha = ORDERS[i];
+			double eps = _rdpSum[i] + Math.log(1.0 / _delta) / (alpha - 1.0);
+			if(eps < gaussianEps)
+				gaussianEps = eps;
+		}
+		return _pureEpsilonSum + Math.max(gaussianEps, 0.0);
+	}
+
+	/** Returns the remaining epsilon budget (negative if the budget is exceeded). */
+	public double remainingBudget() {
+		return _epsilonBudget - totalEpsilonSpent();
+	}
+
+	/** Returns the number of DP releases recorded so far. */
+	public int releaseCount() {
+		return _releaseCount;
+	}
+
+	// -----------------------------------------------------------------------
+	// Private helpers
+	// -----------------------------------------------------------------------
+
+	/**
+	 * Renyi divergence of order alpha for the Gaussian mechanism (Mironov 2017, Proposition 3, example 2):
+	 * D_alpha = alpha * delta_f^2 / (2 * sigma^2)
+	 */
+	private static double rdpGaussian(double alpha, double sensitivity, double sigma) {
+		return alpha * (sensitivity * sensitivity) / (2.0 * sigma * sigma);
+	}
+}
diff --git a/src/test/java/org/apache/sysds/test/component/cp/DPBuiltinCPInstructionTest.java b/src/test/java/org/apache/sysds/test/component/cp/DPBuiltinCPInstructionTest.java
new file mode 100755
index 00000000000..b07e2201dc9
--- /dev/null
+++ b/src/test/java/org/apache/sysds/test/component/cp/DPBuiltinCPInstructionTest.java
@@ -0,0 +1,343 @@
+/*
+ * Licensed to the Apache Software Foundation (ASF) under one
+ * or more contributor license agreements.  See the NOTICE file
+ * distributed with this work for additional information
+ * regarding copyright ownership.  The ASF licenses this file
+ * to you under the Apache License, Version 2.0 (the
+ * "License"); you may not use this file except in compliance
+ * with the License.  You may obtain a copy of the License at
+ *
+ *   http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing,
+ * software distributed under the License is distributed on an
+ * "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
+ * KIND, either express or implied.  See the License for the
+ * specific language governing permissions and limitations
+ * under the License.
+ */
+
+package org.apache.sysds.test.component.cp;
+
+import org.apache.sysds.runtime.DMLRuntimeException;
+import org.apache.sysds.runtime.instructions.cp.DPBuiltinCPInstruction;
+import org.apache.sysds.runtime.matrix.data.MatrixBlock;
+import org.apache.sysds.runtime.privacy.dp.DPBudgetAccountant;
+
+import java.lang.reflect.Method;
+
+import org.junit.Test;
+import org.junit.Assert;
+
+/**
+ * Tests for DPBuiltinCPInstruction and DPBudgetAccountant.
+ *
+ * The tests are grouped into four levels:
+ * - Unit tests on DPBudgetAccountant - verify composition, conversion, and
+ * budget enforcement in isolation, with no dependency on the full SystemDS runtime.
+ * - Noise distribution tests - verify that the noise blocks generated by the Laplace and
+ * Gaussian mechanisms have statistically correct means and variances
+ * (Kolmogorov-Smirnov style sanity checks).
+ * - DPBuiltinCPInstruction structural tests - exercise parseInstruction()'s required-parameter
+ * validation.
+ * - DML integration tests - run complete DML scripts and verify
+ * end-to-end correctness via the existing AutomatedTestBase machinery.
+ * The DML integration tests require a built SystemDS jar and are separated into a companion class
+ * {@link org.apache.sysds.test.functions.privacy.dp.DPBuiltinDMLTest}.
+ */
+public class DPBuiltinCPInstructionTest {
+
+	private static final double EPS = 1e-9;
+
+	// =======================================================================
+	// 1. DPBudgetAccountant unit tests
+	// =======================================================================
+
+	@Test
+	public void testAccountantInitialisesAtZeroCostThenAcceptsSingleReleases() {
+		DPBudgetAccountant acc = new DPBudgetAccountant(1.0, 1e-5);
+		// No releases yet: total cost should be a large negative number
+		// (conversion formula gives -inf when rdpSum = 0 for all orders),
+		// so remainingBudget() should exceed the budget.
+		Assert.assertTrue("No releases should leave budget intact", acc.remainingBudget() > 0);
+		Assert.assertEquals(0, acc.releaseCount());
+
+		// epsilon=0.5, budget=1.0: one Laplace release should consume < budget.
+		acc.compose(0.5, 0.0, 1.0); // Laplace, sensitivity=1
+		Assert.assertEquals(1, acc.releaseCount());
+		Assert.assertTrue("Single release within budget", acc.totalEpsilonSpent() <= 1.0);
+
+		// Likewise, a single Gaussian release on a fresh accountant should stay within budget.
+		DPBudgetAccountant gaussianAcc = new DPBudgetAccountant(1.0, 1e-5);
+		gaussianAcc.compose(0.5, 1e-5, 1.0); // Gaussian
+		Assert.assertEquals(1, gaussianAcc.releaseCount());
+		Assert.assertTrue("Single Gaussian release within budget", gaussianAcc.totalEpsilonSpent() <= 1.0);
+	}
+
+	@Test(expected = DMLRuntimeException.class)
+	public void testBudgetExhaustionThrows() {
+		// Budget = 0.5, but we try to make 10 releases at epsilon=0.1 each.
+		// After enough releases the budget must be exceeded.
+		DPBudgetAccountant acc = new DPBudgetAccountant(0.5, 1e-5);
+		double prevEpsilonSpent = acc.totalEpsilonSpent();
+		double prevRemainingBudget = acc.remainingBudget();
+		for(int i = 0; i < 10; i++) {
+			acc.compose(0.1, 0.0, 1.0); // will throw before the 10th
+			double currentEpsilonSpent = acc.totalEpsilonSpent();
+			double currentRemainingBudget = acc.remainingBudget();
+			int releaseCount = acc.releaseCount();
+			Assert.assertTrue("Epsilon spent must increase with each release", currentEpsilonSpent > prevEpsilonSpent);
+			Assert.assertTrue("Remaining budget must decrease", currentRemainingBudget < prevRemainingBudget);
+			Assert.assertEquals("Release count must match", i+1, releaseCount);
+			prevEpsilonSpent = currentEpsilonSpent;
+			prevRemainingBudget = currentRemainingBudget;
+		}
+	}
+
+	@Test
+	public void testGaussianTighterThanLaplaceForSameEpsilon() {
+		// For the same nominal (epsilon, delta), Gaussian uses RDP composition which
+		// is tighter than Laplace with basic composition. After 5 releases:
+		// Laplace (basic, worst-case): 5epsilon
+		// Gaussian (RDP) : something < 5epsilon
+		double eps = 0.5;
+		double delta = 1e-5;
+
+		DPBudgetAccountant gaussian = new DPBudgetAccountant(100.0, delta);
+		DPBudgetAccountant laplace = new DPBudgetAccountant(100.0, delta);
+
+		for(int i = 0; i < 5; i++) {
+			gaussian.compose(eps, delta, 1.0);
+			laplace.compose(eps, 0.0, 1.0);
+		}
+
+		// After 5 releases, Gaussian RDP bound should be tighter.
+		// (Both may be < 5*eps; the point is Gaussian <= Laplace.)
+		Assert.assertTrue("Gaussian RDP bound should be <= Laplace bound after 5 releases",
+			gaussian.totalEpsilonSpent() <= laplace.totalEpsilonSpent() + 1e-6);
+	}
+
+	@Test
+	public void testHigherEpsilonCostMoreForLaplace() {
+		// For Laplace, the accountant uses basic (pure epsilon-DP) composition: cost = epsilon.
+		// Sensitivity determines noise scale but NOT the budget consumed - that is set
+		// entirely by the caller's epsilon parameter.
+		// A release at epsilon=1.0 costs more budget than one at epsilon=0.5.
+		DPBudgetAccountant acc1 = new DPBudgetAccountant(100.0, 1e-5);
+		DPBudgetAccountant acc2 = new DPBudgetAccountant(100.0, 1e-5);
+		acc1.compose(0.5, 0.0, 1.0); // epsilon=0.5, Laplace
+		acc2.compose(1.0, 0.0, 1.0); // epsilon=1.0, same sensitivity
+
+		Assert.assertTrue("Higher epsilon costs more budget (Laplace basic composition)",
+			acc1.totalEpsilonSpent() < acc2.totalEpsilonSpent());
+	}
+
+	// --- Constructor error paths ------------------------------------
+
+	@Test
+	public void testConstructorRejectsInvalidBudgetParameters() {
+		assertConstructorRejects(0.0, 1e-5); // zero epsilon budget
+		assertConstructorRejects(-0.5, 1e-5); // negative epsilon budget
+		assertConstructorRejects(1.0, 0.0); // delta = 0
+		assertConstructorRejects(1.0, 1.0); // delta = 1
+	}
+
+	private static void assertConstructorRejects(double epsilonBudget, double delta) {
+		try {
+			new DPBudgetAccountant(epsilonBudget, delta);
+			Assert.fail("Expected DMLRuntimeException for epsilonBudget=" + epsilonBudget + ", delta=" + delta);
+		}
+		catch(DMLRuntimeException e) {
+			// expected
+		}
+	}
+
+	@Test(expected = DMLRuntimeException.class)
+	public void testGaussianBudgetExhaustionThrows() {
+		// Budget = 0.5. Each Gaussian release is composed at epsilon=0.3, delta=1e-5,
+		// so the RDP-converted cost accumulates well past the budget within 20 releases.
+		DPBudgetAccountant acc = new DPBudgetAccountant(0.5, 1e-5);
+		for(int i = 0; i < 20; i++) {
+			acc.compose(0.3, 1e-5, 1.0);
+		}
+	}
+
+	@Test
+	public void testMixedCompositionExceedsEitherAlone() {
+		// Compose one Laplace and one Gaussian release. The total cost must
+		// exceed what either mechanism contributes alone, exercising the
+		// _pureEpsilonSum + gaussianEps addition path in totalEpsilonSpent().
+		DPBudgetAccountant mixed = new DPBudgetAccountant(100.0, 1e-5);
+		DPBudgetAccountant lapOnly = new DPBudgetAccountant(100.0, 1e-5);
+		DPBudgetAccountant gauOnly = new DPBudgetAccountant(100.0, 1e-5);
+
+		mixed.compose(0.5, 0.0, 1.0); // Laplace
+		mixed.compose(0.5, 1e-5, 1.0); // Gaussian
+
+		lapOnly.compose(0.5, 0.0, 1.0);
+		gauOnly.compose(0.5, 1e-5, 1.0);
+
+		Assert.assertTrue("Mixed cost must exceed Laplace-only cost",
+			mixed.totalEpsilonSpent() > lapOnly.totalEpsilonSpent());
+		Assert.assertTrue("Mixed cost must exceed Gaussian-only cost",
+			mixed.totalEpsilonSpent() > gauOnly.totalEpsilonSpent());
+	}
+
+	@Test
+	public void testGaussianSensitivityCancelsInRDP() {
+		// For the Gaussian mechanism: sigma = delta_f*sqrt(2*ln(1.25/delta))/epsilon, so
+		// D_alpha = alpha*delta_f^2/(2*sigma^2) = alpha*epsilon^2/(4*ln(1.25/delta)).
+		// Sensitivity cancels. Two accountants with the same (epsilon,delta) but
+		// different sensitivity must report identical totalEpsilonSpent().
+		DPBudgetAccountant acc1 = new DPBudgetAccountant(100.0, 1e-5);
+		DPBudgetAccountant acc2 = new DPBudgetAccountant(100.0, 1e-5);
+		acc1.compose(0.5, 1e-5, 1.0); // sensitivity = 1
+		acc2.compose(0.5, 1e-5, 100.0); // sensitivity = 100, same (epsilon,delta)
+		Assert.assertEquals("Gaussian RDP cost must be independent of sensitivity when (epsilon,delta) are fixed",
+			acc1.totalEpsilonSpent(), acc2.totalEpsilonSpent(), EPS);
+	}
+
+	@Test
+	public void testGaussianLargerEpsilonCostsMoreBudget() {
+		// D_alpha is proportional to epsilon^2, so a release declared at a higher epsilon (less noise, more
+		// privacy loss) must cost more budget than one at a lower epsilon.
+		DPBudgetAccountant lowEps = new DPBudgetAccountant(100.0, 1e-5);
+		DPBudgetAccountant highEps = new DPBudgetAccountant(100.0, 1e-5);
+		lowEps.compose(0.1, 1e-5, 1.0);
+		highEps.compose(0.5, 1e-5, 1.0);
+		Assert.assertTrue("Larger epsilon per Gaussian release must cost more budget",
+			highEps.totalEpsilonSpent() > lowEps.totalEpsilonSpent());
+	}
+
+	// =======================================================================
+	// 2. Noise distribution tests (statistical sanity checks)
+	// =======================================================================
+	// These tests generate many samples and verify that the empirical mean
+	// is near zero and the empirical variance matches the theoretical value
+	// within a reasonable tolerance.
+	//
+
+	@Test
+	public void testLaplaceNoiseDistribution() throws ReflectiveOperationException {
+		// For 10000 samples the empirical mean should be within 5*sigma/sqrt(n) of 0,
+		// and Var[Laplace(0, b)] = 2b^2 should match within 10% relative error.
+		int n = 10_000;
+		double scale = 1.5;
+		double[] samples = sampleLaplace(n, scale);
+
+		double mean = mean(samples);
+		double theoreticalStdErr = scale * Math.sqrt(2.0) / Math.sqrt(n);
+		Assert.assertTrue("Laplace mean should be near 0", Math.abs(mean) < 5 * theoreticalStdErr);
+
+		double variance = variance(samples);
+		double expected = 2.0 * scale * scale;
+		Assert.assertEquals("Laplace variance", expected, variance, 0.1 * expected);
+	}
+
+	@Test
+	public void testGaussianNoiseDistribution() throws ReflectiveOperationException {
+		// Same idea as testLaplaceNoiseDistribution: check mean and variance from one sample set.
+		int n = 10_000;
+		double sigma = 2.0;
+		double[] samples = sampleGaussian(n, sigma);
+
+		double mean = mean(samples);
+		double theoreticalStdErr = sigma / Math.sqrt(n);
+		Assert.assertTrue("Gaussian mean should be near 0", Math.abs(mean) < 5 * theoreticalStdErr);
+
+		double variance = variance(samples);
+		double expected = sigma * sigma;
+		Assert.assertEquals("Gaussian variance", expected, variance, 0.1 * expected);
+	}
+
+	// =======================================================================
+	// 3. DPBuiltinCPInstruction structural tests
+	// =======================================================================
+
+	@Test
+	public void testParseInstructionValidLaplace() {
+		DPBuiltinCPInstruction inst = DPBuiltinCPInstruction
+			.parseInstruction("CP°dp_laplace°target=mVar1°query=colMeans°sensitivity=1.0°epsilon=0.5°_mVar2·MATRIX·FP64");
+		Assert.assertEquals(DPBuiltinCPInstruction.OPCODE_LAPLACE, inst.getOpcode());
+		Assert.assertEquals("mVar1", inst.input1.getName());
+		Assert.assertEquals("_mVar2", inst.getOutput().getName());
+	}
+
+	@Test
+	public void testParseInstructionValidGaussian() {
+		DPBuiltinCPInstruction inst = DPBuiltinCPInstruction.parseInstruction(
+			"CP°dp_gaussian°target=mVar1°query=colMeans°sensitivity=1.0°epsilon=0.5°delta=1e-5°_mVar2·MATRIX·FP64");
+		Assert.assertEquals(DPBuiltinCPInstruction.OPCODE_GAUSSIAN, inst.getOpcode());
+		Assert.assertEquals("mVar1", inst.input1.getName());
+		Assert.assertEquals("_mVar2", inst.getOutput().getName());
+	}
+
+	@Test
+	public void testParseInstructionMissingRequiredKeysThrows() {
+		// Each string below keeps the field COUNT that checkNumFields expects (5 for laplace,
+		// 6 for gaussian), but renames one required key so parseInstruction's own
+		// containsKey(...) checks - not checkNumFields - are what reject it.
+		assertParseInstructionRejects(
+			"CP°dp_laplace°target=mVar1°qry=colMeans°sensitivity=1.0°epsilon=0.5°_mVar2·MATRIX·FP64"); // missing query
+		assertParseInstructionRejects(
+			"CP°dp_laplace°target=mVar1°query=colMeans°sens=1.0°epsilon=0.5°_mVar2·MATRIX·FP64"); // missing sensitivity
+		assertParseInstructionRejects(
+			"CP°dp_laplace°target=mVar1°query=colMeans°sensitivity=1.0°eps=0.5°_mVar2·MATRIX·FP64"); // missing epsilon
+		assertParseInstructionRejects("CP°dp_gaussian°target=mVar1°query=colMeans°sensitivity=1.0°epsilon=0.5°"
+			+ "notdelta=1e-5°_mVar2·MATRIX·FP64"); // missing delta (gaussian only)
+	}
+
+	private static void assertParseInstructionRejects(String instStr) {
+		try {
+			DPBuiltinCPInstruction.parseInstruction(instStr);
+			Assert.fail("Expected DMLRuntimeException for instruction: " + instStr);
+		}
+		catch(DMLRuntimeException e) {
+			// expected
+		}
+	}
+
+	// -----------------------------------------------------------------------
+	// Helpers for noise distribution tests
+	// -----------------------------------------------------------------------
+
+	/** Sample n Laplace(0, scale) values via the fillLaplaceNoise method. */
+	private static double[] sampleLaplace(int n, double scale) throws ReflectiveOperationException {
+		Method m = DPBuiltinCPInstruction.class.getDeclaredMethod("fillLaplaceNoise", int.class, int.class, double.class);
+		m.setAccessible(true);
+		MatrixBlock block = (MatrixBlock) m.invoke(null, n, 1, scale);
+
+		double[] out = new double[n];
+		for(int i = 0; i < n; i++)
+			out[i] = block.get(i, 0);
+		return out;
+	}
+
+	/** Sample n N(0, sigma^2) values via the production fillGaussianNoise method. */
+	private static double[] sampleGaussian(int n, double sigma) throws ReflectiveOperationException {
+		Method m = DPBuiltinCPInstruction.class.getDeclaredMethod("fillGaussianNoise", int.class, int.class,
+			double.class);
+		m.setAccessible(true);
+		MatrixBlock block = (MatrixBlock) m.invoke(null, n, 1, sigma);
+
+		double[] out = new double[n];
+		for(int i = 0; i < n; i++)
+			out[i] = block.get(i, 0);
+		return out;
+	}
+
+	private static double mean(double[] xs) {
+		double s = 0;
+		for(double x : xs)
+			s += x;
+		return s / xs.length;
+	}
+
+	private static double variance(double[] xs) {
+		double m = mean(xs);
+		double s = 0;
+		for(double x : xs)
+			s += (x - m) * (x - m);
+		return s / (xs.length - 1);
+	}
+}
diff --git a/src/test/java/org/apache/sysds/test/functions/privacy/dp/DPBuiltinDMLTest.java b/src/test/java/org/apache/sysds/test/functions/privacy/dp/DPBuiltinDMLTest.java
new file mode 100644
index 00000000000..c6750ceebd9
--- /dev/null
+++ b/src/test/java/org/apache/sysds/test/functions/privacy/dp/DPBuiltinDMLTest.java
@@ -0,0 +1,263 @@
+/*
+ * Licensed to the Apache Software Foundation (ASF) under one
+ * or more contributor license agreements.  See the NOTICE file
+ * distributed with this work for additional information
+ * regarding copyright ownership.  The ASF licenses this file
+ * to you under the Apache License, Version 2.0 (the
+ * "License"); you may not use this file except in compliance
+ * with the License.  You may obtain a copy of the License at
+ *
+ *   http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing,
+ * software distributed under the License is distributed on an
+ * "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
+ * KIND, either express or implied.  See the License for the
+ * specific language governing permissions and limitations
+ * under the License.
+ */
+
+package org.apache.sysds.test.functions.privacy.dp;
+
+import static org.junit.Assert.assertEquals;
+import static org.junit.Assert.assertTrue;
+
+import java.io.File;
+import java.io.IOException;
+import java.nio.file.Files;
+import java.util.HashMap;
+
+import org.apache.sysds.parser.LanguageException;
+import org.apache.sysds.parser.ParseException;
+import org.apache.sysds.runtime.DMLRuntimeException;
+import org.apache.sysds.runtime.matrix.data.MatrixValue.CellIndex;
+import org.apache.sysds.test.AutomatedTestBase;
+import org.apache.sysds.test.TestConfiguration;
+import org.junit.Test;
+
+/*
+ * ==========================================================================
+ * DML integration test
+ * ==========================================================================
+ *
+ * Full integration tests extend AutomatedTestBase and drive the DML runner.
+ * Each test:
+ *   (a) Writes a DML script to a temp file.
+ *   (b) Provides input matrices via TestUtils.
+ *   (c) Calls runTest() and reads the output MatrixBlock.
+ *   (d) Verifies that the noisy result differs from the clean result by a
+ *       statistically plausible amount (not zero, not astronomically large).
+ */
+public class DPBuiltinDMLTest extends AutomatedTestBase {
+
+	private static final String TEST_DIR = "functions/privacy/dp/";
+	private static final String TEST_CLASS = TEST_DIR + DPBuiltinDMLTest.class.getSimpleName() + "/";
+	private static final int ROWS = 100;
+	private static final int COLS = 10;
+
+	private double[][] data;
+
+	private static final String DML_LAPLACE_TEMPLATE = "X = read($1);\n"
+		+ "result = dp_laplace(X, query=\"%s\", sensitivity=1.0, epsilon=$2);\n"
+		+ "write(result, $3, format=\"text\");\n";
+	private static final String DML_GAUSSIAN_TEMPLATE = "X = read($1);\n"
+		+ "result = dp_gaussian(X, query=\"%s\", sensitivity=1.0, epsilon=$2, delta=1e-5);\n"
+		+ "write(result, $3, format=\"text\");\n";
+
+	private static final String DML_LAPLACE = String.format(DML_LAPLACE_TEMPLATE, "colMeans");
+	private static final String DML_GAUSSIAN = String.format(DML_GAUSSIAN_TEMPLATE, "colMeans");
+
+	// dp_set_budget(epsilon, delta) is resolved entirely at compile time (its arguments
+	// must be literals), called via a dummy assignment, then a single dp_laplace release
+	// at $2 records a cost of exactly $2 (Laplace basic composition).
+	private static final String DML_SET_BUDGET_TEMPLATE = "eps = dp_set_budget(%s, 1e-6);\n" + "X = read($1);\n"
+		+ "result = dp_laplace(X, query=\"colMeans\", sensitivity=1.0, epsilon=$2);\n"
+		+ "write(result, $3, format=\"text\");\n";
+
+	// Two dp_set_budget calls in the same script; DMLTranslator must reject this at
+	// compile time (DMLProgram.hasDPBudget()).
+	private static final String DML_SET_BUDGET_TWICE = "eps = dp_set_budget(3.0, 1e-6);\n"
+		+ "eps2 = dp_set_budget(5.0, 1e-6);\n" + "X = read($1);\n"
+		+ "result = dp_laplace(X, query=\"colMeans\", sensitivity=1.0, epsilon=$2);\n"
+		+ "write(result, $3, format=\"text\");\n";
+
+	// A budget argument computed at runtime (not a literal); must be rejected at
+	// compile time by BuiltinFunctionExpression's isConstant() check.
+	private static final String DML_SET_BUDGET_NON_LITERAL = "X = read($1);\n" + "computed = sum(X) / nrow(X);\n"
+		+ "eps = dp_set_budget(computed, 1e-6);\n"
+		+ "result = dp_laplace(X, query=\"colMeans\", sensitivity=1.0, epsilon=$2);\n"
+		+ "write(result, $3, format=\"text\");\n";
+
+	@Override
+	public void setUp() {
+		addTestConfiguration("DPLaplace", new TestConfiguration(TEST_CLASS, "DPLaplace"));
+		addTestConfiguration("DPGaussian", new TestConfiguration(TEST_CLASS, "DPGaussian"));
+		addTestConfiguration("DPSetBudget", new TestConfiguration(TEST_CLASS, "DPSetBudget"));
+		data = this.getRandomMatrix(ROWS, COLS, 0, 1, 1.0, 42);
+	}
+
+	@Test
+	public void testLaplaceOutputDiffersFromCleanMean() {
+		runColMeansDPTest("DPLaplace", DML_LAPLACE, "0.5");
+	}
+
+	@Test
+	public void testGaussianOutputDiffersFromCleanMean() {
+		runColMeansDPTest("DPGaussian", DML_GAUSSIAN, "0.5");
+	}
+
+	@Test
+	public void testLaplaceColSums() {
+		// query="colSums": T is 1 x n filled with 1.0, output is the noisy column-sum row vector.
+		HashMap<CellIndex, Double> result = runAndGetResult("DPLaplace", String.format(DML_LAPLACE_TEMPLATE, "colSums"),
+			"0.5", data);
+		assertShape(result, 1, COLS);
+		double maxDiff = maxAbsDiffFromClean(data, result, DPBuiltinDMLTest::colSum);
+		assertTrue("Result should differ from the clean column sums", maxDiff > 0);
+	}
+
+	@Test
+	public void testGaussianIdentity() {
+		// query="identity": T is the n x n identity, output is a noisy release of X itself.
+		HashMap<CellIndex, Double> result = runAndGetResult("DPGaussian",
+			String.format(DML_GAUSSIAN_TEMPLATE, "identity"), "0.5", data);
+		assertShape(result, ROWS, COLS);
+		// identity releases X row-by-row, so compare cell-by-cell rather than via a per-column reduction.
+		double maxCellDiff = 0;
+		for(int r = 0; r < ROWS; r++) {
+			for(int c = 0; c < COLS; c++) {
+				double noisy = result.get(new CellIndex(r + 1, c + 1));
+				maxCellDiff = Math.max(maxCellDiff, Math.abs(noisy - data[r][c]));
+			}
+		}
+		assertTrue("Result should differ from the clean matrix", maxCellDiff > 0);
+	}
+
+	@Test
+	public void testHighEpsilonIsCloserToTruth() {
+		// Higher epsilon => less noise => result closer to the true mean.
+		// NOTE: the DPBudgetAccountant caps total spend at the default budget
+		// (epsilon = 1.0) regardless of the per-release epsilon requested, so epsilon values
+		// here must stay well under that cap or the release is rejected.
+		double noisyLow = runAndGetMaxAbsColMeansDiffFromClean(data, "DPGaussian", DML_GAUSSIAN, "0.1");
+		double noisyHigh = runAndGetMaxAbsColMeansDiffFromClean(data, "DPGaussian", DML_GAUSSIAN, "0.5");
+		assertTrue("epsilon=0.5 should give less noise than epsilon=0.1", noisyHigh < noisyLow);
+	}
+
+	@Test
+	public void testSetBudgetLiteralAllowsExceedingDefaultBudget() {
+		// Default budget is epsilon=1.0; a single release at epsilon=1.5 would be
+		// rejected unless dp_set_budget(3.0, ...) widens it first.
+		HashMap<CellIndex, Double> result = runAndGetResult("DPSetBudget",
+			String.format(DML_SET_BUDGET_TEMPLATE, "3.0"), "1.5", data);
+		assertShape(result, 1, COLS);
+	}
+
+	@Test
+	public void testSetBudgetNarrowBudgetStillEnforced() {
+		// An explicit narrow budget must still be enforced: epsilon=0.8 exceeds
+		// the explicit budget of 0.5.
+		runExpectingException("DPSetBudget", String.format(DML_SET_BUDGET_TEMPLATE, "0.5"), "0.8", data,
+			DMLRuntimeException.class);
+	}
+
+	@Test
+	public void testSetBudgetCalledTwiceFailsAtCompileTime() {
+		// Thrown from DMLTranslator.processBuiltinFunctionExpression (HOP construction),
+		// which wraps all case-block exceptions in ParseException (see processExpression's
+		// catch-all) - unlike the non-literal check below, which runs during validation
+		// and so surfaces as a bare LanguageException.
+		runExpectingException("DPSetBudget", DML_SET_BUDGET_TWICE, "0.5", data, ParseException.class);
+	}
+
+	@Test
+	public void testSetBudgetRejectsNonLiteralArgs() {
+		runExpectingException("DPSetBudget", DML_SET_BUDGET_NON_LITERAL, "0.5", data, LanguageException.class);
+	}
+
+	private void runColMeansDPTest(String testName, String dml, String epsilonStr) {
+		HashMap<CellIndex, Double> result = runAndGetResult(testName, dml, epsilonStr, data);
+		assertShape(result, 1, COLS);
+		// Must differ from the exact (clean) mean by a non-trivial amount.
+		// (A single-seed exact-equality check is fragile; use range check.)
+		double maxDiff = maxAbsDiffFromClean(data, result, DPBuiltinDMLTest::colMean);
+		assertTrue("Result should differ from the clean mean", maxDiff > 0);
+	}
+
+	private double runAndGetMaxAbsColMeansDiffFromClean(double[][] data, String testName, String dml,
+		String epsilonStr) {
+		HashMap<CellIndex, Double> result = runAndGetResult(testName, dml, epsilonStr, data);
+		return maxAbsDiffFromClean(data, result, DPBuiltinDMLTest::colMean);
+	}
+
+	private static void assertShape(HashMap<CellIndex, Double> result, int expectedRows, int expectedCols) {
+		int maxRow = 0, maxCol = 0;
+		for(CellIndex ci : result.keySet()) {
+			maxRow = Math.max(maxRow, ci.row);
+			maxCol = Math.max(maxCol, ci.column);
+		}
+		assertEquals("Result should have " + expectedRows + " row(s)", expectedRows, maxRow);
+		assertEquals("Result should have " + expectedCols + " column(s)", expectedCols, maxCol);
+	}
+
+	@FunctionalInterface
+	private interface CleanColumnFn {
+		double apply(double[][] data, int col);
+	}
+
+	/** Computes max|noisy(1,c) - clean(data,c)| across the (1 x COLS) row-vector releases. */
+	private static double maxAbsDiffFromClean(double[][] data, HashMap<CellIndex, Double> result,
+		CleanColumnFn cleanFn) {
+		double maxDiff = 0;
+		for(int c = 0; c < COLS; c++) {
+			double clean = cleanFn.apply(data, c);
+			double noisy = result.get(new CellIndex(1, c + 1));
+			maxDiff = Math.max(maxDiff, Math.abs(noisy - clean));
+		}
+		return maxDiff;
+	}
+
+	private static double colMean(double[][] data, int c) {
+		double sum = 0;
+		for(int r = 0; r < ROWS; r++)
+			sum += data[r][c];
+		return sum / ROWS;
+	}
+
+	private static double colSum(double[][] data, int c) {
+		double sum = 0;
+		for(int r = 0; r < ROWS; r++)
+			sum += data[r][c];
+		return sum;
+	}
+
+	private HashMap<CellIndex, Double> runAndGetResult(String testName, String dml, String epsilonStr,
+		double[][] data) {
+		prepareScript(testName, dml, epsilonStr, data);
+		runTest(true, false, null, -1);
+		return readDMLMatrixFromOutputDir("result");
+	}
+
+	private void runExpectingException(String testName, String dml, String epsilonStr, double[][] data,
+		Class<?> expectedException) {
+		prepareScript(testName, dml, epsilonStr, data);
+		runTest(true, true, expectedException, -1);
+	}
+
+	private void prepareScript(String testName, String dml, String epsilonStr, double[][] data) {
+		getAndLoadTestConfiguration(testName);
+		writeInputMatrixWithMTD("X", data, false);
+
+		fullDMLScriptName = getScript();
+		try {
+			File scriptFile = new File(fullDMLScriptName);
+			scriptFile.getParentFile().mkdirs();
+			Files.write(scriptFile.toPath(), dml.getBytes());
+		}
+		catch(IOException e) {
+			throw new RuntimeException(e);
+		}
+
+		programArgs = new String[] {"-args", input("X"), epsilonStr, output("result")};
+	}
+}