"Fossies" - the Fresh Open Source Software Archive

Member "mathmod-branches-r508-trunk/pariso/commun.cpp" (8 Mar 2021, 10508 Bytes) of package /linux/misc/mathmod-11.0-source.zip:


As a special service "Fossies" has tried to format the requested source page into HTML format using (guessed) C and C++ source code syntax highlighting (style: standard) with prefixed line numbers and code folding option. Alternatively you can here view or download the uninterpreted source code file. For more information about "commun.cpp" see the Fossies "Dox" file reference documentation and the latest Fossies "Diffs" side-by-side code changes report: 10.1_vs_11.0.

    1 /***************************************************************************
    2  *   Copyright (C) 2021 by Abderrahman Taha                                *
    3  *                                                                         *
    4  *                                                                         *
    5  *   This program is free software; you can redistribute it and/or modify  *
    6  *   it under the terms of the GNU General Public License as published by  *
    7  *   the Free Software Foundation; either version 2 of the License, or     *
    8  *   (at your option) any later version.                                   *
    9  *                                                                         *
   10  *   This program is distributed in the hope that it will be useful,       *
   11  *   but WITHOUT ANY WARRANTY; without even the implied warranty of        *
   12  *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the         *
   13  *   GNU General Public License for more details.                          *
   14  *                                                                         *
   15  *   You should have received a copy of the GNU General Public License     *
   16  *   along with this program; if not, write to the                         *
   17  *   Free Software Foundation, Inc.,                                       *
   18  *   51 Franklin Street, Fifth Floor,Boston, MA 02110-1301 USA             *
   19  ***************************************************************************/
   20 
   21 #include "commun.h"
   22 
   23 static int p[512];
   24 static int permutation[256] =
   25 {
   26     151, 160, 137, 91,  90,  15,  131, 13,  201, 95,  96,  53,  194, 233, 7,
   27     225, 140, 36,  103, 30,  69,  142, 8,   99,  37,  240, 21,  10,  23,  190,
   28     6,   148, 247, 120, 234, 75,  0,   26,  197, 62,  94,  252, 219, 203, 117,
   29     35,  11,  32,  57,  177, 33,  88,  237, 149, 56,  87,  174, 20,  125, 136,
   30     171, 168, 68,  175, 74,  165, 71,  134, 139, 48,  27,  166, 77,  146, 158,
   31     231, 83,  111, 229, 122, 60,  211, 133, 230, 220, 105, 92,  41,  55,  46,
   32     245, 40,  244, 102, 143, 54,  65,  25,  63,  161, 1,   216, 80,  73,  209,
   33     76,  132, 187, 208, 89,  18,  169, 200, 196, 135, 130, 116, 188, 159, 86,
   34     164, 100, 109, 198, 173, 186, 3,   64,  52,  217, 226, 250, 124, 123, 5,
   35     202, 38,  147, 118, 126, 255, 82,  85,  212, 207, 206, 59,  227, 47,  16,
   36     58,  17,  182, 189, 28,  42,  223, 183, 170, 213, 119, 248, 152, 2,   44,
   37     154, 163, 70,  221, 153, 101, 155, 167, 43,  172, 9,   129, 22,  39,  253,
   38     19,  98,  108, 110, 79,  113, 224, 232, 178, 185, 112, 104, 218, 246, 97,
   39     228, 251, 34,  242, 193, 238, 210, 144, 12,  191, 179, 162, 241, 81,  51,
   40     145, 235, 249, 14,  239, 107, 49,  192, 214, 31,  181, 199, 106, 157, 184,
   41     84,  204, 176, 115, 121, 50,  45,  127, 4,   150, 254, 138, 236, 205, 93,
   42     222, 114, 67,  29,  24,  72,  243, 141, 128, 195, 78,  66,  215, 61,  156,
   43     180
   44 };
   45 float tinyrnd()
   46 {
   47     static unsigned trand = 0;
   48     trand = 1664525u * trand + 1013904223u;
   49     return (float(trand) / 4294967296.0f);
   50 }
   51 float CellNoise::CellNoiseFunc(float x, float y, float z, int seed, int type, int CombineDist)
   52 {
   53     uint lastRandom, numberFeaturePoints;
   54     float randomDiff[4];
   55     int cubeX, cubeY, cubeZ;
   56     float distanceArray[9];
   57     float color = 0;
   58     for (int i = 0; i < 9; i++)
   59         distanceArray[i] = 6666;
   60     int evalCubeX = int(floor(x));
   61     int evalCubeY = int(floor(y));
   62     int evalCubeZ = int(floor(z));
   63 
   64     for (int i = -1; i < 2; ++i)
   65         for (int j = -1; j < 2; ++j)
   66             for (int k = -1; k < 2; ++k)
   67             {
   68                 cubeX = evalCubeX + i;
   69                 cubeY = evalCubeY + j;
   70                 cubeZ = evalCubeZ + k;
   71                 lastRandom = uint(lcgRandom(hash((cubeX + seed), (cubeY), (cubeZ))));
   72                 numberFeaturePoints = uint(probLookup(lastRandom));
   73                 for (uint l = 0; l < numberFeaturePoints; ++l)
   74                 {
   75                     lastRandom = uint(lcgRandom(int(lastRandom)));
   76                     randomDiff[0] = float(lastRandom) / 0x100000000;
   77 
   78                     lastRandom = uint(lcgRandom(int(lastRandom)));
   79                     randomDiff[1] = float(lastRandom) / 0x100000000;
   80 
   81                     lastRandom = uint(lcgRandom(int(lastRandom)));
   82                     randomDiff[2] = float(lastRandom) / 0x100000000;
   83 
   84                     featurePoint[0] = randomDiff[0] + float(cubeX);
   85                     featurePoint[1] = randomDiff[1] + float(cubeY);
   86                     featurePoint[2] = randomDiff[2] + float(cubeZ);
   87                     if (type == 1)
   88                         insert(distanceArray,
   89                                ManhattanDistanceFunc(x, y, z, featurePoint[0],
   90                                                      featurePoint[1], featurePoint[2]));
   91                     else if (type == 2 || type == 4)
   92                         insert(distanceArray,
   93                                EuclidianDistanceFunc(x, y, z, featurePoint[0],
   94                                                      featurePoint[1], featurePoint[2]));
   95                     else if (type == 3 || type == 5)
   96                         insert(distanceArray,
   97                                ChebyshevDistanceFunc(x, y, z, featurePoint[0],
   98                                                      featurePoint[1], featurePoint[2]));
   99                 }
  100             }
  101     if (CombineDist == 1)
  102         color = distanceArray[1] - distanceArray[0];
  103     else if (CombineDist == 2)
  104         color = distanceArray[2] - distanceArray[0];
  105     else
  106         color = distanceArray[0];
  107     if (color < 0)
  108         color = 0;
  109     if (color > 1)
  110         color = 1;
  111     return color;
  112 }
  113 float CellNoise::EuclidianDistanceFunc(float x1, float y1, float z1, float x2, float y2, float z2)
  114 {
  115     return (x1 - x2) * (x1 - x2) + (y1 - y2) * (y1 - y2) + (z1 - z2) * (z1 - z2);
  116 }
  117 float CellNoise::ManhattanDistanceFunc(float x1, float y1, float z1, float x2, float y2, float z2)
  118 {
  119     float tmp = std::abs(x1 - x2) + std::abs(y1 - y2) + std::abs(z1 - z2);
  120     return tmp;
  121 }
  122 float CellNoise::ChebyshevDistanceFunc(float x1, float y1, float z1, float x2, float y2, float z2)
  123 {
  124     float diff[3];
  125     diff[0] = x1 - x2;
  126     diff[1] = y1 - y2;
  127     diff[2] = z1 - z2;
  128     return std::max(std::max(std::abs(diff[0]), std::abs(diff[1])),
  129                     std::abs(diff[2]));
  130 }
  131 int CellNoise::probLookup(uint value)
  132 {
  133     if (value < 393325350U)
  134         return 1;
  135     if (value < 1022645910U)
  136         return 2;
  137     if (value < 1861739990U)
  138         return 3;
  139     if (value < 2700834071U)
  140         return 4;
  141     if (value < 3372109335U)
  142         return 5;
  143     if (value < 3819626178U)
  144         return 6;
  145     if (value < 4075350088U)
  146         return 7;
  147     if (value < 4203212043U)
  148         return 8;
  149     return 9;
  150 }
  151 void CellNoise::insert(float *arr, float value)
  152 {
  153     float temp;
  154     for (int i = 8; i >= 0; i--)
  155     {
  156         if (value > arr[i])
  157             break;
  158         temp = arr[i];
  159         arr[i] = value;
  160         if (i == 0)
  161         {
  162             rd[0] = featurePoint[0];
  163             rd[1] = featurePoint[1];
  164             rd[2] = featurePoint[2];
  165         }
  166         if (i < 8)
  167             arr[i + 1] = temp;
  168     }
  169 }
  170 int CellNoise::lcgRandom(int lastValue)
  171 {
  172     return int(((1103515245u * uint(lastValue) + 12345u) % 0x100000000u));
  173 }
  174 int CellNoise::hash(int i, int j, int k)
  175 {
  176     return ((((((OFFSET_BASIS^i)*FNV_PRIME)^j) * FNV_PRIME)^k)*FNV_PRIME);
  177 }
  178 ImprovedNoise::ImprovedNoise(float xsize, float ysize, float zsize)
  179 {
  180     correction = 1.0f / 100000000.0f;
  181     passes = int(std::log(xsize) / std::log(MAGIC_SCALE) + 0.5f);
  182     passes =
  183         std::max(passes, int(std::log(ysize) / std::log(MAGIC_SCALE) + 0.5f));
  184     passes =
  185         std::max(passes, int(std::log(zsize) / std::log(MAGIC_SCALE) + 0.5f));
  186     float factor = 1.0f;
  187     for (int pass = 0; pass < passes; pass++, factor *= MAGIC_SCALE)
  188         correction += factor * factor;
  189     correction = 1.0f / std::sqrt(correction);
  190 
  191     for (int i = 0; i < 256; i++)
  192         p[256 + i] = p[i] = permutation[i];
  193 
  194     for (int i = 0; i < 256; i++)
  195     {
  196         int k = int((tinyrnd() * (256 - i) + i));
  197 
  198         int l = p[i];
  199 
  200         p[i] = p[k];
  201         p[k] = l;
  202         p[i + 256] = p[i];
  203     }
  204 }
  205 float ImprovedNoise::noise(float x, float y, float z)
  206 {
  207     int X = int(std::floor(x)) & 255, Y = int(std::floor(y)) & 255,
  208         Z = int(std::floor(z)) & 255;
  209     x -= std::floor(x);
  210     y -= std::floor(y);
  211     z -= std::floor(z);
  212     float u = fade(x), v = fade(y), w = fade(z);
  213     int A = p[X] + Y, AA = p[A] + Z, AB = p[A + 1] + Z, B = p[X + 1] + Y,
  214         BA = p[B] + Z, BB = p[B + 1] + Z;
  215 
  216     return lerp(
  217                w,
  218                lerp(v, lerp(u, grad(p[AA], x, y, z), grad(p[BA], x - 1, y, z)),
  219                     lerp(u, grad(p[AB], x, y - 1, z), grad(p[BB], x - 1, y - 1, z))),
  220                lerp(v,
  221                     lerp(u, grad(p[AA + 1], x, y, z - 1),
  222                          grad(p[BA + 1], x - 1, y, z - 1)),
  223                     lerp(u, grad(p[AB + 1], x, y - 1, z - 1),
  224                          grad(p[BB + 1], x - 1, y - 1, z - 1))));
  225 }
  226 float ImprovedNoise::fade(float f)
  227 {
  228     return f*f*f*(f*(f*6-15)+10); // t * t * (3.0 - 2.0 * t);
  229 }
  230 float ImprovedNoise::lerp(float t, float a, float b)
  231 {
  232     return a + t*(b - a);
  233 }
  234 float ImprovedNoise::grad(int hash, float x, float y, float z)
  235 {
  236     int h = hash & 15;       // CONVERT LO 4 BITS OF HASH CODE
  237     float u = h < 8 ? x : y, // INTO 12 GRADIENT DIRECTIONS.
  238           v = h < 4 ? y : h == 12 || h == 14 ? x : z;
  239     return ((h & 1) == 0 ? u : -u) + ((h & 2) == 0 ? v : -v);
  240 }
  241 float ImprovedNoise::FractalNoise3D(float x, float y, float z, int octNum,
  242                                     float lacunarity, float gain)
  243 {
  244     float freq = 1.0, amp = 1.0, sum = 0;
  245 
  246     for (int i = 0; i < octNum; i++)
  247     {
  248         sum += noise(x * freq, y * freq, z * freq) * amp;
  249         freq *= lacunarity;
  250         amp *= gain;
  251     }
  252     return sum;
  253 }
  254 float ImprovedNoise::Marble(float x, float y, float z, int octNum)
  255 {
  256     float t = 0;
  257     float factor = 1.0;
  258     for (int pass = 0; pass < octNum; pass++, factor *= MAGIC_SCALE)
  259     {
  260         float r = 1.0f / factor;
  261         t += noise(x * r, y * r, z * r) * factor;
  262     }
  263     return t * correction;
  264 }
  265 float ImprovedNoise::lookup(float x, float y, float z)
  266 {
  267     float t = 0;
  268     float factor = 1.0;
  269     for (int pass = 0; pass < passes; pass++, factor *= MAGIC_SCALE)
  270     {
  271         float r = 1 / factor;
  272         t += noise(x*r, y*r, z*r)*factor;
  273     }
  274     return t * correction;
  275 }