Groups | Search | Server Info | Keyboard shortcuts | Login | Register [http] [https] [nntp] [nntps]

Groups > comp.lang.c > #42767 > unrolled thread

Request for source code review of simple Ising model

Back to article view | Back to comp.lang.c

Contents

  Request for source code review of simple Ising model Udyant Wig <udyant@panda.goosenet.in> - 2014-04-10 23:38 +0530
    Re: Request for source code review of simple Ising model Malcolm McLean <malcolm.mclean5@btinternet.com> - 2014-04-10 13:30 -0700
      Re: Request for source code review of simple Ising model Udyant Wig <udyantw@gmail.com> - 2014-04-12 00:34 +0530
        Re: Request for source code review of simple Ising model Malcolm McLean <malcolm.mclean5@btinternet.com> - 2014-04-11 12:25 -0700
        Re: Request for source code review of simple Ising model Udyant Wig <udyantw@gmail.com> - 2014-04-16 01:12 +0530
          Re: Request for source code review of simple Ising model Ben Bacarisse <ben.usenet@bsb.me.uk> - 2014-04-16 01:33 +0100
            Re: Request for source code review of simple Ising model Udyant Wig <udyantw@gmail.com> - 2014-04-16 17:57 +0530
              Re: Request for source code review of simple Ising model Malcolm McLean <malcolm.mclean5@btinternet.com> - 2014-04-16 06:00 -0700
              Re: Request for source code review of simple Ising model James Kuyper <jameskuyper@verizon.net> - 2014-04-16 10:58 -0400
                Re: Request for source code review of simple Ising model Keith Thompson <kst-u@mib.org> - 2014-04-16 09:07 -0700
                  Re: Request for source code review of simple Ising model James Kuyper <jameskuyper@verizon.net> - 2014-04-16 12:22 -0400
                    Re: Request for source code review of simple Ising model Keith Thompson <kst-u@mib.org> - 2014-04-16 10:38 -0700
                      Re: Request for source code review of simple Ising model Kaz Kylheku <kaz@kylheku.com> - 2014-04-16 19:37 +0000
                  Re: Request for source code review of simple Ising model Udyant Wig <udyantw@gmail.com> - 2014-04-17 12:00 +0530
                    Re: Request for source code review of simple Ising model glen herrmannsfeldt <gah@ugcs.caltech.edu> - 2014-04-17 09:32 +0000
                      Re: Request for source code review of simple Ising model "BartC" <bc@freeuk.com> - 2014-04-17 12:24 +0100
                        Re: Request for source code review of simple Ising model Keith Thompson <kst-u@mib.org> - 2014-04-17 08:16 -0700
                      Re: Request for source code review of simple Ising model Malcolm McLean <malcolm.mclean5@btinternet.com> - 2014-04-17 04:41 -0700
                        Re: Request for source code review of simple Ising model Les Cargill <lcargill99@comcast.com> - 2014-04-17 07:49 -0500
                    Re: Request for source code review of simple Ising model Keith Thompson <kst-u@mib.org> - 2014-04-17 08:00 -0700
                      Re: Request for source code review of simple Ising model Kaz Kylheku <kaz@kylheku.com> - 2014-04-17 15:06 +0000
                        Re: Request for source code review of simple Ising model Tim Rentsch <txr@alumni.caltech.edu> - 2014-04-18 21:04 -0700
                      Re: Request for source code review of simple Ising model glen herrmannsfeldt <gah@ugcs.caltech.edu> - 2014-04-17 18:38 +0000
                        Re: Request for source code review of simple Ising model James Kuyper <jameskuyper@verizon.net> - 2014-04-17 15:26 -0400
                          Re: Request for source code review of simple Ising model "BartC" <bc@freeuk.com> - 2014-04-17 21:00 +0100
                            Re: Request for source code review of simple Ising model James Kuyper <jameskuyper@verizon.net> - 2014-04-17 16:34 -0400
                            Re: Request for source code review of simple Ising model glen herrmannsfeldt <gah@ugcs.caltech.edu> - 2014-04-18 04:31 +0000
                        Re: Request for source code review of simple Ising model Keith Thompson <kst-u@mib.org> - 2014-04-17 14:05 -0700
                          Re: Request for source code review of simple Ising model Keith Thompson <kst-u@mib.org> - 2014-04-17 17:06 -0700
                          Re: Request for source code review of simple Ising model glen herrmannsfeldt <gah@ugcs.caltech.edu> - 2014-04-18 04:54 +0000
                          Re: Request for source code review of simple Ising model Tim Rentsch <txr@alumni.caltech.edu> - 2014-04-18 11:51 -0700
                Re: Request for source code review of simple Ising model Udyant Wig <udyantw@gmail.com> - 2014-04-17 11:42 +0530
              Re: Request for source code review of simple Ising model Ben Bacarisse <ben.usenet@bsb.me.uk> - 2014-04-16 16:54 +0100
                Re: Request for source code review of simple Ising model Malcolm McLean <malcolm.mclean5@btinternet.com> - 2014-04-16 15:22 -0700
                Re: Request for source code review of simple Ising model Udyant Wig <udyantw@gmail.com> - 2014-04-17 13:03 +0530
                  Re: Request for source code review of simple Ising model Keith Thompson <kst-u@mib.org> - 2014-04-17 08:44 -0700
                    Re: Request for source code review of simple Ising model James Kuyper <jameskuyper@verizon.net> - 2014-04-17 12:12 -0400
                      Re: Request for source code review of simple Ising model Udyant Wig <udyantw@gmail.com> - 2014-04-18 12:38 +0530
                      Re: Request for source code review of simple Ising model Malcolm McLean <malcolm.mclean5@btinternet.com> - 2014-04-18 03:11 -0700
                    Re: Request for source code review of simple Ising model Udyant Wig <udyantw@gmail.com> - 2014-04-18 12:32 +0530
                      Re: Request for source code review of simple Ising model James Kuyper <jameskuyper@verizon.net> - 2014-04-18 10:05 -0400
          Re: Request for source code review of simple Ising model Ike Naar <ike@iceland.freeshell.org> - 2014-04-16 07:10 +0000
          Re: Request for source code review of simple Ising model James Kuyper <jameskuyper@verizon.net> - 2014-04-17 13:10 -0400
            Re: Request for source code review of simple Ising model Udyant Wig <udyantw@gmail.com> - 2014-04-18 18:24 +0530
              Re: Request for source code review of simple Ising model James Kuyper <jameskuyper@verizon.net> - 2014-04-18 09:52 -0400
                Re: Request for source code review of simple Ising model Udyant Wig <udyantw@gmail.com> - 2014-04-18 20:04 +0530
                  Re: Request for source code review of simple Ising model James Kuyper <jameskuyper@verizon.net> - 2014-04-18 11:00 -0400
                    Re: Request for source code review of simple Ising model Udyant Wig <udyantw@gmail.com> - 2014-04-20 01:49 +0530
                      Re: Request for source code review of simple Ising model glen herrmannsfeldt <gah@ugcs.caltech.edu> - 2014-04-19 21:20 +0000
                        Re: Request for source code review of simple Ising model Udyant Wig <udyantw@gmail.com> - 2014-04-20 11:11 +0530
                          Re: Request for source code review of simple Ising model Richard Damon <Richard@Damon-Family.org> - 2014-04-20 09:06 -0400
                        Re: Request for source code review of simple Ising model James Kuyper <jameskuyper@verizon.net> - 2014-04-20 10:50 -0400
                          Re: Request for source code review of simple Ising model Malcolm McLean <malcolm.mclean5@btinternet.com> - 2014-04-21 09:14 -0700
                Re: Request for source code review of simple Ising model Ike Naar <ike@iceland.freeshell.org> - 2014-04-19 05:43 +0000
              Re: Request for source code review of simple Ising model Udyant Wig <udyantw@gmail.com> - 2014-04-18 19:47 +0530
              Re: Request for source code review of simple Ising model Udyant Wig <udyantw@gmail.com> - 2014-04-21 01:21 +0530
                Re: Request for source code review of simple Ising model Keith Thompson <kst-u@mib.org> - 2014-04-20 14:24 -0700
                  Re: Request for source code review of simple Ising model Udyant Wig <udyantw@gmail.com> - 2014-04-21 10:40 +0530
                Re: Request for source code review of simple Ising model Udyant Wig <udyantw@gmail.com> - 2014-04-21 12:29 +0530
            Re: Request for source code review of simple Ising model Keith Thompson <kst-u@mib.org> - 2014-04-18 08:36 -0700
              Re: Request for source code review of simple Ising model James Kuyper <jameskuyper@verizon.net> - 2014-04-18 13:40 -0400
                Re: Request for source code review of simple Ising model Keith Thompson <kst-u@mib.org> - 2014-04-18 11:14 -0700
                Re: Request for source code review of simple Ising model Malcolm McLean <malcolm.mclean5@btinternet.com> - 2014-04-18 13:55 -0700
              Re: Request for source code review of simple Ising model Udyant Wig <udyantw@gmail.com> - 2014-04-19 01:25 +0530
                Re: Request for source code review of simple Ising model Udyant Wig <udyantw@gmail.com> - 2014-04-19 01:44 +0530
                Re: Request for source code review of simple Ising model Kaz Kylheku <kaz@kylheku.com> - 2014-04-18 20:24 +0000
                Re: Request for source code review of simple Ising model Keith Thompson <kst-u@mib.org> - 2014-04-18 13:36 -0700
                  Re: Request for source code review of simple Ising model Udyant Wig <udyantw@gmail.com> - 2014-04-19 20:18 +0530
                Re: Request for source code review of simple Ising model Richard <rgrdev_@gmail.com> - 2014-04-19 14:55 +0100
                  Re: Request for source code review of simple Ising model Udyant Wig <udyantw@gmail.com> - 2014-04-20 02:10 +0530
              Re: Request for source code review of simple Ising model Ike Naar <ike@iceland.freeshell.org> - 2014-04-19 05:40 +0000
                Re: Request for source code review of simple Ising model Keith Thompson <kst-u@mib.org> - 2014-04-19 14:39 -0700
                  Re: Request for source code review of simple Ising model Ian Collins <ian-news@hotmail.com> - 2014-04-20 10:44 +1200
                    Re: Request for source code review of simple Ising model Keith Thompson <kst-u@mib.org> - 2014-04-19 17:11 -0700
                    Re: Request for source code review of simple Ising model Udyant Wig <udyantw@gmail.com> - 2014-04-20 11:24 +0530
                  Re: Request for source code review of simple Ising model Keith Thompson <kst-u@mib.org> - 2014-04-19 17:19 -0700
                    Re: Request for source code review of simple Ising model Ike Naar <ike@iceland.freeshell.org> - 2014-04-20 10:10 +0000
                      Re: Request for source code review of simple Ising model Keith Thompson <kst-u@mib.org> - 2014-04-20 14:21 -0700
                        Re: Request for source code review of simple Ising model Ike Naar <ike@iceland.freeshell.org> - 2014-04-20 22:18 +0000
                          Re: Request for source code review of simple Ising model Keith Thompson <kst-u@mib.org> - 2014-04-20 19:27 -0700
                            Re: Request for source code review of simple Ising model Ike Naar <ike@iceland.freeshell.org> - 2014-04-21 12:38 +0000
                            Re: Request for source code review of simple Ising model Udyant Wig <udyantw@gmail.com> - 2014-04-21 18:12 +0530
                            Re: Request for source code review of simple Ising model Ike Naar <ike@iceland.freeshell.org> - 2014-04-21 13:13 +0000
                              Re: Request for source code review of simple Ising model Udyant Wig <udyantw@gmail.com> - 2014-04-21 19:46 +0530
                              Re: Request for source code review of simple Ising model Keith Thompson <kst-u@mib.org> - 2014-04-21 09:06 -0700
                  Re: Request for source code review of simple Ising model jacob navia <jacob@spamsink.net> - 2014-04-21 00:26 +0200
          Re: Request for source code review of simple Ising model Udyant Wig <udyantw@gmail.com> - 2014-04-18 01:15 +0530
          Re: Request for source code review of simple Ising model Udyant Wig <udyantw@gmail.com> - 2014-04-20 11:46 +0530

