Google Groups no longer supports new Usenet posts or subscriptions. Historical content remains viewable.
Dismiss

ordering integer array

1 view
Skip to first unread message

Ron Ford

unread,
Aug 15, 2008, 7:08:32 PM8/15/08
to

I have a program whose output is an array of integers.

implicit none

integer, parameter:: sides = 6
integer, parameter:: trials = 100

integer, dimension(sides)::A
integer, dimension(trials)::C
integer:: b, clock, seed, ii, i, counter, tab, penultimate
real:: harvest


! seed random num generator
CALL SYSTEM_CLOCK(COUNT=clock)
seed = clock + 37
CALL RANDOM_SEED(PUT = seed)

! prime the pump
call random_number(harvest)
b = 3 + nint(10*harvest)
do i=1,b
call random_number(harvest)
print *, i, harvest
end do


! main control

! outer loop is the number of trials
C = 0
do ii=1,trials
tab = 0
A = 0
! inner loop rolls die till all values attained
do
call random_number(harvest)
b = ceiling(harvest*real(sides))
A(b) = A(b) + 1
! count until all outcomes non-zero
counter = 0
do i = 1, sides
if (real(A(i)) > .5) then
counter = counter + 1
end if
end do

tab = tab + 1
!print *, "counter= ", counter
!print *, A
!print *, "tab= ", tab
if (counter == sides) then
penultimate = tab -1
C(ii) = penultimate
end if

if (counter == sides) exit
end do !inner
end do !outer

print *, "C is", C
end program

1 0.176787
2 0.458366
3 0.228078
4 0.711985
C is 14 26 5 27 12
12
9 15 18
19 12 39 10 11 23
9 16 12
7 11 8 17 18 21
18 15 12
14 18 14 29 26 16
11 21 18
6 7 20 11 12 15
10 7 16
15 8 8 5 11 16
20 11 13
15 13 9 20 13 21
23 8 9
5 12 23 13 12 7
11 10 11
11 10 14 14 13 12
15 7 12
9 7 7 20 9 6
9 6 14
11 11 14 12 6 16
13 8 14
26

Press RETURN to close window . . .


What is a slick way to order C?

Thanks and cheers,

--
When a new source of taxation is found it never means, in practice, that
the old source is abandoned. It merely means that the politicians have two
ways of milking the taxpayer where they had one before. 8
H. L. Mencken

dpb

unread,
Aug 15, 2008, 7:46:58 PM8/15/08
to
Ron Ford wrote:
> I have a program whose output is an array of integers.
>
> implicit none
...
> integer, dimension(trials)::C
...

> What is a slick way to order C?

Need a sort routine...(DOH! :) ).

Depending on the compiler you're using, it may or may not have an
extended intrinsic or compatibility sort module.

CVF has SORTQQ or QSORT; can't speak for others altho the lineage would
suggest Intel would. There are a myriad of sort routines available on
the 'net from Netlib and similar sources. One can also write a
braindead version relatively easily which would be adequate for
small-sized arrays...

--

e p chandler

unread,
Aug 15, 2008, 8:51:58 PM8/15/08
to

If you have a good upper bound on the values in C and it is small,
then you can just throw each value into a bin as you generate it!

Imagine that you want to create a frequency histogram of the values in
the array C. Suppose you know or guess what the largest value in the C
array might be, let's say 50.

integer HIST(50)
....

HIST=0

....

HIST(penultimate) = HIST(penultimate) + 1

....

do ii = 1,50
if (HIST(ii) > 0) print *,ii,HIST(ii)
end do

....

The same idea can produce a "wonder sort" if you have a set of values
in a limited range. See Jon Bentley's "Programming Pearls". Column 1
"Cracking the Oyster".

IMO this trick is worth remembering!

- e

Ron Ford

unread,
Aug 15, 2008, 10:28:54 PM8/15/08
to
On Fri, 15 Aug 2008 18:46:58 -0500, dpb posted:

> Ron Ford wrote:
>> I have a program whose output is an array of integers.
>>
>> implicit none
> ...
>> integer, dimension(trials)::C
> ...
>> What is a slick way to order C?
>
> Need a sort routine...(DOH! :) ).

I'm a bit of an idiot savant with sorting. The savant part is when I
actually do it; the idiot part is when I do it with a syntax.:-)


>
> Depending on the compiler you're using, it may or may not have an
> extended intrinsic or compatibility sort module.
>
> CVF has SORTQQ or QSORT; can't speak for others altho the lineage would
> suggest Intel would. There are a myriad of sort routines available on
> the 'net from Netlib and similar sources. One can also write a
> braindead version relatively easily which would be adequate for
> small-sized arrays...

I was hoping for an intrinsic. minloc and maxloc won't help?

Found a good link searching for this stuff:

http://www.fortran.com/tools.html
--
We must be willing to pay a price for freedom. 4
H. L. Mencken

dpb

unread,
Aug 15, 2008, 10:53:30 PM8/15/08
to
Ron Ford wrote:
...
> I was hoping for an intrinsic. minloc and maxloc won't help?
Hopes dashed, sorry... :)

One could cobble a solution up w/ minloc() and/or maxloc(), sure, but it
would be pretty ugly -- repetitive calling minloc() on the array w/o the
preceding location included either by indexing and/or moving elements or
via MASK optional argument would yield a solution eventually, but it
would be pretty inefficient sort algorithm.

--

Ron Ford

unread,
Aug 15, 2008, 11:11:28 PM8/15/08
to
On Fri, 15 Aug 2008 21:53:30 -0500, dpb posted:

Alright, well I'll do the keystrokes to implement MR&C's selection sort.
Since millions of particles have hit my retinas in the last hours, I can
promise a solition based on typing. It'll take me some time to find the
page again.

Fish, chips and 997 points of light,
--
What men value in this world is not rights but privileges. 7
H. L. Mencken

Ron Ford

unread,
Aug 15, 2008, 11:51:15 PM8/15/08
to
On Fri, 15 Aug 2008 17:51:58 -0700 (PDT), e p chandler posted:


16 bins is what I remember from UPS. You'd be top tan at 48105. It's a
no-looker: one step and behind.


>
> Imagine that you want to create a frequency histogram of the values in
> the array C. Suppose you know or guess what the largest value in the C
> array might be, let's say 50.
>
> integer HIST(50)
> ....
>
> HIST=0
>
> ....
>
> HIST(penultimate) = HIST(penultimate) + 1
>
> ....
>
> do ii = 1,50
> if (HIST(ii) > 0) print *,ii,HIST(ii)
> end do
>
> ....


>
> The same idea can produce a "wonder sort" if you have a set of values
> in a limited range. See Jon Bentley's "Programming Pearls". Column 1
> "Cracking the Oyster".

I don't like shellfish so don't see the motivation.


>
> IMO this trick is worth remembering!
>
> - e

Sorting is not something a person learns from the standard, nor MR&C:
neither selection nor sort appear in the index.

Ron Ford

unread,
Aug 17, 2008, 2:49:58 PM8/17/08
to
On Fri, 15 Aug 2008 17:51:58 -0700 (PDT), e p chandler posted:

> On Aug 15, 7:08 pm, Ron Ford <r...@example.invalid> wrote:

module sort
implicit none
private
public :: selection_sort
integer, parameter :: string_length = 30
type, public :: address
character(len = string_length) :: name, &
street, town, state*2
integer :: zip_code
end type address
contains
recursive subroutine selection_sort (array_arg)
type (address), dimension(:), intent (inout) &
:: array_arg
integer :: current_size
integer :: big
current_size = size (array_arg)
if (current_size > 0) then
big = maxloc(array_arg(:)%zip_code, dim=1)
call swap (big, current_size)
call selection_sort (array_arg(1: current_size -1))
end if
contains
subroutine swap(i,j)
integer, intent (in) :: i,j
type (address) :: temp
temp = array_arg(i)
array_arg(i) = array_arg(j)
array_arg(j) = temp
end subroutine swap
end subroutine selection_sort
end module sort

program zippy
use sort
implicit none
integer, parameter :: array_size = 100
type (address), dimension (array_size) :: data_array
integer :: i,n
do i = 1, array_size
read (*, '(/a/a/a/a2,i8)',end=10) data_array(i)
write (*, '(/a/a/a/a2,i8)') data_array(i)
end do
10 n = i - 1
call selection_sort (data_array(1:n))
write(*, '(//a)') 'after sorting'
do i = 1, n
write (*, '(/a/a/a/a2,i8)') data_array(i)
end do
end program zippy

T. Boone Pickens
350 Swiftboat Lane
Houston, TX 45536

Ev Bayh
1016 Humdity Bath
Indianapolis, IN 47250

Cindy McCain
1000 Richie Rich
Scottsdale, AZ 85250

Diane Feinstein
250 Harvey Milk Circle
San Fransisco, CA 96580

Tim Pawlenty
500 Olson way
St. Paul, MN 56674

Donald Trump
100 Trump Plaza
New York, NY 10080

