       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