242 lines
4 KiB
C
242 lines
4 KiB
C
/***************************************************************************
|
|
|
|
gbx_math.c
|
|
|
|
(c) 2000-2017 Benoît Minisini <benoit.minisini@gambas-basic.org>
|
|
|
|
This program is free software; you can redistribute it and/or modify
|
|
it under the terms of the GNU General Public License as published by
|
|
the Free Software Foundation; either version 2, or (at your option)
|
|
any later version.
|
|
|
|
This program is distributed in the hope that it will be useful,
|
|
but WITHOUT ANY WARRANTY; without even the implied warranty of
|
|
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
|
GNU General Public License for more details.
|
|
|
|
You should have received a copy of the GNU General Public License
|
|
along with this program; if not, write to the Free Software
|
|
Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston,
|
|
MA 02110-1301, USA.
|
|
|
|
***************************************************************************/
|
|
|
|
#define __GBX_MATH_C
|
|
|
|
#include "gb_common.h"
|
|
|
|
#include <math.h>
|
|
#include <time.h>
|
|
#include <sys/time.h>
|
|
|
|
#include "gb_hash.h"
|
|
#include "gbx_math.h"
|
|
|
|
const double MATH_pow10_double[] = {1, 10, 100, 1000, 10000, 100000, 1000000, 10000000, 100000000, 1000000000};
|
|
static const uint _pow10_uint[] = {1, 10, 100, 1000, 10000, 100000, 1000000, 10000000, 100000000, 1000000000};
|
|
|
|
/* This is a twisted generalized feedback shift register
|
|
that generates pseudo-random numbers.
|
|
|
|
Source code from the paper of Yann Guidon
|
|
in March edition of GNU/Linux Magazine France
|
|
*/
|
|
|
|
#define CRC32_STD 0x04C11DB7
|
|
#define GFSR_SIZE 15
|
|
|
|
static uint GFSR_table[GFSR_SIZE];
|
|
static uint GFSR_temp;
|
|
static uint GFSR_index;
|
|
|
|
static void GFSR_init(uint seed)
|
|
{
|
|
int i = 0, j;
|
|
uint t = seed;
|
|
|
|
if (t == 0)
|
|
t = (uint)-1;
|
|
|
|
do
|
|
{
|
|
t ^= (t >> 5) ^ (t << 1);
|
|
j = t >> 31;
|
|
t <<= 1;
|
|
if (j)
|
|
t ^= CRC32_STD;
|
|
GFSR_table[i++] = t;
|
|
}
|
|
while (i < GFSR_SIZE);
|
|
|
|
GFSR_temp = GFSR_table[GFSR_SIZE - 1];
|
|
GFSR_index = 0;
|
|
}
|
|
|
|
static uint GFSR_random(void)
|
|
{
|
|
int t;
|
|
|
|
GFSR_temp ^= GFSR_table[GFSR_index];
|
|
|
|
t = GFSR_temp;
|
|
GFSR_temp <<= 1;
|
|
if (t < 0)
|
|
GFSR_temp ^= CRC32_STD;
|
|
|
|
GFSR_table[GFSR_index] = GFSR_temp;
|
|
|
|
if (++GFSR_index >= GFSR_SIZE)
|
|
GFSR_index = 0;
|
|
|
|
return GFSR_temp;
|
|
}
|
|
|
|
double frac(double x)
|
|
{
|
|
x = fabs(x);
|
|
return x - floor(x);
|
|
}
|
|
|
|
int lsgn(int x)
|
|
{
|
|
return ((x > 0) ? 1 : ((x < 0) ? (-1) : 0));
|
|
}
|
|
|
|
int llsgn(int64_t x)
|
|
{
|
|
return ((x > 0) ? 1 : ((x < 0) ? (-1) : 0));
|
|
}
|
|
|
|
/*int64_t llabs(int64_t x)
|
|
{
|
|
return ((x < 0) ? (-x) : x);
|
|
}*/
|
|
|
|
int fsgn(double x)
|
|
{
|
|
return ((x > 0) ? 1 : ((x < 0) ? (-1) : 0));
|
|
}
|
|
|
|
float fixf(float x)
|
|
{
|
|
if (x >= 0)
|
|
return floorf(x);
|
|
else
|
|
return -floorf(fabsf(x));
|
|
}
|
|
|
|
double fix(double x)
|
|
{
|
|
if (x >= 0)
|
|
return floor(x);
|
|
else
|
|
return -floor(fabs(x));
|
|
}
|
|
|
|
double frexp10(double x, int *exp)
|
|
{
|
|
int p;
|
|
|
|
if (x == 0.0)
|
|
{
|
|
*exp = 0;
|
|
return x;
|
|
}
|
|
|
|
p = (int)log10(x);
|
|
x /= pow10(p);
|
|
|
|
if (x >= 1)
|
|
{
|
|
x /= 10;
|
|
p++;
|
|
}
|
|
|
|
*exp = p;
|
|
return x;
|
|
}
|
|
|
|
|
|
void randomize(bool set, uint seed)
|
|
{
|
|
struct timeval tv;
|
|
|
|
if (!set && gettimeofday(&tv, NULL) == 0)
|
|
seed = 0xD1C2B3A4 + (tv.tv_sec << 20) + tv.tv_usec;
|
|
|
|
GFSR_init(seed);
|
|
}
|
|
|
|
|
|
double rnd(void)
|
|
{
|
|
/*seed = 16807L * (seed % 127773L) - 2836L * (seed / 127773L);
|
|
if (seed <= 0) seed += 2147483647;
|
|
|
|
return (double)seed / 2147483648.0;*/
|
|
|
|
uint64_t val;
|
|
|
|
val = GFSR_random();
|
|
val <<= 32;
|
|
val |= GFSR_random();
|
|
|
|
return (double)val / 18446744073709551616.0; //0xFFFFFFFFFFFFFFFFULL;
|
|
}
|
|
|
|
#ifndef HAVE_EXP10
|
|
double exp10(double x)
|
|
{
|
|
return pow(10, x);
|
|
}
|
|
#endif
|
|
|
|
#ifndef HAVE_LOG2
|
|
double log2(double x)
|
|
{
|
|
return log(x) / M_LN2;
|
|
}
|
|
#endif
|
|
|
|
#ifndef HAVE_EXP2
|
|
double exp2(double x)
|
|
{
|
|
return pow(2, x);
|
|
}
|
|
#endif
|
|
|
|
static int nbits(uint n)
|
|
{
|
|
uint c;
|
|
for (c = 0; n; c++)
|
|
n &= n - 1;
|
|
return c;
|
|
}
|
|
|
|
void MATH_init(void)
|
|
{
|
|
uint seed;
|
|
|
|
randomize(FALSE, 0);
|
|
|
|
// Internet condom
|
|
do
|
|
seed = GFSR_random();
|
|
while ((seed & 1) == 0 || nbits(seed) < 16);
|
|
|
|
HASH_seed = seed;
|
|
}
|
|
|
|
uint64_t pow10_uint64_p(int n)
|
|
{
|
|
uint64_t v = 1;
|
|
|
|
while (n > 8)
|
|
{
|
|
v *= 100000000;
|
|
n -= 8;
|
|
}
|
|
v *= _pow10_uint[n];
|
|
return v;
|
|
}
|
|
|