Page 3 of 5 — ← Prev page 1 2 [3] 4 5  Next page →

#43067

On 04/18/2014 03:02 AM, Udyant Wig wrote:
...
>  But when might this situation arise?  That is, could there be a string
>  some of whose elements were characters < 0?

The basic source and execution character sets are defined as containing
the following characters, which I've listed in the order in which they
are referred to in section 5.2.1:

"\0ABCDEFGHIJKLMNOPQRSTUVWXYZabcdefghijklmnopqrstuvwxyz0123456789!"
"!\"#%&'()*+,-./:;<=>?[\\]^_{|}~ \t\v\f"

In addition, the basic execution character set must also include "\a\b\r\n"

Note: The basic source character set is NOT required to contain '\n'.
Source code files are required to have some method of indicating the end
of a line. That method is treated as if it were a newline character at
the end of each line, regardless of whether or not that is the method
actually used.

Implementations are free to support additional characters, these are
referred to as being members of the extended character set.

"If a member of the basic execution character set is stored in a
char object, its value is guaranteed to be nonnegative. If any other
character is stored in a char object, the resulting value is
implementation-defined but shall be within the range of values that can
be represented in that type." (6.2.5p3)

"The implementation shall define char to have the same range,
representation, and behavior as either signed char or unsigned char."
(6.2.5p15).