This compiles but I get a runtime when I paste in the addresses to the dos
window. I would think that herding these into an external file would be
best. I don't see many examples of this in MR&C but in 10.3 they have
open(2, iostat=ios, err=99, file...

As far as I can tell with what I read of the defaults, I could use
open(2, file='addresses.txt')
or give it better diagnostics with
open(2, file='addresses.txt', iostat=ios, err=99)

The runtime was "invalid character in field."

--
We are here and it is now. Further than that, all human knowledge is
moonshine. 3
H. L. Mencken

Ron Ford

unread,
Aug 18, 2008, 12:05:12 AM8/18/08
to
On Fri, 15 Aug 2008 17:51:58 -0700 (PDT), e p chandler posted:

> On Aug 15, 7:08 pm, Ron Ford <r...@example.invalid> wrote:

>> What is a slick way to order C?
>
> If you have a good upper bound on the values in C and it is small,
> then you can just throw each value into a bin as you generate it!

50 bins would suffice well for this sort. I use maxput to deal with the
outliers:

implicit none

integer, parameter:: sides = 6
integer, parameter:: trials = 100

integer, parameter:: bins = 100
integer, parameter:: percentile = 95

integer, dimension(sides)::A
integer, dimension(trials)::C

integer, dimension(bins)::D
integer:: b, clock, seed, ii, i, counter, &
tab, maxput, penultimate, goal
real:: harvest, tot


! main control

if (counter == sides) then
C(ii) = tab
exit
end if


end do !inner
end do !outer

!print *, "C is", C

! sort into bins
D=0
maxput = bins
do ii =1, bins

if (C(ii) > maxput) then
C(ii) = maxput
end if

D(C(ii))=D(C(ii)) + 1

end do
tot= sum(D)
print *, tot, D

! determine 95th percentile
tot = 0.01 *trials*real(percentile)
print *, "95th % is", tot

counter = 0
do ii = 1, bins
counter = counter + D(ii)
print *, counter

if (real(counter).ge. tot) then
exit
end if
end do
print *, "95th percentile was in bin number", ii


end program

This compiles and behaves, giving the statistic of "how many times do you
have to roll a die until you observe all outcomes, keeping rolling, and
counting the number of times you do it, and when you hit ii, this number is
95% unlikely."

Unfortunately, I get a runtime when the number of bins and the number of
trials is unequal. Right now, I can't step through the main control with
F8'ing in my debugger, ie, stepping over; this might be a good thing to
relegate to a function or subroutine, as I don't need a return value from
the main control.

--
We must respect the other fellow's religion, but only in the sense and to
the extent that we respect his theory that his wife is beautiful and his
children smart. 5
H. L. Mencken

e p chandler

unread,
Aug 18, 2008, 4:35:39 PM8/18/08
to
On Aug 18, 12:05 am, Ron Ford <r...@example.invalid> wrote:
> On Fri, 15 Aug 2008 17:51:58 -0700 (PDT), e p chandler posted:
>
> > On Aug 15, 7:08 pm, Ron Ford <r...@example.invalid> wrote:
> >> What is a slick way to order C?
>
> > If you have a good upper bound on the values in C and it is small,
> > then you can just throw each value into a bin as you generate it!
>
> 50 bins would suffice well for this sort.  I use maxput to deal with the
> outliers:
>
> implicit none
>
> integer, parameter:: sides = 6
> integer, parameter:: trials = 100
> integer, parameter:: bins = 100
> integer, parameter:: percentile = 95
>
> integer, dimension(sides)::A
> integer, dimension(trials)::C
> integer, dimension(bins)::D
> integer:: b, clock, seed, ii, i, counter, &
>           tab, maxput, penultimate, goal
> real:: harvest, tot
>
> ! seed random num generator
> CALL SYSTEM_CLOCK(COUNT=clock)
> seed = clock + 37
> CALL RANDOM_SEED(PUT = seed)

This does not compile with G95. Perhaps it is a compiler error. G95
does not document this call but
under GNU fortran,

integer size
....
CALL RANDOM_SEED(size)

size is the minimum size of seed which is an array.

see http://gcc.gnu.org/onlinedocs/gfortran/RANDOM_005fSEED.html

Perhaps the interface is different for your compiler (Are you using
Silverfrost's FTN95?)

** running with array bounds (-fbounds-check for g95) checking turned
on reveals this error:
** Should be **
do ii =1, trials

See comments above.

[flame on]

For some reason I find it easier to debug a problem program from a
printed listing instead of looking at it on the screen. It's easy to
mark off flow of control, underline used variables, check off
interfaces, etc.. YMMV.

[flame off]

- e


Richard Maine

unread,
Aug 18, 2008, 5:18:02 PM8/18/08
to
e p chandler <ep...@juno.com> wrote:

> On Aug 18, 12:05 am, Ron Ford <r...@example.invalid> wrote:

...


> > integer:: b, clock, seed, ii, i, counter, &
> > tab, maxput, penultimate, goal

...


> > CALL RANDOM_SEED(PUT = seed)
>
> This does not compile with G95. Perhaps it is a compiler error. G95
> does not document this call

No, it is not a compiler error. It is an error in the code. The reason
for the error was discussed previously in this thread - that argument is
required to be a rank 1 array rather than a scalar.

Admitedly, the message I get from g95 is a bit confusing, but that's
another matter. I can sort of see how that message comes about.
(The second message. The first one I got was from the invalid characters
in the continuation line. Hex CA? Anyway, I replaced those by blanks.)

This is a standard intrinsic documented in the Fortran standard. I
wouldn't expect the g95 documentation to duplicate the documentation of
the zillions of standard Fortran intrinsics. That would seem like a
pretty big waste of developer time, which has plenty of other cals on
it. Documenting the compiler-specific bits (such as what random number
generator is used), sure, but what a waste it would be to repeat the
standard's documentation of the standard intrinsic interfaces
(particularly if one goes to the trouble of rewriting all the material
to avoid stepping on ANSI copyright toes instead of literally copying
it.)

--
Richard Maine | Good judgement comes from experience;
email: last name at domain . net | experience comes from bad judgement.
domain: summertriangle | -- Mark Twain

James Giles

unread,
Aug 18, 2008, 6:01:43 PM8/18/08
to
Richard Maine wrote:
> [...] but what a waste it would be to repeat the

> standard's documentation of the standard intrinsic interfaces
> (particularly if one goes to the trouble of rewriting all the material
> to avoid stepping on ANSI copyright toes instead of literally copying
> it.)

Would INCITS or ISO really enforce their copyright if a compiler
manual used some content like: "the standard says: '...'? Wouldn't
that constitute fair use?

--
J. Giles

"I conclude that there are two ways of constructing a software
design: One way is to make it so simple that there are obviously
no deficiencies and the other way is to make it so complicated
that there are no obvious deficiencies." -- C. A. R. Hoare

"Simplicity is prerequisite for reliability" -- E. W. Dijkstra


e p chandler

unread,
Aug 18, 2008, 7:28:16 PM8/18/08
to
On Aug 18, 5:18 pm, nos...@see.signature (Richard Maine) wrote:

> e p chandler <e...@juno.com> wrote:
>
>
>
> > On Aug 18, 12:05 am, Ron Ford <r...@example.invalid> wrote:
> ...
> > > integer:: b, clock, seed, ii, i, counter, &
> > >           tab, maxput, penultimate, goal
> ...
> > > CALL RANDOM_SEED(PUT = seed)
>
> > This does not compile with G95. Perhaps it is a compiler error. G95
> > does not document this call
>
> No, it is not a compiler error. It is an error in the code. The reason
> for the error was discussed previously in this thread - that argument is
> required to be a rank 1 array rather than a scalar.

Now that you explain it the source of the error is clear. I
interpreted the error message that I saw as being about seed not being
a proper positional argument. That led me astray. Also I believe it
was discussed in a _related_ thread - obviously I missed the key
point.

> Admitedly, the message I get from g95 is a bit confusing, but that's
> another matter. I can sort of see how that message comes about.
> (The second message. The first one I got was from the invalid characters
> in the continuation line. Hex CA? Anyway, I replaced those by blanks.)

> This is a standard intrinsic documented in the Fortran standard. I
> wouldn't expect the g95 documentation to duplicate the documentation of
> the zillions of standard Fortran intrinsics.

OK. I should have RTFM. RANDOM_NUMBER was documented BTW.

- e

Richard Maine

unread,
Aug 18, 2008, 8:00:43 PM8/18/08
to
James Giles <james...@worldnet.att.net> wrote:

> Richard Maine wrote:
> > [...] but what a waste it would be to repeat the
> > standard's documentation of the standard intrinsic interfaces
> > (particularly if one goes to the trouble of rewriting all the material
> > to avoid stepping on ANSI copyright toes instead of literally copying
> > it.)
>
> Would INCITS or ISO really enforce their copyright if a compiler
> manual used some content like: "the standard says: '...'? Wouldn't
> that constitute fair use?

If you quoted the entire chapter on intrinsics (or for that matter, the
entire standard) with just the preface that "the standard says", I doubt
that a fair use defense would hold up. And such copying of pretty much
the whole chapter is what I was talking about.

But then I'm most definitely neither an expert on the relevant laws nor
on the other legal ramifications (including such things as the odds of
winning a lawsuit, but losing big bucks on legal fees in the process).
I'm even more so not an expert on what INCITS or ISO would do.

For example, I don't actually think that ISO has a valid copyright
claim. I wrote some of what they claim copyright to. I never assigned
copyright to them or had any relevant employment agreement. (From some
of what I read, they might think they employed me, but assignment of
copyright to an employer requires a written agreement saying such, which
we didn't have. Also, sort of an odd kind of "employment" when I have to
pay them to be "employed" by them). Yet they want me to pay them
royalties for copying some of the material that I wrote myself. Yes,
this *HAS* come up; we decided to rewrite the material in question
instead of paying their "ransom".

I don't think they have a case that would hold up. But I'm also not a
lawyer. Nor am I about to spend the money to hire one to defend myself
in such a case. Even if I won, the lawyer would cost me more than the
royalties in question.

James Giles

unread,
Aug 18, 2008, 10:08:52 PM8/18/08
to
Richard Maine wrote:
> James Giles <james...@worldnet.att.net> wrote:
>
>> Richard Maine wrote:
>>> [...] but what a waste it would be to repeat the
>>> standard's documentation of the standard intrinsic interfaces
>>> (particularly if one goes to the trouble of rewriting all the
>>> material to avoid stepping on ANSI copyright toes instead of
>>> literally copying it.)
>>
>> Would INCITS or ISO really enforce their copyright if a compiler
>> manual used some content like: "the standard says: '...'? Wouldn't
>> that constitute fair use?
>
> If you quoted the entire chapter on intrinsics (or for that matter,
> the entire standard) with just the preface that "the standard says",
> I doubt that a fair use defense would hold up. And such copying of
> pretty much the whole chapter is what I was talking about.

With an assist from John Reid, the wording of three of those
intrinsic function definitons is essentially mine. If I were
writing a book, I couldn't even quote my own work?

Richard Maine

unread,
Aug 18, 2008, 10:42:49 PM8/18/08
to
James Giles <james...@worldnet.att.net> wrote:

> Richard Maine wrote:

> > If you quoted the entire chapter on intrinsics (or for that matter,
> > the entire standard) with just the preface that "the standard says",
> > I doubt that a fair use defense would hold up. And such copying of
> > pretty much the whole chapter is what I was talking about.
>
> With an assist from John Reid, the wording of three of those
> intrinsic function definitons is essentially mine. If I were
> writing a book, I couldn't even quote my own work?

As I said, I'm no expert on the matter. I assume the question is just
rhetorical, because if you are actually asking, you are surely asking in
the wrong place. My opinion is neither expert legal advice (I'd probably
be in trouble with other laws if I pretended to offer that) nor
indicative of ISO's position. You'd need to ask a lawyer for one, and
ISO for the other. Neither of those is likely to give you an answer
here. As has been said on other matters, those who do claim to have
definitive answers on such things are the ones least likely to actually
know the real answers.

Ron Ford

unread,
Aug 18, 2008, 11:28:53 PM8/18/08
to
On Mon, 18 Aug 2008 13:35:39 -0700 (PDT), e p chandler posted:

> On Aug 18, 12:05 am, Ron Ford <r...@example.invalid> wrote:


>> ! seed random num generator
>> CALL SYSTEM_CLOCK(COUNT=clock)
>> seed = clock + 37
>> CALL RANDOM_SEED(PUT = seed)
>
> This does not compile with G95. Perhaps it is a compiler error. G95
> does not document this call but
> under GNU fortran,
>
> integer size
> ....
> CALL RANDOM_SEED(size)
>
> size is the minimum size of seed which is an array.
>
> see http://gcc.gnu.org/onlinedocs/gfortran/RANDOM_005fSEED.html
>
> Perhaps the interface is different for your compiler (Are you using
> Silverfrost's FTN95?)

Silverfrost does accept the above (is that non-standard?). I'd been
wracking my brain how to get an appropriate seed for gfortran. I'll try
the code at that link and see what happens.

As a work-around, I just created an integer array of one hundred dimensions
and assigned the clock to each element. What I have now compiles with
gfortran.


>> ! sort into bins
>> D=0
>> maxput = bins
>> do ii =1, bins
>
> ** running with array bounds (-fbounds-check for g95) checking turned
> on reveals this error:
> ** Should be **
> do ii =1, trials

Thanks, elliot. Nice comments from a compiler, I'm impressed, in
particular because that gets rid of the runtime.


> [flame on]
>
> For some reason I find it easier to debug a problem program from a
> printed listing instead of looking at it on the screen. It's easy to
> mark off flow of control, underline used variables, check off
> interfaces, etc.. YMMV.
>
> [flame off]
>
> - e

I used to do that more, and I was a more-careful coder because of it. With
bins at 50 and trials at ten thousand, this is the output:

1 0.887136
2 0.989418
3 0.847468
4 0.820818
5 0.369006
10000.0 0 0 0 0
0
157 414 579
750 808 829 840 759 685
607 507 469
438 364 308 258 208 159
140 118 102
97 70 57 36 32 32
26 23 18
13 17 17 6 9 11
10 4 4
1 1 2 4 1 10
95th % is 9500.00
0
0
0
0
0
157
571
1150
1900
2708
3537
4377
5136
5821
6428
6935
7404
7842
8206
8514
8772
8980
9139
9279
9397
9499
9596
95th percentile was in bin number 27

Ron Ford

unread,
Aug 18, 2008, 11:39:05 PM8/18/08
to
On Mon, 18 Aug 2008 14:18:02 -0700, Richard Maine posted:

>> This does not compile with G95. Perhaps it is a compiler error. G95
>> does not document this call
>
> No, it is not a compiler error. It is an error in the code. The reason
> for the error was discussed previously in this thread - that argument is
> required to be a rank 1 array rather than a scalar.

If we're thinking of the same instance, it was actually gfortran that
insisted that the seed be an array. I'd never made an array out of the
clock before, but this seems to fit the bill for gfortran:

integer, dimension(100)::E


! seed random num generator

do i=1,100
CALL SYSTEM_CLOCK(COUNT=clock)
E(i)=clock

end do
!print *, E

CALL RANDOM_SEED(PUT = E)

Apparently, gfortran wants a vector larger than six integers, too. Do I
take it that a rank zero seed is non-standard?


>
> Admitedly, the message I get from g95 is a bit confusing, but that's
> another matter. I can sort of see how that message comes about.
> (The second message. The first one I got was from the invalid characters
> in the continuation line. Hex CA? Anyway, I replaced those by blanks.)

Hex CA sounds like a leftover from the zipcode sort. (?)


--
Wealth - any income that is at least one hundred dollars more a year than
the income of one's wife's sister's husband. 6
H. L. Mencken

Ron Ford

unread,
Aug 19, 2008, 12:09:21 AM8/19/08
to
On Mon, 18 Aug 2008 16:28:16 -0700 (PDT), e p chandler posted:

I took a gander at random_number in 04-007.pdf, and I didn't follow the
point it was making about successive puts and gets. When I get what I just
put, it gives me the same value. I was hoping when I "got" the seed, it
would give me a hint as to how my compiler represents it.

e p chandler

unread,
Aug 19, 2008, 12:24:41 AM8/19/08
to

Sorry, no expreience with this yet. But I did write a little test
program

x.f95:

implicit none

integer array_size

call random_seed(array_size)
print *,'minimum array size =',array_size

end

This reports the minimum size of the seed array.

C:\Users\epc\temp>g95 x.f95

C:\Users\epc\temp>a
minimum array size = 4

C:\Users\epc\temp>gfortran x.f95

C:\Users\epc\temp>a
minimum array size = 8

C:\Users\epc\temp>g95 -v
...
gcc version 4.1.2 (g95 0.92!) Aug 14 2008

C:\Users\epc\temp>gfortran -v
...
gcc version 4.4.0 20080603 (experimental) [trunk revision 136333]
(GCC)

C:\Users\epc\temp>

[Time to pull out my corrected copy of MR&C.]

- e


Richard Maine

unread,
Aug 19, 2008, 12:33:22 AM8/19/08
to
Ron Ford <r...@example.invalid> wrote:

> On Mon, 18 Aug 2008 14:18:02 -0700, Richard Maine posted:
>
> >> This does not compile with G95. Perhaps it is a compiler error. G95
> >> does not document this call
> >
> > No, it is not a compiler error. It is an error in the code. The reason
> > for the error was discussed previously in this thread - that argument is
> > required to be a rank 1 array rather than a scalar.
>
> If we're thinking of the same instance, it was actually gfortran that
> insisted that the seed be an array.

No, I wasn't talking about gfortran.

I was referring to exactly what was in the post I replied to (while I
don't always do that, it is reasonably common), where e p chandler said
that it didn't compile with g95. When I comented about the messages I
also got from g95, I also meant exactly what I said.

Since the standard requires the seed be an array, it would not be at all
surprising if other compilers also stuck to the standard.

Steven G. Kargl

unread,
Aug 19, 2008, 1:14:52 AM8/19/08
to
In article <yd5my4btmwut$.dlg@example.invalid>,

Ron Ford <r...@example.invalid> writes:
>
> Silverfrost does accept the above (is that non-standard?). I'd been
> wracking my brain how to get an appropriate seed for gfortran. I'll try
> the code at that link and see what happens.

Have you read the gfortran manual? The description of random_seed()
gives you a portable method.


SUBROUTINE init_random_seed()
INTEGER :: i, n, clock
INTEGER, DIMENSION(:), ALLOCATABLE :: seed

CALL RANDOM_SEED(size = n)
ALLOCATE(seed(n))

CALL SYSTEM_CLOCK(COUNT=clock)

seed = clock + 37 * (/ (i - 1, i = 1, n) /)
CALL RANDOM_SEED(PUT = seed)

DEALLOCATE(seed)
END SUBROUTINE

Of course, I suspect you've read either the manual or a post
based on the manual because 37 appears elsewhere in this
thread (e.g., your initial post). Yes, I know the origin
of 37 above.

> As a work-around, I just created an integer array of one hundred dimensions
> and assigned the clock to each element.

Why? If you actually read the Standard or the gfortran manual,
you'll see that you can query random_seed() for the size of
the array.

--
Steve
http://troutmask.apl.washington.edu/~kargl/

Kurt Kallblad

unread,
Aug 19, 2008, 4:09:05 AM8/19/08
to

"Richard Maine" <nos...@see.signature> wrote in message
news:1ilw0ka.ztp2fc1p84pboN%nos...@see.signature...
>
> <snip>

>
> Since the standard requires the seed be an array, it would not
> be at all
> surprising if other compilers also stuck to the standard.
>
> --
> Richard Maine | Good judgement comes from
> experience;
> email: last name at domain . net | experience comes from bad
> judgement.
> domain: summertriangle | -- Mark Twain

The old CVF6.6 gives

C:\Kurt\CLF\random_seed.f90(23) : Error: The data types of the
argument(s) are invalid. [RANDOM_SEED]
call random_seed(put = x)
-------------------------^
C:\Kurt\CLF\random_seed.f90(23) : Error: One-dimensional
array-valued arguments are required in this context.
[RANDOM_SEED]
call random_seed(put = x)

Kurt

Kurt Kallblad

unread,
Aug 19, 2008, 5:10:11 AM8/19/08
to
"e p chandler" <ep...@juno.com> wrote in message
news:52f53a9f-78cf-4a03...@s50g2000hsb.googlegroups.com...

>On Aug 19, 12:09 am, Ron Ford <r...@example.invalid> wrote:
>> On Mon, 18 Aug 2008 16:28:16 -0700 (PDT), e p chandler posted:
>
><snip>

>>
>> I took a gander at random_number in 04-007.pdf, and I didn't
>> follow the
>> point it was making about successive puts and gets. When I get
>> what I just
>> put, it gives me the same value. I was hoping when I "got" the
>> seed, it
>> would give me a hint as to how my compiler represents it.
>
>Sorry, no expreience with this yet. But I did write a little
>test
>program
>
><snip>
>

The random sequence given by the random_number intrinsic is (for
a given platform) always the same when started from a given seed
(at least in the compilers I have used). Furthermore, the seed is
well defind before each random number if given.

Thus, the "got" is usful if, especially when testing a new
routine, one want to restart the random sequence from a given
point as illustrated below.

Kurt

PS. I don't like ;-langauges but short news.
-------------------------------
program test

real :: x(6); integer :: i
integer, allocatable :: old1(:), old2(:)

! get the needed size
call random_seed(size = i)
allocate(old1(i)); allocate(old2(i))

! get the compiler dependent seed
call random_seed()
! save the seed for later use
call random_seed(get = old1)
! get the first 6 compiler & seed dependent random numbers
call random_number(x); write(*, '(6f10.5)') x

! save the current seed for later use
call random_seed(get = old2)
call random_number(x); write(*, '(6f10.5)') x

! repeate the last 6 numbers once again
call random_seed(put = old2)
call random_number(x); write(*, '(/6f10.5)') x

! repeate the first 6 numbers once again
call random_seed(put = old1)
call random_number(x); write(*, '(/6f10.5)') x

end program test

-------------------------------

Ron Shepard

unread,
Aug 19, 2008, 10:28:11 AM8/19/08
to
In article
<bjmqk.168070$102....@bgtnsc05-news.ops.worldnet.att.net>,
"James Giles" <james...@worldnet.att.net> wrote:

> Would INCITS or ISO really enforce their copyright if a compiler
> manual used some content like: "the standard says: '...'? Wouldn't
> that constitute fair use?

I used a compiler once that had this kind of documentation. It was
the first version of the FPS-164 compiler in 1981 based on the f77
standard. Large sections of it were taken directly from the ANSI
standard. ANSI, or someone, complained, and the subsequent
documentation came in two parts, the ANSI standard document along
with a separate, smaller, user's guide.

$.02 -Ron Shepard

glen herrmannsfeldt

unread,
Aug 19, 2008, 3:25:12 PM8/19/08
to
Steven G. Kargl wrote:
(snip)

> Have you read the gfortran manual? The description
> of random_seed() gives you a portable method.

> SUBROUTINE init_random_seed()
> INTEGER :: i, n, clock
> INTEGER, DIMENSION(:), ALLOCATABLE :: seed
>
> CALL RANDOM_SEED(size = n)
> ALLOCATE(seed(n))
>
> CALL SYSTEM_CLOCK(COUNT=clock)
>
> seed = clock + 37 * (/ (i - 1, i = 1, n) /)
> CALL RANDOM_SEED(PUT = seed)
>
> DEALLOCATE(seed)
> END SUBROUTINE

> Of course, I suspect you've read either the manual or a post
> based on the manual because 37 appears elsewhere in this
> thread (e.g., your initial post). Yes, I know the origin
> of 37 above.

Well, 37 is also a convenient prime number. It is hard
to say what might make a high quality seed, though.

>>As a work-around, I just created an integer array of one hundred dimensions
>>and assigned the clock to each element.

Dynamically allocating to the exact size seems to be
the suggested way, though it is a fair amount of work.

100 sounds a little small, but probably big enough for most.
I might have suggested 1000 or maybe 2048. If 2048 bits
(64 32 bit INTEGERs) is enough for cryptographic applications,
it is likely enough for other random number applications.

> Why? If you actually read the Standard or the gfortran manual,
> you'll see that you can query random_seed() for the size of
> the array.

The standard also says that the array size must be greater
than or equal to the returned size.

-- glen

Steven G. Kargl

unread,
Aug 19, 2008, 2:50:49 PM8/19/08
to
In article <KoOdnfCeL85xkjbV...@comcast.com>,

glen herrmannsfeldt <g...@ugcs.caltech.edu> writes:
> Steven G. Kargl wrote:
> (snip)
>
>> Have you read the gfortran manual? The description
>> of random_seed() gives you a portable method.
>
>> SUBROUTINE init_random_seed()
>> INTEGER :: i, n, clock
>> INTEGER, DIMENSION(:), ALLOCATABLE :: seed
>> CALL RANDOM_SEED(size = n)
>> ALLOCATE(seed(n))
>> CALL SYSTEM_CLOCK(COUNT=clock)
>> seed = clock + 37 * (/ (i - 1, i = 1, n) /)
>> CALL RANDOM_SEED(PUT = seed)
>> DEALLOCATE(seed)
>> END SUBROUTINE
>
>> Of course, I suspect you've read either the manual or a post
>> based on the manual because 37 appears elsewhere in this
>> thread (e.g., your initial post). Yes, I know the origin
>> of 37 above.
>
> Well, 37 is also a convenient prime number. It is hard
> to say what might make a high quality seed, though.

I never made any claims about the quality of the seeds. The point
is OP must have seen either the gfortran manual or a post involving
the above code because of the number 37. Thus, OP should have
been able to avoid his problems with (1) PUT taking a rank 1 array
and (2) the required size of that rank 1 array.


>
>>>As a work-around, I just created an integer array of one hundred dimensions
>>>and assigned the clock to each element.
>
> Dynamically allocating to the exact size seems to be
> the suggested way, though it is a fair amount of work.

3 lines of code is a fair amount of work?

> 100 sounds a little small, but probably big enough for most.
> I might have suggested 1000 or maybe 2048. If 2048 bits
> (64 32 bit INTEGERs) is enough for cryptographic applications,
> it is likely enough for other random number applications.
>
>> Why? If you actually read the Standard or the gfortran manual,
>> you'll see that you can query random_seed() for the size of
>> the array.
>
> The standard also says that the array size must be greater
> than or equal to the returned size.

It's a rank 1 array. It the array is larger than needed, those
elements are simply ignored.

--
Steve
http://troutmask.apl.washington.edu/~kargl/

glen herrmannsfeldt

unread,
Aug 19, 2008, 4:47:05 PM8/19/08
to
Steven G. Kargl wrote:
(snip, I wrote)

>>Dynamically allocating to the exact size seems to be
>>the suggested way, though it is a fair amount of work.

> 3 lines of code is a fair amount of work?

No, the implementation of allocatable arrays is usually
much more work than implementing static arrays.

The number of lines of Fortran code is often not
directly related to the work required to implement
them.

-- glen

Richard Maine

unread,
Aug 19, 2008, 4:22:03 PM8/19/08
to
glen herrmannsfeldt <g...@ugcs.caltech.edu> wrote:

That seems rather irrelevant to any practical consideration in the
current discussion to me.

The extra work involved in making the compiler support allocatable
arrays is already done and has nothing in particular to do with this one
isolated case. Even if one wants to argue for using f77, you aren't
going to find a compiler that implements the f90 random number
intrinsics, complete with keyword arguments, but doesn't already support
allocatables. (And this is a trivially simple simple case of
allocatables - not the TR stuff). So that line of argument is a
non-starter.

And the extra runtime work done by the application is going to be
completely trivial. Allocating an array for the seed of a random number
generator is the kind of thing done once in an application, or perhaps a
few times if you have a few independent random sequences. If you are
using it in such a way that the work is worth even the most cursory of
attention, then you are doing something very wrong more than just
allocating the array.

The only measure of relevance here is the one that Steven alludes to -
the 3 lines (or so) of code. I'd further maintain that any "reasonable"
alternative would actually be more work by such relevant measures. I
don't consider it "reasonable" to just declare a large static array and
hope that it is large enough. Rationalizations of why some particular
size ought to be large enough don't interest me (and besides, it would
be a lot more work to deduce a "safe" size than to just do the
allocation). I don't do fixed-size arrays without checking somewhere
that the size is adequate. By the time you throw in such a check, along
with appropriate code for when the check fails (even if only to print a
message and stop), you've just added quite a lot more code than it would
have taken to just do the allocation.

In short (I know, I blew that already), no allocating the array is not a
"fair amount of work" by any relevant measure. I don't buy it. I'd say
that avoiding the allocation was significantly more work.

Ron Ford

unread,
Aug 19, 2008, 7:46:58 PM8/19/08
to
On Tue, 19 Aug 2008 10:09:05 +0200, Kurt Kallblad posted:

> "Richard Maine" <nos...@see.signature> wrote in message
> news:1ilw0ka.ztp2fc1p84pboN%nos...@see.signature...
>>
>> <snip>
>>
>> Since the standard requires the seed be an array, it would not
>> be at all
>> surprising if other compilers also stuck to the standard.

Richard,

Is silverfrost conforming with the output below?

> The old CVF6.6 gives
>
> C:\Kurt\CLF\random_seed.f90(23) : Error: The data types of the
> argument(s) are invalid. [RANDOM_SEED]
> call random_seed(put = x)
> -------------------------^
> C:\Kurt\CLF\random_seed.f90(23) : Error: One-dimensional
> array-valued arguments are required in this context.
> [RANDOM_SEED]
> call random_seed(put = x)
>
> Kurt

Thanks for your reply, kurt. I'm gonna fire up the source you posted as a
reply to elliot next.

Steve posted the source that also was at the link that elliot posted:

http://gcc.gnu.org/onlinedocs/gfortran/RANDOM_005fSEED.html

Now that I've gotten my head around the whole business of creating an array
to seed an RNG, it makes a lot more sense. In C, we had a single word to
seed, and until you *have* randomness, you can't well create it in the form
of guessing how big an appropriate array should be.

That, of course, is the role of the size argument. Furthermore, you don't
just have one call to random_seed to do it properly. You need at least 2.

program seed1

call init_random_seed


contains


SUBROUTINE init_random_seed()
INTEGER :: i, n, clock
INTEGER, DIMENSION(:), ALLOCATABLE :: seed

CALL RANDOM_SEED(size = n)
print *, "n=", n
ALLOCATE(seed(n))

CALL SYSTEM_CLOCK(COUNT=clock)
print *, "clock=", clock

seed = clock + 37 * (/ (i - 1, i = 1, n) /)
CALL RANDOM_SEED(PUT = seed)

print *, "seed= ", seed

DEALLOCATE(seed)
END SUBROUTINE

end program

n= 1
clock= 17164366
seed= 17164366

Press RETURN to close window . . .

--
War will never cease until babies begin to come into the world with larger
cerebrums and smaller adrenal glands. 2
H. L. Mencken

Richard Maine

unread,
Aug 19, 2008, 8:24:16 PM8/19/08
to
Ah well. I was trying to ignore this whole thread, since I consider the
question of "how does one sort things" to be off topic. But I see I
haven't suceeded very well, I suppose partly because the thread has
actually drifted from an off-topic start into material that is actually
about Fortran. That's the opposite of usual thread drift.

Ron Ford <r...@example.invalid> wrote:

> Is silverfrost conforming with the output below?

I didn't notice anything "below" that appeared to be from Silverfrost.
But speculating on what you probably meant based on earlier posts, I
thought that someone else had already answered, and correctly so.

In 90% of the cases (a statistic I just made up, but I'm confident is as
good as most others) questions about whether a compiler is conforming
all have the same answer:

The standard is primarily about whether user code conforms - not about
whether compilers conform. The standard says so itself (in slightly
different, but simillar terms).

There is material about whether compilers conform, but by far the
biggest part of that is just about doing the right thing to user code
that conforms. There are a few other things, but they pale in
comparison.

With limitted exceptions (this isn't one of them, and I don't feel it
worth elaborating) a compiler can do anything with code that doesn't
conform. That doesn't say anything about the conformance of the
compiler. The code in question doesn't conform (at least if I've guessed
correctly about what you were asking). Q.E.D.

Steven G. Kargl

unread,
Aug 19, 2008, 8:32:53 PM8/19/08
to
In article <dpzmlr04...@example.invalid>,

Ron Ford <r...@example.invalid> writes:
> On Tue, 19 Aug 2008 10:09:05 +0200, Kurt Kallblad posted:
>
>> "Richard Maine" <nos...@see.signature> wrote in message
>> news:1ilw0ka.ztp2fc1p84pboN%nos...@see.signature...
>>>
>>> <snip>
>>>
>>> Since the standard requires the seed be an array, it would not
>>> be at all
>>> surprising if other compilers also stuck to the standard.
>
> Richard,
>
> Is silverfrost conforming with the output below?

I'm not Richard, but yes, it's conforming.

> n= 1
> clock= 17164366
> seed= 17164366
>

Apparently, Silverfrost's implementation only requires
a single seed.

--
Steve
http://troutmask.apl.washington.edu/~kargl/

Richard Maine

unread,
Aug 19, 2008, 8:48:26 PM8/19/08
to

Ah, and your answer is better than mine. I appear to have been answering
a different question. If this was the question, then yes, it is fine. An
array of size 1 is still an array. There is nothing restricting the size
from being 1.

Steven G. Kargl

unread,
Aug 19, 2008, 9:10:58 PM8/19/08
to
In article <1ilxlba.1fqzxdojwzw8cN%nos...@see.signature>,

nos...@see.signature (Richard Maine) writes:
> Steven G. Kargl <ka...@troutmask.apl.washington.edu> wrote:
>
>> In article <dpzmlr04...@example.invalid>,
>> Ron Ford <r...@example.invalid> writes:
>
>> > Is silverfrost conforming with the output below?
>>
>> I'm not Richard, but yes, it's conforming.
>>
>> > n= 1
>> > clock= 17164366
>> > seed= 17164366
>> >
>>
>> Apparently, Silverfrost's implementation only requires
>> a single seed.
>
> Ah, and your answer is better than mine. I appear to have been answering
> a different question. If this was the question, then yes, it is fine. An
> array of size 1 is still an array. There is nothing restricting the size
> from being 1.
>

Your previous response is quite understandable. Ron Ford chose
to respond to 3 different posts in one go.

--
Steve
http://troutmask.apl.washington.edu/~kargl/

Ron Ford

unread,
Aug 19, 2008, 10:47:30 PM8/19/08
to
On Tue, 19 Aug 2008 09:28:11 -0500, Ron Shepard posted:

Secret documentation is just becoming hilarious. My sysadmin buddy has
been telling me lately of some of the nuances in copyright and patent.

Novelle just embraced open source.

Please contact ANSI to pay for a gideon's bible.

Ron Ford

unread,
Aug 19, 2008, 11:15:00 PM8/19/08
to
On Tue, 19 Aug 2008 13:22:03 -0700, Richard Maine posted:

> In short (I know, I blew that already), no allocating the array is not a
> "fair amount of work" by any relevant measure. I don't buy it. I'd say
> that avoiding the allocation was significantly more work.

As the person who demonstrably *did* the work to understand fortran's
implementation of random_seed, I can assure you that it was work and that
it is not yet behind me. No one is selling.

Some of the work was in the metasynatctics. With put, get, and seed as
arguments, mutually exclusive, the subtulties are many.

What does a person need to do, in terms of a random_seed call, to invoke
the Mersenne twister?
--
We must be willing to pay a price for freedom. 4
H. L. Mencken

Steven G. Kargl

unread,
Aug 20, 2008, 12:52:18 AM8/20/08
to
In article <ubfx6ay2wghf$.dlg@example.invalid>,

Ron Ford <r...@example.invalid> writes:
> On Tue, 19 Aug 2008 13:22:03 -0700, Richard Maine posted:
>
>> In short (I know, I blew that already), no allocating the array is not a
>> "fair amount of work" by any relevant measure. I don't buy it. I'd say
>> that avoiding the allocation was significantly more work.
>
> As the person who demonstrably *did* the work to understand fortran's
> implementation of random_seed, I can assure you that it was work and that
> it is not yet behind me. No one is selling.

well, if you actually read the documentation that more than one person
pointed you to, understanding randomw_seed() would have been trivial.


> What does a person need to do, in terms of a random_seed call, to invoke
> the Mersenne twister?

Well, you'ld need a compiler that used MT as it PRNG. Of course,
you could overload random_seed() for MT if your compiler did not
use MT.

For the record, I ripped gfortran's implementation of MT out years
ago. For the reasons for its removal see

http://gcc.gnu.org/viewcvs/trunk/libgfortran/intrinsics/random.c?revision=138078&view=markup

--
Steve
http://troutmask.apl.washington.edu/~kargl/

Ron Ford

unread,
Aug 20, 2008, 6:45:13 PM8/20/08
to
On Tue, 19 Aug 2008 11:10:11 +0200, Kurt Kallblad posted:

> The random sequence given by the random_number intrinsic is (for
> a given platform) always the same when started from a given seed
> (at least in the compilers I have used). Furthermore, the seed is
> well defind before each random number if given.
>
> Thus, the "got" is usful if, especially when testing a new
> routine, one want to restart the random sequence from a given
> point as illustrated below.
>
> Kurt

I modified your source slightly and tried it on silverfrost and gfortran.
What follows is the source, the output for gfortran and then the output for
silverfrost.

program test

real :: x(6); integer :: i
integer, allocatable :: old1(:), old2(:)

! get the needed size
call random_seed(size = i)

print *, "size=", i
allocate(old1(i)); allocate(old2(i))

! get the compiler dependent seed
call random_seed()

! save the seed for later use
call random_seed(get = old1)

print *, "old1=", old1


! get the first 6 compiler & seed dependent random numbers
call random_number(x); write(*, '(6f10.5)') x

! save the current seed for later use
call random_seed(get = old2)
call random_number(x); write(*, '(6f10.5)') x

! repeate the last 6 numbers once again
call random_seed(put = old2)
call random_number(x); write(*, '(/6f10.5)') x

! repeate the first 6 numbers once again
call random_seed(put = old1)
call random_number(x); write(*, '(/6f10.5)') x

end program test
! gfortran -o prob freeformat35.f95


F:\gfortran\source>gfortran -o prob freeformat35.f95

F:\gfortran\source>prob
size= 8
old1= 437395160 1404128605 572505362 -1187264075 454383258
525702629
973594203 1758310677
0.99756 0.56682 0.96592 0.74793 0.36739 0.48064
0.07375 0.00536 0.34708 0.34224 0.21795 0.13316

0.07375 0.00536 0.34708 0.34224 0.21795 0.13316

0.99756 0.56682 0.96592 0.74793 0.36739 0.48064

F:\gfortran\source>prob
size= 8
old1= 437395160 1404128605 572505362 -1187264075 454383258
525702629
973594203 1758310677
0.99756 0.56682 0.96592 0.74793 0.36739 0.48064
0.07375 0.00536 0.34708 0.34224 0.21795 0.13316

0.07375 0.00536 0.34708 0.34224 0.21795 0.13316

0.99756 0.56682 0.96592 0.74793 0.36739 0.48064

F:\gfortran\source>

This shows that the default seed for gfortran does not vary.

For silverfrost:

size= 1
old1= -2074922541
0.38205 0.70350 0.95220 0.53890 0.66540 0.52014
0.93099 0.97529 0.78729 0.55684 0.02209 0.80699

0.38205 0.70350 0.95220 0.53890 0.66540 0.52014

0.38205 0.70350 0.95220 0.53890 0.66540 0.52014

Press RETURN to close window . . .

size= 1
old1= 609462396
0.67998 0.13839 0.54638 0.19263 0.37788 0.81034
0.80125 0.11799 0.45571 0.42196 0.89159 0.56965

0.67998 0.13839 0.54638 0.19263 0.37788 0.81034

0.67998 0.13839 0.54638 0.19263 0.37788 0.81034

Press RETURN to close window . . .

This shows that the default behavior for silverfrost *does* vary the seed.
Both compilers show the get and put behavior that the standard requires.

Ron Ford

unread,
Aug 20, 2008, 9:20:02 PM8/20/08
to
On Wed, 20 Aug 2008 04:52:18 +0000 (UTC), Steven G. Kargl posted:

> In article <ubfx6ay2wghf$.dlg@example.invalid>,
> Ron Ford <r...@example.invalid> writes:
>> On Tue, 19 Aug 2008 13:22:03 -0700, Richard Maine posted:
>>
>>> In short (I know, I blew that already), no allocating the array is not a
>>> "fair amount of work" by any relevant measure. I don't buy it. I'd say
>>> that avoiding the allocation was significantly more work.
>>
>> As the person who demonstrably *did* the work to understand fortran's
>> implementation of random_seed, I can assure you that it was work and that
>> it is not yet behind me. No one is selling.
>
> well, if you actually read the documentation that more than one person
> pointed you to, understanding randomw_seed() would have been trivial.

You must play different trivia games than do I. I've never hit the science
and math category and been asked "how many random_seed calls do you have to
make for an effective fortran pseudorandom implementation?"


>
>
>> What does a person need to do, in terms of a random_seed call, to invoke
>> the Mersenne twister?
>
> Well, you'ld need a compiler that used MT as it PRNG. Of course,
> you could overload random_seed() for MT if your compiler did not
> use MT.
>
> For the record, I ripped gfortran's implementation of MT out years
> ago. For the reasons for its removal see
>
> http://gcc.gnu.org/viewcvs/trunk/libgfortran/intrinsics/random.c?revision=138078&view=markup

I see in here some of the same information as was in gfortran.pdf, to wit:

libgfortran currently uses George Marsaglia's KISS (Keep It Simple
Stupid)
random number generator. This PRNG combines:

(1) The congruential generator x(n)=69069*x(n-1)+1327217885 with a
period
of 2^32,
(2) A 3-shift shift-register generator with a period of 2^32-1,
(3) Two 16-bit multiply-with-carry generators with a period of
597273182964842497 > 2^59.

The overall period exceeds 2^123.

What I don't see here is what makes fortran's implementation of
pseudorandoms different than C's, at least gcc's. Fortran's is, by
reputation, good, while C's--at least my vanilla rand call with
Microsoft--was famously clunky. Is it all the same stuff stuff with
different window-dressing?
--
What men value in this world is not rights but privileges. 7
H. L. Mencken

Steven G. Kargl

unread,
Aug 20, 2008, 11:14:34 PM8/20/08
to
In article <jvt0njth...@example.invalid>,

Ron Ford <r...@example.invalid> writes:
> On Wed, 20 Aug 2008 04:52:18 +0000 (UTC), Steven G. Kargl posted:
>
>> In article <ubfx6ay2wghf$.dlg@example.invalid>,
>> Ron Ford <r...@example.invalid> writes:
>
>>> What does a person need to do, in terms of a random_seed call, to invoke
>>> the Mersenne twister?
>>
>> Well, you'ld need a compiler that used MT as it PRNG. Of course,
>> you could overload random_seed() for MT if your compiler did not
>> use MT.
>>
>> For the record, I ripped gfortran's implementation of MT out years
>> ago. For the reasons for its removal see
>>
>> http://gcc.gnu.org/viewcvs/trunk/libgfortran/intrinsics/random.c?revision=138078&view=markup
>
> I see in here some of the same information as was in gfortran.pdf,

The gfortran developers don't expect an ordinary user to read the
source code. Heck, at least one former gfortran developer has come
to realize that ordinary user typically don't read the manual, too. :)

> to wit:
>
> libgfortran currently uses George Marsaglia's KISS (Keep It Simple
> Stupid) random number generator. This PRNG combines:
>
> (1) The congruential generator x(n)=69069*x(n-1)+1327217885 with a
> period of 2^32,
> (2) A 3-shift shift-register generator with a period of 2^32-1,
> (3) Two 16-bit multiply-with-carry generators with a period of
> 597273182964842497 > 2^59.
>
> The overall period exceeds 2^123.

