• Welcome to the new Internet Infidels Discussion Board, formerly Talk Freethought.

The Programming Thread

Found a better way to create fractional random numbers using rand() on the net.

First floating point numbers. For a 64 bit number there are a fixed number of digits. Spread between the integer and fractional part of the number. If you run the first code you will see the frictional digiits decrease as the integer digits increase. There is plenty of information on floating point numbers on the net.

I never figured out exactly how print() works. It looks like there is a buffer the end of which can be printed. A simple test is to assign a sequence of decimal digits to a a number and see how many are correctly, but that does not give a complete picture. The smallest number is 1/nbits the lsb. Fractional numbers will be exact only when in multiples of the lsb. So when entering fractional numbers and then displaying it may not exactly match. Quantization errors.

So, when creating float random numbers you may want to set the number of fractional digits, The second code rounds at a specified number of digits. When doing chain calculations sometimes ti helps to round intermediate values to avoid underflow. You can end up with meaningless results at times. Carrying more digits than the digits of data often has little value. A false sense of accuracy and precision.

Code:
#include <iostream>
#include "math.h"
#include "stdlib.h"
#include "stdio.h"
#include<iomanip>
using namespace std;


#define FLOOR 1
#define CEIL 2
#define ROUND 3

    double rand_float(int flag,int scale_factor, int dec_places,double lo_num ) {
        //random float between lo_num and scale_fctor

        double rfrac = 0.;
        int rint = 0,max_iterations = 100000;

        while(max_iterations--){
            rint =(rand() % scale_factor) + 1;
            if(rint > lo_num && rint < scale_factor)break;
        }//while
        if(!max_iterations){cout<<"failed to find number"<<endl;return 0;}

        rfrac = double(rand()%RAND_MAX)/double(RAND_MAX);
        rfrac = rfrac * pow(10.,dec_places);
        switch(flag){
            case ROUND:rfrac = double(round(rfrac));break;
            case CEIL:rfrac = double(ceil(rfrac));break;
            case FLOOR:rfrac = double(floor(rfrac));break;
        }
        return(rint + (rfrac/pow(10.,dec_places)));
    }//rand_float()

double dec_round(double x, int flag,int dec_places){
        // rounds the decimal part of a double
        x = x * pow(10.,dec_places);
        switch(flag){
            case ROUND:x = double(round(x));break;
            case CEIL:x = double(ceil(x));break;
            case FLOOR:x = double(floor(x));break;
        }
        return(x/pow(10.,dec_places));
}//dec_round()



int main() {
    double x;
    int i,rint;;
    unsigned seed = unsigned (time(NULL));
    srand(seed);

    //rounding
    x = 123.098765432109876543210987654321;
    printf("%.20f\n",x);
    x = dec_round(x,FLOOR,4);
    printf("X Rounded  %f\n\n",x);

   //random numbers
   printf("random numbers\n");
   for(i = 0;i<10;i++){
        x =  rand_float(FLOOR,10,3,5) ;
        printf("%f\n",x);
        }

   //loss of digits
    x = .098765432109876543210987654321;
    printf("\n%.21f\n",x);
    double z,y;
    for(i = 0;i < 51;i++){
        z = 0.;
        y = 0.;
        z =  double(pow(2,i));
        y = x + z;
        printf("%d   %.17f\n",i,y);
    }

    return(0);
}//main()
 
You never know when you might need a randomized probability distribution.


First randomizing and sorting with code to test i

Code:
const int MAXITER = 200000;
const int UP = 0;
const int DOWN = 1;

int randomize(int n,double *y,double *yr){
    int flag = 0, j = 0,i = 0, x = 0,  cnt = 0;
    int *seq = new int [n];
    for(i=0;i<n;i++)seq[i] = 0;
    for(i=0;i<n;i++){
    cnt = MAXITER;
    // create rand nunbers 1 to n with no duplication
    while(cnt--){
        x = int(rand()%n)+1;
        if(i == 0)break;
        flag = 0;
        for(j = 0;j < n;j++)if(seq[j] == x)flag = 1;
        if(!flag)break;
    }//while
    if(cnt <= 0)return(999);
    seq[i] = x;
    yr[i] = y[x-1];//arrays start at index 0
    }//i
    delete[] seq;
    return(0);;
}//randimize()