This wasn't an arbitrary decision by the committee - at the time that
the standard was written, there were many implementations where plain
char was signed, an many others where it was unsigned.

If an implementation chooses for char to have the same range as signed
char, then for any member of the extended execution character set, if
that character is stored in a char object, the implementation is free to
define it to have a negative value. This is normally the case for some
of those characters on any system where char is signed that uses just
about any character encoding other than 7-bit ASCII.
-- 
James Kuyper

[toc] | [prev] | [next] | [standalone]

#43005

On 2014-04-15, Udyant Wig <udyantw@gmail.com> wrote:
> bool is_all_zeroes (char *string)
> {
>     char *sp;
>
>     for (sp = string; *sp != '\0'; sp++) {
> 	if (*sp != '0') {
> 	    return false;
> 	}
>     }
>
>     return true;
> }    

For operations like these, the strspn function comes in handy:

bool is_all_zeroes(char const *string)
{
  return string[strspn(string,"0")] == '\0';
}

[toc] | [prev] | [next] | [standalone]

#43046

On 04/15/2014 03:42 PM, Udyant Wig wrote:
...
> 	    *cell |= ((byte) (rand () % 2)) << 1;

The C standard doesn't impose any significant restrictions on the
quality of implementation of rand(), and on many implementations it's
not very good. It might be good enough, if your needs are simple enough.
However, if you're going to rely upon rand() despite the fact that it
might not be very good, you should be aware of the fact that, in poor
quality implementations, lower-order bits are likely to be significantly
less random than higher-order ones. RAND_MAX is required to be >= 32767,
so 0x4000 is the highest bit that guaranteed to be <RAND_MAX. If you're
only going extract one bit from rand(), I'd recommend extracting that
one: ((rand()%0x4000) == 0x4000 ? 2 : 1).

If, instead, you have access to high-quality random number generator,
it's a waste to extract only one bit from each number: you should use
each of the bits to initialize a different cell of your lattice.

And that brings to mind another issue. Once you're sure that your
program is working, and you've reached the point where it's reasonable
to worry about performance, you might consider an alternative data
storage strategy. You're using one byte to store one cell representing a
spin 1/2 object: it therefore has only two quantized orientations, up
and down, so you're storing only one bit of information for that cell.
It might be better to use, for example, a 32-bit integer to store the
values of 32 consecutive cells. Not only will this use up a lot less
memory, but with (quite) a bit of ingenuity, you can use bit-wise
operations to handle all 32 cells at the same time, for a naive speed-up
by a factor of 32. In practice, the speed up will actually be less than
that, but should still be substantial.
On the other hand, if you implement this idea your code will be a lot
more complicated, and much harder to understand. You'll have to decide
whether that's a acceptable cost for the increased processing speed and
decreased memory requirements.

[toc] | [prev] | [next] | [standalone]

#43064

James Kuyper <jameskuyper@verizon.net> writes:
| On 04/15/2014 03:42 PM, Udyant Wig wrote:
| ...
|> 	    *cell |= ((byte) (rand () % 2)) << 1;
|
| The C standard doesn't impose any significant restrictions on the
| quality of implementation of rand(), and on many implementations it's
| not very good. It might be good enough, if your needs are simple enough.
| However, if you're going to rely upon rand() despite the fact that it
| might not be very good, you should be aware of the fact that, in poor
| quality implementations, lower-order bits are likely to be significantly
| less random than higher-order ones. RAND_MAX is required to be >= 32767,
| so 0x4000 is the highest bit that guaranteed to be <RAND_MAX. If you're
| only going extract one bit from rand(), I'd recommend extracting that
| one: ((rand()%0x4000) == 0x4000 ? 2 : 1).

 How does this work?  Is it ever the case that

                     x mod 16384 = 16384?
                     
