added my programs for "symbolic computation", used to generate

some of the test data files


git-svn-id: svn://scm.gforge.inria.fr/svn/mpc/trunk@306 211d60ee-9f03-0410-a15a-8952a2c7a4e4

10 files changed, 695 insertions, 0 deletions

diff --git a/tools/README.testgen b/tools/README.testgen
new file mode 100644
index 0000000..281ee17
--- /dev/null
+++ b/tools/README.testgen
@@ -0,0 +1,3 @@
+The subdirectory "testgen" contains some C programs for computing with
+special values that are represented as strings "+inf", "nan" and so on.
+These have been used to obtain data files for some tests.
diff --git a/tools/testgen/makefile b/tools/testgen/makefile
new file mode 100644
index 0000000..5f94463
--- /dev/null
+++ b/tools/testgen/makefile
@@ -0,0 +1,30 @@
+OPT = -pedantic -Wall -Werror -g --std=c99
+CC = gcc
+DIR = .
+
+write_test_frsub : write_test_frsub.o
+	$(CC) -o $(DIR)/write_test_frsub write_test_frsub.o $(OPT)
+
+write_test_subfr : write_test_subfr.o
+	$(CC) -o $(DIR)/write_test_subfr write_test_subfr.o $(OPT)
+
+write_test_frdiv : write_test_frdiv.o
+	$(CC) -o $(DIR)/write_test_frdiv write_test_frdiv.o $(OPT)
+
+write_test_div : write_test_div.o
+	$(CC) -o $(DIR)/write_test_div write_test_div.o $(OPT)
+
+write_test_divfr : write_test_divfr.o
+	$(CC) -o $(DIR)/write_test_divfr write_test_divfr.o $(OPT)
+
+write_test_norm : write_test_norm.o
+	$(CC) -o $(DIR)/write_test_norm write_test_norm.o $(OPT)
+
+write_test_mulfr : write_test_mulfr.o
+	$(CC) -o $(DIR)/write_test_mulfr write_test_mulfr.o $(OPT)
+
+%.o : %.c write_test.h
+	$(CC) $< -c -o $@ $(OPT)
+
+clean :
+	rm -f *.o *~
diff --git a/tools/testgen/write_test.h b/tools/testgen/write_test.h
new file mode 100644
index 0000000..4c3f391
--- /dev/null
+++ b/tools/testgen/write_test.h
@@ -0,0 +1,231 @@
+#include <stdio.h>
+#include <string.h>
+
+char* licence =
+"# This file is part of the MPC Library.\n"
+"#\n"
+"# The MPC Library is free software; you can redistribute it and/or modify\n"
+"# it under the terms of the GNU Lesser General Public License as published by\n"
+"# the Free Software Foundation; either version 2.1 of the License, or (at your\n"
+"# option) any later version.\n"
+"#\n"
+"# The MPC Library is distributed in the hope that it will be useful, but\n"
+"# WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY\n"
+"# or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public\n"
+"# License for more details.\n"
+"#\n"
+"# You should have received a copy of the GNU Lesser General Public License\n"
+"# along with the MPC Library; see the file COPYING.LIB.  If not, write to\n"
+"# the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,\n"
+"# MA 02111-1307, USA.\n"
+"#\n"
+;
+
+char* values [] = {"-inf", "-1", "-0", "+0", "1", "+inf", "nan"};
+int valueslen = 7;
+char* infvalues [] = {"-inf", "+inf", "nan"};
+int infvalueslen = 3;
+
+int equals (char* a, char* b) {
+   return !strcmp (a, b);
+}
+
+int index (char** table, int length, char* value) {
+   // returns the index of value in the table,
+   // or -1 if it is not found
+   int i;
+   for (i = 0; i < length; i++)
+      if (equals (table [i], value))
+         return i;
+   return -1;
+}
+
+char* mul_fr_fr (char *op1, char *op2) {
+   // returns op1*op2
+   if (equals (op1, "nan") || equals (op2, "nan"))
+      return "nan";
+   // none is nan
+   else if (equals (op1, "-inf")) {
+      // res = op1 * values
+      char* res [] = {"+inf", "+inf", "nan", "nan", "-inf", "-inf", ""};
+      return res [index (values, valueslen, op2)];
+   }
+   else if (equals (op1, "-1")) {
+      char* res [] = {"+inf", "1", "+0", "-0", "-1", "-inf", ""};
+      return res [index (values, valueslen, op2)];
+   }
+   else if (equals (op1, "-0")) {
+      char* res [] = {"nan", "+0", "+0", "-0", "-0", "nan", ""};
+      return res [index (values, valueslen, op2)];
+   }
+   else if (equals (op1, "+0")) {
+      char* res [] = {"nan", "-0", "-0", "+0", "+0", "nan", ""};
+      return res [index (values, valueslen, op2)];
+   }
+   else if (equals (op1, "1")) {
+      char* res [] = {"-inf", "-1", "-0", "+0", "1", "+inf", ""};
+      return res [index (values, valueslen, op2)];
+   }
+   else if (equals (op1, "+inf")) {
+      char* res [] = {"-inf", "-inf", "nan", "nan", "+inf", "+inf", ""};
+      return res [index (values, valueslen, op2)];
+   }
+   else
+      return "?";
+}
+
+char* div_fr_fr (char *op1, char *op2) {
+   // returns op1/op2
+   char *r;
+//   printf (" %s %s", op1, op2);
+   if (equals (op1, "nan") || equals (op2, "nan"))
+      r = "nan";
+   else if (equals (op2, "-inf")) {
+      char* res [] = {"nan", "+0", "+0", "-0", "-0", "nan", ""};
+      r = res [index (values, valueslen, op1)];
+   }
+   else if (equals (op2, "-2")) {
+      char* res [] = {"+inf", "0.5", "+0", "-0", "-0.5", "-inf", ""};
+      if (equals (op1, "-2"))
+         r = "1";
+      else if (equals (op1, "2"))
+         r =  "-1";
+      else if (equals (op1, "0"))
+         r = "0";
+      else
+         r = res [index (values, valueslen, op1)];
+   }
+   else if (equals (op2, "-1")) {
+      char* res [] = {"+inf", "1", "+0", "-0", "-1", "-inf", ""};
+      r = res [index (values, valueslen, op1)];
+   }
+   else if (equals (op2, "-0")) {
+      char* res [] = {"+inf", "+inf", "nan", "nan", "-inf", "-inf", ""};
+      if (equals (op1, "0"))
+         r = "nan";
+      else
+         r = res [index (values, valueslen, op1)];
+   }
+   else if (equals (op2, "+0")) {
+      char* res [] = {"-inf", "-inf", "nan", "nan", "+inf", "+inf", ""};
+      if (equals (op1, "0"))
+         r = "nan";
+      else
+         r = res [index (values, valueslen, op1)];
+   }
+   else if (equals (op2, "1")) {
+      char* res [] = {"-inf", "-1", "-0", "+0", "1", "+inf", ""};
+      r = res [index (values, valueslen, op1)];
+   }
+   else if (equals (op2, "2")) {
+      char* res [] = {"-inf", "-0.5", "-0", "+0", "0.5", "+inf", ""};
+      if (equals (op1, "-2"))
+         r = "-1";
+      else if (equals (op1, "2"))
+         r = "1";
+      else if (equals (op1, "0"))
+         r = "0";
+      else
+         r = res [index (values, valueslen, op1)];
+   }
+   else if (equals (op2, "+inf")) {
+      char* res [] = {"nan", "-0", "-0", "+0", "+0", "nan", ""};
+      r = res [index (values, valueslen, op1)];
+   }
+   else
+      r = "?";
+//   printf ("  %s\n", r);
+   return r;
+}
+
+char* add_fr_fr (char *op1, char *op2) {
+   // returns op1+op2
+   if (equals (op1, "nan") || equals (op2, "nan"))
+      return "nan";
+   else if (equals (op1, "-inf")) {
+      char* res [] = {"-inf", "-inf", "-inf", "-inf", "-inf", "nan", ""};
+      return res [index (values, valueslen, op2)];
+   }
+   else if (equals (op1, "-1")) {
+      char* res [] = {"-inf", "-2", "-1", "-1", "0", "+inf", ""};
+      return res [index (values, valueslen, op2)];
+   }
+   else if (equals (op1, "-0")) {
+      char* res [] = {"-inf", "-1", "-0", "0", "1", "+inf", ""};
+      return res [index (values, valueslen, op2)];
+   }
+   else if (equals (op1, "+0")) {
+      char* res [] = {"-inf", "-1", "0", "+0", "1", "+inf", ""};
+      return res [index (values, valueslen, op2)];
+   }
+   else if (equals (op1, "1")) {
+      char* res [] = {"-inf", "0", "1", "1", "2", "+inf", ""};
+      return res [index (values, valueslen, op2)];
+   }
+   else if (equals (op1, "+inf")) {
+      char* res [] = {"nan", "+inf", "+inf", "+inf", "+inf", "+inf", ""};
+      return res [index (values, valueslen, op2)];
+   }
+   else
+      return "?";
+}
+
+char* neg_fr (char *op) {
+   // returns -op
+   if (index (values, valueslen, op) != -1) {
+      char* res [] = {"+inf", "1", "+0", "-0", "-1", "-inf", "nan"};
+      return res [index (values, valueslen, op)];
+   }
+   else if (equals (op, "0"))
+      return "0";
+   else
+      return "?";
+}
+
+char* sub_fr_fr (char *op1, char *op2) {
+   // returns op1-op2
+   return add_fr_fr (op1, neg_fr (op2));
+}
+
+char* norm (char *opre, char *opim) {
+   // returns opre^2+opim^2
+   if (equals (opre, "nan") || equals (opim, "nan"))
+      if (   (equals (opre, "-inf") || equals (opre, "+inf"))
+          || (equals (opim, "-inf") || equals (opim, "+inf")))
+         return "+inf";
+      else
+         return "nan";
+   else
+      return add_fr_fr (mul_fr_fr (opre, opre), mul_fr_fr (opim, opim));
+}
+
+char* div_re (char *op1re, char *op1im, char *op2re, char *op2im) {
+   // returns real part of complex division
+   if (   equals (op1re, "nan") || equals (op1im, "nan")
+       || equals (op2re, "nan") || equals (op2im, "nan"))
+       return "nan";
+   else if (equals (op2im, "-0") || equals (op2im, "+0"))
+      return div_fr_fr (op1re, op2re);
+   else if (equals (op2re, "-0") || equals (op2re, "+0"))
+      return div_fr_fr (op1im, op2im);
+   else
+      return div_fr_fr (
+         add_fr_fr (mul_fr_fr (op1re, op2re), mul_fr_fr (op1im, op2im)),
+         norm (op2re, op2im));
+}
+
+char* div_im (char *op1re, char *op1im, char *op2re, char *op2im) {
+   // returns imaginary part of complex division
+   if (   equals (op1re, "nan") || equals (op1im, "nan")
+       || equals (op2re, "nan") || equals (op2im, "nan"))
+       return "nan";
+   else if (equals (op2im, "-0") || equals (op2im, "+0"))
+      return div_fr_fr (op1im, op2re);
+   else if (equals (op2re, "-0") || equals (op2re, "+0"))
+      return neg_fr (div_fr_fr (op1re, op2im));
+   else
+      return div_fr_fr (
+         sub_fr_fr (mul_fr_fr (op1im, op2re), mul_fr_fr (op1re, op2im)),
+         norm (op2re, op2im));
+}
+
diff --git a/tools/testgen/write_test_div.c b/tools/testgen/write_test_div.c
new file mode 100644
index 0000000..c00f577
--- /dev/null
+++ b/tools/testgen/write_test_div.c
@@ -0,0 +1,68 @@
+#include "write_test.h"
+
+char* header1 =
+"# Data file for mpc_div.\n"
+"#\n"
+"# Copyright (C) 2008 Andreas Enge\n"
+"#\n"
+;
+char* header2 =
+"# For explanations on the file format, see add.dat.\n\n"
+;
+
+int main () {
+   int op1rei, op1imi, op2rei, op2imi;
+      // index in values table
+   char *op1re, *op1im, *op2re, *op2im;
+
+   printf ("%s%s%s", header1, licence, header2);
+
+   printf ("# special values\n");
+   for (op1rei = 0; op1rei < valueslen; op1rei++)
+   for (op1imi = 0; op1imi < valueslen; op1imi++)
+   for (op2rei = 0; op2rei < valueslen; op2rei++)
+   for (op2imi = 0; op2imi < valueslen; op2imi++)
+   {
+      op1re = values [op1rei];
+      op1im = values [op1imi];
+      op2re = values [op2rei];
+      op2im = values [op2imi];
+      if (   (index (infvalues, infvalueslen, op1re) != -1)
+          || (index (infvalues, infvalueslen, op1im) != -1)
+          || (index (infvalues, infvalueslen, op2re) != -1)
+          || (index (infvalues, infvalueslen, op2im) != -1))
+      {
+         // at least one infinite value
+         printf ("7 %4s 7 %4s   ",
+            div_re (op1re, op1im, op2re, op2im),
+            div_im (op1re, op1im, op2re, op2im));
+         printf ("7 %4s 7 %4s  7 %4s 7 %4s N N\n",
+            op1re, op1im, op2re, op2im);
+      }
+   }
+
+   printf ("\n# finite values\n");
+   for (op1rei = 0; op1rei < valueslen; op1rei++)
+   for (op1imi = 0; op1imi < valueslen; op1imi++)
+   for (op2rei = 0; op2rei < valueslen; op2rei++)
+   for (op2imi = 0; op2imi < valueslen; op2imi++)
+   {
+      op1re = values [op1rei];
+      op1im = values [op1imi];
+      op2re = values [op2rei];
+      op2im = values [op2imi];
+      if (   (index (infvalues, infvalueslen, op1re) == -1)
+          && (index (infvalues, infvalueslen, op1im) == -1)
+          && (index (infvalues, infvalueslen, op2re) == -1)
+          && (index (infvalues, infvalueslen, op2im) == -1))
+      {
+         printf ("7 %4s 7 %4s   ",
+            div_re (op1re, op1im, op2re, op2im),
+            div_im (op1re, op1im, op2re, op2im));
+         printf ("7 %4s 7 %4s  7 %4s 7 %4s N N\n",
+            op1re, op1im, op2re, op2im);
+      }
+   }
+
+   return 1;
+}
diff --git a/tools/testgen/write_test_divfr.c b/tools/testgen/write_test_divfr.c
new file mode 100644
index 0000000..60d182f
--- /dev/null
+++ b/tools/testgen/write_test_divfr.c
@@ -0,0 +1,60 @@
+#include "write_test.h"
+
+char* header1 =
+"# Data file for mpc_div_fr.\n"
+"#\n"
+"# Copyright (C) 2008 Andreas Enge\n"
+"#\n"
+;
+char* header2 =
+"# For explanations on the file format, see add_fr.dat.\n\n"
+;
+
+int main () {
+   int op1rei, op1imi, op2i;
+      // index in values table
+   char *op1re, *op1im, *op2;
+
+   printf ("%s%s%s", header1, licence, header2);
+
+   printf ("# special values\n");
+   for (op1rei = 0; op1rei < valueslen; op1rei++)
+   for (op1imi = 0; op1imi < valueslen; op1imi++)
+   for (op2i = 0; op2i < valueslen; op2i++)
+   {
+      op1re = values [op1rei];
+      op1im = values [op1imi];
+      op2 = values [op2i];
+      if (   (index (infvalues, infvalueslen, op1re) != -1)
+          || (index (infvalues, infvalueslen, op1im) != -1)
+          || (index (infvalues, infvalueslen, op2) != -1))
+      {
+         // at least one infinite value
+         printf ("7 %4s 7 %4s   ",
+            div_fr_fr (op1re, op2), div_fr_fr (op1im, op2));
+         printf ("7 %4s 7 %4s  7 %4s  N N\n",
+            op1re, op1im, op2);
+      }
+   }
+
+   printf ("\n# values with only 1, -1, +0, -0\n");
+   for (op1rei = 0; op1rei < valueslen; op1rei++)
+   for (op1imi = 0; op1imi < valueslen; op1imi++)
+   for (op2i = 0; op2i < valueslen; op2i++)
+   {
+      op1re = values [op1rei];
+      op1im = values [op1imi];
+      op2 = values [op2i];
+      if (   (index (infvalues, infvalueslen, op1re) == -1)
+          && (index (infvalues, infvalueslen, op1im) == -1)
+          && (index (infvalues, infvalueslen, op2) == -1))
+      {
+         printf ("7 %4s 7 %4s   ",
+            div_fr_fr (op1re, op2), div_fr_fr (op1im, op2));
+         printf ("7 %4s 7 %4s  7 %4s  N N\n",
+            op1re, op1im, op2);
+      }
+   }
+
+   return 1;
+}
diff --git a/tools/testgen/write_test_frdiv.c b/tools/testgen/write_test_frdiv.c
new file mode 100644
index 0000000..e7d9b38
--- /dev/null
+++ b/tools/testgen/write_test_frdiv.c
@@ -0,0 +1,63 @@
+#include "write_test.h"
+
+char* header1 =
+"# Data file for mpc_fr_div.\n"
+"#\n"
+"# Copyright (C) 2008 Andreas Enge\n"
+"#\n"
+;
+char* header2 =
+"# For explanations on the file format, see add.dat, and take into\n"
+"# account the necessary modifications for one real argument.\n\n"
+;
+
+int main () {
+   int op1i, op2rei, op2imi;
+      // index in values table
+   char *op1, *op2re, *op2im;
+
+   printf ("%s%s%s", header1, licence, header2);
+
+   printf ("# special values\n");
+   for (op1i = 0; op1i < valueslen; op1i++)
+   for (op2rei = 0; op2rei < valueslen; op2rei++)
+   for (op2imi = 0; op2imi < valueslen; op2imi++)
+   {
+      op1 = values [op1i];
+      op2re = values [op2rei];
+      op2im = values [op2imi];
+      if (   (index (infvalues, infvalueslen, op1) != -1)
+          || (index (infvalues, infvalueslen, op2re) != -1)
+          || (index (infvalues, infvalueslen, op2im) != -1))
+      {
+         // at least one infinite value
+         printf ("7 %4s 7 %4s   ",
+            div_re (op1, "+0", op2re, op2im),
+            div_im (op1, "+0", op2re, op2im));
+         printf ("7 %4s  7 %4s 7 %4s N N\n",
+            op1, op2re, op2im);
+      }
+   }
+
+   printf ("\n# finite values\n");
+   for (op1i = 0; op1i < valueslen; op1i++)
+   for (op2rei = 0; op2rei < valueslen; op2rei++)
+   for (op2imi = 0; op2imi < valueslen; op2imi++)
+   {
+      op1 = values [op1i];
+      op2re = values [op2rei];
+      op2im = values [op2imi];
+      if (   (index (infvalues, infvalueslen, op1) == -1)
+          && (index (infvalues, infvalueslen, op2re) == -1)
+          && (index (infvalues, infvalueslen, op2im) == -1))
+      {
+         printf ("7 %4s 7 %4s   ",
+            div_re (op1, "+0", op2re, op2im),
+            div_im (op1, "+0", op2re, op2im));
+         printf ("7 %4s  7 %4s 7 %4s N N\n",
+            op1, op2re, op2im);
+      }
+   }
+
+   return 1;
+}
diff --git a/tools/testgen/write_test_frsub.c b/tools/testgen/write_test_frsub.c
new file mode 100644
index 0000000..e201ad6
--- /dev/null
+++ b/tools/testgen/write_test_frsub.c
@@ -0,0 +1,60 @@
+#include "write_test.h"
+
+char* header1 =
+"# Data file for mpc_fr_sub.\n"
+"#\n"
+"# Copyright (C) 2008 Andreas Enge\n"
+"#\n"
+;
+char* header2 =
+"# For explanations on the file format, see add_fr.dat.\n\n"
+;
+
+int main () {
+   int op1i, op2rei, op2imi;
+      // index in values table
+   char *op1, *op2re, *op2im;
+
+   printf ("%s%s%s", header1, licence, header2);
+
+   printf ("# special values\n");
+   for (op1i = 0; op1i < valueslen; op1i++)
+   for (op2rei = 0; op2rei < valueslen; op2rei++)
+   for (op2imi = 0; op2imi < valueslen; op2imi++)
+   {
+      op1 = values [op1i];
+      op2re = values [op2rei];
+      op2im = values [op2imi];
+      if (   (index (infvalues, infvalueslen, op1) != -1)
+          || (index (infvalues, infvalueslen, op2re) != -1)
+          || (index (infvalues, infvalueslen, op2im) != -1))
+      {
+         // at least one infinite value
+         printf ("7 %4s 7 %4s   ",
+            sub_fr_fr (op1, op2re), neg_fr (op2im));
+         printf ("7 %4s 7 %4s  7 %4s  N N\n",
+            op1, op2re, op2im);
+      }
+   }
+
+   printf ("\n# values with only 1, -1, +0, -0\n");
+   for (op1i = 0; op1i < valueslen; op1i++)
+   for (op2rei = 0; op2rei < valueslen; op2rei++)
+   for (op2imi = 0; op2imi < valueslen; op2imi++)
+   {
+      op1 = values [op1i];
+      op2re = values [op2rei];
+      op2im = values [op2imi];
+      if (   (index (infvalues, infvalueslen, op1) == -1)
+          && (index (infvalues, infvalueslen, op2re) == -1)
+          && (index (infvalues, infvalueslen, op2im) == -1))
+      {
+         printf ("7 %4s 7 %4s   ",
+            sub_fr_fr (op1, op2re), neg_fr (op2im));
+         printf ("7 %4s 7 %4s  7 %4s  N N\n",
+            op1, op2re, op2im);
+      }
+   }
+
+   return 1;
+}
diff --git a/tools/testgen/write_test_mulfr.c b/tools/testgen/write_test_mulfr.c
new file mode 100644
index 0000000..1dfc577
--- /dev/null
+++ b/tools/testgen/write_test_mulfr.c
@@ -0,0 +1,60 @@
+#include "write_test.h"
+
+char* header1 =
+"# Data file for mpc_mul_fr.\n"
+"#\n"
+"# Copyright (C) 2008 Andreas Enge\n"
+"#\n"
+;
+char* header2 =
+"# For explanations on the file format, see add_fr.dat.\n\n"
+;
+
+int main () {
+   int op1rei, op1imi, op2i;
+      // index in values table
+   char *op1re, *op1im, *op2;
+
+   printf ("%s%s%s", header1, licence, header2);
+
+   printf ("# special values\n");
+   for (op1rei = 0; op1rei < valueslen; op1rei++)
+   for (op1imi = 0; op1imi < valueslen; op1imi++)
+   for (op2i = 0; op2i < valueslen; op2i++)
+   {
+      op1re = values [op1rei];
+      op1im = values [op1imi];
+      op2 = values [op2i];
+      if (   (index (infvalues, infvalueslen, op1re) != -1)
+          || (index (infvalues, infvalueslen, op1im) != -1)
+          || (index (infvalues, infvalueslen, op2) != -1))
+      {
+         // at least one infinite value
+         printf ("7 %4s 7 %4s   ",
+            mul_fr_fr (op1re, op2), mul_fr_fr (op1im, op2));
+         printf ("7 %4s 7 %4s  7 %4s  N N\n",
+            op1re, op1im, op2);
+      }
+   }
+
+   printf ("\n# values with only 1, -1, +0, -0\n");
+   for (op1rei = 0; op1rei < valueslen; op1rei++)
+   for (op1imi = 0; op1imi < valueslen; op1imi++)
+   for (op2i = 0; op2i < valueslen; op2i++)
+   {
+      op1re = values [op1rei];
+      op1im = values [op1imi];
+      op2 = values [op2i];
+      if (   (index (infvalues, infvalueslen, op1re) == -1)
+          && (index (infvalues, infvalueslen, op1im) == -1)
+          && (index (infvalues, infvalueslen, op2) == -1))
+      {
+         printf ("7 %4s 7 %4s   ",
+            mul_fr_fr (op1re, op2), mul_fr_fr (op1im, op2));
+         printf ("7 %4s 7 %4s  7 %4s  N N\n",
+            op1re, op1im, op2);
+      }
+   }
+
+   return 1;
+}
diff --git a/tools/testgen/write_test_norm.c b/tools/testgen/write_test_norm.c
new file mode 100644
index 0000000..2532e8e
--- /dev/null
+++ b/tools/testgen/write_test_norm.c
@@ -0,0 +1,49 @@
+#include "write_test.h"
+
+char* header1 =
+"# Data file for mpc_norm.\n"
+"#\n"
+"# Copyright (C) 2008 Andreas Enge\n"
+"#\n"
+;
+char* header2 =
+"# For explanations on the file format, see abs.dat.\n\n"
+;
+
+int main () {
+   int oprei, opimi;
+      // index in values table
+   char *opre, *opim;
+
+   printf ("%s%s%s", header1, licence, header2);
+
+   printf ("# special values\n");
+   for (oprei = 0; oprei < valueslen; oprei++)
+   for (opimi = 0; opimi < valueslen; opimi++)
+   {
+      opre = values [oprei];
+      opim = values [opimi];
+      if (   (index (infvalues, infvalueslen, opre) != -1)
+          || (index (infvalues, infvalueslen, opim) != -1))
+      {
+         printf ("7 %4s   ", norm (opre, opim));
+         printf ("7 %4s 7 %4s N\n", opre, opim);
+      }
+   }
+
+   printf ("\n# values with only 1, -1, +0, -0\n");
+   for (oprei = 0; oprei < valueslen; oprei++)
+   for (opimi = 0; opimi < valueslen; opimi++)
+   {
+      opre = values [oprei];
+      opim = values [opimi];
+      if (   (index (infvalues, infvalueslen, opre) == -1)
+          && (index (infvalues, infvalueslen, opim) == -1))
+      {
+         printf ("7 %4s   ", norm (opre, opim));
+         printf ("7 %4s 7 %4s N\n", opre, opim);
+      }
+   }
+
+   return 1;
+}
diff --git a/tools/testgen/write_test_subfr.c b/tools/testgen/write_test_subfr.c
new file mode 100644
index 0000000..393355c
--- /dev/null
+++ b/tools/testgen/write_test_subfr.c
@@ -0,0 +1,71 @@
+#include "write_test.h"
+
+char* header1 =
+"# Data file for mpc_sub_fr.\n"
+"#\n"
+"# Copyright (C) 2008 Paul Zimmermann, Andreas Enge\n"
+"#\n"
+;
+char* header2 =
+"# For explanations on the file format, see add_fr.dat.\n\n"
+;
+
+int main () {
+   int op1rei, op1imi, op2i;
+      // index in values table
+   char *op1re, *op1im, *op2;
+
+   printf ("%s%s%s", header1, licence, header2);
+
+   printf ("# special values\n");
+   for (op1rei = 0; op1rei < valueslen; op1rei++)
+   for (op1imi = 0; op1imi < valueslen; op1imi++)
+   for (op2i = 0; op2i < valueslen; op2i++)
+   {
+      op1re = values [op1rei];
+      op1im = values [op1imi];
+      op2 = values [op2i];
+      if (   (index (infvalues, infvalueslen, op1re) != -1)
+          || (index (infvalues, infvalueslen, op1im) != -1)
+          || (index (infvalues, infvalueslen, op2) != -1))
+      {
+         // at least one infinite value
+         printf ("7 %4s 7 %4s   ",
+            sub_fr_fr (op1re, op2), op1im);
+         printf ("7 %4s 7 %4s  7 %4s  N N\n",
+            op1re, op1im, op2);
+      }
+   }
+
+   printf ("\n# values with only 1, -1, +0, -0\n");
+   for (op1rei = 0; op1rei < valueslen; op1rei++)
+   for (op1imi = 0; op1imi < valueslen; op1imi++)
+   for (op2i = 0; op2i < valueslen; op2i++)
+   {
+      op1re = values [op1rei];
+      op1im = values [op1imi];
+      op2 = values [op2i];
+      if (   (index (infvalues, infvalueslen, op1re) == -1)
+          && (index (infvalues, infvalueslen, op1im) == -1)
+          && (index (infvalues, infvalueslen, op2) == -1))
+      {
+         printf ("7 %4s 7 %4s   ",
+            sub_fr_fr (op1re, op2), op1im);
+         printf ("7 %4s 7 %4s  7 %4s  N N\n",
+            op1re, op1im, op2);
+      }
+   }
+
+   printf (
+      "\n# also check the imaginary part is correctly rounded\n"
+      "2 nan 2 16     2 nan 5 17     2 nan   N N\n"
+      "2 nan 2 24     2 nan 5 17     2 nan   N U\n"
+      "2 nan 2 16     2 nan 5 17     2 nan   N D\n"
+      "2 nan 2 16     2 nan 5 17     2 nan   N Z\n"
+      "2 nan 2 -16    2 nan 5 -17    2 nan   N N\n"
+      "2 nan 2 -16    2 nan 5 -17    2 nan   N U\n"
+      "2 nan 2 -24    2 nan 5 -17    2 nan   N D\n"
+      "2 nan 2 -16    2 nan 5 -17    2 nan   N Z\n");
+
+   return 1;
+}