int comp_doub_descend (const void * a, const void * b) {
   return ( *(double*)b - *(double*)a );
}
int comp_doub_ascend (const void * a, const void * b) {
   return ( *(double*)a - *(double*)b );
}

int comp_int_descend (const void * a, const void * b) {
   return ( *(int*)b - *(int*)a );
}
int comp_int_ascend (const void * a, const void * b) {
   return ( *(int*)a - *(int*)b );
}

void sort_doub(double *yin,int flag,int n){
    switch(flag){
    case DOWN:
    qsort(yin, n, sizeof(double), comp_doub_descend);
    break;
    case UP:
    qsort(yin, n, sizeof(double), comp_doub_ascend);
    break;
    }
}//sort_doub()

void sort_int(int *yin,int flag,int n){
    switch(flag){
    case DOWN:
    qsort(yin, n, sizeof(int), comp_int_descend);
    break;
    case UP:
    qsort(yin, n, sizeof(int), comp_int_ascend);
    break;
    }
}//sort_int()
int main() {
    int I;
    double x[10];
    double xr[10];
    for(i=0;i<10    x[9] = 6.;
    int g = randomize(10,x,xr);
    cout<<endl;
    for(i=0;i<10;i++)printf("%f   %f\n",x[i],xr[i]);
    cout<<"DOWN"<<endl;
    sort_doub(xr,DOWN,10);
    for(i=0;i<10;i++)printf("%f   %f\n",x[i],xr[i]);
    cout<<"up"<<endl;
    sort_doub(xr,UP,10);
    for(i=0;i<10;i++)printf("%f   %f\n",x[i],xr[i]);
return(0);
}//main()
 
Adding a distribution with a random mea value, average and standard deviation, and cumulative distribution that finds the median value, the 50% point.

The cumulative distribution takes sorted data and integrates or sums point by point and computes the percent of the total summation of the data.


Code:
int exp_dist(int n,double *t,double *y, int ulo,int uhi,
                    double dec_places,double *ux){
    //find random mean
    double rfrac = 0,u;
    int i, rint = 0,cnt = MAXITERATIONS;
    while(cnt--){
        rint =(rand() % uhi) + 1;
        if(rint >= ulo && rint <= uhi)break;
    }//while
    if(cnt <= 0)return(1);//failed to find a number
    rfrac = double(rand()%RAND_MAX)/double(RAND_MAX);
    rfrac = rfrac * pow(10.,dec_places);
    rfrac = double(round(rfrac));
    rfrac = rfrac / pow(10.,dec_places);
    u = double(rint) + rfrac;

    //generate distribution and ranomize
    double k = 5.*u;
    double T = 5.*u, dt = T/double(n);
    u = u *T/k;
    *ux = u;
    double tsum = 0.;
    for(i = 0; i < n; i++){
       t[i] = tsum;
       tsum += dt;
       y[i] = k*exp(-t[i]/u);
   }//i
    return(0);
}//exp_dist()

void cum_dist(int n,double *t,double *y,double *cd, double *med){
    int i, flag = 1;
    double cum_sum = 0.;
    double ysum = 0.;
    for(i = 0;i < n;i++)cum_sum += y[i];
    for(i = 0;i < n;i++){
        ysum += y[i];
        cd[i] = 100.* ysum/cum_sum;
        if(flag && (cd[i] >= 50.)){
           flag = 0;
            *med = t[i];
        }//if

    }//for

}//cum_dist()

int  save(int n,double *x,double*y,double *z){
    FILE *ptr = fopen("data.txt","w");
    if(!ptr)return(1);
    fprintf(ptr,"%d\n",n);
    for(int i = 0;i <n;i++)
        fprintf(ptr,"%.4f    %.4f    %.4f\n",x[i],y[i],z[i]);
    fclose(ptr);
    cout<<"SAVE DONE"<<endl;
    return(0);
}//save()

void ave(int n,double *y,double *av,double *std){
    int i;
    double asum = 0, ssum = 0;
    for(i=0;i<n;i++)asum += y[i];
    *av = asum/double(n);
    for(i=0;i<n;i++)ssum += pow(*av-y[i],2);
    *std = sqrt(ssum/double(n));
}//ave()
int main() {
    unsigned seed = unsigned (time(NULL));
    srand(seed);
    int Npoints = int(1e3),k,i;
    double avg = 0, std = 0, med = 0;
    double u = 0.;
    double *y = new double[Npoints];
    double *yr = new double[Npoints];
    double *ys = new double[Npoints];
    double *t = new double[Npoints];
    double *cd = new double[Npoints];

    k = exp_dist(Npoints,t,y,180,220,4,&u);
    if(k) return(999);
    k = randomize(Npoints,y,yr);
    if(k) return(999);
    sort_doub(yr,UP,Npoints);
    cum_dist(Npoints,t,y,cd,&med);
    ave(Npoints,yr,&avg,&std);
    printf("u %.4f\n",u);
    printf("avg %.4f  std  %.4f\n",avg,std);
    printf("Median  Actaul %.4f Calc %.4f\n",med,u*log(2));
    k = save(Npoints,t,yr,cd);
    if(k) return(888)
return(0);
}//main()
 
I ws palying around with random numbers and wondered how to get around the 32 bit limit for the C library rand() function.

After deep meditation and contemplation of the nature of the universe it came to me.

Build up the random number digit by digit using rnd) to crte numbers 0-9. With this random numbers can be created for 64 bit ad 128 bit flots and ints.