| If, instead, you have access to high-quality random number generator,
| it's a waste to extract only one bit from each number: you should use
| each of the bits to initialize a different cell of your lattice.

 It took me a while to work this out.  I used the exact integers
 facility of C99 via <stdint.h>.  I allocated an array of int16_t and
 then, for as many cells in the lattice, I took bits from a 16-bit
 integer in the array until I had run through all bits in the one,
 whereupon I went to the next one.

 This necessitated changes to the allocators and initializers:

/* 16-bit integer bit source */
int16_t *allocate_bitsource (size_t size)
{
    int16_t *bitsource;
    bitsource = malloc (size);
    if (bitsource == NULL) {
	fprintf (stderr, "allocate_bitsource: %s\n", strerror (errno));
	exit (EXIT_FAILURE);
    }

    return bitsource;
}

void initialize_bitsource (int16_t *bitsource)
{
    int i;
    int int_count = minimum_multiple_16 (dimension * dimension);

    for (i = 0; i < int_count; i++) {
	*bitsource = (int16_t)(rand () % RAND_MAX);
    }
}


void initialize_lattice (byte *lattice, int16_t *bitsource)
{
    int16_t *ip = bitsource;
    byte *bp = lattice, *end = bp + dimension * dimension;
    int i;
    
    for (i = 0; bp < end; i++) {
	if (i == 16) {
	    ip++;
	    i = 0;
	}

	*bp = ((byte)(*ip >> i)) & 0x01;
	bp++;
    }
}

/* defined in utilities.c */
int minimum_multiple_16 (int n)
{
    int mm16 = n >> 4;
    if (n % 16 != 0) {
	mm16++;
    }
    return mm16;
}    


| And that brings to mind another issue. Once you're sure that your
| program is working, and you've reached the point where it's reasonable
| to worry about performance, you might consider an alternative data
| storage strategy. You're using one byte to store one cell representing
| a spin 1/2 object: it therefore has only two quantized orientations,
| up and down, so you're storing only one bit of information for that
| cell.  It might be better to use, for example, a 32-bit integer to
| store the values of 32 consecutive cells. Not only will this use up a
| lot less memory, but with (quite) a bit of ingenuity, you can use
| bit-wise operations to handle all 32 cells at the same time, for a
| naive speed-up by a factor of 32. In practice, the speed up will
| actually be less than that, but should still be substantial.  On the
| other hand, if you implement this idea your code will be a lot more
| complicated, and much harder to understand. You'll have to decide
| whether that's a acceptable cost for the increased processing speed
| and decreased memory requirements.

 This would be another rewrite of the program, as much of the existing
 source would have to be modified for this data representation.

 It should be worth looking into, because, as I said elsewhere, a given
 run of a 10x10 lattice can take upwards of three hours and not
 complete, while a 6x6 one finishes in mere minutes.

[toc] | [prev] | [next] | [standalone]

#43065

On 04/18/2014 08:54 AM, Udyant Wig wrote:
> James Kuyper <jameskuyper@verizon.net> writes:
...
> | so 0x4000 is the highest bit that guaranteed to be <RAND_MAX. If you're
> | only going extract one bit from rand(), I'd recommend extracting that
> | one: ((rand()%0x4000) == 0x4000 ? 2 : 1).
> 
>  How does this work?  Is it ever the case that
> 
>                      x mod 16384 = 16384?
It doesn't - that was a typo. It should have been

	(rand()&0x4000 == 0x4000) ? 2 : 1

Sorry!
-- 
James Kuyper

[toc] | [prev] | [next] | [standalone]

#43070

James Kuyper <jameskuyper@verizon.net> writes:
| 	(rand()&0x4000 == 0x4000) ? 2 : 1

 Is this what this does: if RAND_MAX is in fact 32767, then this piece
 of code produces 2 or 1 with roughly equal probability?
 
| Sorry!

 No harm done.  The original piece was puzzling, is all.

[toc] | [prev] | [next] | [standalone]

#43071

On 04/18/2014 10:34 AM, Udyant Wig wrote:
> James Kuyper <jameskuyper@verizon.net> writes:
> | 	(rand()&0x4000 == 0x4000) ? 2 : 1
> 
>  Is this what this does: if RAND_MAX is in fact 32767, then this piece
>  of code produces 2 or 1 with roughly equal probability?

Yes. If RAND_MAX is not 2^n-1 for some integer n (which is permitted),
then the equality would be much rougher. For best results, there's no
substitute for paying close attention to the value of RAND_MAX, and
adjusting your method for using the values returned by rand()
accordingly, but if you need to achieve a precisely specified
probability, that's not easy to do. For instance, your code uses

double random_value = (double) (rand () % 100) / 100.0;

It is guaranteed that 0 <= random_value && random_value < 1.0, and
(after allowing for floating point round-off), random_value will be an
integer multiple of 0.01.  However, it will not produce all 100
different possible values with equal probability (not even if rand()
were an ideal random number generator). See if you can figure out why.
Hint: if RAND_MAX were 32799, it would work.

For this purpose, I'd normally use rand()/(RAND_MAX+1.0), which would
produce each possible value with equal probability if rand() were an
ideal random number generator. However, if you need to implement a
precise rational probability m/n, where m and n are small integers, it's
a very real issue.
-- 
James Kuyper

[toc] | [prev] | [next] | [standalone]

#43132

James Kuyper <jameskuyper@verizon.net> writes:
| For instance, your code uses
|
| double random_value = (double) (rand () % 100) / 100.0;
|
| It is guaranteed that 0 <= random_value && random_value < 1.0, and
| (after allowing for floating point round-off), random_value will be an
| integer multiple of 0.01.  However, it will not produce all 100
| different possible values with equal probability (not even if rand()
| were an ideal random number generator). See if you can figure out why.
| Hint: if RAND_MAX were 32799, it would work.

 Qualitatively, if RAND_MAX is 32767, then rand () % 100 will produce
 values between 0 and 67 with a higher probability than it will those
 between 68 and 99.

