1 files changed, 104 insertions, 0 deletions
diff --git a/libm/double/cpmul.c b/libm/double/cpmul.c
new file mode 100644
index 000000000..3880ac5a1
--- /dev/null
+++ b/libm/double/cpmul.c
@@ -0,0 +1,104 @@
+/*							cpmul.c
+ *
+ *	Multiply two polynomials with complex coefficients
+ *
+ *
+ *
+ * SYNOPSIS:
+ *
+ * typedef struct
+ *		{
+ *		double r;
+ *		double i;
+ *		}cmplx;
+ *
+ * cmplx a[], b[], c[];
+ * int da, db, dc;
+ *
+ * cpmul( a, da, b, db, c, &dc );
+ *
+ *
+ *
+ * DESCRIPTION:
+ *
+ * The two argument polynomials are multiplied together, and
+ * their product is placed in c.
+ *
+ * Each polynomial is represented by its coefficients stored
+ * as an array of complex number structures (see the typedef).
+ * The degree of a is da, which must be passed to the routine
+ * as an argument; similarly the degree db of b is an argument.
+ * Array a has da + 1 elements and array b has db + 1 elements.
+ * Array c must have storage allocated for at least da + db + 1
+ * elements.  The value da + db is returned in dc; this is
+ * the degree of the product polynomial.
+ *
+ * Polynomial coefficients are stored in ascending order; i.e.,
+ * a(x) = a[0]*x**0 + a[1]*x**1 + ... + a[da]*x**da.
+ *
+ *
+ * If desired, c may be the same as either a or b, in which
+ * case the input argument array is replaced by the product
+ * array (but only up to terms of degree da + db).
+ *
+ */
+
+/*							cpmul	*/
+
+typedef struct
+	{
+	double r;
+	double i;
+	}cmplx;
+
+int cpmul( a, da, b, db, c, dc )
+cmplx *a, *b, *c;
+int da, db;
+int *dc;
+{
+int i, j, k;
+cmplx y;
+register cmplx *pa, *pb, *pc;
+
+if( da > db )	/* Know which polynomial has higher degree */
+	{
+	i = da;	/* Swapping is OK because args are on the stack */
+	da = db;
+	db = i;
+	pa = a;
+	a = b;
+	b = pa;
+	}
+	
+k = da + db;
+*dc = k;		/* Output the degree of the product */
+pc = &c[db+1];
+for( i=db+1; i<=k; i++ )	/* Clear high order terms of output */
+	{
+	pc->r = 0;
+	pc->i = 0;
+	pc++;
+	}
+/* To permit replacement of input, work backward from highest degree */
+pb = &b[db];
+for( j=0; j<=db; j++ )
+	{
+	pa = &a[da];
+	pc = &c[k-j];
+	for( i=0; i<da; i++ )
+		{
+		y.r = pa->r * pb->r  -  pa->i * pb->i;	/* cmpx multiply */
+		y.i = pa->r * pb->i  +  pa->i * pb->r;
+		pc->r += y.r;	/* accumulate partial product */
+		pc->i += y.i;
+		pa--;
+		pc--;
+		}
+	y.r = pa->r * pb->r  -  pa->i * pb->i;	/* replace last term,	*/
+	y.i = pa->r * pb->i  +  pa->i * pb->r;	/* ...do not accumulate	*/
+	pc->r = y.r;
+	pc->i = y.i;
+	pb--;
+	}
+  return 0;
+}