diff --git a/CHANGELOG.md b/CHANGELOG.md
index b40089fb..94f6b6f0 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -9,6 +9,7 @@ To install or update LODA, please follow the [installation instructions](https:/
 
 ### Enhancements
 
+* Support `prg` in `unfold` command
 * Internal `replace` command
 * Remove deprecated dot export
 
diff --git a/src/cmd/test.cpp b/src/cmd/test.cpp
index 14d70ece..fe63df3a 100644
--- a/src/cmd/test.cpp
+++ b/src/cmd/test.cpp
@@ -59,6 +59,7 @@ void Test::fast() {
   config();
   steps();
   blocks();
+  fold();
   unfold();
   incEval();
   linearMatcher();
@@ -701,10 +702,57 @@ void Test::knownPrograms() {
                     60480, 60480, 60480, 120960, 120960, 241920, 1209600}));
 }
 
+void Test::fold() {
+  auto tests =
+      loadInOutTests(std::string("tests") + FILE_SEP + "fold" + FILE_SEP + "F");
+  size_t i = 1;
+  Parser parser;
+  Evaluator evaluator(settings, false);
+  std::map<int64_t, int64_t> cellMap;
+  for (const auto& t : tests) {
+    Log::get().info("Testing fold " + std::to_string(i));
+    Operation seqOrPrg;
+    for (const auto& op : t.second.ops) {
+      if (op.type == Operation::Type::SEQ || op.type == Operation::Type::PRG) {
+        seqOrPrg = op;
+        break;
+      }
+    }
+    if (seqOrPrg.type == Operation::Type::NOP) {
+      Log::get().error("No seq or prg in output found", true);
+    }
+    int64_t subId = seqOrPrg.source.value.asInt();
+    std::string path;
+    if (seqOrPrg.type == Operation::Type::SEQ) {
+      path = ProgramUtil::getProgramPath(subId);
+    } else {
+      path = ProgramUtil::getProgramPath(subId, "prg", "P");
+    }
+    auto sub = parser.parse(path);
+    auto p = t.first;
+    cellMap.clear();
+    if (!Subprogram::fold(p, sub, subId, cellMap)) {
+      Log::get().error("Folding not supported", true);
+    }
+    if (p != t.second) {
+      ProgramUtil::print(p, std::cout);
+      Log::get().error("Unexpected program", true);
+    }
+    Sequence expected, got;
+    evaluator.eval(t.first, expected, 20);
+    evaluator.eval(t.second, got, 20);
+    if (expected != got) {
+      Log::get().error("Unexpected sequence", true);
+    }
+    i++;
+  }
+}
+
 void Test::unfold() {
   auto tests = loadInOutTests(std::string("tests") + FILE_SEP + "unfold" +
                               FILE_SEP + "U");
   size_t i = 1;
+  Evaluator evaluator(settings, false);
   for (const auto& t : tests) {
     Log::get().info("Testing unfold " + std::to_string(i));
     auto p = t.first;
@@ -713,7 +761,13 @@ void Test::unfold() {
     }
     if (p != t.second) {
       ProgramUtil::print(p, std::cout);
-      Log::get().error("Unexpected result", true);
+      Log::get().error("Unexpected program", true);
+    }
+    Sequence expected, got;
+    evaluator.eval(t.first, expected, 20);
+    evaluator.eval(t.second, got, 20);
+    if (expected != got) {
+      Log::get().error("Unexpected sequence", true);
     }
     i++;
   }
diff --git a/src/cmd/test.hpp b/src/cmd/test.hpp
index 4e9bfd38..47830b2b 100644
--- a/src/cmd/test.hpp
+++ b/src/cmd/test.hpp
@@ -34,6 +34,8 @@ class Test {
 
   void knownPrograms();
 
+  void fold();
+
   void unfold();
 
   void incEval();
diff --git a/src/lang/program_util.cpp b/src/lang/program_util.cpp
index e6998c14..d017f20f 100644
--- a/src/lang/program_util.cpp
+++ b/src/lang/program_util.cpp
@@ -162,6 +162,12 @@ bool ProgramUtil::isWritingRegion(Operation::Type t) {
           t == Operation::Type::SRT);
 }
 
+bool ProgramUtil::hasRegionOperation(const Program &p) {
+  return std::any_of(p.ops.begin(), p.ops.end(), [](const Operation &op) {
+    return isWritingRegion(op.type);
+  });
+}
+
 bool ProgramUtil::hasIndirectOperand(const Operation &op) {
   const auto num_ops = Operation::Metadata::get(op.type).num_operands;
   return (num_ops > 0 && op.target.type == Operand::Type::INDIRECT) ||
diff --git a/src/lang/program_util.hpp b/src/lang/program_util.hpp
index cffc573d..ac919df2 100644
--- a/src/lang/program_util.hpp
+++ b/src/lang/program_util.hpp
@@ -39,6 +39,8 @@ class ProgramUtil {
 
   static bool isWritingRegion(Operation::Type t);
 
+  static bool hasRegionOperation(const Program &p);
+
   static bool hasIndirectOperand(const Operation &op);
 
   static bool hasIndirectOperand(const Program &p);
diff --git a/src/lang/subprogram.cpp b/src/lang/subprogram.cpp
index f5e47489..61c8a191 100644
--- a/src/lang/subprogram.cpp
+++ b/src/lang/subprogram.cpp
@@ -92,26 +92,35 @@ int64_t Subprogram::search(const Program &main, const Program &sub,
   return main_pos - sub_pos;
 }
 
-void updateOperand(Operand &op, int64_t target, int64_t largest_used) {
+void updateOperand(Operand &op, int64_t start, int64_t shared_region_length,
+                   int64_t largest_used) {
   if (op.type != Operand::Type::DIRECT) {
     return;
   }
-  if (op.value.asInt() == Program::INPUT_CELL) {
-    op.value = target;
+  auto v = op.value.asInt();
+  if (v < shared_region_length) {
+    op.value += start;
   } else {
     op.value += largest_used;
   }
 }
 
-bool prepareEmbedding(int64_t id, Program &sub) {
+bool prepareEmbedding(int64_t id, Program &sub, Operation::Type embeddingType) {
   // load and check program to be embedded
+  std::string path;
+  if (embeddingType == Operation::Type::SEQ) {
+    path = ProgramUtil::getProgramPath(id);
+  } else if (embeddingType == Operation::Type::PRG) {
+    path = ProgramUtil::getProgramPath(id, "prg", "P");
+  } else {
+    throw std::runtime_error("Unsupported embedding type");
+  }
   Parser parser;
-  sub = parser.parse(ProgramUtil::getProgramPath(id));
+  sub = parser.parse(path);
   if (ProgramUtil::hasIndirectOperand(sub)) {
     return false;
   }
-  // prepare program for embedding
-  // remove nops and comments
+  // prepare program for embedding: remove nops and comments
   ProgramUtil::removeOps(sub, Operation::Type::NOP);
   for (auto &op : sub.ops) {
     if (op.type != Operation::Type::SEQ) {
@@ -120,7 +129,16 @@ bool prepareEmbedding(int64_t id, Program &sub) {
   }
   // find cells that are read and uninitialized
   std::set<int64_t> initialized, uninitialized;
-  initialized.insert(Program::INPUT_CELL);
+  if (embeddingType == Operation::Type::SEQ) {
+    initialized.insert(Program::INPUT_CELL);
+  } else if (embeddingType == Operation::Type::PRG) {
+    auto numInputs = sub.getDirective("inputs");
+    for (int64_t i = 0; i < numInputs; i++) {
+      initialized.insert(i);
+    }
+  } else {
+    throw std::runtime_error("Unsupported embedding type");
+  }
   ProgramUtil::getUsedUninitializedCells(sub, initialized, uninitialized);
   // initialize cells that are read and were uninitialized
   for (auto cell : uninitialized) {
@@ -136,34 +154,40 @@ bool Subprogram::unfold(Program &main) {
   if (ProgramUtil::hasIndirectOperand(main)) {
     return false;
   }
-  // TODO unfold prg instructions, too
-  // find first seq operation
-  int64_t seq_index = -1;
+  // find first seq or prg operation
+  int64_t pos = -1;
   for (size_t i = 0; i < main.ops.size(); i++) {
-    if (main.ops[i].type == Operation::Type::SEQ) {
-      seq_index = i;
+    if (main.ops[i].type == Operation::Type::SEQ ||
+        main.ops[i].type == Operation::Type::PRG) {
+      pos = i;
       break;
     }
   }
-  if (seq_index < 0) {
+  if (pos < 0) {
     return false;
   }
-  auto sub_id = main.ops[seq_index].source.value.asInt();
+  const auto &emb_op = main.ops[pos];
+  auto sub_id = emb_op.source.value.asInt();
   Program sub;
-  if (!prepareEmbedding(sub_id, sub)) {
+  if (!prepareEmbedding(sub_id, sub, emb_op.type)) {
     return false;
   }
   // shift used operands
-  auto target = main.ops[seq_index].target.value.asInt();
-  auto largest_used = ProgramUtil::getLargestDirectMemoryCell(main);
+  int64_t start = main.ops[pos].target.value.asInt();
+  int64_t shared_region_length = 1;
+  if (emb_op.type == Operation::Type::PRG) {
+    shared_region_length = std::max<int64_t>(sub.getDirective("inputs"),
+                                             sub.getDirective("outputs"));
+  }
+  int64_t largest_used = ProgramUtil::getLargestDirectMemoryCell(main);
   for (auto &op : sub.ops) {
-    updateOperand(op.target, target, largest_used);
-    updateOperand(op.source, target, largest_used);
+    updateOperand(op.target, start, shared_region_length, largest_used);
+    updateOperand(op.source, start, shared_region_length, largest_used);
   }
   // delete seq operation
-  main.ops.erase(main.ops.begin() + seq_index);
+  main.ops.erase(main.ops.begin() + pos);
   // embed program
-  main.ops.insert(main.ops.begin() + seq_index, sub.ops.begin(), sub.ops.end());
+  main.ops.insert(main.ops.begin() + pos, sub.ops.begin(), sub.ops.end());
   return true;
 }
 
@@ -189,7 +213,7 @@ bool Subprogram::autoUnfold(Program &main) {
 bool Subprogram::shouldFold(const Program &main) {
   int64_t level = 0;
   int64_t numLoops = 0;
-  bool hasRootSeq = false;
+  bool hasRootRef = false;
   for (const auto &op : main.ops) {
     switch (op.type) {
       case Operation::Type::LPB:
@@ -199,51 +223,54 @@ bool Subprogram::shouldFold(const Program &main) {
       case Operation::Type::LPE:
         level--;
         break;
+      case Operation::Type::PRG:
       case Operation::Type::SEQ:
         if (level == 0) {
-          hasRootSeq = true;
+          hasRootRef = true;
         }
         break;
       default:
         break;
     }
   }
-  return (numLoops > 1) || (numLoops > 0 && hasRootSeq);
+  return (numLoops > 1) || (numLoops > 0 && hasRootRef);
 }
 
 bool Subprogram::fold(Program &main, Program sub, size_t subId,
-                      std::map<int64_t, int64_t> &cell_map) {
+                      std::map<int64_t, int64_t> &cellMap) {
   // prepare and check subprogram
   ProgramUtil::removeOps(sub, Operation::Type::NOP);
   if (sub.ops.empty()) {
     return false;
   }
-  int64_t output_cell = Program::OUTPUT_CELL;
+  int64_t outputCell = Program::OUTPUT_CELL;
   auto last = sub.ops.back();
   if (sub.ops.size() > 1 && last.type == Operation::Type::MOV &&
       last.target == Operand(Operand::Type::DIRECT, Program::OUTPUT_CELL) &&
       last.source.type == Operand::Type::DIRECT) {
-    output_cell = last.source.value.asInt();
+    outputCell = last.source.value.asInt();
     sub.ops.pop_back();
   }
-  auto main_pos = Subprogram::search(main, sub, cell_map);
-  // no indirect ops allowed
+  auto mainPos = Subprogram::search(main, sub, cellMap);
+  // no indirect operands  allowed
   if (ProgramUtil::hasIndirectOperand(main) ||
       ProgramUtil::hasIndirectOperand(sub)) {
     return false;
   }
   // perform folding on main program
-  const Number mapped_input(cell_map.at(Program::INPUT_CELL));
-  const Number mapped_output(cell_map.at(output_cell));
-  main.ops.erase(main.ops.begin() + main_pos,
-                 main.ops.begin() + main_pos + sub.ops.size());
-  main.ops.insert(main.ops.begin() + main_pos,
+  const Number mappedInput(cellMap.at(Program::INPUT_CELL));
+  const Number mappedOutput(cellMap.at(outputCell));
+  main.ops.erase(main.ops.begin() + mainPos,
+                 main.ops.begin() + mainPos + sub.ops.size());
+  main.ops.insert(main.ops.begin() + mainPos,
                   Operation(Operation::Type::SEQ,
-                            Operand(Operand::Type::DIRECT, mapped_input),
+                            Operand(Operand::Type::DIRECT, mappedInput),
                             Operand(Operand::Type::CONSTANT, Number(subId))));
-  main.ops.insert(main.ops.begin() + main_pos + 1,
-                  Operation(Operation::Type::MOV,
-                            Operand(Operand::Type::DIRECT, mapped_output),
-                            Operand(Operand::Type::DIRECT, mapped_input)));
+  if (mappedOutput != mappedInput) {
+    main.ops.insert(main.ops.begin() + mainPos + 1,
+                    Operation(Operation::Type::MOV,
+                              Operand(Operand::Type::DIRECT, mappedOutput),
+                              Operand(Operand::Type::DIRECT, mappedInput)));
+  }
   return true;
 }
diff --git a/tests/fold/F001.asm b/tests/fold/F001.asm
new file mode 100644
index 00000000..fed4d9da
--- /dev/null
+++ b/tests/fold/F001.asm
@@ -0,0 +1,12 @@
+; in
+mul $0,5
+mov $1,$0
+lpb $1
+  div $0,10
+  sub $1,$0
+lpe
+div $0,3
+; out
+mul $0,5
+seq $0,30
+div $0,3
diff --git a/tests/fold/F002.asm b/tests/fold/F002.asm
new file mode 100644
index 00000000..708b079f
--- /dev/null
+++ b/tests/fold/F002.asm
@@ -0,0 +1,16 @@
+; in
+mul $0,5
+mov $3,$0
+mov $2,$3
+lpb $2
+  div $3,10
+  sub $2,$3
+lpe
+add $0,$3
+div $0,3
+; out
+mul $0,5
+mov $3,$0
+seq $3,30
+add $0,$3
+div $0,3
diff --git a/tests/programs/prg/000/P000002.asm b/tests/programs/prg/000/P000002.asm
new file mode 100644
index 00000000..fd189bec
--- /dev/null
+++ b/tests/programs/prg/000/P000002.asm
@@ -0,0 +1,9 @@
+; Test internal memory
+
+#inputs 1
+#outputs 1
+
+mov $1,$0
+add $1,2
+bin $1,$0
+mov $0,$1
diff --git a/tests/unfold/U003.asm b/tests/unfold/U003.asm
new file mode 100644
index 00000000..7fc7718b
--- /dev/null
+++ b/tests/unfold/U003.asm
@@ -0,0 +1,22 @@
+; in
+mov $1,$0
+add $1,2
+bin $1,$0
+mov $3,5
+prg $1,1
+add $1,$2
+add $1,$3
+mov $0,$1
+; out
+mov $1,$0
+add $1,2
+bin $1,$0
+mov $3,5
+mov $2,0
+lpb $1
+  add $2,1
+  sub $1,$2
+lpe
+add $1,$2
+add $1,$3
+mov $0,$1
diff --git a/tests/unfold/U004.asm b/tests/unfold/U004.asm
new file mode 100644
index 00000000..d535c3fc
--- /dev/null
+++ b/tests/unfold/U004.asm
@@ -0,0 +1,13 @@
+; in
+mov $1,$0
+mov $2,4
+prg $1,2
+add $1,$2
+; out
+mov $1,$0
+mov $2,4
+mov $3,$1
+add $3,2
+bin $3,$1
+mov $1,$3
+add $1,$2