cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c Program to construct Golomb rulers from projective planes
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c IBM SOFTWARE DISCLAIMER
c
c conpp.f (version 1.1)
c Copyright (1998,1986)
c International Business Machines Corporation
c
c Permission to use, copy, modify and distribute this software for
c any purpose and without fee is hereby granted, provided that this
c copyright and permission notice appear on all copies of the software.
c The name of the IBM corporation may not be used in any advertising or
c publicity pertaining to the use of the software. IBM makes no
c warranty or representations about the suitability of the software
c for any purpose. It is provided "AS IS" without any express or
c implied warranty, including the implied warranties of merchantability,
c fitness for a particular purpose and non-infringement. IBM shall not
c be liable for any direct, indirect, special or consequential damages
c resulting from the loss of use, data or projects, whether in action
c of contract or tort, arising out or in the connection with the use or
c performance of this software.
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c Author: James B. Shearer
c email: jbs@watson.ibm.com
c website: http://www.research.ibm.com/people/s/shearer/
c date: 1998 (based on code written in 1986)
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c Version 1.1 (12/11/98) - Renamed variables to conform with
c exhaustive search routines.
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c Requires: ESSL library or portable versions of ESSL routines
c durand, isort (see essl.f)
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c This program constructs good but not necessarily optimal
c Golomb rulers from finte projective planes.
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c Theory
c
c Suppose p is a prime power. Consider the Galois field GF(p) and
c the extension field GF(p**3). Let x be a generator of the cyclic
c multiplicative group of GF(p**3). Then the elements GF(p**3) can
c be represented in the form a+b*x+c*x**2 where a,b,c are elements of
c GF(p). (Note in particular x**3 can be so written, so we may take
c x to be the root of a cubic polynomial over GF(p).) Let two non-
c zero elements, {y,z}, of GF(p**3) be equivalent if one is a scalar
c of the other (ie y/z is an element of the base field GF(p)). This
c partitions the p**3-1 nonzero elements of GF(p**3) into p**2+p+1
c classes (of size p-1). As is well known these classes can be
c thought of as the points of a finite projective plane. Consider
c such a point consisting of the class {y1, y2 ...}. Let y1=x**n1,
c y2=x**n2 ... . We claim n1=n2 mod (p**2+p+1). (Because the
c elements of the base field are generated by x**(p**2+p+1).) Hence
c it is easy to see that we can associate each point of the plane
c with an unique residue mod p**2+p+1. Consider the residues
c associated with the p+1 points on a line in the projective plane.
c We claim these p+1 residues form a distinct difference set mod
c p**2+p+1. Consider for example the points (a+b*x+c*x**2) with
c third coordinate (c) zero. There are p+1 such points which we
c can take to be a+x (a in GF(p)) and 1. Suppose the associated
c residues are not a modular distinct difference set. Then we
c would have for example (a+x)/(b+x)=e*(c+x)/(d+x) (a,b,c,d,e in
c GF(p)). But then x**2+(a+d)*x+a*d=e*(x**2+(b+c)*x+b*c). Or
c equating powers of x, e=1, a+d=b+c, a*d=b*c. So {a,d}={b,c}
c (since they are roots of the same quadratic polynomial). The
c claim follows by contradiction. The other cases involving the
c point 1 are similar.
c Modular distinct difference sets can be unwound and truncated
c to form Golomb rulers. Note we may multiply a modular distinct
c difference set by anything prime to the modulus to obtain another
c modular distinct difference set. The program below tests all
c possibilities to obtain the shortest Golomb rulers.
c The modular difference set construction is due to Singer [2],
c the application to Golomb rulers to Robinson and Bernstein [1].
c The program below finds the best Golomb rulers using this
c construction for prime powers up to maxn-1. It will start to
c fail as maxn**4 overflows integer*4 arithmetic (loop 230). A
c program which just handled primes and not prime powers would be
c simpler.
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c References
c
c 1. J. P. Robinson and A. J. Bernstein, "A class of binary recurrent
c codes with limited error propagation", IEEE Transactions on
c Information Theory, IT-13(1967), p. 106-113.
c 2. J. Singer, "A theorem in finite projective geometry and some
c applications to number theory", Transactions American
c Mathematical Society, 43(1938), p.377-385.
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
Parameter (maxn=160,maxpow=10)
integer*4 len(maxn),nval(maxn),mrec(maxn,maxn)
integer*4 ids(maxn)
integer*4 mw(2*maxn)
integer*4 ipc(3*maxpow),iit(3*maxpow),itemp(3*maxpow)
real*8 buf(3*maxpow)
integer*4 ibase(3*maxpow,2*maxpow),iperp(3*maxpow,maxpow)
integer*4 irep(maxn),ipiv(maxn)
integer*8 array_len, array_size
integer(4), external :: cmp_function
* .. External Subroutines ..
EXTERNAL DLARNV
* .. Local Arrays ..
INTEGER ISEED( 4 )
*
* Initialize seed for random number generator DLARNV.
*
DO 11 I = 1, 4
ISEED( I ) = 1
11 CONTINUE
c initialize best rulers so far
do 5 j=2,maxn
len(j)=maxn*maxn
nval(j)=0
5 continue
c loop over n
do 10 n=2,maxn-1
c check if n is a prime power
c first find smallest prime divisor
do 20 j=2,n
if(mod(n,j).eq.0)go to 30
20 continue
stop "error 20"
30 np=j
npow=1
nprod=j
c next check if n is a power of the smallest prime divisor
40 if(nprod.eq.n)go to 50
nprod=nprod*np
npow=npow+1
if(nprod.le.n)go to 40
c n is not a prime power, go to next n
go to 10
c n is a prime power, construct GF(n**3)=GF(np**ndeg)
50 ndeg=3*npow
c generate random coefficients for monic polynomial P,
c with degree ndeg over GF(np)
60 call irand(ipc,np,ndeg,buf,iseed)
c check if constant term is 0, if so generate another polynomial
if(ipc(1).eq.0)go to 60
c check if x is a multiplicative generator mod P
c initialize iit (x**0) to the unit vector
do 70 j=1,ndeg
iit(j)=0
70 continue
iit(1)=1
c generate powers of x
do 80 j=1,n*n*n-2
c multiply iit by x
itemp(1)=ipc(1)*iit(ndeg)
do 90 i=2,ndeg
itemp(i)=iit(i-1)+ipc(i)*iit(ndeg)
90 continue
do 100 i=1,ndeg
iit(i)=mod(itemp(i),np)
100 continue
c check if power of x is 1 prematurely
if(iit(1).ne.1)go to 80
do 110 i=2,ndeg
if(iit(i).ne.0)go to 80
110 continue
c this polynomial no good, go generate another
go to 60
80 continue
c All powers of x ok. We now have a representation of GF(n**3) as
c arithmetic mod a polynomial over GF(np)
c The field GF(n**3) is an extension of the field GF(n). We need
c a way of identifying elements of the subspace generated by GF(n)
c and x. This is easy when n is a prime rather than a prime power,
c these are the elements with x**2 term 0. The prime power case is
c more complicated, we find a linear basis for the space and then
c for the perpendicular space. We can then test by looking at
c inner products with the basis of the perpendicular space.
c First npow powers of x**(n**2+n+1) span GF(n)
c Set first element of GF(n) basis to 1
do 85 j=1,ndeg
ibase(j,1)=0
85 continue
ibase(1,1)=1
c find remaining elements of basis
ibp=n**2+n+1
c initialize iit (x**0) to the unit vector
do 120 j=1,ndeg
iit(j)=0
120 continue
iit(1)=1
c generate powers of x
do 125 ja=2,npow
do 130 j=1,ibp
c multiply iit by x
itemp(1)=ipc(1)*iit(ndeg)
do 140 i=2,ndeg
itemp(i)=iit(i-1)+ipc(i)*iit(ndeg)
140 continue
do 150 i=1,ndeg
iit(i)=mod(itemp(i),np)
150 continue
130 continue
c add to basis vectors
do 160 j=1,ndeg
ibase(j,ja)=iit(j)
160 continue
125 continue
c now find basis of space generated by GF(n) and x.
c generate n additional basis vectors by multiplying first n vectors by x
do 170 ja=1,npow
itemp(1)=ipc(1)*ibase(ndeg,ja)
do 180 i=2,ndeg
itemp(i)=ibase(i-1,ja)+ipc(i)*ibase(ndeg,ja)
180 continue
do 190 i=1,ndeg
ibase(i,npow+ja)=mod(itemp(i),np)
190 continue
170 continue
c we now generate a table of inverses to help us do arithmetic in GF(np)
do 200 j=1,np-1
do 210 i=1,np-1
if(mod(i*j,np).ne.1)go to 210
irep(j)=i
go to 200
210 continue
stop "error 210"
200 continue
c put basis in a normal form.
do 220 j=1,2*npow
c find first non-zero in column j
do 230 i=1,ndeg
if(ibase(i,j).ne.0)go to 240
230 continue
stop "error 230"
c record position of first non-zero
240 ipiv(j)=i
c multiply column so non-zero becomes 1
imult=irep(ibase(i,j))
do 250 i=1,ndeg
ibase(i,j)=mod(ibase(i,j)*imult,np)
250 continue
c zero remaining elements in row
do 260 ja=1,2*npow
if(ja.eq.j)go to 260
imult=ibase(ipiv(j),ja)
do 270 i=1,ndeg
ibase(i,ja)=mod(ibase(i,ja)-imult*ibase(i,j),np)
if(ibase(i,ja).lt.0)ibase(i,ja)=ibase(i,ja)+np
270 continue
260 continue
220 continue
c now construct perpendicular basis
jp=0
do 280 j=1,ndeg
c look for row not among the 2*npow recorded in ipiv
do 290 i=1,2*npow
if(ipiv(i).eq.j)go to 280
290 continue
jp=jp+1
c zero jp'th vector in perpendicular basis
do 300 i=1,ndeg
iperp(i,jp)=0
300 continue
c fill in elements in ipiv rows so as to make inner products 0
do 310 ja=1,2*npow
iperp(ipiv(ja),jp)=ibase(j,ja)
310 continue
c note -1 = np - 1 mod np
iperp(j,jp)=np-1
280 continue
c construct perfect difference set mod n**2+n+1
c look for powers of x in space spanned by GF(n), x
c put 0 in set
ids(1)=0
nc=1
c initialize iit
do 320 j=1,ndeg
iit(j)=0
320 continue
iit(1)=1
c look for powers of x in space spanned by GF(n), x
do 330 j=1,n*n+n
c multiply by x
itemp(1)=ipc(1)*iit(ndeg)
do 340 i=2,ndeg
itemp(i)=iit(i-1)+ipc(i)*iit(ndeg)
340 continue
do 350 i=1,ndeg
iit(i)=mod(itemp(i),np)
350 continue
c check inner products
do 360 ja=1,npow
iip=0
do 370 i=1,ndeg
iip=iip+iperp(i,ja)*iit(i)
370 continue
iip=mod(iip,np)
if(iip.ne.0)go to 330
360 continue
nc=nc+1
ids(nc)=j
330 continue
c check that difference set is right size
if(nc.ne.n+1)stop "error 330"
c output difference set
c write(6,1000)n,(ids(j),j=1,n+1)
c1000 format(1x,i5,5x,(10i5))
c check for better rulers
c cycle over multipliers, don't need to try j and -j mod (n*n+n+1)
c do 420 j=1,n*n+n
do 420 j=1,(n*n+n+1)/2
c check if multiplier prime to modulus
if(igcd(j,n*n+n+1).ne.1)go to 420
c multiply difference set
do 430 i=1,n+1
mw(i)=mod(ids(i)*j,n*n+n+1)
430 continue
c sort new difference set
array_len = n + 1
array_size = 4
CALL qsort(mw,array_len, array_size,cmp_function)
c unwrap difference set
do 440 i=1,n+1
mw(i+n+1)=mw(i)+n*n+n+1
440 continue
c check for new records
do 450 ia=1,n+1
do 460 ib=1,n
if(mw(ia+ib)-mw(ia).ge.len(ib+1))go to 460
c new record ruler
len(ib+1)=mw(ia+ib)-mw(ia)
nval(ib+1)=n
do 470 ja=1,ib+1
mrec(ja,ib+1)=mw(ia+ja-1)-mw(ia)
470 continue
460 continue
450 continue
420 continue
10 continue
c output maximum rulers
do 500 j=2,maxn
if(nval(j).eq.0)go to 500
c put ruler in standard form (flip if needed)
if(mrec((j+1)/2,j)+mrec((j+2)/2,j).lt.len(j))go to 520
c flip ruler
do 510 i=1,(j+1)/2
mtemp=mrec(i,j)
mrec(i,j)=len(j)-mrec(j+1-i,j)
mrec(j+1-i,j)=len(j)-mtemp
510 continue
520 continue
c
c write results for j marks
c
c j,length and prime power to unit 6 (terminal)
c j,length and prime power to unit 1 (disk)
c ruler to unit 1 (disk)
c
write(6,1010)j,len(j),nval(j)
write(1,1020)j,len(j),nval(j)
write(1,1030)(mrec(i,j),i=1,j)
1010 format(1x,3i10)
1020 format(3i10)
1030 format(10i6)
500 continue
c mark end of disk file
write(1,1020)0,0,0
stop
end
c generate random vector of n integers 0,...,np-1
subroutine irand(l,np,n,x,iseed)
integer*4 l(*)
real*8 x(*)
real*8 dseed/1.d0/
save dseed
c generate random 0-1 real*8 vector
CALL DLARNV( 1, ISEED, N, X )
c convert to integer 0,...,np-1
do 10 j=1,n
l(j)=np*x(j)
10 continue
return
end
c find gcd of ia,ib using Euler's method
function igcd(ia,ib)
ja=ia
jb=ib
1 jc=mod(ja,jb)
ja=jb
jb=jc
if(jb.ne.0)go to 1
igcd=ja
return
end
integer(4) function cmp_function(a1, a2)
integer(4) a1, a2
cmp_function=a1-a2
end function
16.12.2015, 10:32 
лишних элементов (голомбовость при этом сохраняется).