1 // ****************************************************************************** 2 // 3 // Title: Force Field X. 4 // Description: Force Field X - Software for Molecular Biophysics. 5 // Copyright: Copyright (c) Michael J. Schnieders 2001-2024. 6 // 7 // This file is part of Force Field X. 8 // 9 // Force Field X is free software; you can redistribute it and/or modify it 10 // under the terms of the GNU General Public License version 3 as published by 11 // the Free Software Foundation. 12 // 13 // Force Field X is distributed in the hope that it will be useful, but WITHOUT 14 // ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS 15 // FOR A PARTICULAR PURPOSE. See the GNU General Public License for more 16 // details. 17 // 18 // You should have received a copy of the GNU General Public License along with 19 // Force Field X; if not, write to the Free Software Foundation, Inc., 59 Temple 20 // Place, Suite 330, Boston, MA 02111-1307 USA 21 // 22 // Linking this library statically or dynamically with other modules is making a 23 // combined work based on this library. Thus, the terms and conditions of the 24 // GNU General Public License cover the whole combination. 25 // 26 // As a special exception, the copyright holders of this library give you 27 // permission to link this library with independent modules to produce an 28 // executable, regardless of the license terms of these independent modules, and 29 // to copy and distribute the resulting executable under terms of your choice, 30 // provided that you also meet, for each linked independent module, the terms 31 // and conditions of the license of that module. An independent module is a 32 // module which is not derived from or based on this library. If you modify this 33 // library, you may extend this exception to your version of the library, but 34 // you are not obligated to do so. If you do not wish to do so, delete this 35 // exception statement from your version. 36 // 37 // ****************************************************************************** 38 package ffx.numerics.math; 39 40 import static org.apache.commons.math3.util.FastMath.PI; 41 import static org.apache.commons.math3.util.FastMath.exp; 42 43 /** 44 * The ScalarMath class is a simple math library that operates on single variables 45 * 46 * @author Jacob M. Litman 47 * @author Michael J. Schnieders 48 * @since 1.0 49 */ 50 public class ScalarMath { 51 52 private static final double eightPi2 = 8.0 * PI * PI; 53 54 /** 55 * b2u 56 * 57 * @param b a double. 58 * @return a double. 59 */ 60 public static double b2u(double b) { 61 return b / eightPi2; 62 } 63 64 /** 65 * binomial 66 * 67 * @param n a long. 68 * @param k a long. 69 * @return Returns the binomial of (n,k). 70 */ 71 public static long binomial(long n, long k) { 72 return factorial(n) / (factorial(n - k) * factorial(k)); 73 } 74 75 /** 76 * Returns n!! Precondition: n .GE. -1 Returning 1 for -1 input is analogous to Maple behavior. 77 * 78 * @param n long. 79 * @return Returns the n!!. 80 */ 81 public static long doubleFactorial(long n) { 82 if (n < -1) { 83 throw new RuntimeException("Underflow error in doubleFactorial"); 84 } else if (n == -1 || n == 0 || n == 1) { 85 return 1; 86 } else { 87 return n * doubleFactorial(n - 2); 88 } 89 } 90 91 /** 92 * Returns n! <br> Precondition: n .GE. 0 and n .LE. 20 <br> Max long = 9223372036854775807 <br> 93 * 20! = 2432902008176640000 is ok. <br> 21! returns an overflow: -4249290049419214848 94 * 95 * @param n long. 96 * @return Returns n!. 97 */ 98 public static long factorial(long n) { 99 if (n < 0) { 100 throw new RuntimeException("Underflow error in factorial"); 101 } else if (n > 20) { 102 throw new RuntimeException("Overflow error in factorial"); 103 } else if (n == 0) { 104 return 1; 105 } else { 106 return n * factorial(n - 1); 107 } 108 } 109 110 /** 111 * Compute 1.0 / (1.0 + exp(x)). 112 * 113 * @param x Input. 114 * @return Returns 1.0 / (1.0 + exp(x)). 115 */ 116 public static double fermiFunction(double x) { 117 return 1.0 / (1.0 + exp(x)); 118 } 119 120 /** 121 * This is an atypical mod function used by crystallography methods. 122 * 123 * <p>mod 124 * 125 * @param a Value to mod. 126 * @param b Value to mod by. 127 * @return Positive a % b. 128 */ 129 public static double mod(double a, double b) { 130 var res = a % b; 131 if (res < 0.0) { 132 res += b; 133 } 134 return res; 135 } 136 137 /** 138 * This is an atypical mod function used by crystallography methods. 139 * 140 * @param a an int. 141 * @param b an int. 142 * @return an int. 143 */ 144 public static int mod(int a, int b) { 145 int res = a % b; 146 if (res < 0) { 147 res += b; 148 } 149 return res; 150 } 151 152 /** 153 * Atypical mod function used to move a value into the range lb <= value < ub, assuming the 154 * domain is periodic with a period of (ub - lb). 155 * 156 * @param value Value to move between bounds. 157 * @param lb Lower bound. 158 * @param ub Upper bound. 159 * @return Returns periodic copy of value, in the range lb <= value < ub. 160 */ 161 public static double modToRange(double value, double lb, double ub) { 162 value -= lb; 163 var range = ub - lb; 164 value = mod(value, range); 165 value += lb; 166 return value; 167 } 168 169 /** 170 * u2b 171 * 172 * @param u a double. 173 * @return a double. 174 */ 175 public static double u2b(double u) { 176 return u * eightPi2; 177 } 178 179 /** 180 * quadForm 181 * 182 * @param v an array of double. 183 * @param mat an array of double. 184 * @return a double. 185 */ 186 public static double quadForm(double[] v, double[][] mat) { 187 return (v[0] * (v[0] * mat[0][0] + 2 * (v[1] * mat[0][1] + v[2] * mat[0][2])) 188 + v[1] * (v[1] * mat[1][1] + 2 * (v[2] * mat[1][2])) 189 + v[2] * v[2] * mat[2][2]); 190 } 191 192 /** 193 * Reflect proposed angles to be within 0 and 180 degrees. 194 * 195 * @param angle Proposed angle in radians. 196 * @return value of mirrored angle in radians. 197 */ 198 public static double mirrorRadians(double angle) { 199 double angleDegrees = angle * 180 / PI; 200 return mirrorDegrees(angleDegrees) * PI / 180; 201 } 202 203 /** 204 * Reflect proposed angles to be within 0 and 180 degrees. 205 * 206 * @param angle Proposed angle in radians. 207 * @return value of mirrored angle in radians. 208 */ 209 public static double mirrorDegrees(double angle) { 210 if (angle > 180.0) { 211 angle = 180.0 - (angle - 180.0); 212 } else if (angle < 0.0) { 213 angle = 0.0 - angle; 214 } 215 return angle; 216 } 217 }