0001: /*************************************************************************
0002: * *
0003: * 1) This source code file, in unmodified form, and compiled classes *
0004: * derived from it can be used and distributed without restriction, *
0005: * including for commercial use. (Attribution is not required *
0006: * but is appreciated.) *
0007: * *
0008: * 2) Modified versions of this file can be made and distributed *
0009: * provided: the modified versions are put into a Java package *
0010: * different from the original package, edu.hws; modified *
0011: * versions are distributed under the same terms as the original; *
0012: * and the modifications are documented in comments. (Modification *
0013: * here does not include simply making subclasses that belong to *
0014: * a package other than edu.hws, which can be done without any *
0015: * restriction.) *
0016: * *
0017: * David J. Eck *
0018: * Department of Mathematics and Computer Science *
0019: * Hobart and William Smith Colleges *
0020: * Geneva, New York 14456, USA *
0021: * Email: eck@hws.edu WWW: http://math.hws.edu/eck/ *
0022: * *
0023: *************************************************************************/package edu.hws.jcm.data;
0024:
0025: /**
0026: * An ExprssionProgram represents a mathematical expression such as "3" or "sin(x^2)", stored
0027: * in the form of a program for a stack machine. The program consists of a sequence of
0028: * commands that, when executed, will compute the value of the expression.
0029: *
0030: * <p>Each command is encoded as an integer. There are three types of commands that can
0031: * occur: (1) A negative integer must be one of the 36 constant values PLUS, MINUS,..., CUBERT.
0032: * These constants represent unary and binary operators and standard functions. (2) An integer in the range
0033: * 0 <= n < 0x3FFFFFFF encodes an operation of the form "push a constant onto the stack".
0034: * The constant that is being pushed is encoded as an index in the array "constant",
0035: * which is a private member of this class that holds all the constants that
0036: * occur in this ExpressionProgram. (3) An integer >= 0x3FFFFFFF represents
0037: * an ExpressionCommand object. When 0x3FFFFFFF is subtracted from the integer, the
0038: * result is an index into the array "command", which is a private member of this
0039: * class that holds all the ExpressionCommands that occur in this
0040: * ExpressionProgram.
0041: */
0042: public class ExpressionProgram implements Expression {
0043:
0044: /**
0045: * Code for a unary or binary operator or a standard function.
0046: */
0047: public static final int
0048:
0049: PLUS = -1, MINUS = -2, TIMES = -3, DIVIDE = -4, POWER = -5,
0050: EQ = -6, NE = -7, LT = -8, GT = -9, LE = -10, GE = -11,
0051: AND = -12, OR = -13, NOT = -14, UNARY_MINUS = -15,
0052: FACTORIAL = -16, SIN = -17, COS = -18, TAN = -19,
0053: COT = -20, SEC = -21, CSC = -22, ARCSIN = -23,
0054: ARCCOS = -24, ARCTAN = -25, ABS = -26, SQRT = -27,
0055: EXP = -28, LN = -29, LOG2 = -30, LOG10 = -31, TRUNC = -32,
0056: ROUND = -33, FLOOR = -34, CEILING = -35, CUBERT = -36;
0057:
0058: /**
0059: * If this is non-null, it is used as the print string
0060: * for this expression in the toString() method. (When an
0061: * expression is created by a Parser by parsing a string,
0062: * the parse stores that string in this variable.)
0063: */
0064: public String sourceString;
0065:
0066: private int[] prog = new int[1]; // Contains the commands that make up this program.
0067: // (Constants and ExpressionCommands are encoded as
0068: // positive numbers, while negative numbers are
0069: // opCodes from the above list of constants.)
0070: // This array expands as necessary.
0071:
0072: private int progCt; // The number of commands in prog.
0073:
0074: //--------------------- Methods for creating the program ------------------------------
0075:
0076: /**
0077: * Default constructor creates an initially empty program.
0078: */
0079: public ExpressionProgram() {
0080: }
0081:
0082: /**
0083: * Adds com as the next command in the program. Among other things, for example, com can
0084: * be a Variable or Constant. In that case, the meaning of the command is the stack
0085: * operation "push (value of com)".
0086: *
0087: * @param com added as next command in the program.
0088: */
0089: public void addCommandObject(ExpressionCommand com) {
0090: int loc = findCommand(com);
0091: addCommandCode(loc + 0x3FFFFFFF);
0092: }
0093:
0094: /**
0095: * Add the number d as the next command in the program. The meaning of this command is
0096: * actually the stack operation "push d".
0097: *
0098: * @param d added as next command in program.
0099: */
0100: public void addConstant(double d) {
0101: int loc = findConstant(d);
0102: addCommandCode(loc);
0103: }
0104:
0105: /**
0106: * Add a command code to the program, where code is one of the opCode constants
0107: * that are public final members of this class, from CUBERT to PLUS. Each code
0108: * represents either a binary or unary operator or a standard function that operates on the stack by
0109: * poping its argument(s) from the stack, perfroming the operation, and pushing
0110: * the result back onto the stack.
0111: */
0112: public void addCommand(int code) {
0113: if (code >= 0 || code < CUBERT)
0114: throw new IllegalArgumentException(
0115: "Internal Error. Illegal command code.");
0116: addCommandCode(code);
0117: }
0118:
0119: /**
0120: * To save space, cut the arrays that holds the program data down to the actual
0121: * amount of data that they contain. This should be called after the complete
0122: * program has been generated.
0123: */
0124: public void trim() {
0125: if (progCt != prog.length) {
0126: int[] temp = new int[progCt];
0127: System.arraycopy(prog, 0, temp, 0, progCt);
0128: prog = temp;
0129: }
0130: if (commandCt != command.length) {
0131: ExpressionCommand[] temp = new ExpressionCommand[commandCt];
0132: System.arraycopy(command, 0, temp, 0, commandCt);
0133: command = temp;
0134: }
0135: if (constantCt != constant.length) {
0136: double[] temp = new double[constantCt];
0137: System.arraycopy(constant, 0, temp, 0, constantCt);
0138: constant = temp;
0139: }
0140: }
0141:
0142: private void addCommandCode(int code) { // System.out.println("Add code " + code + " " + progCt);
0143: if (progCt == prog.length) {
0144: int[] temp = new int[prog.length * 2];
0145: System.arraycopy(prog, 0, temp, 0, progCt);
0146: prog = temp;
0147: }
0148: prog[progCt++] = code;
0149: }
0150:
0151: //------------------- Methods for evaluating the program -----------------------
0152:
0153: private Cases cases; // If this is non-null while the expression is being evaluated, then
0154: // information about "cases" is recorded in it as the evaluation is
0155: // done. See the Cases class for more information.
0156:
0157: private StackOfDouble stack = new StackOfDouble();
0158:
0159: // An ExpressionProgram is a program for a stack machien. This is the
0160: // stack that is used when the program is executed.
0161:
0162: /**
0163: * Run the ExprssionProgram and return the value that it computes.
0164: */
0165: synchronized public double getVal() {
0166: cases = null;
0167: return basicGetValue();
0168: }
0169:
0170: /**
0171: * Run the ExprssionProgram and return the value that it computes.
0172: * If the Cases object, c, is non-null, then information about "cases" is recorded in c.
0173: * This information can be used to help detect possible "discontinuities"
0174: * between two evaluations. See the Cases class for more information.
0175: */
0176: synchronized public double getValueWithCases(Cases c) {
0177: cases = c;
0178: double d = basicGetValue();
0179: cases = null;
0180: return d;
0181: }
0182:
0183: private double basicGetValue() {
0184: // Run the program and return its value.
0185: stack.makeEmpty();
0186: for (int pc = 0; pc < progCt; pc++) {
0187: int code = prog[pc];
0188: if (code < 0) {
0189: double ans = applyCommandCode(code);
0190: if (Double.isNaN(ans) || Double.isInfinite(ans)) {
0191: if (cases != null)
0192: cases.addCase(0);
0193: return Double.NaN;
0194: }
0195: stack.push(ans);
0196: } else if (code < 0x3FFFFFFF) {
0197: stack.push(constant[code]);
0198: } else {
0199: command[code - 0x3FFFFFFF].apply(stack, cases);
0200: }
0201: }
0202: if (stack.size() != 1)
0203: throw new IllegalArgumentException(
0204: "Internal Error: Improper stack state after expression evaluation.");
0205: double val = stack.pop();
0206: if (cases != null)
0207: cases.addCase(Double.isNaN(val) ? 0 : 1);
0208: return val;
0209: }
0210:
0211: private void addCase(int c) {
0212: // If the member variable cases is not null, record c as the next item of case information in cases.
0213: if (cases != null)
0214: cases.addCase(c);
0215: }
0216:
0217: /**
0218: * Apply the stack operation represented by code (a number < 0) to the stack.
0219: */
0220: protected double applyCommandCode(int code) {
0221: double ans;
0222: if (code < OR)
0223: ans = eval(code, stack.pop());
0224: else
0225: ans = eval(code, stack.pop(), stack.pop());
0226: return ans;
0227: }
0228:
0229: private double eval(int commandCode, double x) {
0230: // Compute the value of the unary operator represented by commandCode, when
0231: // applied to the value x.
0232: switch (commandCode) {
0233: case NOT:
0234: return (x == 0) ? 1 : 0;
0235: case UNARY_MINUS:
0236: return -x;
0237: case FACTORIAL:
0238: return factorial(x);
0239: case SIN:
0240: return Math.sin(x);
0241: case COS:
0242: return Math.cos(x);
0243: case TAN:
0244: addCase((int) Math.floor((x - Math.PI / 2.0) / Math.PI));
0245: return Math.tan(x);
0246: case COT:
0247: addCase((int) Math.floor(x / Math.PI));
0248: return Math.cos(x) / Math.sin(x);
0249: case SEC:
0250: addCase((int) Math.floor((x - Math.PI / 2.0) / Math.PI));
0251: return 1 / Math.cos(x);
0252: case CSC:
0253: addCase((int) Math.floor(x / Math.PI));
0254: return 1 / Math.sin(x);
0255: case ARCSIN:
0256: return Math.asin(x);
0257: case ARCCOS:
0258: return Math.acos(x);
0259: case ARCTAN:
0260: return Math.atan(x);
0261: case ABS:
0262: addCase((x > 0) ? 1 : ((x < 0) ? -1 : 0));
0263: return Math.abs(x);
0264: case SQRT:
0265: return (x < 0) ? Double.NaN : Math.sqrt(x);
0266: case EXP:
0267: return Math.exp(x);
0268: case LN:
0269: return (x <= 0) ? Double.NaN : Math.log(x);
0270: case LOG2:
0271: return (x <= 0) ? Double.NaN : Math.log(x) / Math.log(2);
0272: case LOG10:
0273: return (x <= 0) ? Double.NaN : Math.log(x) / Math.log(10);
0274: case TRUNC:
0275: addCase((int) x);
0276: return (long) x;
0277: case ROUND:
0278: addCase((int) Math.floor(x + 0.5));
0279: return Math.floor(x + 0.5);
0280: case FLOOR:
0281: addCase((int) Math.floor(x));
0282: return Math.floor(x);
0283: case CEILING:
0284: addCase((int) Math.floor(x));
0285: return Math.ceil(x);
0286: case CUBERT:
0287: addCase((x > 0) ? 1 : -1);
0288: return (x >= 0) ? Math.pow(x, 1.0 / 3.0) : -Math.pow(-x,
0289: 1.0 / 3.0);
0290: default:
0291: return Double.NaN;
0292: }
0293: }
0294:
0295: private double factorial(double x) {
0296: // Compute x!. x is rounded to the nearest integer. If x > 170, then the
0297: // answer is too big to represent in a value of type double, so the value
0298: // is given as Double.NaN.
0299: if (x <= -0.5 || x > 170.5) {
0300: addCase(-1);
0301: return Double.NaN;
0302: }
0303: int n = (int) x;
0304: addCase(n);
0305: double ans = 1;
0306: for (int i = 1; i <= n; i++)
0307: ans *= i;
0308: return ans;
0309: }
0310:
0311: private double eval(int commandCode, double y, double x) {
0312: // Compute the value of the unary operator represented by commandCode, when
0313: // applied to the values x and y. (Note that the second operand comes
0314: // first in the parameter list, since that is the way they were popped
0315: // off the stack.)
0316: switch (commandCode) {
0317: case PLUS:
0318: return x + y;
0319: case MINUS:
0320: return x - y;
0321: case TIMES:
0322: return x * y;
0323: case DIVIDE:
0324: addCase((y > 0) ? 1 : ((y < 0) ? -1 : 0));
0325: return x / y;
0326: case POWER:
0327: return Math.pow(x, y);
0328: case EQ:
0329: return (x == y) ? 1 : 0;
0330: case NE:
0331: return (x != y) ? 1 : 0;
0332: case GT:
0333: return (x > y) ? 1 : 0;
0334: case LT:
0335: return (x < y) ? 1 : 0;
0336: case GE:
0337: return (x >= y) ? 1 : 0;
0338: case LE:
0339: return (x <= y) ? 1 : 0;
0340: case AND:
0341: return ((x != 0) && (y != 0)) ? 1 : 0;
0342: case OR:
0343: return ((x != 0) || (y != 0)) ? 1 : 0;
0344: default:
0345: return Double.NaN;
0346: }
0347: }
0348:
0349: //------------- Getting the print string of this expression ------------------------
0350:
0351: /**
0352: * If a source string has been saved, use it as the print string. (When a Parser creates
0353: * an expression by parsing a string, it saves the source string in the ExpressionProgram.)
0354: * Otherwise, construct the print string based on the commands in the program.
0355: */
0356: public String toString() {
0357: if (sourceString != null)
0358: return sourceString;
0359: else {
0360: StringBuffer buffer = new StringBuffer();
0361: appendOutputString(progCt - 1, buffer);
0362: return buffer.toString();
0363: }
0364: }
0365:
0366: /**
0367: * Add a string representing part of the expression to the output buffer.
0368: * You probably are only interested in this if you write a ParserExtension
0369: * or ExpressionCommand.
0370: * (The command at position index in the program represents a
0371: * subexpression. It could be a constant or a variable, for example,
0372: * which is complete subexpression in itself. Or it could be the
0373: * final operator in a larger subexpression. In that case, the operands,
0374: * which are located in lower positions in the program, are considered to
0375: * be part of the expression. This routine appends a print string
0376: * for the entire subexpression to the buffer. When this is called
0377: * with index = progCt-1 (the last command in the program), it processes
0378: * the entire program.
0379: * Note that the hard part here is deciding when to put in parentheses.
0380: * This is done based on the precedence of the operators. The result is not always pretty.
0381: */
0382: public void appendOutputString(int index, StringBuffer buffer) {
0383: if (prog[index] >= 0x3FFFFFFF) {
0384: command[prog[index] - 0x3FFFFFFF].appendOutputString(this ,
0385: index, buffer);
0386: } else if (prog[index] >= 0) {
0387: buffer.append(NumUtils.realToString(constant[prog[index]]));
0388: } else if (prog[index] >= OR) {
0389: int op2 = index - 1; // Position of command representing the second operand.
0390: int op1 = op2 - extent(op2); // Position of command representing the first operand.
0391: if (precedence(prog[op1]) < precedence(prog[index])
0392: || (prog[index] == POWER && precedence(prog[op1]) == precedence(prog[index]))) {
0393: buffer.append('(');
0394: appendOutputString(op1, buffer);
0395: buffer.append(')');
0396: } else
0397: appendOutputString(op1, buffer);
0398: switch (prog[index]) {
0399: case PLUS:
0400: buffer.append(" + ");
0401: break;
0402: case MINUS:
0403: buffer.append(" - ");
0404: break;
0405: case TIMES:
0406: buffer.append("*");
0407: break;
0408: case DIVIDE:
0409: buffer.append("/");
0410: break;
0411: case POWER:
0412: buffer.append("^");
0413: break;
0414: case AND:
0415: buffer.append(" AND ");
0416: break;
0417: case OR:
0418: buffer.append(" OR ");
0419: break;
0420: case EQ:
0421: buffer.append(" = ");
0422: break;
0423: case NE:
0424: buffer.append(" <> ");
0425: break;
0426: case GE:
0427: buffer.append(" >= ");
0428: break;
0429: case LE:
0430: buffer.append(" <= ");
0431: break;
0432: case GT:
0433: buffer.append(" > ");
0434: break;
0435: case LT:
0436: buffer.append(" < ");
0437: break;
0438: }
0439: if (prog[op2] == UNARY_MINUS
0440: || precedence(prog[op2]) < precedence(prog[index])
0441: || ((prog[index] == MINUS || prog[index] == DIVIDE) && precedence(prog[op2]) == precedence(prog[index]))) {
0442: buffer.append('(');
0443: appendOutputString(op2, buffer);
0444: buffer.append(')');
0445: } else
0446: appendOutputString(op2, buffer);
0447: } else if (prog[index] <= SIN) {
0448: buffer.append(StandardFunction
0449: .standardFunctionName(prog[index]));
0450: buffer.append('(');
0451: appendOutputString(index - 1, buffer);
0452: buffer.append(')');
0453: } else if (prog[index] == UNARY_MINUS) {
0454: buffer.append('-');
0455: if (precedence(prog[index - 1]) <= precedence(UNARY_MINUS)) {
0456: buffer.append('(');
0457: appendOutputString(index - 1, buffer);
0458: buffer.append(')');
0459: } else
0460: appendOutputString(index - 1, buffer);
0461: } else if (prog[index] == NOT) {
0462: buffer.append("NOT (");
0463: appendOutputString(index - 1, buffer);
0464: buffer.append(')');
0465: } else if (prog[index] == FACTORIAL) {
0466: if (prog[index - 1] >= 0
0467: && (prog[index - 1] < 0x3FFFFFFF
0468: || command[prog[index - 1] - 0x3FFFFFFF] instanceof Variable || command[prog[index - 1] - 0x3FFFFFFF] instanceof Constant)) {
0469: appendOutputString(index - 1, buffer);
0470: } else {
0471: buffer.append('(');
0472: appendOutputString(index - 1, buffer);
0473: buffer.append(')');
0474: }
0475: buffer.append('!');
0476: }
0477: }
0478:
0479: private int precedence(int opcode) {
0480: // Return the precedence of the operator. This is used
0481: // by the prceding method to decide when parentheses are needed.
0482: if (opcode >= 0) {
0483: if (opcode >= 0x3FFFFFFF
0484: && (command[opcode - 0x3FFFFFFF] instanceof ConditionalExpression))
0485: return 0;
0486: else
0487: return 7;
0488: } else
0489: switch (opcode) {
0490: case FACTORIAL:
0491: case OR:
0492: case AND:
0493: return 1;
0494: case GE:
0495: case LE:
0496: case GT:
0497: case EQ:
0498: case NE:
0499: case LT:
0500: return 2;
0501: case PLUS:
0502: case MINUS:
0503: case UNARY_MINUS:
0504: return 3;
0505: case TIMES:
0506: case DIVIDE:
0507: return 4;
0508: case POWER:
0509: return 6;
0510: default:
0511: return 7;
0512: }
0513: }
0514:
0515: //----------------------- Methods for computing the derivative -------------------------
0516:
0517: /**
0518: * Compute the derivative of this expression with respect to the Variable wrt.
0519: * The value returned is actually an ExpressionProgram.
0520: */
0521: public Expression derivative(Variable wrt) {
0522: ExpressionProgram deriv = new ExpressionProgram();
0523: compileDerivative(progCt - 1, deriv, wrt);
0524: deriv.trim();
0525: return deriv;
0526: }
0527:
0528: /**
0529: * The command at position index in the program represents a subexpression of
0530: * the whole expression. This routine adds commands to deriv for computing the
0531: * derivative of that subexpression with respect to the variable wrt.
0532: * You probably are not interested in this unless you write a ParserExtension
0533: * or an ExpressionCommand.
0534: */
0535: public void compileDerivative(int index, ExpressionProgram deriv,
0536: Variable wrt) {
0537: if (!dependsOn(index, wrt))
0538: deriv.addConstant(0);
0539: else if (prog[index] >= 0x3FFFFFFF)
0540: command[prog[index] - 0x3FFFFFFF].compileDerivative(this ,
0541: index, deriv, wrt);
0542: else if (prog[index] >= 0)
0543: deriv.addConstant(0);
0544: else if (prog[index] >= POWER) {
0545: int indexOp2 = index - 1;
0546: int indexOp1 = indexOp2 - extent(indexOp2);
0547: doBinDeriv(prog[index], indexOp1, indexOp2, deriv, wrt);
0548: } else if (prog[index] <= SIN) {
0549: doFuncDeriv(prog[index], index - 1, deriv, wrt);
0550: } else if (prog[index] == UNARY_MINUS) {
0551: compileDerivative(index - 1, deriv, wrt);
0552: deriv.addCommand(UNARY_MINUS);
0553: } else if (prog[index] == FACTORIAL) {
0554: deriv.addConstant(Double.NaN);
0555: } else if (prog[index] >= NOT)
0556: throw new IllegalArgumentException(
0557: "Internal Error: Attempt to take the derivative of a logical-valued expression.");
0558: else
0559: throw new IllegalArgumentException(
0560: "Internal Error: Unknown opcode.");
0561: }
0562:
0563: /**
0564: * The command at position index in the program represents a subexpression of
0565: * the whole expression. This routine finds and returns the number of commands
0566: * in the program that are part of that subexpression. That is, the subexpresssion
0567: * occupies the part of the program between index - extent + 1 and index.
0568: * You probably are not interested in this unless you write a ParserExtension
0569: * or an ExpressionCommand.
0570: */
0571: public int extent(int index) {
0572: if (prog[index] <= NOT)
0573: return 1 + extent(index - 1);
0574: else if (prog[index] < 0) {
0575: int extentOp1 = extent(index - 1); // Extent of second operand (which is on top in the program).
0576: int extentOp2 = extent(index - 1 - extentOp1); // Extent of second operand.
0577: return extentOp1 + extentOp2 + 1;
0578: } else if (prog[index] < 0x3FFFFFFF)
0579: return 1;
0580: else
0581: return command[prog[index] - 0x3FFFFFFF]
0582: .extent(this , index);
0583: }
0584:
0585: /**
0586: * The command at position index in the program represents a subexpression of
0587: * the whole expression. This routine copies the commands for the entire
0588: * subexpression to the destination program.
0589: * You probably are not interested in this unless you write a ParserExtension
0590: * or an ExpressionCommand.
0591: */
0592: public void copyExpression(int index, ExpressionProgram destination) {
0593: int size = extent(index);
0594: for (int i = index - size + 1; i <= index; i++) {
0595: if (prog[i] < 0)
0596: destination.addCommand(prog[i]);
0597: else if (prog[i] >= 0x3FFFFFFF)
0598: destination
0599: .addCommandObject(command[prog[i] - 0x3FFFFFFF]);
0600: else
0601: destination.addConstant(constant[prog[i]]);
0602: }
0603: }
0604:
0605: /**
0606: * The command at position index in the program represents a subexpression of
0607: * the whole expression. If that subexpression includes some dependence on
0608: * the variable x, then true is returned. If the subexpression is constant
0609: * with respect to the variable x, then false is returned.
0610: * You probably are not interested in this unless you write a ParserExtension
0611: * or an ExpressionCommand.
0612: */
0613: public boolean dependsOn(int index, Variable x) {
0614: int size = extent(index);
0615: for (int i = index - size + 1; i <= index; i++) {
0616: if (prog[i] >= 0x3FFFFFFF) {
0617: ExpressionCommand c = command[prog[i] - 0x3FFFFFFF];
0618: if (c == x || c.dependsOn(x))
0619: return true;
0620: }
0621: }
0622: return false;
0623: }
0624:
0625: /**
0626: * Checks whether the expression as a whole has any dependence on the variable x.
0627: */
0628: public boolean dependsOn(Variable x) {
0629: return dependsOn(progCt - 1, x);
0630: }
0631:
0632: private void doBinDeriv(int opCode, int op1, int op2,
0633: ExpressionProgram deriv, Variable wrt) {
0634: // Add commands to deriv to compute the derivative of a subexpression of this program
0635: // with respect to the variable wrt, where the main operation in the subexpression is
0636: // the binary operator opCode, the position of the first operand of the operator is
0637: // at the index op1 in the program, and the position of second operand is op2.
0638: switch (opCode) {
0639: case PLUS:
0640: if (!dependsOn(op1, wrt))
0641: compileDerivative(op2, deriv, wrt);
0642: else if (!dependsOn(op2, wrt))
0643: compileDerivative(op1, deriv, wrt);
0644: else {
0645: compileDerivative(op1, deriv, wrt);
0646: compileDerivative(op2, deriv, wrt);
0647: deriv.addCommand(PLUS);
0648: }
0649: break;
0650: case MINUS:
0651: if (!dependsOn(op1, wrt)) {
0652: compileDerivative(op2, deriv, wrt);
0653: deriv.addCommand(UNARY_MINUS);
0654: } else if (!dependsOn(op2, wrt))
0655: compileDerivative(op1, deriv, wrt);
0656: else {
0657: compileDerivative(op1, deriv, wrt);
0658: compileDerivative(op2, deriv, wrt);
0659: deriv.addCommand(MINUS);
0660: }
0661: break;
0662: case TIMES: {
0663: int cases = 0;
0664: if (dependsOn(op2, wrt)) {
0665: copyExpression(op1, deriv);
0666: if (prog[op2] < 0x3FFFFFFF
0667: || command[prog[op2] - 0x3FFFFFFF] != wrt) {
0668: compileDerivative(op2, deriv, wrt);
0669: deriv.addCommand(TIMES);
0670: }
0671: cases++;
0672: }
0673: if (dependsOn(op1, wrt)) {
0674: copyExpression(op2, deriv);
0675: if (prog[op1] < 0x3FFFFFFF
0676: || command[prog[op1] - 0x3FFFFFFF] != wrt) {
0677: compileDerivative(op1, deriv, wrt);
0678: deriv.addCommand(TIMES);
0679: }
0680: cases++;
0681: }
0682: if (cases == 2)
0683: deriv.addCommand(PLUS);
0684: }
0685: break;
0686: case DIVIDE:
0687: if (!dependsOn(op2, wrt)) {
0688: compileDerivative(op1, deriv, wrt);
0689: copyExpression(op2, deriv);
0690: deriv.addCommand(DIVIDE);
0691: } else if (!dependsOn(op1, wrt)) {
0692: copyExpression(op1, deriv);
0693: deriv.addCommand(UNARY_MINUS);
0694: copyExpression(op2, deriv);
0695: deriv.addConstant(2);
0696: deriv.addCommand(POWER);
0697: deriv.addCommand(DIVIDE);
0698: if (prog[op2] < 0x3FFFFFFF
0699: || command[prog[op2] - 0x3FFFFFFF] != wrt) {
0700: compileDerivative(op2, deriv, wrt);
0701: deriv.addCommand(TIMES);
0702: }
0703: } else {
0704: copyExpression(op2, deriv);
0705: if (prog[op1] < 0x3FFFFFFF
0706: || command[prog[op1] - 0x3FFFFFFF] != wrt) {
0707: compileDerivative(op1, deriv, wrt);
0708: deriv.addCommand(TIMES);
0709: }
0710: copyExpression(op1, deriv);
0711: if (prog[op2] < 0x3FFFFFFF
0712: || command[prog[op2] - 0x3FFFFFFF] != wrt) {
0713: compileDerivative(op2, deriv, wrt);
0714: deriv.addCommand(TIMES);
0715: }
0716: deriv.addCommand(MINUS);
0717: copyExpression(op2, deriv);
0718: deriv.addConstant(2);
0719: deriv.addCommand(POWER);
0720: deriv.addCommand(DIVIDE);
0721: }
0722: break;
0723: case POWER:
0724: if (!dependsOn(op2, wrt)) {
0725: copyExpression(op2, deriv);
0726: copyExpression(op1, deriv);
0727: if (prog[op2] >= 0 && prog[op2] < 0x3FFFFFFF) {
0728: if (constant[prog[op2]] != 2) {
0729: deriv.addConstant(constant[prog[op2]] - 1);
0730: deriv.addCommand(POWER);
0731: }
0732: } else if (prog[op2] == UNARY_MINUS
0733: && prog[op2 - 1] >= 0
0734: && prog[op2 - 1] < 0x3FFFFFFF) {
0735: deriv.addConstant(constant[prog[op2 - 1]] + 1);
0736: deriv.addCommand(UNARY_MINUS);
0737: deriv.addCommand(POWER);
0738: } else {
0739: copyExpression(op2, deriv);
0740: deriv.addConstant(1);
0741: deriv.addCommand(MINUS);
0742: deriv.addCommand(POWER);
0743: }
0744: deriv.addCommand(TIMES);
0745: if (prog[op1] < 0x3FFFFFFF
0746: || command[prog[op1] - 0x3FFFFFFF] != wrt) {
0747: compileDerivative(op1, deriv, wrt);
0748: deriv.addCommand(TIMES);
0749: }
0750: } else if (!dependsOn(op1, wrt)) {
0751: copyExpression(op1, deriv);
0752: copyExpression(op2, deriv);
0753: deriv.addCommand(POWER);
0754: if (!(prog[op1] >= 0x3FFFFFFF
0755: && command[prog[op1] - 0x3FFFFFFF] instanceof Constant && ((Constant) command[prog[op1] - 0x3FFFFFFF])
0756: .getVal() == Math.E)) {
0757: copyExpression(op1, deriv);
0758: deriv.addCommand(LN);
0759: deriv.addCommand(TIMES);
0760: }
0761: if (prog[op2] < 0x3FFFFFFF
0762: || command[prog[op2] - 0x3FFFFFFF] != wrt) {
0763: compileDerivative(op2, deriv, wrt);
0764: deriv.addCommand(TIMES);
0765: }
0766: } else {
0767: copyExpression(op1, deriv);
0768: copyExpression(op2, deriv);
0769: deriv.addCommand(POWER);
0770: boolean eq = true;
0771: int ext1 = extent(op1);
0772: int ext2 = extent(op2);
0773: if (ext1 != ext2)
0774: eq = false;
0775: else
0776: for (int i = 0; i < ext1; i++)
0777: if (prog[op1 - i] != prog[op2 - i]) {
0778: eq = false;
0779: break;
0780: }
0781: if (eq)
0782: deriv.addConstant(1);
0783: else {
0784: copyExpression(op2, deriv);
0785: copyExpression(op1, deriv);
0786: deriv.addCommand(DIVIDE);
0787: }
0788: if (prog[op1] < 0x3FFFFFFF
0789: || command[prog[op1] - 0x3FFFFFFF] != wrt) {
0790: compileDerivative(op1, deriv, wrt);
0791: deriv.addCommand(TIMES);
0792: }
0793: copyExpression(op1, deriv);
0794: deriv.addCommand(LN);
0795: if (prog[op2] < 0x3FFFFFFF
0796: || command[prog[op2] - 0x3FFFFFFF] != wrt) {
0797: compileDerivative(op2, deriv, wrt);
0798: deriv.addCommand(TIMES);
0799: }
0800: deriv.addCommand(PLUS);
0801: deriv.addCommand(TIMES);
0802: }
0803: break;
0804: }
0805: }
0806:
0807: private void doFuncDeriv(int opCode, int op,
0808: ExpressionProgram deriv, Variable wrt) {
0809: // Add commands to deriv to compute the derivative of a subexpression of this program
0810: // with respect to the variable wrt, where the main operation in the subexpression is
0811: // the standard function represented by opCode, and the position of the operand of the
0812: // standard function is at the index op in the program.
0813: /* if (opCode == TRUNC || opCode == ROUND || opCode == FLOOR || opCode == CEILING) {
0814: // I'm pretending that the derivative is zero, even though it's undefined at
0815: // certain values.
0816: deriv.addConstant(0);
0817: return;
0818: }
0819: */
0820: switch (opCode) {
0821: case FLOOR:
0822: case CEILING:
0823: case TRUNC:
0824: case ROUND:
0825: copyExpression(op, deriv);
0826: if (opCode == ROUND) {
0827: deriv.addConstant(0.5);
0828: deriv.addCommand(PLUS);
0829: }
0830: deriv.addCommand(ROUND);
0831: copyExpression(op, deriv);
0832: if (opCode == ROUND) {
0833: deriv.addConstant(0.5);
0834: deriv.addCommand(PLUS);
0835: }
0836: deriv.addCommand(NE);
0837: if (opCode == TRUNC) {
0838: copyExpression(op, deriv);
0839: deriv.addConstant(0);
0840: deriv.addCommand(EQ);
0841: deriv.addCommand(OR);
0842: }
0843: ExpressionProgram zero = new ExpressionProgram();
0844: zero.addConstant(0);
0845: deriv
0846: .addCommandObject(new ConditionalExpression(zero,
0847: null));
0848: return;
0849: case SIN:
0850: copyExpression(op, deriv);
0851: deriv.addCommand(COS);
0852: break;
0853: case COS:
0854: copyExpression(op, deriv);
0855: deriv.addCommand(SIN);
0856: deriv.addCommand(UNARY_MINUS);
0857: break;
0858: case TAN:
0859: copyExpression(op, deriv);
0860: deriv.addCommand(SEC);
0861: deriv.addConstant(2);
0862: deriv.addCommand(POWER);
0863: break;
0864: case COT:
0865: copyExpression(op, deriv);
0866: deriv.addCommand(CSC);
0867: deriv.addConstant(2);
0868: deriv.addCommand(POWER);
0869: deriv.addCommand(UNARY_MINUS);
0870: break;
0871: case SEC:
0872: copyExpression(op, deriv);
0873: deriv.addCommand(SEC);
0874: copyExpression(op, deriv);
0875: deriv.addCommand(TAN);
0876: deriv.addCommand(TIMES);
0877: break;
0878: case CSC:
0879: copyExpression(op, deriv);
0880: deriv.addCommand(CSC);
0881: copyExpression(op, deriv);
0882: deriv.addCommand(COT);
0883: deriv.addCommand(TIMES);
0884: deriv.addCommand(UNARY_MINUS);
0885: break;
0886: case ARCSIN:
0887: case ARCCOS:
0888: deriv.addConstant(1);
0889: if (opCode == ARCCOS)
0890: deriv.addCommand(UNARY_MINUS);
0891: deriv.addConstant(1);
0892: copyExpression(op, deriv);
0893: deriv.addConstant(2);
0894: deriv.addCommand(POWER);
0895: deriv.addCommand(MINUS);
0896: deriv.addCommand(SQRT);
0897: deriv.addCommand(DIVIDE);
0898: break;
0899: case ARCTAN:
0900: deriv.addConstant(1);
0901: deriv.addConstant(1);
0902: copyExpression(op, deriv);
0903: deriv.addConstant(2);
0904: deriv.addCommand(POWER);
0905: deriv.addCommand(PLUS);
0906: deriv.addCommand(DIVIDE);
0907: break;
0908: case ABS: {
0909: ExpressionProgram pos = new ExpressionProgram();
0910: ExpressionProgram neg = new ExpressionProgram();
0911: compileDerivative(op, pos, wrt);
0912: compileDerivative(op, neg, wrt);
0913: neg.addCommand(UNARY_MINUS);
0914: ExpressionProgram negTest = new ExpressionProgram();
0915: copyExpression(op, negTest);
0916: negTest.addConstant(0);
0917: negTest.addCommand(LT);
0918: negTest.addCommandObject(new ConditionalExpression(neg,
0919: null));
0920: copyExpression(op, deriv);
0921: deriv.addConstant(0);
0922: deriv.addCommand(GT);
0923: deriv.addCommandObject(new ConditionalExpression(pos,
0924: negTest));
0925: }
0926: return;
0927: case SQRT:
0928: deriv.addConstant(1);
0929: deriv.addConstant(2);
0930: copyExpression(op, deriv);
0931: deriv.addCommand(SQRT);
0932: deriv.addCommand(TIMES);
0933: deriv.addCommand(DIVIDE);
0934: break;
0935: case EXP:
0936: copyExpression(op, deriv);
0937: deriv.addCommand(EXP);
0938: break;
0939: case LN:
0940: case LOG2:
0941: case LOG10:
0942: ExpressionProgram d = new ExpressionProgram();
0943: d.addConstant(1);
0944: copyExpression(op, d);
0945: d.addCommand(DIVIDE);
0946: if (opCode != LN) {
0947: d.addConstant((opCode == LOG2) ? 2 : 10);
0948: d.addCommand(LN);
0949: d.addCommand(DIVIDE);
0950: }
0951: copyExpression(op, deriv);
0952: deriv.addConstant(0);
0953: deriv.addCommand(GT);
0954: deriv.addCommandObject(new ConditionalExpression(d, null));
0955: break;
0956: case CUBERT:
0957: deriv.addConstant(1);
0958: deriv.addConstant(3);
0959: copyExpression(op, deriv);
0960: deriv.addConstant(2);
0961: deriv.addCommand(POWER);
0962: deriv.addCommand(CUBERT);
0963: deriv.addCommand(TIMES);
0964: deriv.addCommand(DIVIDE);
0965: break;
0966: }
0967: if (prog[op] < 0x3FFFFFFF
0968: || command[prog[op] - 0x3FFFFFFF] != wrt) {
0969: compileDerivative(op, deriv, wrt);
0970: deriv.addCommand(TIMES);
0971: }
0972: }
0973:
0974: //--------- private stuff for storing constants and ExpressionCommands ----------
0975:
0976: private double[] constant = new double[1]; // Holds all the constants that have been added
0977: // to this ExpressionProgram. In a program, a
0978: // constant is represented by an index into this
0979: // array.
0980:
0981: private int constantCt; // The number of values in the constant array.
0982:
0983: private ExpressionCommand[] command = new ExpressionCommand[1]; // Holds all the ExpressionCommands that have been added
0984: // to this ExpressionProgram. In a program, an
0985: // ExpressionCommand is represented by an index into this
0986: // array, with 0x3FFFFFFF added to the index to give
0987: // a number >= 0x3FFFFFFF.
0988:
0989: private int commandCt; // The number of items in the command array.
0990:
0991: private int findConstant(double d) {
0992: // Find the index of d in the constant array, adding it if it is not already there.
0993: // The array will be expanded if necessary.
0994: for (int i = 0; i < constantCt; i++)
0995: if (constant[i] == d) {
0996: return i;
0997: }
0998: if (constantCt == constant.length) {
0999: double[] temp = new double[constant.length * 2];
1000: System.arraycopy(constant, 0, temp, 0, constantCt);
1001: constant = temp;
1002: }
1003: constant[constantCt++] = d;
1004: return constantCt - 1;
1005: }
1006:
1007: private int findCommand(ExpressionCommand com) {
1008: // Find the index of com in the command array, adding it if it is not already there.
1009: // The array will be expanded if necessary.
1010: for (int i = 0; i < commandCt; i++)
1011: if (command[i] == com)
1012: return i;
1013: if (commandCt == command.length) {
1014: ExpressionCommand[] temp = new ExpressionCommand[command.length * 2];
1015: System.arraycopy(command, 0, temp, 0, commandCt);
1016: command = temp;
1017: }
1018: command[commandCt++] = com;
1019: return commandCt - 1;
1020: }
1021: } // end class ExpressionProgram
|