[toc] | [prev] | [next] | [standalone]

#43136

Udyant Wig <udyantw@gmail.com> wrote:

(snip)

> Qualitatively, if RAND_MAX is 32767, then rand () % 100 will produce
> values between 0 and 67 with a higher probability than it will those
> between 68 and 99.

And how many calls will it take for the difference to be
statistically significant?

-- glen

[toc] | [prev] | [next] | [standalone]

#43164

ram@zedat.fu-berlin.de (Stefan Ram) writes:
|   This program shows the relative deviations in per thousand
|   when using »rand() %« in one single case:
|
| #include <stdio.h>  /* printf */
| #include <stdlib.h> /* RAND_MAX, rand() */
| #define N 5e6
| #define COVER 10
| int a[ COVER + 2 ];
| int main( void )
| { for( int i = 0; i < N; ++i )++a[ 1 + rand() % COVER ];
|   for( int i = 0; i < COVER; ++i )
|   printf( "%d:%+f\n", i, 1000*( a[ i + 1 ]/( N / COVER )- 1 )); }
|
| 0:+0.396000
| 1:-1.370000
| 2:+0.070000
| 3:-0.572000
| 4:-0.154000
| 5:-1.236000
| 6:+0.004000
| 7:+4.234000
| 8:-0.470000
| 9:-0.902000

 I get this:
 
0:-2.388000
1:-0.510000
2:-1.028000
3:-1.276000
4:+3.250000
5:-0.390000
6:-1.664000
7:+1.018000
8:+2.154000
9:+0.834000

 There is a greater number of 4's and 8's.

[toc] | [prev] | [next] | [standalone]

#43174

On 4/20/14, 1:41 AM, Udyant Wig wrote:
> 
>  There is a greater number of 4's and 8's.
> 

The first question to ask yourself, what answer are you really expecting.

If you are expecting that you will get EXACTLY the same number of each
result, that is wrong. You expect that there will be some (small)
variation in each number, and that variation is expected to have certain
statistics itself.

One thing to note is that the expected variation will go down by the
square root of the number of samples taken, so to a rough order, taking
1e6 samples expects errors on the order of 1e-3 (there are a few scaling
factors, but this gives you a rough order of magnitude).

The big test, is if you run it with a different seed, do different
numbers come out ahead. That is a strong indication of bias, as would be
running a LOT more samples. What you can do from this sample is set a
likely limit to how biased the numbers will be.

[toc] | [prev] | [next] | [standalone]

#43182

On 04/19/2014 05:20 PM, glen herrmannsfeldt wrote:
> Udyant Wig <udyantw@gmail.com> wrote:
> 
> (snip)
> 
>> Qualitatively, if RAND_MAX is 32767, then rand () % 100 will produce
>> values between 0 and 67 with a higher probability than it will those
>> between 68 and 99.
> 
> And how many calls will it take for the difference to be
> statistically significant?

A relatively tiny number, for any program involving serious number
crunching. It's a difference of 1/327, which is fairly large number for
certain purposes, and a completely negligible one for other purposes -
the key point is making sure you know which category your particular
case involves.
-- 
James Kuyper

[toc] | [prev] | [next] | [standalone]

#43230

On Sunday, April 20, 2014 3:50:32 PM UTC+1, James Kuyper wrote:
> On 04/19/2014 05:20 PM, glen herrmannsfeldt wrote:
> 
> 
> A relatively tiny number, for any program involving serious number
> crunching. It's a difference of 1/327, which is fairly large number for
> certain purposes, and a completely negligible one for other purposes -
> the key point is making sure you know which category your particular
> case involves.
> 
As you take the temperatue down slowly enough, an Ising model will always
settle in one ground state (all blue or all white) or the other. Without 
doing  the maths, which is largely beyond me, a slight bias to blue or 
white will mean that it always settles in that state and never the other.

[toc] | [prev] | [next] | [standalone]

#43114

On 2014-04-18, James Kuyper <jameskuyper@verizon.net> wrote:
> On 04/18/2014 08:54 AM, Udyant Wig wrote:
>> James Kuyper <jameskuyper@verizon.net> writes:
> ...
>> | so 0x4000 is the highest bit that guaranteed to be <RAND_MAX. If you're
>> | only going extract one bit from rand(), I'd recommend extracting that
>> | one: ((rand()%0x4000) == 0x4000 ? 2 : 1).
>> 
>>  How does this work?  Is it ever the case that
>> 
>>                      x mod 16384 = 16384?
> It doesn't - that was a typo. It should have been
>
> 	(rand()&0x4000 == 0x4000) ? 2 : 1

This produces 1 or 2, depending on the output of rand().

The original expression was

  ((byte) (rand () % 2)) << 1

which produces 0 or 2, depending on the output of rand().

[toc] | [prev] | [next] | [standalone]

#43068

Udyant Wig <udyantw@gmail.com> writes:
| void initialize_bitsource (int16_t *bitsource)
| {
|     int i;
|     int int_count = minimum_multiple_16 (dimension * dimension);
|
|     for (i = 0; i < int_count; i++) {
| 	*bitsource = (int16_t)(rand () % RAND_MAX);
|     }
| }

 Note that this code does not work as would be expected (it keeps
 setting the same integer.)  Either add

bitsource++;

 to the loop body, or use this rewrite (using the loop form suggested by
 Ben Bacarisse):

void initialize_bitsource (int16_t *bitsource)
{
    int int_count = minimum_multiple_16 (dimension * dimension);

    int16_t *ip = bitsource, *end = ip + int_count * sizeof *ip;
    while (ip < end) {
    	*ip = (int16_t)(rand () % RAND_MAX);
    	ip++;
    }
}    