Actually, gfortran's random_number() will use up to 3 independent
KISS generators. It also will scramble user supplied seeds.

> What I don't see here is what makes fortran's implementation of
> pseudorandoms different than C's, at least gcc's. Fortran's is, by
> reputation, good, while C's--at least my vanilla rand call with
> Microsoft--was famously clunky. Is it all the same stuff stuff with
> different window-dressing?

Neither C nor Fortran specify the underlying PRNG algorithm. One
should not be surprised to find that the PRNG is in fact the same
on a given system for both languages. With gfortran, you can call
use RAND() or RAN() if you want a rather poor PRNG.

One other item to consider is that gfortran provides an implementation
of a PRNG to avoid potential problems with PRNGs on various operating
systems. If a bug is reported against random_number(), the gfortran
developers don't have to ask about the OS and the quality of its
PRNG.

PS: If I really want a random stream of bits, I do

program testing
implicit none
integer i,j
open(unit=10, file='/dev/random', access='stream', status='old')
do i = 1, 10
read(10) j
print *, j
end do
end program testing

--
Steve
http://troutmask.apl.washington.edu/~kargl/

Ron Ford

unread,
Aug 21, 2008, 7:21:06 PM8/21/08
to
On Fri, 15 Aug 2008 17:08:32 -0600, Ron Ford posted:

