c RSA6.f  Jun 14 2008
c Find how close can match N from above
c by a multiple of e, e odd

c The program looks at the cyclotomic cosets mod e
c Cf. MacWilliams and Sloane, Error-Correcting Codes, p. 104.

	implicit integer(a-z)
	integer coset(130,10000)
	integer mum(10000), length(130)
	integer nblis(10000),height(130),noad(1000),lisad(130,1000)
	integer seen(130),ans(10000)
c mum(i) tells which coset i belongs to
c coset(i,*) lists ith coset

c max number of cosets
	Morb=130
c max e
	Me=10000
c set e
	lim=1000
	do 300 ie=1,lim
	e=2*ie+1
	if(e.le.Me)goto 30
	write(*,*)"ERROR e too big"
	stop
30	continue
	write(*,*)" "
	write(*,*)"Starting e= ",e
	do 4 i=1,e
4	mum(i)=0
	do 44 i=1,Morb
44	seen(i)=0
c get base coset
c norb = coset number
c s1 = current size of coset
c t1 = current member of coset
	norb=1
	t1=1
	s1=1
	coset(norb,s1)=t1
	mum(t1)=norb
	do 1 i=1,e
	t1=mod(2*t1,e)
	if(t1.eq.1)goto 2
	s1=s1+1
	coset(norb,s1)=t1
	mum(t1)=norb
1	continue
2	continue
	length(norb)=s1

c get remaining cosets, if any
	do 3 i=3,e-2,2
	if(mum(i).gt.0)goto 3
c have new multiplier, i
	norb=norb+1
	if(norb.le.Morb)goto 7
	write(*,*)" ERROR - too many cosets"
	stop
7	continue
	coset(norb,1)=i
	len=1
	mum(i)=norb
	do 5 j=2,s1
	t1=mod(i*coset(1,j),e)
	if(t1.eq.i)goto 6
	len=len+1
	coset(norb,len)=t1
	mum(t1)=norb
5	continue
6	continue
	length(norb)=len
3	continue
c now have all cosets
	nblis(ie)=norb
	write(*,*)"Number of cosets = ",norb
	write(*,*)"lengths of cosets"
	write(*,100)(length(i),i=1,norb)
	do 8 i=1,norb
	write(*,*)"coset ", i
	write(*,99)(coset(i,j),j=1,length(i))
8	continue
99	format(12i5)
100	format(10i6)
	write(*,*)"mum:"
	write(*,99)(i,i=1,e-1)
	write(*,99)(mum(i),i=1,e-1)

c get heights of cosets

c noad(h) = number of cosets at height h
c lisad(h,*) = list of cosets at height h
	noad(1)=1
	lisad(1,1)=1
	height(1)=1
	seen(1)=1
	maxht=1
c is there only one coset?
	if(norb.eq.1)goto 299
c there is more than one coset
c add coset 1 to itself
	noad(2)=0
	do 10 i1=1,s1-1
	do 10 i2=i1+1,s1
	i3=mod(coset(1,i1)+coset(1,i2),e) 
	if(i3.eq.0)goto 10
	i4=mum(i3)
c	write(*,99)i1,i2,i3,i4,seen(i4)
	if(seen(i4).eq.1)goto 10
	height(i4)=2
	seen(i4)=1
	noad(2)=noad(2)+1
	lisad(2,noad(2))=i4
10	continue
	maxht=2
c were there only 2 cosets?
	if(norb.eq.2)go to 299
c	write(*,*)"seen"
c	write(*,99)(seen(i),i=1,norb)

290	continue
c	write(*,*)"Point 290, maxht, newht= ",maxht,newht
	do 122 i=1,maxht
	write(*,*)noad(i)
	write(*,99)(lisad(i,j),j=1,noad(i))
122	continue
	write(*,*)"seen"
	write(*,99)(seen(i),i=1,norb)
c have we seen all the cosets?
	do 21 i=1,norb
	if(seen(i).ne.1)goto 22
21	continue
c yes done
	goto 299
22	continue
c no, there are unvisited cosets
c add base coset to all cosets at height maxht
	newht=maxht+1
	noad(newht)=0
	do 23 io=1,noad(maxht)
c	write(*,*)"pt 23",maxht,newht,noad(maxht),io,i1
	il=lisad(maxht,io)
c need to add base coset to coset il 
        do 24 i1=1,s1
        do 24 i2=1,length(il)
        i3=mod(coset(1,i1)+coset(il,i2),e)
        if(i3.eq.0)goto 24
        i4=mum(i3)
c       write(*,99)i1,i2,i3,i4,seen(i4)
        if(seen(i4).eq.1)goto 24
        height(i4)=newht
        seen(i4)=1
        noad(newht)=noad(newht)+1
        lisad(newht,noad(newht))=i4
24      continue
23	continue
        maxht=newht
	goto 290
299	continue
	ans(ie)=maxht
	write(*,*)"maxht = ",maxht
	write(*,*)"List of cosets by height"
	do 12 i=1,maxht
	write(*,*)noad(i)
	write(*,99)(lisad(i,j),j=1,noad(i))
12	continue
300	continue

c print summary
	write(*,*)"All done, print summary"
	write(*,99)(2*ie+1,ie=1,lim)
	write(*,*)"List of max heights"
	write(*,99)(ans(ie),ie=1,lim)
	write(*,*)"List of numbers of cosets:"
	write(*,99)(nblis(ie),ie=1,lim)
	stop
	end