[toc] | [prev] | [next] | [standalone]

#43187

Udyant Wig <udyantw@gmail.com> writes:
| James Kuyper <jameskuyper@verizon.net| writes:
| | And that brings to mind another issue. Once you're sure that your
| | program is working, and you've reached the point where it's
| | reasonable to worry about performance, you might consider an
| | alternative data storage strategy. You're using one byte to store
| | one cell representing a spin 1/2 object: it therefore has only two
| | quantized orientations, up and down, so you're storing only one bit
| | of information for that cell.  It might be better to use, for
| | example, a 32-bit integer to store the values of 32 consecutive
| | cells. Not only will this use up a lot less memory, but with (quite)
| | a bit of ingenuity, you can use bit-wise operations to handle all 32
| | cells at the same time, for a naive speed-up by a factor of 32. In
| | practice, the speed up will actually be less than that, but should
| | still be substantial.  On the other hand, if you implement this idea
| | your code will be a lot more complicated, and much harder to
| | understand. You'll have to decide whether that's a acceptable cost
| | for the increased processing speed and decreased memory
| | requirements.
|
|  This would be another rewrite of the program, as much of the existing
|  source would have to be modified for this data representation.
|
|  It should be worth looking into, because, as I said elsewhere, a
|  given run of a 10x10 lattice can take upwards of three hours and not
|  complete, while a 6x6 one finishes in mere minutes.

 I found this to be somewhat difficult initially.  After thinking and
 writing, and some more thinking and rewriting, I developed the appended
 source.  Though I have not yet found out a way to handle 32 bits at
 once, I did manage to use as many 32-bit integers as are required to
 house dimension * dimension bits.  So, a single int32_t can contain
 upto 5x5 bits, while two can contain upto 8x8, etc.

 Testing showed that this is faster than the earlier version wherein I
 used the integers as mere sources of bits for the lattice.

 The general structure remains the same.

 
/* common.h begins */

#ifndef COMMON_H
#define COMMON_H

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <errno.h>
#include <math.h>
#include <stdbool.h>
#include <time.h>
#include <ctype.h>
#include <stdint.h>

#endif

/* common.h ends */


/* utilities.h begins */

#ifndef UTILITIES_H
#define UTILITIES_H

#include "common.h"

bool is_all_zeroes (char *string);
bool is_positive_integer (char *string);
int  how_many (const char *s, int c);
char *copy_substring (char *substring, const char *string, size_t start, size_t end);
int  minimum_multiple_32 (int n);

#endif

/* utilities.h ends */


/* utilities.c begins */

#include "utilities.h"

bool is_all_zeroes (char *string)
{
    return string [strspn (string, "0")] == '\0';
}    

bool is_positive_integer (char *string)
{
    char *sp;

    for (sp = string; *sp != '\0'; sp++) {
	if (!isdigit (*sp)) {
	    return false;
	}
    }

    return !is_all_zeroes (string);
}    

/* From /C: A Reference Manual/, 5th ed., by Harbison and Steele
 * page 352
 */
int how_many (const char *s, int c)
{
    int n = 0;
    if (c == 0) return 0;
    while (s) {
	s = strchr (s, c);
	if (s) n++, s++;
    }
    return n;
}    

char *copy_substring (char *substring, const char *string, size_t start, size_t end)
{
    size_t length = end - start + 1;

    errno = 0;
    substring = malloc (1 + length);
    if (substring == NULL) {
	fprintf (stderr, "copy_substring: %s\n", strerror (errno));
	exit (1);
    }

    memset (substring, 0, 1 + length);
    strncpy (substring, string + start, length - 1);
    substring [length] = '\0';
    
    return substring;
}

int minimum_multiple_32 (int n)
{
    int mm32 = n >> 5;
    if (n % 32 != 0) {
	mm32++;
    }
    return mm32;
}    

/* utilities.c ends */


/* bitsing.h begins */

#ifndef BITSING_H
#define BITSING_H

/* 32-bit integer bit source that doubles as lattice */
int32_t *allocate_lattice (size_t size);
void    initialize_lattice (int32_t *lattice);
void    print_lattice (int32_t *lattice);

/* Energies vector */
int  *allocate_energies (size_t size);
void initialize_energies (int *energies);

/* Core */
int  check_and_add_neighbor (int32_t *lattice, int32_t cell, int col, int row);
void iterate (int32_t *lattice, int *energies);

/* Miscellaneous */
int  count_ones (int32_t *lattice);
void print_iteration (int32_t *lattice, const char *format_string, int count);

/* Error checking */
void minority_error (void);
void beta_error (void);
void dimension_error (void);

#endif

/* bitsing.h ends */


/* bitsing.c begins */

#include "common.h"
#include "bitsing.h"
#include "utilities.h"

#define INTSRC_BITS 32

static int    dimension		= 3;
static double beta		= 0.5;
static int    minority_size	= 1;

/* 32-bit integer bit source that doubles as lattice */
int32_t *allocate_lattice (size_t size)
{
    int32_t *lattice;
    lattice = malloc (size);
    if (lattice == NULL) {
	fprintf (stderr, "allocate_lattice: %s\n", strerror (errno));
	exit (2);
    }
    return lattice;
}

void initialize_lattice (int32_t *lattice)
{
    int int_count = minimum_multiple_32 (dimension * dimension);
    int32_t *ip = lattice, *end = ip + int_count;
    while (ip < end) {
	*ip = (int32_t) (rand () % RAND_MAX);
	ip++;
    }
}    