> I have a program whose output is an array of integers.
>
> implicit none
>
> integer, parameter:: sides = 6
> integer, parameter:: trials = 100
>
> integer, dimension(sides)::A
> integer, dimension(trials)::C
> integer:: b, clock, seed, ii, i, counter, tab, penultimate
> real:: harvest


>
>
> ! seed random num generator

> CALL SYSTEM_CLOCK(COUNT=clock)


> seed = clock + 37

> CALL RANDOM_SEED(PUT = seed)
>
> ! prime the pump
> call random_number(harvest)
> b = 3 + nint(10*harvest)
> do i=1,b
> call random_number(harvest)
> print *, i, harvest
> end do
>
>
> ! main control
>
> ! outer loop is the number of trials
> C = 0
> do ii=1,trials
> tab = 0
> A = 0
> ! inner loop rolls die till all values attained
> do
> call random_number(harvest)
> b = ceiling(harvest*real(sides))
> A(b) = A(b) + 1
> ! count until all outcomes non-zero
> counter = 0
> do i = 1, sides
> if (real(A(i)) > .5) then
> counter = counter + 1
> end if
> end do
>
> tab = tab + 1
> !print *, "counter= ", counter
> !print *, A
> !print *, "tab= ", tab
> if (counter == sides) then
> penultimate = tab -1
> C(ii) = penultimate
> end if
>
> if (counter == sides) exit
> end do !inner
> end do !outer
>
> print *, "C is", C
> end program
>
> 1 0.176787
> 2 0.458366
> 3 0.228078
> 4 0.711985
> C is 14 26 5 27 12
> 12
> 9 15 18
> 19 12 39 10 11 23
> 9 16 12
> 7 11 8 17 18 21
> 18 15 12
> 14 18 14 29 26 16
> 11 21 18
> 6 7 20 11 12 15
> 10 7 16
> 15 8 8 5 11 16
> 20 11 13
> 15 13 9 20 13 21
> 23 8 9
> 5 12 23 13 12 7
> 11 10 11
> 11 10 14 14 13 12
> 15 7 12
> 9 7 7 20 9 6
> 9 6 14
> 11 11 14 12 6 16
> 13 8 14
> 26


