/* This file is part of the source code for 3D-XplorMath-J, Version 1.0 (January 2008). * Copyright (c) 2008 The 3D-XplorMath Consortium (http://3d-xplormath.org). * This source code is released under a BSD License, which allows redistribution * in source and binary form, with or without modification, provided copyright * and license information are included, and with no warranty or guarantee of * any kind. For details, see http://3d-xplormath.org/j/source/BSDLicense.txt */ package vmm.functions; import java.util.HashMap; import vmm.core.Complex; /** * An object of type EvalStack is used in the evaluation of an expression or function * (as defined, for example, in the classes {@link Expression}, {@link Function}, * and {@link ComplexFunction}). The value() methods in such expressions * and fuctions can, optionally, take a parameter of type EvalStack. In fact, the * only thing that you can do with an EvalStack, outside of the package vmm.functions, * is construct it and pass it as a parameter to one of these value() * methods; the only reason to do this is to evaluate functions and expressions in * a thread-safe way. */ public final class EvalStack { private double[] stack; private int top; private int addressBase; // The base of the "activation record" that holds the // arguments of a ProgFunction that is being evaluaed. // This value changes in startFunctionCall(), endRealValuedFunction() // and endComplexValuedFunction(). private static HashMap perThreadStacks = new HashMap(); /** * Create an EvalStack with an initial size of 30. */ public EvalStack() { this(30); } /** * Create an EvalStack with a specified initial size. The size will grow, if necessary, * so it is generally sufficient to use the default constructor. * @param initialSize the initial size of the stack. If a size less than 1 is specified, * then the initial size will be 1. */ public EvalStack(int initialSize) { if (initialSize < 1) initialSize = 1; stack = new double[initialSize]; } /** * Returns an EvalStack that is unique to the current thread. This * method simply calls {@link #perThread(Thread)}, with Thread.currentThread() * as its parameter. */ final public static EvalStack perThread() { return perThread(Thread.currentThread()); } /** * The first time this method is called by a given thread, it creates and returns * a new EvalStack; subsequent calls with return the same EvalStack rather than create * a new one. This meant to support thread-safe evalution of functions; this capability * is used in evaluation methods such as {@link Expression#value()}, {@link Function#value(double[])} * and {@link ComplexFunction#value(Complex[])}. The method {@link #perThreadRelease(Thread)} * can be called to discard the EvalStack that has been created for a given thread. * @param thread the Thread whose associated EvalStack is to be returned. If thread is null, * the return value is also null. */ final public static EvalStack perThread(Thread thread) { EvalStack stack = perThreadStacks.get(thread); if (stack == null && thread != null) { stack = new EvalStack(); perThreadStacks.put(thread,stack); } return stack; } /** * Discard the EvalStack, if any, that has been created for the specified thread by * {@link #perThread()} or {@link #perThread(Thread)}. * @param thread The thread whose EvalStack is to be discarded. If the value is * null or if there is not EvalStack for the specified stack, nothing is done; * this is not an error. */ public static void perThreadRelease(Thread thread) { perThreadStacks.remove(thread); } void reset() { top = addressBase = 0; } boolean isEmpty() { return top == 0; } void startFunctionCall(int sizeOfArgBlock) { // Called before applying a ProgFunction push(addressBase); addressBase = top-1-sizeOfArgBlock; } void endRealValuedFunction(int sizeOfArgBlock) { // Called after applying a ProgFunction double x = pop(); addressBase = (int)pop(); top -= sizeOfArgBlock; push(x); } void endComplexValuedFunction(int sizeOfArgBlock) { // Called afer applying a ProgFunction double a = pop(); double b = pop(); addressBase = (int)pop(); top -= sizeOfArgBlock; push(b); push(a); } void push(double x) { if (stack.length == top) { double[] newStack = new double[top*2]; System.arraycopy(stack, 0, newStack, 0, top); stack = newStack; } stack[top++] = x; } void push(double re, double im) { push(re); push(im); } void push(Complex c) { push(c.re); push(c.im); } double pop() { return stack[--top]; } Complex popComplex() { double im = pop(); double re = pop(); return new Complex(re,im); } void popComplex(Complex c) { c.im = pop(); c.re = pop(); } void fetch(int address) { // Used to copy a real-valued function argument from its position lower on the stack. push(stack[address+addressBase]); } void fetchComplex(int address) { // Used to copy a complex-valued function argument from its position lower on the stack. push(stack[address+addressBase]); push(stack[address+addressBase+1]); } void apply(StackOp op) { double x, y; double xr,xi, yr,yi; switch(op) { case COMPLEX_TO_REAL: pop(); break; case REAL_TO_COMPLEX: push(0); break; case FIRST_OP_TO_COMPLEX: x = pop(); y = pop(); push(0); push(y); push(x); break; case IMAGINARY_PART: xi = pop(); pop(); push(xi); break; case PLUS: push(pop()+pop()); break; case MINUS: y = pop(); push(pop() - y); break; case TIMES: push(pop()*pop()); break; case DIVIDE: y = pop(); push(pop() / y); break; case POWER: y = pop(); push(Math.pow(pop(), y)); break; case C_PLUS: yi = pop(); yr = pop(); xi = pop(); xr = pop(); push(xr+yr, xi+yi); break; case C_MINUS: yi = pop(); yr = pop(); xi = pop(); xr = pop(); push(xr-yr, xi-yi); break; case C_TIMES: yi = pop(); yr = pop(); xi = pop(); xr = pop(); push(xr*yr - xi*yi, xr*yi + yr*xi); break; case C_DIVIDE: { yi = pop(); yr = pop(); xi = pop(); xr = pop(); double denom = yr*yr + yi*yi; push( (xr*yr+xi*yi)/denom, (xi*yr-xr*yi)/denom ); } break; case C_POWER: { yi = pop(); yr = pop(); xi = pop(); xr = pop(); double modulus = Math.sqrt(xr*xr+xi*xi); double modulus_log = Math.log(modulus); double arg = Math.atan2(xi,xr); double modulus_ans = Math.exp(yr*modulus_log - yi*arg); double arg_ans = yi*modulus_log + yr*arg; push(modulus_ans*Math.cos(arg_ans), modulus_ans*Math.sin(arg_ans)); } break; case C_REAL_POWER: { x = pop(); xi = pop(); xr = pop(); double modulus = Math.sqrt(xi*xi + xr*xr); double arg = Math.atan2(xi,xr); double log_re = Math.log(modulus); double log_im = arg; double x_log_re = x * log_re; double x_log_im = x * log_im; double modulus_ans = Math.exp(x_log_re); push(modulus_ans*Math.cos(x_log_im), modulus_ans*Math.sin(x_log_im)); } break; case EQ: push( (pop() == pop())? 1 : 0 ); break; case NE: push( (pop() != pop())? 1 : 0 ); break; case GE: push( (pop() <= pop())? 1 : 0 ); break; case GT: push( (pop() < pop())? 1 : 0 ); break; case LE: push( (pop() >= pop())? 1 : 0 ); break; case LT: push( (pop() > pop())? 1 : 0 ); break; case C_EQ: yi = pop(); yr = pop(); xi = pop(); xr = pop(); push ( (xr==yr && xi==yi)? 1 : 0 ); break; case C_NE: yi = pop(); yr = pop(); xi = pop(); xr = pop(); push ( (xr==yr && xi==yi)? 0 : 1 ); break; case AND: y = pop(); x = pop(); push( (x != 0 && y != 0)? 1 : 0); break; case OR: y = pop(); x = pop(); push( (x != 0 || y != 0)? 1 : 0); break; case NOT: push( (pop() != 0)? 0 : 1); break; case UNARY_MINUS: push(-pop()); break; case C_UNARY_MINUS: xi = pop(); xr = pop(); push(-xr,-xi); break; case ABS: push(Math.abs(pop())); break; case SQRT: push(Math.sqrt(pop())); break; case CUBERT: push(Math.cbrt(pop())); break; case EXP: push(Math.exp(pop())); break; case LOG: push(Math.log(pop())); break; case LOG2: push(Math.log(pop())/Math.log(2)); break; case LOG10: push(Math.log10(pop())); break; case C_ABS: xi = pop(); xr = pop(); push(Math.sqrt(xr*xr+xi*xi)); break; case C_SQRT: push(0.5); apply(StackOp.C_REAL_POWER); break; case C_CUBERT: push(1.0/3.0); apply(StackOp.C_REAL_POWER); break; case C_EXP: { xi = pop(); xr = pop(); double e = Math.exp(xr); push( e*Math.cos(xi), e*Math.sin(xi) ); } break; case C_LOG: xi = pop(); xr = pop(); push( 0.5*Math.log(xr*xr+xi*xi), Math.atan2(xi,xr) ); break; case C_LOG2: apply(StackOp.C_LOG); push(Math.log(2),0); apply(StackOp.C_DIVIDE); break; case C_LOG10: apply(StackOp.C_LOG); push(Math.log(10),0); apply(StackOp.C_DIVIDE); break; case ARG: xi = pop(); xr = pop(); push(Math.atan2(xi,xr)); break; case TRUNC: push( (int)pop() ); break; case ROUND: push( (int)Math.round(pop()) ); break; case CEILING: push( Math.ceil(pop()) ); break; case FLOOR: push( Math.floor(pop()) ); break; case SIGNUM: push( Math.signum(pop()) ); break; case SIN: push( Math.sin(pop()) ); break; case COS: push( Math.cos(pop()) ); break; case TAN: push( Math.tan(pop()) ); break; case SEC: push( 1.0 / Math.cos(pop()) ); break; case COT: x = pop(); push( Math.cos(x) / Math.sin(x) ); break; case CSC: push( 1.0 / Math.sin(pop()) ); break; case SINH: push( Math.sinh(pop()) ); break; case COSH: push( Math.cosh(pop()) ); break; case TANH: push( Math.tanh(pop()) ); break; case SECH: push( 1.0 / Math.cosh(pop()) ); break; case COTH: x = pop(); push( Math.cosh(x) / Math.sinh(x) ); break; case CSCH: push( 1.0 / Math.sinh(pop()) ); break; case ARCSIN: push( Math.asin(pop()) ); break; case ARCCOS: push( Math.acos(pop()) ); break; case ARCTAN: push( Math.atan(pop()) ); break; case ARCSINH: x = pop(); push( Math.log(Math.sqrt(x*x+1) + x)); break; case ARCCOSH: x = pop(); push( x < 1 ? Double.NaN : Math.log(Math.sqrt(x*x - 1) + x)); break; case ARCTANH: x = pop(); push( Math.abs(x) >= 1 ? Double.NaN : Math.log( (1+x)/(1-x) )/2 ); break; case C_SIN: { xi = pop(); xr = pop(); double e2 = Math.exp(xi); double e1 = 1.0 / e2; push( (e1+e2)*Math.sin(xr)/2, (e2-e1)*Math.cos(xr)/2 ); } break; case C_COS: { xi = pop(); xr = pop(); double e2 = Math.exp(xi); double e1 = 1.0 / e2; push( (e1+e2)*Math.cos(xr)/2, (e1-e2)*Math.sin(xr)/2 ); } break; case C_TAN: { xi = pop(); xr = pop(); double e2 = Math.exp(xi); double e1 = 1.0 / e2; double a = (e1+e2)*Math.sin(xr)/2; // sin double b = (e2-e1)*Math.cos(xr)/2; double c = (e1+e2)*Math.cos(xr)/2; // cos double d = (e1-e2)*Math.sin(xr)/2; double denom = c*c + d*d; push( (a*c+b*d)/denom, (b*c-a*d)/denom ); } break; case C_SEC: { xi = pop(); xr = pop(); double e2 = Math.exp(xi); double e1 = 1.0 / e2; double a = (e1+e2)*Math.cos(xr)/2; // cos double b = (e1-e2)*Math.sin(xr)/2; double denom = a*a + b*b; push( a/denom, -b/denom ); } break; case C_COT: { xi = pop(); xr = pop(); double e2 = Math.exp(xi); double e1 = 1.0 / e2; double a = (e1+e2)*Math.cos(xr)/2; // cos double b = (e1-e2)*Math.sin(xr)/2; double c = (e1+e2)*Math.sin(xr)/2; // sin double d = (e2-e1)*Math.cos(xr)/2; double denom = c*c + d*d; push( (a*c+b*d)/denom, (b*c-a*d)/denom ); } break; case C_CSC: { xi = pop(); xr = pop(); double e2 = Math.exp(xi); double e1 = 1.0 / e2; double a = (e1+e2)*Math.sin(xr)/2; // sin double b = (e2-e1)*Math.cos(xr)/2; double denom = a*a + b*b; push( a/denom, -b/denom ); } break; case C_SINH: { xi = pop(); xr = pop(); double e2 = Math.exp(xr); double e1 = 1.0 / e2; push( (e2-e1)*Math.cos(xi)/2, (e1+e2)*Math.sin(xi)/2 ); } break; case C_COSH: { xi = pop(); xr = pop(); double e2 = Math.exp(xr); double e1 = 1.0 / e2; push( (e1+e2)*Math.cos(xi)/2, (e2-e1)*Math.sin(xi)/2 ); } break; case C_TANH: { xi = pop(); xr = pop(); double e2 = Math.exp(xr); double e1 = 1.0 / e2; double a = (e2-e1)*Math.cos(xi)/2; // sinh double b = (e1+e2)*Math.sin(xi)/2; double c = (e1+e2)*Math.cos(xi)/2; // cosh double d = (e2-e1)*Math.sin(xi)/2; double denom = c*c + d*d; push( (a*c+b*d)/denom, (b*c-a*d)/denom ); } break; case C_SECH: { xi = pop(); xr = pop(); double e2 = Math.exp(xr); double e1 = 1.0 / e2; double a = (e1+e2)*Math.cos(xi)/2; // cosh double b = (e2-e1)*Math.sin(xi)/2; double denom = a*a + b*b; push( a/denom, -b/denom ); } break; case C_COTH: { xi = pop(); xr = pop(); double e2 = Math.exp(xr); double e1 = 1.0 / e2; double a = (e1+e2)*Math.cos(xi)/2; // cosh double b = (e2-e1)*Math.sin(xi)/2; double c = (e2-e1)*Math.cos(xi)/2; // sinh double d = (e1+e2)*Math.sin(xi)/2; double denom = c*c + d*d; push( (a*c+b*d)/denom, (b*c-a*d)/denom ); } break; case C_CSCH: { xi = pop(); xr = pop(); double e2 = Math.exp(xr); double e1 = 1.0 / e2; double a = (e2-e1)*Math.cos(xi)/2; // sinh double b = (e1+e2)*Math.sin(xi)/2; double denom = a*a + b*b; push( a/denom, -b/denom ); } break; case C_ARCSIN: // -i log[iz + (1-z^2)^0.5] xi = pop(); xr = pop(); push(-xi,xr); // iz push(1,0); push(xr,xi); push(2); apply(StackOp.C_REAL_POWER); apply(StackOp.C_MINUS); push(0.5); apply(StackOp.C_REAL_POWER); apply(StackOp.C_PLUS); apply(StackOp.C_LOG); xi = pop(); xr = pop(); push(xi,-xr); // *(-i) break; case C_ARCCOS: // -i log[z+i(1-z^2)^0.5] xi = pop(); xr = pop(); push(xr,xi); // z push(1,0); push(xr,xi); push(2); apply(StackOp.C_REAL_POWER); apply(StackOp.C_MINUS); push(0.5); apply(StackOp.C_REAL_POWER); xi = pop(); xr = pop(); push(-xi,xr); // *i apply(StackOp.C_PLUS); apply(StackOp.C_LOG); xi = pop(); xr = pop(); push(xi,-xr); // *(-i) break; case C_ARCTAN: // i/2 * log[ (i+z)/(i-z) ] xi = pop(); xr = pop(); push(xr,xi+1); // i+z push(-xr,1-xi); // i-z apply(StackOp.C_DIVIDE); apply(StackOp.C_LOG); xi = pop(); xr = pop(); push(-xi/2,xr/2); break; case C_ARCSINH: // log(z+(z^2+1)^0.5) xi = pop(); xr = pop(); push(xr,xi); push(xr,xi); push(2); apply(StackOp.C_REAL_POWER); push(1,0); apply(StackOp.C_PLUS); push(0.5); apply(StackOp.C_REAL_POWER); apply(StackOp.C_PLUS); apply(StackOp.C_LOG); break; case C_ARCCOSH: // log(a+(x^2-1)^0.5) xi = pop(); xr = pop(); push(xr,xi); push(xr,xi); push(2); apply(StackOp.C_REAL_POWER); push(1,0); apply(StackOp.C_MINUS); push(0.5); apply(StackOp.C_REAL_POWER); apply(StackOp.C_PLUS); apply(StackOp.C_LOG); break; case C_ARCTANH: // 1/2 log[ (1+z)/(1-z) ] xi = pop(); xr = pop(); push(xr+1,xi); push(1-xr,-xi); apply(StackOp.C_DIVIDE); apply(StackOp.C_LOG); xi = pop(); xr = pop(); push(xr/2,xi/2); break; case CONJ: xi = pop(); push(-xi); break; } } }