void print_lattice (int32_t *lattice)
{
    int row, col;
    int position;

    for (col = 0; col < dimension; col++) {
	for (row = 0; row < dimension; row++) {
	    position = col * dimension + row;
	    printf ("%d", (lattice [position / INTSRC_BITS] >> (position % INTSRC_BITS)) & 1);
	}
	putchar ('\n');
    }
}


/* Energies vector */
int *allocate_energies (size_t size)
{
    int *energies;
    energies = malloc (size);
    if (energies == NULL) {
	fprintf (stderr, "allocate_energies: %s\n", strerror (errno));
	exit (2);
    }
    return energies;
}

void initialize_energies (int *energies)
{
    int *ep = energies, *end = ep + dimension * dimension;
    while (ep < end) {
	*ep = 0;
	ep++;
    }
}


/* Core */
int check_and_add_neighbor (int32_t *lattice, int cell_bit, int col, int row)
{
    int neighbor_bit;
    int position;
    if (0 <= col && 0 <= row && col < dimension && row < dimension) {
	position = col * dimension + row;
	neighbor_bit = (lattice [position / INTSRC_BITS] >> (position % INTSRC_BITS)) & 1;
	return cell_bit ^ neighbor_bit;
    }
    return 0;
}

void iterate (int32_t *lattice, int *energies)
{
    int rrow, rcol;
    int position;
    int rcell_bit, ccell_bit;
    int old_energy;
    int new_energy;
    int delta;

    rrow = rand () % dimension;
    rcol = rand () % dimension;
    position = rcol * dimension + rrow;
    rcell_bit =  (lattice [position / INTSRC_BITS] >> (position % INTSRC_BITS)) & 1;
    ccell_bit = rcell_bit ^ 1;	/* comlemented cell */
    old_energy = energies [position];

    new_energy = check_and_add_neighbor (lattice, ccell_bit, rrow - 1, rcol);
    new_energy = check_and_add_neighbor (lattice, ccell_bit, rrow + 1, rcol);
    new_energy = check_and_add_neighbor (lattice, ccell_bit, rrow, rcol - 1);
    new_energy = check_and_add_neighbor (lattice, ccell_bit, rrow, rcol + 1);
	
    delta = new_energy - old_energy;
    if (delta < 0) {
	lattice [position / INTSRC_BITS] ^= (1 << (position % INTSRC_BITS));
	energies [position] = new_energy;
    }
    else {
	double ising_window = exp (-(beta * delta));
	double random_value = (double) (rand () % 100) / 100.0;

	if (random_value < ising_window) {
	    lattice [position / INTSRC_BITS] ^= (1 << (position % INTSRC_BITS));
	    energies [position] = new_energy;
	}
    }
}    


/* Miscellaneous */
int count_ones (int32_t *lattice)
{
    int count = 0;
    int int_count = minimum_multiple_32 (dimension * dimension);
    int32_t *ip = lattice, *end = ip + int_count;
    int field;
    int i = 0;
    
    while (ip < end) {
	for (field = *ip; field; field >>= 1) {
	    if (i == dimension * dimension) {
		return count;
	    }
	    count += field & 1;
	    i++;
	}	
	ip++;
    }

    return count;
}

static const char initial_format [] = \
    "Initial configuration %6d\n----------------------------\n";
static const char iteration_format [] = \
    "Iteration %6d\n----------------\n";
static const char final_format [] = \
    "Final configuration %6d\n--------------------------\n";

void print_iteration (int32_t *lattice, const char *format_string, int count)
{
    putchar ('\n');
    printf (format_string, count);
    print_lattice (lattice);
}


/* Error checking */
void minority_error (void)
{
    exit (2);
}

void beta_error (void)
{
    exit (2);
}

void dimension_error (void)
{
    exit (2);
}    


/* Return codes:
 * 0 -- success
 * 1 -- memory allocation failure
 * 2 -- invalid command line argument
 */
int main (int argc, char *argv [])
{
    if (argc > 3) {
	if (is_positive_integer (argv [3]) || is_all_zeroes (argv [3])) {
	    minority_size = atoi (argv [3]);
	}
	else {
	    minority_error ();
	}
	
	if (minority_size < 0) {
	    minority_size = 0;
	}
	else if (minority_size > dimension * dimension) {
	    minority_size = dimension * dimension;
	}
	else if (minority_size > (dimension * dimension) / 2) {
	    minority_size = (dimension * dimension) - minority_size;
	}
    }
    if (argc > 2) {
	char *beta_string = argv [2];
	size_t length = strlen (beta_string);
	
	if (how_many (beta_string, '.') == 1) {
	    size_t position_decimal = strcspn (beta_string, ".");
	    char *integral_part = NULL;
	    char *fractional_part = NULL;

	    if (position_decimal == 0) {
		beta_error ();
	    }
	    if (position_decimal == strlen (beta_string) - 1) {
		beta_error ();
	    }

	    integral_part = copy_substring (integral_part, \
					    beta_string, \
					    0, \
					    position_decimal - 1);
	    if (is_positive_integer (integral_part) || \
		is_all_zeroes (integral_part)) {
		fractional_part = copy_substring (fractional_part, \
						  beta_string, \
						  position_decimal + 1, \
						  length - 1);
		if (is_positive_integer (fractional_part) || \
		    is_all_zeroes (fractional_part)) {
		    beta = atof (argv [2]);
		}
		else {
		    beta_error ();
		}

		free (fractional_part);
	    }
	    else {
		beta_error ();
	    }

	    free (integral_part);
	}
	else {
	    beta_error ();
	}
    }
    if (argc > 1) {
	if (is_positive_integer (argv [1]))
	    dimension = atoi (argv [1]);
	else {
	    dimension_error ();
	}
    }

    srand (time (0));

    int32_t *lattice;
    lattice = allocate_lattice (minimum_multiple_32 (dimension * dimension) \
				* sizeof *lattice);
    initialize_lattice (lattice);

    int *energies;
    energies = allocate_energies (dimension * dimension * sizeof *energies);
    initialize_energies (energies);

    print_iteration (lattice, initial_format, 0);

    long count = 1;
    while (true) {
    	iterate (lattice, energies);

    	int ones_count = count_ones (lattice);
    	if ((ones_count == minority_size) || \
    	    ((dimension * dimension - ones_count) == minority_size)) {
    	    break;
    	}
	
    	print_iteration (lattice, iteration_format, count++);
    }

    print_iteration (lattice, final_format, count);
    free (energies);
    free (lattice);

    exit (0);
}    