>
> Press RETURN to close window . . .
>
>

> What is a slick way to order C?
>


I've got all the ingredients to order C now; I just have to put it all
together. I modified MR&C's zipcode sort as follows:

module sort3
implicit none
private
public :: selection_sort
integer, parameter :: string_length = 30
type, public :: address

integer :: zip_code
end type address
contains
recursive subroutine selection_sort (array_arg)
type (address), dimension(:), intent (inout) &
:: array_arg
integer :: current_size
integer :: big
current_size = size (array_arg)
if (current_size > 0) then
big = maxloc(array_arg(:)%zip_code, dim=1)
call swap (big, current_size)
call selection_sort (array_arg(1: current_size -1))
end if
contains
subroutine swap(i,j)
integer, intent (in) :: i,j
type (address) :: temp
temp = array_arg(i)
array_arg(i) = array_arg(j)
array_arg(j) = temp
end subroutine swap
end subroutine selection_sort
end module sort3

program zippy
use sort3
implicit none
integer, parameter :: array_size = 1000
type (address), dimension (array_size) :: data_array
! temp variable - attempt to read here
type (address) :: input_record
integer :: i, n, ios

! no more statement 99
open(2, file='zips5.txt', status='old', iostat=ios)
if (ios /= 0) then
write(*,*) 'error ',ios,' opening input file'
stop
end if