A problem with __int128 numbers is neither printf or cout with print them. I found code on the net to do it.

Found the original Urand algorithm paper and am working on porting it from FORTRAN to C.


Code:
long double rand_dec128(int ndec,int ndig,double sign){
    unsigned seed = unsigned (time(NULL));
    srand(seed);

    long double val = 0;
    long double x = 0;
    int i;
    for(i=0;i<ndig;i++){
        x = (long double)(rand()%10);
        val += x*pow(10.,i);
        }
    return(sign*val/pow(10.,ndec));
}//rand_dec128()

unsigned __int128 rand_int128(int ndig,int sign){
    unsigned seed = unsigned (time(NULL));
    srand(seed);
    unsigned __int128 val = 0,x = 0;
    int i;
    for(i=0;i<ndig;i++){
        x = (unsigned __int128)(rand()%10);
        val += x*(unsigned __int128)pow(10.,i);
    }
    return(val*sign);;
}//rand_int128()

unsigned long long  rand_int64(int ndig,int sign){
    unsigned seed = unsigned (time(NULL));
    srand(seed);
    unsigned long long  val = 0, x = 0;
    int i;
    for(i=0;i<ndig;i++){
        x = (unsigned long long)(rand()%10);
        val += x*(unsigned long long)pow(10.,i);
        }
    return(sign*val);;
}//rand_int64()

void print_int128(__int128 x) {
    // requires external new line
    // source -- https://codeforces.com/blog/entry/75044

    if (x < 0) {
        putchar('-');
        x = -x;
    }
    if (x > 9) print_int128(x / 10);
    putchar(x % 10 + '0');

}//print_int128()


int main()
{

    long double val = rand_dec128(4,10,-1);
    printf("%Lf\n",val);
    cout<<fixed<<setprecision(6);
    cout<<val<<endl;

    unsigned __int128  a = rand_int128(10,-1);
    print_int128(a);
    cout<<endl;

    cout<<endl
    <<"double  "<<sizeof(double)*8<<endl
    <<"long double  "<<sizeof(long double)*8<<endl
    <<"long int  "<<sizeof(long int)*8<<endl
    <<"long long  "<<sizeof(long long)*8<<endl
    <<"__int128  "<<sizeof(__int128)*8<<endl;
    return 0;
}
 
This recent RNG thing has me thinking about a system I whipped up for making a cryptographic PRNG that allows the RNG to produce deterministic results even for asynchronously executed threads by using an addressing method.

The idea is that you start with a key and 2 64 bit values. You set the first 64 bit value, the "seed", and then generate using the second 64 bit value, the "address", so you can index from zero.