/* bitsing.c ends */

[toc] | [prev] | [next] | [standalone]

#43191

Udyant Wig <udyantw@gmail.com> writes:
[...]
> bool is_positive_integer (char *string)
> {
>     char *sp;
>
>     for (sp = string; *sp != '\0'; sp++) {
> 	if (!isdigit (*sp)) {
> 	    return false;
> 	}
>     }
>
>     return !is_all_zeroes (string);
> }    

I see you're still not casting the argument to isdigit().  That causes
undefined behavior if plain char is signed and your string contains a
character with a negative value (which probably represents a character
greater than 127).

[...]

-- 
Keith Thompson (The_Other_Keith) kst-u@mib.org  <http://www.ghoti.net/~kst>
Working, but not speaking, for JetHead Development, Inc.
"We must do something.  This is something.  Therefore, we must do this."
    -- Antony Jay and Jonathan Lynn, "Yes Minister"

[toc] | [prev] | [next] | [standalone]

#43211

Keith Thompson <kst-u@mib.org> writes:

| Udyant Wig <udyantw@gmail.com> writes:
| [...]
>> bool is_positive_integer (char *string)
>> {
>>     char *sp;
>>
>>     for (sp = string; *sp != '\0'; sp++) {
>> 	if (!isdigit (*sp)) {
>> 	    return false;
>> 	}
>>     }
>>
>>     return !is_all_zeroes (string);
>> }    
>
| I see you're still not casting the argument to isdigit().  That causes
| undefined behavior if plain char is signed and your string contains a
| character with a negative value (which probably represents a character
| greater than 127).
>
| [...]

 That correction managed to evade commission until now.

bool is_positive_integer (char *string)
{
    char *sp;

    for (sp = string; *sp != '\0'; sp++) {
	if (!isdigit ((unsigned char)*sp)) {
	    return false;
	}
    }

    return !is_all_zeroes (string);
}    
 

[toc] | [prev] | [next] | [standalone]

#43215

Udyant Wig <udyantw@gmail.com> writes:
| int count_ones (int32_t *lattice)
| {
|     int count = 0;
|     int int_count = minimum_multiple_32 (dimension * dimension);
|     int32_t *ip = lattice, *end = ip + int_count;
|     int field;
|     int i = 0;
|     
|     while (ip < end) {
| 	for (field = *ip; field; field >>= 1) {
| 	    if (i == dimension * dimension) {
| 		return count;
| 	    }
| 	    count += field & 1;
| 	    i++;
| 	}	
| 	ip++;
|     }
|
|     return count;
| }

int count_ones (int32_t *lattice)
{
    int count = 0;
    int int_count = minimum_multiple_32 (dimension * dimension);
    int32_t *ip = lattice, *end = ip + int_count;
    int32_t field; /* <-- type change */
    int i = 0;
    
    while (ip < end) {
	for (field = *ip; field; field >>= 1) {
	    if (i == dimension * dimension) {
		return count;
	    }
	    count += field & 1;
	    i++;
	}	
	ip++;
    }

    return count;
}

[toc] | [prev] | [next] | [standalone]

#43074

James Kuyper <jameskuyper@verizon.net> writes:
> On 04/15/2014 03:42 PM, Udyant Wig wrote:
> ...
>> 	    *cell |= ((byte) (rand () % 2)) << 1;
>
> The C standard doesn't impose any significant restrictions on the
> quality of implementation of rand(), and on many implementations it's
> not very good. It might be good enough, if your needs are simple enough.
> However, if you're going to rely upon rand() despite the fact that it
> might not be very good, you should be aware of the fact that, in poor
> quality implementations, lower-order bits are likely to be significantly
> less random than higher-order ones. RAND_MAX is required to be >= 32767,
> so 0x4000 is the highest bit that guaranteed to be <RAND_MAX. If you're
> only going extract one bit from rand(), I'd recommend extracting that
> one: ((rand()%0x4000) == 0x4000 ? 2 : 1).

(As discussed downthread, that "%" was meant to be "&".)

Is it still the case that the low-order bits are commonly less random
than the high-order bits?  Question 13.18 of the FAQ
<http://www.c-faq.com/> says so, but I suspect it's out of date.

I think I've seen PRNGs where the low-order bit actually alternates
0, 1, 0, 1, 0, 1, ..., but the sample implementation in the standard
doesn't do that.

Here's a test program if anyone wants to try it.  If it prints either
"01010101..." or "10101010...", please let us know what implementation
you're using.

#include <stdio.h>
#include <stdlib.h>

int main(void) {
    for (int i = 0; i < 64; i ++) {
        printf("%d", rand() % 2);
    }
    putchar('\n');
}

On one of my systems, the output is

1011110011010110000010110001111000111010111101001010100100011101

-- 
Keith Thompson (The_Other_Keith) kst-u@mib.org  <http://www.ghoti.net/~kst>
Working, but not speaking, for JetHead Development, Inc.
"We must do something.  This is something.  Therefore, we must do this."
    -- Antony Jay and Jonathan Lynn, "Yes Minister"

[toc] | [prev] | [next] | [standalone]

Page 3 of 5 — ← Prev page 1 2 [3] 4 5  Next page →

Back to top | Article view | comp.lang.c

csiph-web