open(3, file='output5.txt', status='replace')

n = 0

do
read (2, '(/i8)',iostat=ios) input_record
! quit on end of file or other errors
if (ios /= 0) exit
! feel free to check the value of ios and do something with it here

write (3, '(/i8)') input_record
! ok, we got another one
! who cares if we copy it to the array?
n = n + 1
data_array(n) = input_record
if (n == array_size) exit
end do

call selection_sort (data_array(1:n))

write(3, '(//a)') 'after sorting'
do i = 1, n
write (3, '(/i8)') data_array(i)
end do

end program zippy

zips5.txt reads:

47506

12345

54321

56789

76543

45678

98786

45678

87654
and output5.txt reads:

47506

12345

54321

56789

76543

45678

98786

45678

87654


after sorting

12345

45678

45678

47506

54321

56789

76543

87654

98786

So I've got a way to order integers but need to take out a bunch of spare
parts like:
integer, parameter :: string_length = 30
, and reading from a file, when I've already got the array.

As to topicality, there are many aspects of sorting that have nothing to do
with the fortran syntax. However, my game is not sufficiently advanced to
address them. I don't have a specific question right now but am fishing
for tips to achieve the objective of this thread.

> Thanks and cheers,
--
We are here and it is now. Further than that, all human knowledge is
moonshine. 3
H. L. Mencken

