1 package MathUtil; 2 # 3 # $RCSfile: MathUtil.pm,v $ 4 # $Date: 2015/02/28 20:47:17 $ 5 # $Revision: 1.28 $ 6 # 7 # Author: Manish Sud <msud@san.rr.com> 8 # 9 # Copyright (C) 2015 Manish Sud. All rights reserved. 10 # 11 # This file is part of MayaChemTools. 12 # 13 # MayaChemTools is free software; you can redistribute it and/or modify it under 14 # the terms of the GNU Lesser General Public License as published by the Free 15 # Software Foundation; either version 3 of the License, or (at your option) any 16 # later version. 17 # 18 # MayaChemTools is distributed in the hope that it will be useful, but without 19 # any warranty; without even the implied warranty of merchantability of fitness 20 # for a particular purpose. See the GNU Lesser General Public License for more 21 # details. 22 # 23 # You should have received a copy of the GNU Lesser General Public License 24 # along with MayaChemTools; if not, see <http://www.gnu.org/licenses/> or 25 # write to the Free Software Foundation Inc., 59 Temple Place, Suite 330, 26 # Boston, MA, 02111-1307, USA. 27 # 28 29 use strict; 30 use Exporter; 31 use Constants; 32 use Math::Trig (); 33 use POSIX (); 34 35 use vars qw(@ISA @EXPORT @EXPORT_OK %EXPORT_TAGS); 36 37 @ISA = qw(Exporter); 38 @EXPORT = qw(acos asin atan tan ceil floor log10 min max srandom random round GeneratePrimeNumbersUpToLimit GeneratePrimeNumbersUpToCount); 39 @EXPORT_OK = qw(); 40 41 %EXPORT_TAGS = (all => [@EXPORT, @EXPORT_OK] 42 ); 43 44 45 # Return next largest integer... 46 sub ceil ($) { 47 my($Value) = @_; 48 49 return POSIX::ceil($Value); 50 } 51 52 # Return previous smallest integer... 53 sub floor ($) { 54 my($Value) = @_; 55 56 return POSIX::floor($Value); 57 } 58 59 # Calculate log value using base 10... 60 sub log10 ($) { 61 my($Value) = @_; 62 63 return CORE::log($Value)/CORE::log(10); 64 } 65 66 # Return the smaller of two numbers... 67 sub min ($$) { 68 my($Value1, $Value2) = @_; 69 70 return ($Value1 <= $Value2) ? $Value1 : $Value2; 71 } 72 73 # Return the larger of two numbers... 74 sub max ($$) { 75 my($Value1, $Value2) = @_; 76 77 return ($Value1 >= $Value2) ? $Value1 : $Value2; 78 } 79 80 # The random number generator implemented in MayaChemTools is a variant of linear 81 # congruential generator (LCG) as described by Miller et al. [ Ref 120 ]. It is 82 # also referred to as Lehmer random number generator or Park-Miller random number 83 # generator. 84 # 85 # Unlike Perl's core random number generator function rand, the random number 86 # generator implemented in MayaChemTools generates consistent random values 87 # across different platforms - Windows, CygWin, Linux, Unix - for a specific random 88 # seed. 89 # 90 91 # $RandomModulus = 2**31 - 1; 92 # $RandomMultiplier = 16807; 93 # $RandomQuotient = $RandomModulus / $RandomMultiplier; 94 # $RandomRemainder = $RandomModulus % $RandomMultiplier 95 # 96 # $MaxRandomSeed = 2*31 -2 97 # 98 my($MaxRandomSeed, $RandomSeed, $RandomModulus, $RandomMultiplier, $RandomQuotient, $RandomRemainder); 99 100 $MaxRandomSeed = 2147483646; 101 $RandomSeed = 123456789; 102 103 $RandomModulus = 2147483647; 104 $RandomMultiplier = 16807; 105 $RandomQuotient = 127773; 106 $RandomRemainder = 2836; 107 108 # Set random number seed... 109 # 110 # The intial value of random number seed is recommeded to be an integer between 1 111 # and 2**31 - 2 [Ref 120] which translates to be 1 and 2147483646 112 # 113 sub srandom ($) { 114 my($Seed) = @_; 115 116 if ($Seed <= 0 ) { 117 die "Error: srandom: Specified seed value must be greater than 0..."; 118 } 119 120 $RandomSeed = ($Seed > $MaxRandomSeed) ? ($Seed % $MaxRandomSeed) : $Seed; 121 122 return $RandomSeed; 123 } 124 125 # Retrun a random number between 0 and less than 1 or specified size... 126 # 127 sub random (;$) { 128 my($Size) = @_; 129 my($Value, $LowValue, $HighValue); 130 131 $Size = defined $Size ? $Size : 1.0; 132 133 $HighValue = $RandomSeed / $RandomQuotient; 134 $LowValue = $RandomSeed % $RandomQuotient; 135 136 $Value = $RandomMultiplier * $LowValue - $RandomRemainder * $HighValue; 137 138 $RandomSeed = ($Value > 0) ? $Value : ($Value + $RandomModulus); 139 140 return ($RandomSeed / $RandomModulus) * $Size; 141 } 142 143 # Round a integer/real number to: 144 # . A nearest integer 145 # . Specified number of decimal places 146 # 147 sub round ($;$) { 148 my($Value, $DecimalPlaces) = @_; 149 my($RoundedValue); 150 151 if (defined($DecimalPlaces) && $DecimalPlaces > 0) { 152 $RoundedValue = sprintf "%.${DecimalPlaces}f", $Value; 153 } 154 else { 155 if ($Value < 0) { 156 $RoundedValue = int($Value - 0.5); 157 } 158 else { 159 $RoundedValue = int($Value + 0.5); 160 } 161 } 162 return $RoundedValue; 163 } 164 165 # Return tangent of an angle expressed in radians. 166 sub tan { 167 my($Value) = @_; 168 169 return (CORE::sin($Value)/CORE::cos($Value)); 170 } 171 172 # Return inverse sine of an angle expressed in radians. 173 # 174 # For a right angle triangle defined by sides X and Y in a unit circle, Pythagorean theorem implies 175 # X**2 + Y**2 = 1 and sin value corresponds to Y. So asin is equivalent to atan2(Y, sqrt(1-Y**2)). 176 # However, taking sqrt of negative numbers is problematic; Math::Trig::asin handles it using complex 177 # numbers. 178 # 179 sub asin ($) { 180 my($Value) = @_; 181 182 return Math::Trig::asin($Value); 183 } 184 185 # Return inverse cosine of an angle expressed in radians. 186 # 187 # For a right angle triangle defined by sides X and Y in a unit circle, Pythagorean theorem implies 188 # X**2 + Y**2 = 1 and cos value corresponds to X. So asin is equivalent to atan2(sqrt(1-X**2), X) 189 # However, taking sqrt of negative numbers is problematic; Math::Trig::acos handles it using complex 190 # numbers. 191 # 192 sub acos ($) { 193 my($Value) = @_; 194 195 return Math::Trig::acos($Value); 196 } 197 198 # Generate prime numbers up to a specified limit and return a reference to an 199 # array containing the prime numbers. 200 # 201 # By default, the first 1000 prime numbers are generated. The 1000th prime 202 # number is 7919 and that's why default limit is set to 7920. 203 # 204 sub GeneratePrimeNumbersUpToLimit (;$) { 205 my($Limit) = @_; 206 207 $Limit = defined $Limit ? $Limit : 7920; 208 209 return _GeneratePrimeNumbers('ByLimit', $Limit) 210 } 211 212 # Generate prime numbers up to specified count of prime numbers and return a 213 # reference to an array containing the prime numbers. 214 # 215 # By default, the first 1000 prime numbers are generated. The 1000th prime 216 # number is 7919. 217 # 218 sub GeneratePrimeNumbersUpToCount (;$) { 219 my($Count) = @_; 220 221 $Count = defined $Count ? $Count : 1000; 222 223 return _GeneratePrimeNumbers('ByCount', $Count) 224 } 225 226 # Generate prime numbers up to specified limit or count and return a reference 227 # to an array containing the prime numbers. 228 # 229 # The algorithm to generate prime numbers is a modification of Sieve of Erastothenes 230 # prime number generator. 231 # 232 sub _GeneratePrimeNumbers { 233 my($Mode, $Value) = @_; 234 my($ByLimit, $PrimeNumber, $Number, $SqrtOfNumber, $NumberIsPrime, @PrimeNumbers); 235 236 $ByLimit = ($Mode =~ /^ByLimit$/i) ? 1 : 0; 237 238 @PrimeNumbers = (2, 3); 239 $Number = 3; 240 241 # while ($Number <= $Limit) { 242 while ($ByLimit ? ($Number < $Value) : (@PrimeNumbers < $Value)) { 243 $Number += 2; 244 $SqrtOfNumber = sqrt $Number; 245 246 $NumberIsPrime = 1; 247 PRIMENUMBER: for $PrimeNumber (@PrimeNumbers) { 248 if ($PrimeNumber > $SqrtOfNumber) { 249 last PRIMENUMBER; 250 } 251 if (!($Number % $PrimeNumber)) { 252 $NumberIsPrime = 0; 253 last PRIMENUMBER; 254 } 255 } 256 if ($NumberIsPrime) { 257 push @PrimeNumbers, $Number; 258 } 259 } 260 return \@PrimeNumbers; 261 } 262