src/felib/femath.cpp

   1 /*
   2  *
   3  *  Iter Vehemens ad Necem (IVAN)
   4  *  Copyright (C) Timo Kiviluoto
   5  *  Released under the GNU General
   6  *  Public License
   7  *
   8  *  See LICENSING which should be included
   9  *  along with this file for more details
  10  *
  11  */
  12 #include <cmath>
  13
  14 #include "femath.h"
  15 #include "error.h"
  16 #include "fesave.h"
  17
  18
  19 cint basequadricontroller::OrigoDeltaX[4] = { 0, 1, 0, 1 };
  20 cint basequadricontroller::OrigoDeltaY[4] = { 0, 0, 1, 1 };
  21 int basequadricontroller::OrigoX, basequadricontroller::OrigoY;
  22 int basequadricontroller::StartX, basequadricontroller::StartY;
  23 int basequadricontroller::XSize, basequadricontroller::YSize;
  24 int basequadricontroller::RadiusSquare;
  25 truth basequadricontroller::SectorCompletelyClear;
  26
  27
  28 /* A C-program for MT19937: Integer     version                   */
  29 /*  genrand() generates one pseudorandom unsigned integer (32bit) */
  30 /* which is uniformly distributed among 0 to 2^32-1  for each     */
  31 /* call. sgenrand(seed) set initial values to the working area    */
  32 /* of 624 words. Before genrand(), sgenrand(seed) must be         */
  33 /* called once. (seed is any 32-bit integer except for 0).        */
  34 /*   Coded by Takuji Nishimura, considering the suggestions by    */
  35 /* Topher Cooper and Marc Rieffel in July-Aug. 1997.              */
  36
  37 /* This library is free software; you can redistribute it and/or   */
  38 /* modify it under the terms of the GNU Library General Public     */
  39 /* License as published by the Free Software Foundation; either    */
  40 /* version 2 of the License, or (at your option) any later         */
  41 /* version.                                                        */
  42 /* This library is distributed in the hope that it will be useful, */
  43 /* but WITHOUT ANY WARRANTY; without even the implied warranty of  */
  44 /* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.            */
  45 /* See the GNU Library General Public License for more details.    */
  46 /* You should have received a copy of the GNU Library General      */
  47 /* Public License along with this library; if not, write to the    */
  48 /* Free Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA   */
  49 /* 02111-1307  USA                                                 */
  50
  51 /* Copyright (C) 1997 Makoto Matsumoto and Takuji Nishimura.       */
  52 /* Any feedback is very welcome. For any question, comments,       */
  53 /* see http://www.math.keio.ac.jp/matumoto/emt.html or email       */
  54 /* matumoto@math.keio.ac.jp                                        */
  55
  56 /* Period parameters */
  57 #define N1 624
  58 #define M 397
  59 #define MATRIX_A 0x9908b0df   /* constant vector a */
  60 #define UPPER_MASK 0x80000000 /* most significant w-r bits */
  61 #define LOWER_MASK 0x7fffffff /* least significant r bits */
  62
  63 /* Tempering parameters */
  64 #define TEMPERING_MASK_B 0x9d2c5680
  65 #define TEMPERING_MASK_C 0xefc60000
  66 #define TEMPERING_SHIFT_U(y) (y >> 11)
  67 #define TEMPERING_SHIFT_S(y) (y << 7)
  68 #define TEMPERING_SHIFT_T(y) (y << 15)
  69 #define TEMPERING_SHIFT_L(y) (y >> 18)
  70
  71 uLong femath::mt[N1]; /* the array for the state vector  */
  72 sLong femath::mti = N1+1; /* mti==N+1 means mt[N] is not initialized */
  73
  74 /* backups */
  75
  76 uLong femath::mtb[N1];
  77 sLong femath::mtib;
  78
  79 void femath::SetSeed (uLong Seed) {
  80   /* setting initial seeds to mt[N] using         */
  81   /* the generator Line 25 of Table 1 in          */
  82   /* [KNUTH 1981, The Art of Computer Programming */
  83   /*    Vol. 2 (2nd Ed.), pp102]                  */
  84   mt[0] = Seed & 0xffffffff;
  85   for (mti=1; mti < N1; mti++) mt[mti] = (69069*mt[mti-1])&0xffffffff;
  86 }
  87
  88
  89 sLong femath::Rand () {
  90   uLong y;
  91   static uLong mag01[2]={0x0, MATRIX_A};
  92
  93   /* mag01[x] = x * MATRIX_A  for x=0,1 */
  94
  95   if (mti >= N1) { /* generate N words at one time */
  96     int kk;
  97
  98     if (mti == N1+1) /* if sgenrand() has not been called, */
  99       SetSeed(4357); /* a default initial seed is used   */
 100
 101     for (kk=0;kk<N1-M;kk++) {
 102       y = (mt[kk]&UPPER_MASK)|(mt[kk+1]&LOWER_MASK);
 103       mt[kk] = mt[kk+M] ^ (y >> 1) ^ mag01[y & 0x1];
 104     }
 105     for (;kk<N1-1;kk++) {
 106       y = (mt[kk]&UPPER_MASK)|(mt[kk+1]&LOWER_MASK);
 107       mt[kk] = mt[kk+(M-N1)] ^ (y >> 1) ^ mag01[y & 0x1];
 108     }
 109     y = (mt[N1-1]&UPPER_MASK)|(mt[0]&LOWER_MASK);
 110     mt[N1-1] = mt[M-1] ^ (y >> 1) ^ mag01[y & 0x1];
 111
 112     mti = 0;
 113   }
 114
 115   y = mt[mti++];
 116   y ^= TEMPERING_SHIFT_U(y);
 117   y ^= TEMPERING_SHIFT_S(y) & TEMPERING_MASK_B;
 118   y ^= TEMPERING_SHIFT_T(y) & TEMPERING_MASK_C;
 119   y ^= TEMPERING_SHIFT_L(y);
 120
 121   return y & 0x7FFFFFFF;
 122 }
 123
 124
 125 int femath::WeightedRand (sLong *Possibility, sLong TotalPossibility) {
 126   sLong Rand = RAND()%TotalPossibility, PartialSum = 0;
 127   for (int c = 0; ; ++c) {
 128     PartialSum += Possibility[c];
 129     if (PartialSum > Rand) return c;
 130   }
 131 }
 132
 133
 134 int femath::WeightedRand (const std::vector<sLong> &Possibility, sLong TotalPossibility) {
 135   sLong Rand = RAND()%TotalPossibility, PartialSum = 0;
 136   for (int c = 0;; ++c) {
 137     PartialSum += Possibility[c];
 138     if (PartialSum > Rand) return c;
 139   }
 140 }
 141
 142
 143 double femath::CalculateAngle(v2 Direction) {
 144   if (Direction.X < 0) return atan(double(Direction.Y)/Direction.X)+FPI;
 145   if (Direction.X > 0) {
 146     if (Direction.Y < 0) return atan(double(Direction.Y)/Direction.X)+2*FPI;
 147     return atan(double(Direction.Y)/Direction.X);
 148   }
 149   if (Direction.Y < 0) return 3*FPI/2;
 150   if (Direction.Y > 0) return FPI/2;
 151   ABORT("Illegal direction (0, 0) passed to femath::CalculateAngle()!");
 152   return 0;
 153 }
 154
 155
 156 void femath::CalculateEnvironmentRectangle (rect &Rect, const rect &MotherRect, v2 Origo, int Radius) {
 157   Rect.X1 = Origo.X-Radius;
 158   Rect.Y1 = Origo.Y-Radius;
 159   Rect.X2 = Origo.X+Radius;
 160   Rect.Y2 = Origo.Y+Radius;
 161   if (Rect.X1 < MotherRect.X1) Rect.X1 = MotherRect.X1;
 162   if (Rect.Y1 < MotherRect.Y1) Rect.Y1 = MotherRect.Y1;
 163   if (Rect.X2 > MotherRect.X2) Rect.X2 = MotherRect.X2;
 164   if (Rect.Y2 > MotherRect.Y2) Rect.Y2 = MotherRect.Y2;
 165 }
 166
 167
 168 truth femath::Clip (int &SourceX, int &SourceY, int &DestX, int &DestY, int &Width, int &Height, int XSize, int YSize, int DestXSize, int DestYSize) {
 169   /* This sentence is usually true */
 170   if (SourceX >= 0 && SourceY >= 0 && DestX >= 0 && DestY >= 0 &&
 171       SourceX+Width <= XSize && SourceY+Height <= YSize &&
 172       DestX+Width <= DestXSize && DestY+Height <= DestYSize) return true;
 173   if (SourceX < 0) {
 174     Width += SourceX;
 175     DestX -= SourceX;
 176     SourceX = 0;
 177   }
 178   if (SourceY < 0) {
 179     Height += SourceY;
 180     DestY -= SourceY;
 181     SourceY = 0;
 182   }
 183   if (DestX < 0) {
 184     Width += DestX;
 185     SourceX -= DestX;
 186     DestX = 0;
 187   }
 188   if (DestY < 0) {
 189     Height += DestY;
 190     SourceY -= DestY;
 191     DestY = 0;
 192   }
 193   if (SourceX+Width > XSize) Width = XSize-SourceX;
 194   if (SourceY+Height > YSize) Height = YSize-SourceY;
 195   if (DestX+Width > DestXSize) Width = DestXSize-DestX;
 196   if (DestY+Height > DestYSize) Height = DestYSize-DestY;
 197   return Width>0 && Height>0;
 198 }
 199
 200
 201 void femath::SaveSeed () {
 202   mtib = mti;
 203   for (int c = 0; c < N1; ++c) mtb[c] = mt[c];
 204 }
 205
 206
 207 void femath::LoadSeed () {
 208   mti = mtib;
 209   for (int c = 0; c < N1; ++c) mt[c] = mtb[c];
 210 }
 211
 212
 213 void ReadData(interval &I, inputfile &SaveFile) {
 214   I.Min = SaveFile.ReadNumber(HIGHEST, true);
 215   festring Word;
 216   SaveFile.ReadWord(Word);
 217   if (Word == ";" || Word == ",") I.Max = I.Min;
 218   else if (Word == ":") I.Max = Max(SaveFile.ReadNumber(), I.Min);
 219   else ABORT("Odd interval terminator %s detected, file %s line %d!", Word.CStr(), SaveFile.GetFileName().CStr(), SaveFile.TellLine());
 220 }
 221
 222
 223 void ReadData (region &R, inputfile &SaveFile) {
 224   ReadData(R.X, SaveFile);
 225   ReadData(R.Y, SaveFile);
 226 }
 227
 228
 229 outputfile &operator << (outputfile &SaveFile, const interval &I) {
 230   SaveFile.Write(reinterpret_cast<cchar *>(&I), sizeof(I));
 231   return SaveFile;
 232 }
 233
 234
 235 inputfile &operator >> (inputfile &SaveFile, interval &I) {
 236   SaveFile.Read(reinterpret_cast<char *>(&I), sizeof(I));
 237   return SaveFile;
 238 }
 239
 240
 241 outputfile &operator << (outputfile &SaveFile, const region &R) {
 242   SaveFile.Write(reinterpret_cast<cchar *>(&R), sizeof(R));
 243   return SaveFile;
 244 }
 245
 246
 247 inputfile &operator >> (inputfile &SaveFile, region &R) {
 248   SaveFile.Read(reinterpret_cast<char *>(&R), sizeof(R));
 249   return SaveFile;
 250 }
 251
 252
 253 sLong femath::SumArray (const fearray<sLong> &Vector) {
 254   sLong Sum = 0;
 255   for (uInt c = 0; c < Vector.Size; ++c) Sum += Vector.Data[c];
 256   return Sum;
 257 }
 258
 259
 260 void femath::GenerateFractalMap (int **Map, int Side, int StartStep, int Randomness) {
 261   cint Limit = Side-1;
 262   Map[0][0] = 0;
 263   Map[0][Limit] = 0;
 264   Map[Limit][0] = 0;
 265   Map[Limit][Limit] = 0;
 266   for (int Step = StartStep, HalfStep = Step>>1; HalfStep;
 267     Step = HalfStep, HalfStep>>=1, Randomness = ((Randomness<<3)-Randomness)>>3)
 268   {
 269     int x, y, RandMod = (Randomness<<1)+1;
 270
 271     for (x = HalfStep; x < Side; x += Step)
 272       for (y = HalfStep; y < Side; y += Step)
 273         Map[x][y] =
 274           ((Map[x-HalfStep][y-HalfStep]+
 275             Map[x-HalfStep][y+HalfStep]+
 276             Map[x+HalfStep][y-HalfStep]+
 277             Map[x+HalfStep][y+HalfStep])>>2)-Randomness+RAND()%RandMod;
 278
 279     for (x = HalfStep; x < Side; x += Step) {
 280       for (y = 0; y < Side; y += Step) {
 281         int HeightSum = Map[x-HalfStep][y]+Map[x+HalfStep][y];
 282         int Neighbours = 2;
 283         if (y) {
 284           HeightSum += Map[x][y-HalfStep];
 285           ++Neighbours;
 286         }
 287         if (y != Limit) {
 288           HeightSum += Map[x][y+HalfStep];
 289           ++Neighbours;
 290         }
 291         if (Neighbours == 4) HeightSum >>= 2; else HeightSum /= Neighbours;
 292         Map[x][y] = HeightSum-Randomness+RAND()%RandMod;
 293       }
 294     }
 295
 296     for (x = 0; x < Side; x += Step) {
 297       for (y = HalfStep; y < Side; y += Step) {
 298         int HeightSum = Map[x][y-HalfStep]+Map[x][y+HalfStep];
 299         int Neighbours = 2;
 300         if (x) {
 301           HeightSum += Map[x-HalfStep][y];
 302           ++Neighbours;
 303         }
 304         if (x != Limit) {
 305           HeightSum += Map[x+HalfStep][y];
 306           ++Neighbours;
 307         }
 308         if (Neighbours == 4) HeightSum >>= 2; else HeightSum /= Neighbours;
 309         Map[x][y] = HeightSum-Randomness+RAND()%RandMod;
 310       }
 311     }
 312   }
 313 }