Ron Ford

unread,
Aug 21, 2008, 10:53:55 PM8/21/08
to
On Thu, 21 Aug 2008 17:21:06 -0600, Ron Ford posted:

> On Fri, 15 Aug 2008 17:08:32 -0600, Ron Ford posted:
>
>> I have a program whose output is an array of integers.

>> What is a slick way to order C?
> I've got all the ingredients to order C now; I just have to put it all
> together.

I'm stuck. This compiled until I tried to populate data_array with C.
I've got it remarked where this is:

module sort3
implicit none
private
public :: selection_sort

! type definition includes only an integer

use sort3


implicit none

integer, parameter:: sides = 8
! next 2 params must be equal


integer, parameter:: trials = 100

integer, parameter :: array_size = 100
integer, parameter:: bins = 50
integer, parameter:: percentile = 95


integer, dimension(sides)::A
integer, dimension(trials)::C

integer, dimension(bins)::D

!! from zippy


type (address), dimension (array_size) :: data_array


! integer :: i, n, ios
integer :: n


integer:: b, ii, i, counter, &
tab, maxput
real:: harvest, tot

! seed random num generator

CALL init_random_seed


! main control

if (counter == sides) then
C(ii) = tab
exit
end if


end do !inner
end do !outer

!print *, "C is", C

! sort into bins
D=0
maxput = bins

do ii =1, trials

if (C(ii) > maxput) then
C(ii) = maxput
end if

D(C(ii))=D(C(ii)) + 1

end do
tot= sum(D)
print *, "total trials=",tot

! determine 95th percentile
tot = 0.01 *trials*real(percentile)
print *, "95th % is", tot

counter = 0
do ii = 1, bins
counter = counter + D(ii)

if (real(counter).ge. tot) then
exit
end if
end do
print *, "95th percentile was in bin number", ii
!! paste in zippy

! fill data_array with C

do i=1,array_size
data_array%zip_code(i)=C(i) !this doesn't compile
end do


call selection_sort (data_array(1:n))

! output
write(*, '(//a)') 'after sorting'


do i = 1, n

write (*, '(/i8)') data_array(i)
end do


contains
SUBROUTINE init_random_seed()
INTEGER :: i, n, clock
INTEGER, DIMENSION(:), ALLOCATABLE :: seed

CALL RANDOM_SEED(size = n)
print *, "n=", n
ALLOCATE(seed(n))

CALL SYSTEM_CLOCK(COUNT=clock)
print *, "clock=", clock

seed = clock + 37 * (/ (i - 1, i = 1, n) /)
CALL RANDOM_SEED(PUT = seed)
print *, "seed= ", seed

DEALLOCATE(seed)
END SUBROUTINE


end program

Grateful for any tips,

e p chandler

unread,
Aug 22, 2008, 12:55:51 AM8/22/08
to

You have an array of elements of derived type (an array of
"structures" or an array of UDT = "user defined type"). Try

data_array(i)%zip_code=C(i)

instead.

> end do
>
The value of n is WILD (undefined) here.

> call selection_sort (data_array(1:n))

Try

call selection_sort (data_array(1:array_size))

instead

>
> ! output
>   write(*, '(//a)') 'after sorting'
>   do i = 1, n

Try
do i = 1, array_size

instead

>     write (*,  '(/i8)') data_array(i)
>   end do
>
> contains
>   SUBROUTINE init_random_seed()
>             INTEGER :: i, n, clock
>             INTEGER, DIMENSION(:), ALLOCATABLE :: seed
>
>             CALL RANDOM_SEED(size = n)
>             print *, "n=", n
>             ALLOCATE(seed(n))
>
>             CALL SYSTEM_CLOCK(COUNT=clock)
>             print *, "clock=", clock
>
>             seed = clock + 37 * (/ (i - 1, i = 1, n) /)
>             CALL RANDOM_SEED(PUT = seed)
>             print *, "seed= ", seed
>
>             DEALLOCATE(seed)
>           END SUBROUTINE
>
> end program
>
> Grateful for any tips,
>

