1 // ****************************************************************************** 2 // 3 // Title: Force Field X. 4 // Description: Force Field X - Software for Molecular Biophysics. 5 // Copyright: Copyright (c) Michael J. Schnieders 2001-2025. 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 private ScalarMath() { 55 // Prevent instantiation. 56 } 57 58 /** 59 * b2u 60 * 61 * @param b a double. 62 * @return a double. 63 */ 64 public static double b2u(double b) { 65 return b / eightPi2; 66 } 67 68 /** 69 * binomial 70 * 71 * @param n a long. 72 * @param k a long. 73 * @return Returns the binomial of (n,k). 74 */ 75 public static long binomial(long n, long k) { 76 return factorial(n) / (factorial(n - k) * factorial(k)); 77 } 78 79 /** 80 * Returns n!! Precondition: n .GE. -1 Returning 1 for -1 input is analogous to Maple behavior. 81 * 82 * @param n long. 83 * @return Returns the n!!. 84 */ 85 public static long doubleFactorial(long n) { 86 if (n < -1) { 87 throw new RuntimeException("Underflow error in doubleFactorial"); 88 } else if (n == -1 || n == 0 || n == 1) { 89 return 1; 90 } else { 91 return n * doubleFactorial(n - 2); 92 } 93 } 94 95 /** 96 * Returns n! <br> Precondition: n .GE. 0 and n .LE. 20 <br> Max long = 9223372036854775807 <br> 97 * 20! = 2432902008176640000 is ok. <br> 21! returns an overflow: -4249290049419214848 98 * 99 * @param n long. 100 * @return Returns n!. 101 */ 102 public static long factorial(long n) { 103 if (n < 0) { 104 throw new RuntimeException("Underflow error in factorial"); 105 } else if (n > 20) { 106 throw new RuntimeException("Overflow error in factorial"); 107 } else if (n == 0) { 108 return 1; 109 } else { 110 return n * factorial(n - 1); 111 } 112 } 113 114 /** 115 * Compute 1.0 / (1.0 + exp(x)). 116 * 117 * @param x Input. 118 * @return Returns 1.0 / (1.0 + exp(x)). 119 */ 120 public static double fermiFunction(double x) { 121 return 1.0 / (1.0 + exp(x)); 122 } 123 124 /** 125 * This is an atypical mod function used by crystallography methods. 126 * 127 * <p>mod 128 * 129 * @param a Value to mod. 130 * @param b Value to mod by. 131 * @return Positive a % b. 132 */ 133 public static double mod(double a, double b) { 134 var res = a % b; 135 if (res < 0.0) { 136 res += b; 137 } 138 return res; 139 } 140 141 /** 142 * This is an atypical mod function used by crystallography methods. 143 * 144 * @param a an int. 145 * @param b an int. 146 * @return an int. 147 */ 148 public static int mod(int a, int b) { 149 int res = a % b; 150 if (res < 0) { 151 res += b; 152 } 153 return res; 154 } 155 156 /** 157 * Atypical mod function used to move a value into the range lb <= value < ub, assuming the 158 * domain is periodic with a period of (ub - lb). 159 * 160 * @param value Value to move between bounds. 161 * @param lb Lower bound. 162 * @param ub Upper bound. 163 * @return Returns periodic copy of value, in the range lb <= value < ub. 164 */ 165 public static double modToRange(double value, double lb, double ub) { 166 value -= lb; 167 var range = ub - lb; 168 value = mod(value, range); 169 value += lb; 170 return value; 171 } 172 173 /** 174 * u2b 175 * 176 * @param u a double. 177 * @return a double. 178 */ 179 public static double u2b(double u) { 180 return u * eightPi2; 181 } 182 183 /** 184 * quadForm 185 * 186 * @param v an array of double. 187 * @param mat an array of double. 188 * @return a double. 189 */ 190 public static double quadForm(double[] v, double[][] mat) { 191 return (v[0] * (v[0] * mat[0][0] + 2 * (v[1] * mat[0][1] + v[2] * mat[0][2])) 192 + v[1] * (v[1] * mat[1][1] + 2 * (v[2] * mat[1][2])) 193 + v[2] * v[2] * mat[2][2]); 194 } 195 196 /** 197 * Reflect proposed angles to be within 0 and 180 degrees. 198 * 199 * @param angle Proposed angle in radians. 200 * @return value of mirrored angle in radians. 201 */ 202 public static double mirrorRadians(double angle) { 203 double angleDegrees = angle * 180 / PI; 204 return mirrorDegrees(angleDegrees) * PI / 180; 205 } 206 207 /** 208 * Reflect proposed angles to be within 0 and 180 degrees. 209 * 210 * @param angle Proposed angle in radians. 211 * @return value of mirrored angle in radians. 212 */ 213 public static double mirrorDegrees(double angle) { 214 if (angle > 180.0) { 215 angle = 180.0 - (angle - 180.0); 216 } else if (angle < 0.0) { 217 angle = 0.0 - angle; 218 } 219 return angle; 220 } 221 }