bet@ecsvax.UUCP (Bennett E. Todd III) (10/15/85)
#!/bin/sh
# Cut above the preceeding line, or cut here if you must.
# This is a shar archive. Extract with sh, not csh.
# The rest of this file will extract:
# README complex.h cubic.c complex.c
sed 's/^X//' > README << '/*EOF'
XThis is a program to compute the roots to a cubic polynomial in real
Xcoefficients using the closed form solution as described in CRC's
XStandard Math Tables. It is written in C, which might have been a
Xmistake, but was fun anyhow. Everything is reasonably straightforward
Xand documented in comments, until you get to the complex math. C doesn't
Xdo complex math. C doesn't even like to be forced to do complex math.
XWithout operator overloading you cannot express the equations in complex
Xvariables using infix notation; you are left with a functional notation
Xnot (sufficiently, for my tastes) unlike publication LISP syntax. Worse
Xstill, I really truly despise structure passing and return. Therefore I
Xdon't use them; rather I allocate my own structures and pass pointers to
Xthem. Unfortunately, nested procedure calls for intricate expressions
Xwould result in intermediate structures being dropped on the floor.
XTherefore I have the called routine free its arguments. And therefore,
Xif you wish to pass it a program variable (which you presumably don't
Xwant freed) then you must invoke a utility routine to make a copy of
Xthat variable. I have solicited suggestions for names for this
Xas-far-as-I-know unique parameter passing mechanism, and the current
Xwinner (by default) is "call by copy-and-destroy" by D. Gary Grady
X(dgary@ecsvax.UUCP and dpiggy@tucc.BITNET).
X
X-Bennett
/*EOF
ls -l README
sed 's/^X//' > complex.h << '/*EOF'
X/*
X * complex math routines' header file
X */
X/*
X * Written by Bennett Todd
X * Date: 10/15/85 -- initial release
X * This code is released to the public domain; do with it as you wish.
X * If you use it, credit would be appreciated.
X * This is written for and in DeSmet C for the IBM-PC;
X * as far as I know it should be exceedingly portable.
X * Tabsets are at every 4 columns.
X * Please send any bug fixes, enhancements, or general comments to
X * Bennett Todd
X * Duke Computation Center
X * Durham, NC 27706-7756
X * +1 919 684 3695
X * UUCP: ...{decvax,seismo,philabs,ihnp4,akgua}!mcnc!ecsvax!duccpc!bet
X * BITNET: dbtodd@tucc
X */
X
Xtypedef struct complex_struct_type {
X double float real, imaginary;
X struct complex_struct_type *next;
X} COMPLEX;
X
Xextern COMPLEX *Complex(), /* returns a complex made from the real args */
X *C_copy(), /* returns a copy of its argument */
X *C_negate(), /* negate -- monadic */
X *C_sq(), /* square -- monadic */
X *C_pow(), /* raise to real power -- mixed dyadic */
X *C_cb(), /* cube -- monadic */
X *C_sqrt(), /* square root -- monadic */
X *C_cbrt(), /* cube root -- monadic */
X *C_times(), /* product -- dyadic */
X *C_div(), /* ratio -- dyadic */
X *C_plus(), /* sum -- dyadic */
X *C_minus(); /* difference -- dyadic */
/*EOF
ls -l complex.h
sed 's/^X//' > cubic.c << '/*EOF'
X#include <stdio.h>
X#include <math.h>
X#include "complex.h"
X/*
X * cubic -- apply closed-form solution to an arbitrary cubic
X * in real coefficients
X * based on the exposition of the algorithm in the CRC Standard
X * Mathematical Tables, 27'th edition, CRC Press, Boca Raton FL,
X * 1984, on page 9.
X * syntax: cubic o*x^3 + p*x^2 + q*x + r = s
X * Where o, p, q, r, and s are arbitrary constants in format suitable
X * for conversion using "%f" format in scanf().
X * Obviously, a grossly inefficient command-line syntax. Sorry.
X */
X/*
X * Written by Bennett Todd
X * Date: 10/15/85 -- initial release
X * This code is released to the public domain; do with it as you wish.
X * If you use it, credit would be appreciated.
X * This is written for and in DeSmet C for the IBM-PC;
X * as far as I know it should be exceedingly portable.
X * Tabsets are at every 4 columns.
X * Please send any bug fixes, enhancements, or general comments to
X * Bennett Todd
X * Duke Computation Center
X * Durham, NC 27706-7756
X * +1 919 684 3695
X * UUCP: ...{decvax,seismo,philabs,ihnp4,akgua}!mcnc!ecsvax!duccpc!bet
X * BITNET: dbtodd@tucc
X */
X
X#define STR_NEQ(s1, s2) strcmp(s1, s2) /* take THAT you stylists! */
X
Xmain(argc, argv)
Xint argc;
Xchar **argv;
X{
X /* Sorry about the variable names; they are chosen to match CRC */
X double o, p, q, r, s, a, b;
X COMPLEX *A, *B, *root_1, *root_2, *root_3;
X
X /* (attempt to) parse command line */
X if (argc != 10) {
X fprintf(stderr, "syntax: cubic o*x^3 + p*x^2 + q*x + r = s\n");
X fprintf(stderr, "Indicate missing terms with zeroes.\n");
X exit(1);
X }
X if (sscanf(argv[1], "%lf*x^3", &o) != 1) {
X fprintf(stderr, "cubic: cannot parse %s\n", argv[1]);
X exit(1);
X }
X if (sscanf(argv[3], "%lf*x^2", &p) != 1) {
X fprintf(stderr, "cubic: cannot parse %s\n", argv[3]);
X exit(1);
X }
X if (sscanf(argv[5], "%lf*x", &q) != 1) {
X fprintf(stderr, "cubic: cannot parse %s\n", argv[5]);
X exit(1);
X }
X if (sscanf(argv[7], "%lf", &r) != 1) {
X fprintf(stderr, "cubic: cannot parse %s\n", argv[7]);
X exit(1);
X }
X if (sscanf(argv[9], "%lf", &s) != 1) {
X fprintf(stderr, "cubic: cannot parse %s\n", argv[9]);
X exit(1);
X }
X if (STR_NEQ(argv[2], "+") || STR_NEQ(argv[4], "+") ||
X STR_NEQ(argv[6], "+") || STR_NEQ(argv[8], "=")) {
X fprintf(stderr, "syntax: cubic o*x^3 + p*x^2 + q*x + r = s\n");
X fprintf(stderr, "Indicate missing terms with zeroes.\n");
X exit(1);
X }
X printf("\t%s + %s + %s + %s = %s\n",
X argv[1], argv[3], argv[5], argv[7], argv[9]);
X
X /*
X * normalize "ox^3 + px^2 + qx +r = s" into
X * "x^3 + p'x^2 + q'x + r' = 0" by subtracting "s" from both sides
X * and dividing through by "o"
X */
X r = r-s;
X /* sanity check */
X if (o == 0) {
X fprintf(stderr, "cubic: sorry buddy, that's a quadratic.\n");
X exit(1);
X }
X p /= o;
X q /= o;
X r /= o;
X printf("transforms to\n\tx^3 + %lfx^2 + %lfx + %lf = 0\n", p, q, r);
X
X /*
X * According to CRC, the following assignments implement a change
X * of variable from "x^3 + px^2 + qx + r = 0" to
X * "x'^3 + ax' + b = 0" by propogating the substitution
X * "x = x' - p/3". Let's try and remember to convert back...
X */
X a = (3*q - p*p)/3;
X b = (2*p*p*p - 9*p*q + 27*r)/27;
X printf("which in turn transforms to\n\tx^3 + %lfx + %lf = 0\n", a, b);
X printf("via the substitution\n\tx ==> x - %lf\n", p/(float) 3);
X
X /*
X * No motivating explanation given in CRC; it seems that the
X * form of the ultimate solution suggests the following intermediate
X * values to simplify the algebra
X *
X * A = cube_root(sqrt(b*b/(float)4 + a*a*a/(float)27) - b/(float)2);
X * B = - cube_root(sqrt(b*b/(float)4 + a*a*a/(float)27) + b/(float)2);
X *
X * then the roots of the transformed equation "x'^3 + ax' + b = 0"
X * are
X *
X * A+B, -(A+B)/2 + ((A-B)/2)*sqrt(-3), and
X * -(A+B)/2 - ((A-B)/2)*sqrt(-3)
X *
X * However, since CRC immediately goes on to say
X *
X * "If b^2 / 4 + a^3 / 27 > 0 there will be one real root and two
X * complex conjugate roots;
X * If b^2 / 4 + a^3 / 27 = 0 there will be three real roots at least
X * two of which are equal
X * If b^2 / 4 + a^3 / 27 < 0 there will be three real unequal roots"
X *
X * we'd better go ahead and handle complex numbers here. Friends and
X * neighbors, if you are inextricably wedded to infix notation you had
X * better skip the next part. In fact, the next part undoubtledly
X * shouldn't be written in C. What can I say. I'm stubborn.
X *
X * Note: to avoid dropping intermediate values on the floor, without
X * using structure passing and return (which I despise), the complex
X * routines have the following semantics: they take one or two complex
X * arguments as pointers to the structures, and return a pointer to a
X * new structure containing the result. The routines free their complex
X * arguments -- if you wish to pass a program variable (which you don't
X * want freed) to a complex subroutine, pass in C_copy(variable) which
X * doesn't free its argument, and returns a pointer to a copy. What,
X * ladies and gentlemen, is this to be called? Call by copy-and-destroy?
X */
X A = C_cbrt(
X C_plus(
X C_sqrt(
X C_plus(
X C_div(
X C_sq(Complex(b, 0.0)),
X Complex(4.0, 0.0)
X ),
X C_div(
X C_cb(Complex(a, 0.0)),
X Complex(27.0, 0.0)
X )
X )
X ),
X C_negate(
X C_div(
X Complex(b, 0.0),
X Complex(2.0, 0.0)
X )
X )
X )
X );
X B = C_negate(
X C_cbrt(
X C_plus(
X C_sqrt(
X C_plus(
X C_div(
X C_sq(Complex(b, 0.0)),
X Complex(4.0, 0.0)
X ),
X C_div(
X C_cb(Complex(a, 0.0)),
X Complex(27.0, 0.0)
X )
X )
X ),
X C_div(
X Complex(b, 0.0),
X Complex(2.0, 0.0)
X )
X )
X )
X );
X
X root_1 = C_plus(C_copy(A), C_copy(B));
X root_2 = C_plus(
X C_negate(
X C_div(
X C_copy(root_1),
X Complex(2.0, 0.0)
X )
X ),
X C_times(
X C_div(
X C_minus(C_copy(A), C_copy(B)),
X Complex(2.0, 0.0)
X ),
X Complex(0.0, sqrt(3))
X )
X );
X root_3 = C_negate(
X C_plus(
X C_div(
X C_copy(root_1),
X Complex(2.0, 0.0)
X ),
X C_times(
X C_div(
X C_minus(C_copy(A), C_copy(B)),
X Complex(2.0, 0.0)
X ),
X Complex(0.0, sqrt(3))
X )
X )
X );
X
X printf("The original equation has the roots\n");
X printf("\t %lf", root_1->real - p/(float)3);
X if (root_1->imaginary != 0)
X printf("+%lfi", root_1->imaginary);
X printf(", %lf", root_2->real - p/(float)3);
X if (root_2->imaginary != 0)
X printf("+%lfi", root_2->imaginary);
X printf(", and %lf", root_3->real - p/(float)3);
X if (root_3->imaginary != 0)
X printf("+%lfi", root_3->imaginary);
X printf("\n");
X}
/*EOF
ls -l cubic.c
sed 's/^X//' > complex.c << '/*EOF'
X#include <stdio.h>
X#include <math.h>
X#include "complex.h"
X/*
X * Implementation of complex math routines.
X * Routines expect one or two complex arguments, passed as
X * pointers to COMPLEX structs, and return a pointer to
X * a new COMPLEX struct containing the result. So that calls
X * to the complex routines can be nested in intricate expressions
X * without the structs which were created for return values being
X * lost and dropped on the floor, these routines FREE their arguments.
X * You should therefore use the service routine C_copy() to make
X * copies of any constants or program variables you wish to pass into
X * a complex routine.
X */
X
X/*
X * Written by Bennett Todd
X * Date: 10/15/85 -- initial release
X * This code is released to the public domain; do with it as you wish.
X * If you use it, credit would be appreciated.
X * This is written for and in DeSmet C for the IBM-PC;
X * as far as I know it should be exceedingly portable.
X * Tabsets are at every 4 columns.
X * Please send any bug fixes, enhancements, or general comments to
X * Bennett Todd
X * Duke Computation Center
X * Durham, NC 27706-7756
X * +1 919 684 3695
X * UUCP: ...{decvax,seismo,philabs,ihnp4,akgua}!mcnc!ecsvax!duccpc!bet
X * BITNET: dbtodd@tucc
X */
X
X#define PI 3.14159265359
X
Xstatic COMPLEX *complex_free_head = NULL;
X
Xstatic COMPLEX *newcomplex()
X{
X COMPLEX *r;
X
X if (complex_free_head == NULL)
X if ((r = (COMPLEX *) malloc(sizeof(COMPLEX))) == NULL) {
X fprintf(stderr, "complex library malloc() failure\n");
X exit(1);
X } else
X return(r);
X r = complex_free_head;
X complex_free_head = complex_free_head->next;
X return(r);
X}
X
Xstatic freecomplex(c)
XCOMPLEX *c;
X{
X c->next = complex_free_head;
X complex_free_head = c;
X}
X
XCOMPLEX *Complex(x, y)
Xdouble float x, y;
X{
X COMPLEX *r;
X
X r = newcomplex();
X r->real = x;
X r->imaginary = y;
X return(r);
X}
X
XCOMPLEX *C_copy(c)
XCOMPLEX *c;
X{
X COMPLEX *r;
X
X r = newcomplex();
X r->real = c->real;
X r->imaginary = c->imaginary;
X return(r);
X}
X
XCOMPLEX *C_negate(c)
XCOMPLEX *c;
X{
X c->real = -c->real;
X c->imaginary = -c->imaginary;
X return(c);
X}
X
XCOMPLEX *C_sq(c)
XCOMPLEX *c;
X{
X COMPLEX *r;
X
X r = newcomplex();
X r->real = c->real*c->real - c->imaginary*c->imaginary;
X r->imaginary = 2*c->real*c->imaginary;
X freecomplex(c);
X return(r);
X}
X
XCOMPLEX *C_pow(c, p)
XCOMPLEX *c;
Xdouble float p;
X{
X double float magnitude, angle;
X
X magnitude = sqrt(c->real*c->real + c->imaginary*c->imaginary);
X if (c->real != 0)
X angle = atan(c->imaginary/c->real);
X else if (c->imaginary > 0)
X angle = PI/(float)2;
X else if (c->imaginary < 0)
X angle = -PI/(float)2;
X else /* What is the angle of 0+0i? */
X angle = 0.0;
X if (c->real < 0)
X angle += PI;
X c->real = pow(magnitude, p) * cos(p*angle);
X c->imaginary = pow(magnitude, p) * sin(p*angle);
X return(c);
X}
X
XCOMPLEX *C_cb(c)
XCOMPLEX *c;
X{
X return(C_pow(c, 3.0));
X}
X
XCOMPLEX *C_sqrt(c)
XCOMPLEX *c;
X{
X return(C_pow(c, 1/(float)2));
X}
X
XCOMPLEX *C_cbrt(c)
XCOMPLEX *c;
X{
X return(C_pow(c, 1/(float)3));
X}
X
XCOMPLEX *C_times(a, b)
XCOMPLEX *a, *b;
X{
X COMPLEX *r;
X
X r = newcomplex();
X r->real = a->real*b->real - a->imaginary*b->imaginary;
X r->imaginary = a->real*b->imaginary + a->imaginary*b->real;
X freecomplex(a);
X freecomplex(b);
X return(r);
X}
X
XCOMPLEX *C_div(a, b)
XCOMPLEX *a, *b;
X{
X COMPLEX *r;
X double float tmp;
X
X r = newcomplex();
X tmp = b->real*b->real + b->imaginary*b->imaginary;
X if (tmp == 0) {
X fprintf(stderr, "complex divide overflow\n");
X exit(1);
X }
X r->real = (a->real*b->real + a->imaginary*b->imaginary)/tmp;
X r->imaginary = (b->real*a->imaginary - a->real*b->imaginary)/tmp;
X freecomplex(a);
X freecomplex(b);
X return(r);
X}
X
XCOMPLEX *C_plus(a, b)
XCOMPLEX *a, *b;
X{
X a->real = a->real + b->real;
X a->imaginary = a->imaginary + b->imaginary;
X freecomplex(b);
X return(a);
X}
X
XCOMPLEX *C_minus(a, b)
XCOMPLEX *a, *b;
X{
X a->real = a->real - b->real;
X a->imaginary = a->imaginary - b->imaginary;
X freecomplex(b);
X return(a);
X}
/*EOF
ls -l complex.c
exit
--
"Hypocrisy is the vaseline of social intercourse." (Who said that?)
Bennett Todd -- Duke Computation Center, Durham, NC 27706-7756; (919) 684-3695
UUCP: ...{decvax,seismo,philabs,ihnp4,akgua}!mcnc!ecsvax!duccpc!bet