With these changes, the program runs and produces output:

n= 4
clock= 0
seed= 0 37 74 111
1 0.9665488
2 0.8606074
3 0.6702148
4 0.1392371
5 0.5079537
6 0.33468115
7 0.8950028
8 0.3026237
9 0.00034302473
10 0.020161092
11 0.3237766
total trials= 100.
95th % is 95.
95th percentile was in bin number 34


after sorting

9

10

11

12

12

12

[data snipped]

29

30

33

34

35

35

37

39

41

This is G95 (MinGW-32) so your output will not necessarily be the
same. [IIRC you are using Silverfrost.]

HTH

I'm happy to see that you copied data from C(i) to data_array(i)
%zip_code [as corrected] element by element. I don't think you can
assume that an array of "structures" is necessarily the same as an
array if the "structure" (derived type) has only a single element.
[Sorry for the C-ism "structure" :-(.]
Perhaps I am wrong?

- e

Ron Ford

unread,
Aug 22, 2008, 4:29:17 AM8/22/08
to
On Thu, 21 Aug 2008 21:55:51 -0700 (PDT), e p chandler posted:

> This is G95 (MinGW-32) so your output will not necessarily be the
> same. [IIRC you are using Silverfrost.]

It's looking sexy, e:

! tja

use sort3
implicit none

integer, parameter:: trials = 50000
integer, parameter :: array_size = trials


integer, parameter:: bins = 50
integer, parameter:: percentile = 95


integer, dimension(sides)::A
integer, dimension(trials)::C
integer, dimension(bins)::D

type (address), dimension (array_size) :: data_array

integer:: b, ii, i, counter, &

CALL init_random_seed


! main control

! end main control


!print *, "C is", C

! sort into bins
D=0
maxput = bins
do ii =1, trials
if (C(ii) > maxput) then
C(ii) = maxput
end if
D(C(ii))=D(C(ii)) + 1
end do
tot= sum(D)
print *, "total trials=",tot

! determine 95th percentile
tot = 0.01 *trials*real(percentile)
print *, "95th % is", tot

counter = 0
do ii = 1, bins
counter = counter + D(ii)

if (real(counter).ge. tot) then
exit
end if
end do
print *, "95th percentile was in bin number", ii

! equate C with data_array

do i=1,array_size
data_array(i)%zip_code=C(i)
end do

call selection_sort (data_array(1:array_size))


! output
write(*, '(//a)') 'after sorting'

do i = 1, array_size
!write (*, '(/i8)') data_array(i)
end do

b = nint(tot)
print *, "b=", b
print *, data_array(b)

! end output


contains
SUBROUTINE init_random_seed()
INTEGER :: i, n, clock
INTEGER, DIMENSION(:), ALLOCATABLE :: seed

CALL RANDOM_SEED(size = n)
print *, "n=", n
ALLOCATE(seed(n))

CALL SYSTEM_CLOCK(COUNT=clock)
print *, "clock=", clock

seed = clock + 37 * (/ (i - 1, i = 1, n) /)
CALL RANDOM_SEED(PUT = seed)
print *, "seed= ", seed

DEALLOCATE(seed)
END SUBROUTINE
end program

n= 1
clock= 12506972
seed= 12506972
1 0.523549
2 0.547325
3 8.567679E-02
4 0.597775
5 0.978123
6 0.780981
7 0.827251
total trials= 50000.0
95th % is 47500.0
95th percentile was in bin number 38


after sorting
b= 47500
38

Press RETURN to close window . . .

>
> HTH

I'm happy that you helped.


>
> I'm happy to see that you copied data from C(i) to data_array(i)
> %zip_code [as corrected] element by element. I don't think you can
> assume that an array of "structures" is necessarily the same as an
> array if the "structure" (derived type) has only a single element.
> [Sorry for the C-ism "structure" :-(.]
> Perhaps I am wrong?

Wrong, schmong, you've got results. As it is, it took about 15 seconds for
my machine. I'm curious how others' mileage fares as one ratchets up
trials by succesive orders of magnitude.

One of my proudest moments in life was calling volleyballs out against
Michigan that clearly hit the line at my feet, and doing so serially. I
can't explain things I saw on the tele tonight without thinking the
officiating paradigm were the same.

glen herrmannsfeldt

unread,
Aug 22, 2008, 2:47:15 PM8/22/08
to
Ron Ford wrote:
(snip)

> What I don't see here is what makes fortran's implementation of
> pseudorandoms different than C's, at least gcc's. Fortran's is, by
> reputation, good, while C's--at least my vanilla rand call with
> Microsoft--was famously clunky. Is it all the same stuff stuff with
> different window-dressing?

The C standard rand() is often considered a poor random
number generator, with RAND_MAX required to be at least 32767.

Many include rand48() which generates 48 random bits, and uses
a 48 bit seed. That is much better, but not good enough for
all uses.

-- glen

Ron Ford

unread,
Aug 23, 2008, 11:12:37 PM8/23/08
to
On Fri, 22 Aug 2008 10:47:15 -0800, glen herrmannsfeldt posted:

Usually they were interested in an integer on a closed interval smaller
than RAND_MAX, not unlike rolling an eight-sided die.

The problem that I saw in code other than mine was outcomes getting unequal
weight. If RAND_MAX is 32,005, and you take the rand() return modulo 8,
you'll have 5 outcomes heavy or three outcomes light (+- 1).

They talk of the Mersenne Twister as if it would be something relatively
easy to implement in C. That sounds like the Tasmanian devil showing up to
a tea party. That would be a good one for Bugs. How many lumps?

I've got Looney Tunes in my netflix. Woo-hoo.

Ron Ford

unread,
Aug 24, 2008, 12:02:42 AM8/24/08
to
On Wed, 20 Aug 2008 01:10:58 +0000 (UTC), Steven G. Kargl posted:

http://www.apl.washington.edu/people/profile.php?last=Kargl&first=Steve

In sports, people always post against me, and in numbers greater than one.

Your bio has Dayton and '85 in it, so it wouldn't be long until we had an
airdish number of 2.

Do you remember that the Falcons beat the nation that year in college
hockey?

Steven G. Kargl

unread,
Aug 24, 2008, 1:26:08 AM8/24/08
to
In article <1j4l47xr...@example.invalid>,

Ron Ford <r...@example.invalid> writes:
>
> http://www.apl.washington.edu/people/profile.php?last=Kargl&first=Steve
>
> In sports, people always post against me, and in numbers greater than one.
>
> Your bio has Dayton and '85 in it, so it wouldn't be long until we had an
> airdish number of 2.
>
> Do you remember that the Falcons beat the nation that year in college
> hockey?

I have no idea what you're talking about. First, I don't
follow hockey. Second, UD is known as a Flyer. Third,
according to the NCAA, Rensselaer won the hockey champion
in 1985 with a record of 35-2-1.

--
Steve
http://troutmask.apl.washington.edu/~kargl/

Jan Vorbrüggen

unread,
Aug 25, 2008, 2:20:46 AM8/25/08
to
> I don't think they have a case that would hold up. But I'm also not a
> lawyer. Nor am I about to spend the money to hire one to defend myself
> in such a case. Even if I won, the lawyer would cost me more than the
> royalties in question.

ISO is chartered in Switzerland. I believe you would be able to choose
Switzerland as the venue of litigation. If so, Swiss law would likely
apply, and in that case, if you win, ISO pays _all_ the bills, including
yours.

Jan

Ron Ford

unread,
Aug 25, 2008, 4:57:06 PM8/25/08
to
On Sun, 24 Aug 2008 05:26:08 +0000 (UTC), Steven G. Kargl posted:

The Falcons are BG, and they won in 84-85.

We would go down and play Dayton in club volleyball. Unlike against
Michigan, who beat us despite me giving my club two points from being in
plain clothes as a ref, Dayton was not a challenge for us, and they let me
play.

My teammates would laugh at me because I would go block, even if I was in
the back row, but nothing got past me, so they just laughed harder.
--
When a new source of taxation is found it never means, in practice, that
the old source is abandoned. It merely means that the politicians have two
ways of milking the taxpayer where they had one before. 8
H. L. Mencken

Ron Ford

unread,
Aug 26, 2008, 11:41:59 PM8/26/08
to
On Sun, 24 Aug 2008 05:26:08 +0000 (UTC), Steven G. Kargl posted:

>> Do you remember that the Falcons beat the nation that year in college
>> hockey?
>
> I have no idea what you're talking about. First, I don't
> follow hockey. Second, UD is known as a Flyer. Third,
> according to the NCAA, Rensselaer won the hockey champion
> in 1985 with a record of 35-2-1.

BG, with Tom Shirkey and a Hobey Baker award winner--to say nothing of the
only student athletic manager who was bigger than the largest defenseman on
skates-- won the state and the nation in 83-84.

By '85, I was in the middle of a large, jelly-filled donut.

Dayton becomes relevant soon thereafter for its promotion of peace,
consistent with a democratic ohio.

0 new messages