When you get a result, it is guaranteed to be a unique 128 bit value, so to guarantee at least some collisions, you can XOR the top 64 bit integer of the result against the bottom 64 bit integer (I suppose this isn't necessary; you could just pick half, but I wanted their fold...).

This gives the ability to assign RNG interactions as if assigning address space, so the thread running the calculations for an asynchronously launched event can always produce the same result, and there is no contention on the PRNG.

If you want 128 bit numbers, re-write for 256 bit encryption.

I asked chatGPT to throw a class together that accomplishes this behavior and this is what it spit out:

C++:
#include <openssl/aes.h>
#include <iostream>

class RandomNumberGenerator {
public:
    RandomNumberGenerator(const unsigned char* key, const uint64_t addr_base = 0)
        : m_key(key), m_addr_base(addr_base) {}

    // Generate a random 64-bit number using a given address offset and base
    uint64_t generate(uint64_t addr_offset, uint64_t addr_base) {
        uint64_t counter[2] = {htobe64(addr_offset), htobe64(addr_base)};
        AES_KEY aes_key;
        AES_set_encrypt_key(m_key, 128, &aes_key);
        AES_encrypt((const unsigned char*)counter, (unsigned char*)m_buffer, &aes_key);
        return be64toh(((uint64_t*)m_buffer)[0] ^ ((uint64_t*)m_buffer)[1]);
    }

    // Generate a random 64-bit number using the internal address base and a given address offset
    uint64_t generate(uint64_t addr_offset) {
        return generate(addr_offset, m_addr_base);
    }

    // Get the current address base
    uint64_t get_addr_base() const {
        return m_addr_base;
    }

    // Set the address base
    void set_addr_base(uint64_t addr_base) {
        m_addr_base = addr_base;
    }

private:
    const unsigned char* m_key;
    uint64_t m_addr_base;
    uint8_t m_buffer[16];
};

The inspiration here is I wanted to be able to perhaps write a device that instead of addressing values, addressed positions in a crazy manifold to generate its answers, without having to actually roll anything until the access is needed.

A better version would allow getting whole ranges of values.
 
Last edited:
In this code, how is __int128 defined?

One can do it with a plain-C struct, or else a C++ class with operator-overload functions for doing arithmetic with this type.

ETA: also with C++, one can define a template class for an arbitrary length of fixed-length integer, complete with operator-overloaded arithmetic.
 
Last edited:
In this code, how is __int128 defined?

One can do it with a plain-C struct, or else a C++ class with operator-overload functions for doing arithmetic with this type.

ETA: also with C++, one can define a template class for an arbitrary length of fixed-length integer, complete with operator-overloaded arithmetic.
It is defined in an include file. Somewhere in the gcc distribution.

I think there is a standard for it.

I was not aware of the method you mentioned.

You can create inters of any length you want, the way it was done going back to the 8 bit days. At the machine code level there are integer add and multiply functions with and without carry and without sign. You can add and mutply in sections of the word length of a processor, it just takes more time.

Same with floating point numbers.

On a 64 bit machine __int128 arithmetic probably takes more time than 32 or 64 bit numbers.
 
So it's implmented behind the scenes by gcc, and all one sees of it is a super long integer.
 
So it's implmented behind the scenes by gcc, and all one sees of it is a super long integer.
C is a high level abstraction.

When the compiler sees an arithmetic operation like addition it looks at the varibles and genates assembly code to do it. The code will vary depending on whether the variables are floating point or integer. Or complex numbers.

Integer addition and multiplication are done in hardware in the processor and the limitation is the processor word size. Floating point arithmetic is the same as integer, the difference is in the meaning.

All arithmetic is binary integer .

When you multiply a float by an integer in C thenthe compiler has to make adjusrmets to accomplish it. The int has to be converted to a folat in memory for the operation to work. That is what hap[ens when you cast a varianle like double(x) where x is an int.

In untyped languages it all happens behind the scene and you can freely mix numbers. There are positives and negatives to that, sometimes it causes numercal problems.

If I say x =1. with a decimal point and x = 1 with no decimal point Scilab sees it differently. One is a float and one is an integer. If I leave off the decimal point sometimes the answer is off.
 
Wrote my first code in the late 70s.

I've made a living at it for about 30 years now. Professionally, Quickbasic, Delphi and C#.
You've been round for a long time. I began in 1965, using Autocoder for a 1401 and Basic Assembler language for the 360/40 and am long retired. Then a few years using COBOL on mainframes. When last active in the field, i was using Visual Basic. No longer active. IOW, I'm a dinosaur. I was a math major and had no idea I would spend my life working with computers.
 
I looked in the likely places for the __int128 definition like defs.h. math.h, limits.h. and cmath.h and don't see it. There are include files called inside include files.

It is an extension so it may be a linked .obj in a library or something like that. I am using Code Blocks IDE so I'd have to look at the compiler and linker commands it is sending to gcc..

I know the code I posted works so its in there somewhere.
 
I compiled with GNU C++ __int128 x = 1; with no include files and there were no errors. Looks like it is integral to the compiler. When I switched to ISO C++ 17 it fails.


Extensions can lead to portability problems.



As an extension the integer scalar type __int128 is supported for targets which have an integer mode wide enough to hold 128 bits. Simply write __int128 for a signed 128-bit integer, or unsigned __int128 for an unsigned 128-bit integer. There is no support in GCC for expressing an integer constant of type __int128 for targets with long long integer less than 128 bits wide.
Looks like it may be coming
What is __ int128 in C++?


The int128 type defines a signed 128-bit integer. The API is meant to mimic an intrinsic integer type as closely as possible, so that any forthcoming int128_t can be a drop-in replacement. The int128 type supports the following: Implicit conversion from signed integral types and unsigned types narrower than 128 bits.
 
The joys of low level programming, learn something new everyday.

Looks like the target machine has to be 64 bits. Looks like _int128 requires two 64 bit ints.

The gcc command


g++.exe -Wall -std=c++17 -m64 -g -fexceptions -Wall -g -g -s -c C:\Users\admin\Documents
 
Universal Random Number Generator



The Scilab help files say their random number functions are based in URAND which put me onto the paper, dated around 1971.

Once the constants are calculated the actual sequence generation is simple.

In the FORTRAN code the code at line 10 appears to find the size of the int. Back then there wrer o data type standards and it varied fron system to system.

The DATA statement in the FORTRA code seems to be similar to a C static variable.

I don't know how the GNU rand() function is coded. One way is a static variable where a function like C srand() used to initialize the generator and another function to generate the sequence. It looks like that is what the C library does.

Another way is to encapsulate it in a class object. That way you can have multiple generators with different seeds.

The C code uses the variable names form the paper.

URAND appears to be better than rand(). You can compare them side by side over sections of the counts.

So, now I know how C rand() probably works.


Code:
static unsigned int  IY = 1;
void usrand32(unsigned int seed){
    IY = seed;
}//usrand()

unsigned int urand32(void){
    int n = (sizeof(int)*4) - 1;
    unsigned int m = pow(2,n) - 1; //matches RAND_MAX
    unsigned int  IC = 2*floor(m*(.5 - sqrt(3)/6.)) + 1;
    unsigned int  IA = 8*m*ceil(atan(1.)/8.) + 5;
    return ((IY = IY * IA + IC) % ((unsigned int)m + 1));
}

class Urand32{

    private:

    unsigned int IY = 1;
    unsigned int m,IC,IA;
    int n;

   public:

    Urand32(){
     //constructor
        n = (sizeof(int)*4) - 1;
        m = pow(2,n) - 1; //matches RAND_MAX
        IC = 2*floor(m*(.5 - sqrt(3)/6.)) + 1;
        IA = 8*m*ceil(atan(1.)/8.) + 5;
    }

    void usrand(unsigned int seed){
        IY = seed;
    }

    unsigned int urand(void){
        return((IY = IY*IA + IC)%((unsigned int)m + 1));
    }
};


int main(){

    Urand32 ur;
    unsigned int rn1,rn2;
    int i;
    int *cnt1 = new int[RAND_MAX];
    int *cnt2 = new int[RAND_MAX];
    for(i = 0;i<RAND_MAX-1;i++){cnt1[i] = 0;cnt2[i] = 0;}
    unsigned seed = unsigned (time(NULL));
    srand(seed);
    ur.usrand(seed);
    int n = int(1e9);
    for(int i = 0; i < n;i++){
        rn1 = ur.urand();
        rn2 = rand();
        cnt1[rn1] += 1;
        cnt2[rn2] += 1;
    }

    for(int i=0;i <100;i++)
     printf("%d  %d  %d\n",i,cnt1[i],cnt2[i]);

return 0;
}
 
As I've suggested before, 20th-century random number generators are about as relevant or useful these days as phonograph players or dial telephones. The PRNGs in use today, e.g. the Mersenne Twister from 1997, are far better. They are available as public-domain source code; and any worthy code library (including BSD circa 1980) will have a PRNG much better than the simple PRNGs seen in old textbooks.

These new cryptographically-secure PRNGs are overkill for many applications, but So what? If the time spent computing a random number is a bottleneck, one can use a library routine that reduces bit consumption:
* Generate a random integer 0 ≤ x < 33 --> this needs only 5.0444 random bits on average
* Generate three binary switches with p = 0.83--> this needs only 1.973 random bits on average.
(Writing the library code which recovers the unneeded bits is an interesting exercise!)

More fun, perhaps, is manipulating a random number into some other random form. Here are some one-liners that do that:

double RUF(void); // Return a uniform random variate 0 < x < 1

int rr_binary(double likely)
{ return RUF() < likely; }

int rr_small(int siz)
{ return RUF() * siz; }

int rr_geomd(double p)
{ return log(RUF()) / log(1 - p) + 1; }

double rr_expon(void)
{ return - log(RUF()); }

struct twodgaussian { double r, th; };
void rr_gaussian(struct twodgaussian *p)
{ p->th = sqrt(rr_expon() * 2); p->r = RUF() * 2 * 3.141592653589; }

double rr_cauchy()
{ return tan(RUF() * 3.141592653589); }

int rr_hispick(struct hphelper *p, int psiz)
{ double t = RUF() * psiz; return p[t]->ix[p[t]->thresh > t]; }
 
Today I had to fix a problem that arose because a dev in another team decided to persist a date (with no time) as a UTC datetime value in an application that accepts localised inputs.

It reminded me of this:



You don't even need to watch the video. Just look at the thumbnail.

Programmers have great datetime libraries in every widely-used programming language, which make datetime manipulations much easier, but none of that will save us from devs who don't understand how their system converts to/from UTC.
 
I'm probably going to stop posting here. This may be my last post.

I discover a HUGE anamoly in the Gospel genealogy, of which the many Christian "scholars" on view with Google seem oblivious! It is an enigma or a puzzle, and I asked Infidels for help resolving it. Only one response, and it implied that any solution would just be a fiction anyway, and therefore of little interest. (Briefly, James is one of three principal disciples in BOTH Paul's writing AND the synoptic Gospels. BUT the latter is James ben Zebedee brother of John, while the former is ostensibly the Lord's brother AND also James ben Alphaeus. This latter James is sometimes called James the Just and sometimes James the Less. The synoptic Gospels repeatedly refer to "another Mary" attending the crucifixion. John calls her Jesus' mother but the synoptics call her mother of James and Joses ben Alphaeus. If Jesus' father really was Joseph (or Yahweh!) then James and Joses are his HALF-brothers. ... Another possible solution is that James ben Alphaeus was "the disciple Jesus loved" that "John" writes about, and became a stepbrother with the Lord's dying words.


I post in other threads as well. Call me childish but I smile when someone clicks Like, and am less happy when someone quotes me without a Like.

I was flabbergasted by the response I got from Politesse in a recent thread. If I do continue on, I'll avoid any further dialog with [whatever his pronoun].

In this thread, I thought it fun to post some one-liners that convert a uniform variate to some other type of random variable. As usual, no response. Perhaps you all were waiting for my corrections! :)

Here [corrected]:
#define MASK (MAXRAND >> 2)
double RUF(void)
{ return (random() & MASK | 1) / (double)(MASK+1); }

struct twodgaussian { double r, th; };
void rr_2dgauss(struct twodgaussian *p)
{ p->r = sqrt(log(RUF() * -2); p->th = RUF() * 2 * 3.141592653589; }

int rr_hispick(struct hphelper p[], int psiz)
{ double t = RUF() * psiz; return p[t]->ix[p[t]->thresh > t]; }

For the math nerd out there: WHY is it convenient (natural?) to return a variate over an actual bell-shape (2D manifold) rather than over the "bell curve", i.e. the cross-section of such a bell with a piece of paper?

rr_hispick() actually implements a very useful and oft-seen random selection. Huge kudos if you if you can write code to setup the hphelper array by reverse-engineering the one-liner.
 
Back
Top Bottom