############################################################################
#
#	File:     genrfncs.icn
#
#	Subject:  Procedures to generate sequences
#
#	Author:   Ralph E. Griswold
#
#	Date:     March 4, 2003
#
############################################################################
#
#   This file is in the public domain.
#
############################################################################
#  
#  These procedures generate sequences of results.
#
#  arandseq(i, j)	arithmetic sequence starting at i with randomly
#			chosen increment between 1 and j
#  
#  arithseq(i, j)	arithmetic sequence starting at i with increment j
#
#  beatty1seq()		Beatty's first sequence i * &phi
#
#  beatty2seq()		Beatty's second sequence i * &phi ^ 2
#
#  catlnseq(i)		sequence of generalized Catalan numbers
#
#  cfseq(i, j)		continued-fraction sequence for i / j
#
#  chaosseq()		chaotic sequence
#
#  chexmorphseq()	sequence of centered hexamorphic numbers
#
#  connellseq(p)	generalized Connell sequence
#
#  dietzseq(s)		Dietz sequence for polynomial
#
#  dressseq(i)		dress sequence with increment i, default 1 (Schroeder)
#
#  eisseq(i)		EIS A sequence for i
#
#  factseq()		factorial sequence
#
#  fareyseq(i, k)	Farey fraction sequence; k = 0, the default, produces
#			numerator sequence; k = 1 produces denominator
#			sequence
#
#  fibseq(i, j, k, m)	generalized Fibonacci sequence (Lucas sequence)
#			with initial values i and j and additive constant
#			k.  If m is supplied, the results are produced
#			mod m.
#
#  figurseq(i)		series of ith figurate number
#
#  fileseq(s, i)	generate from file s; if i is null, lines are generated.
#			Otherwise characters, except line terminators.
#
#  friendseq(k)		generate random friendly sequence from k values, 1 to k
#			(in a friendly sequence, successive terms differ by 1).
#
#
#  geomseq(i, j)	geometric sequence starting at i with multiplier j
#
#  hailseq(i)		hailstone sequence starting at i
#
#  irepl(i, j)		j instances of i
#
#  lindseq(f, i)	generate symbols from L-system in file f; i if
#			present overrides the number of generations specified
#			in the L-system.
#
#  logmapseq(k, x)	logistic map
#
#  lrrcseq(L1, L2)
#			generalized linear recurrence with constant
#			coefficients; L1 is a list of initial terms,
#			L2 is a list of coefficients for n previous values,
#			where n = *L2
#
#  meanderseq(s, n)	sequences of all characters that contain all n-tuples
#			of characters from s
#	
#  mthueseq()		Morse-Thue sequence
#
#  mthuegseq(i)		Morse-Thue sequence for base i
#
#  multiseq(i, j, k)	sequence of (i * j + k) i's
#
#  ngonalseq(i)		sequence of the ith polygonal number
#
#  nibonacciseq(values[])
#			generalized Fibonacci sequence that sums the
#			previous n terms, where n = *values.
#
#  partitseq(i, j, k)	sequence of integer partitions of i with minimum j
#			and maximum k
#
#  pellseq(i, j, k)	generalized Pell's sequence starting with i, j and
#			using multiplier k
#
#  perrinseq()		Perrin sequence
#
#  polyseq(coeff[])	polynomial in x evaluated for x := seq()
#	
#  primeseq()		the sequence of prime numbers
#
#  powerseq(i)		sequence n ^ i, n = 1, 2, 3, 4, ...
#
#  powersofseq(i)	sequence i ^ n, n = 1, 2, 3, 4, ...n
#
#  rabbitseq()		rabbit sequence
#
#  ratsseq(i)		versumseq() with sort
#
#  signaseq(r)		signature sequence of r
#
#  spectseq(r)		spectral sequence integer(i * r), i - 1, 2, 3, ...
#
#  srpseq(n, m)		palindromic part of the continued-fraction sequence
#			for sqrt(n^2+m)
#
#  versumseq(i, j)	generalized sequence of added reversed integers with
#			seed i (default 196) and increment j (default 0)
#
#  versumopseq(i, p)	procedure p (default 1) applied to versumseq(i)
#
#  vishwanathseq()	random variation on Fibonacci sequence
#
#  zebra(values[])	zebra colors, alternating 2 and 1, for number of
#			times given by successive values
#  
############################################################################
#
#  Requires:  co-expressions
#
############################################################################
#
#  Links:  convert, fastfncs, io, partit, numbers, rational, xcode
#	   polynom, strings
#
############################################################################

link convert
link lists
link fastfncs
link io
link numbers
link partit
link polynom
link rational
link xcode
link periodic
link factors
link strings

procedure arandseq(i, j)	#: arithmetic sequence with random intervals

   /i := 1
   /j := 1

   suspend seq(i) + ?j

end

procedure arithseq(i, j)	#: arithmetic sequence

   /i := 1
   /j := 0

   suspend seq(i) + j

end

procedure beatty1seq(r)		#: Beatty sequence 1

   /r := &phi

   suspend integer(seq() * r)

end

procedure beatty2seq(r)		#: Beatty sequence 2

   /r := &phi

   suspend integer(seq() * (r / (r - 1)))

end

procedure catlnseq(i)		#: generalized Catalan sequence
   local k

   /i := 1

   suspend (i := 1, k := seq(), i *:= 4 * k + 2, i /:= k + 2)

end

procedure chaosseq()		#: Hofstadter's chaotic sequence

   suspend q(seq())

end

#  The generalization here is to allow a generating procedure, p to
#  be specified.  The default is seq().  Arguments are given in args.

procedure connellseq(p, args[])	#: generalized Connell sequence
   local i, j, count, parity, parity2, C

   C := create (\p | seq) ! args

   count := 0
   parity := 0
   parity2 := 1

   repeat {
      count +:= 1
      parity :=: parity2
      j := 0
      repeat {
         i := @C | fail
         if i % 2 = parity then {
            suspend i
            j +:= 1
            if j = count then break
            }
         }
      }

end
   
procedure chexmorphseq()	#: sequence of centered hexamorphic numbers
   local i, j

   every (i := seq(), j := 3 * i * (i - 1) + 1, j ? {
      tab(-*i)
      if =i then suspend j
      })

end

procedure cfseq(i, j)	#: continued-fraction sequence
   local r

   until j = 0 do {
      suspend integer(i / j)
      r := i % j
      i := j
      j := r
      }

end

procedure dietzseq(str)

   suspend !poly2profile(peval(str))

end

procedure dressseq(i)
   local seq, seq1, n

   /i := 1

   seq := [0]

   suspend seq[1]

   repeat {
      seq1 := copy(seq)
      every n := !seq + i do {
         suspend n
         put(seq1, n)
         }
      seq := seq1
      }

end

procedure eisseq(i)		#: EIS A sequence
   local input, seq
   static lst

   initial {
      input := dopen("eis.seq") | fail
      lst := xdecode(input) | fail
      close(input)
      }

   seq := \lst[integer(i)] | fail

   suspend !seq

end

procedure factseq()		#: factorial sequence
   local i

   i := 1

   suspend i *:= seq()

end

record farey(magnitude, n, d)

procedure fareyseq(i, k)		#: Farey fraction sequence
   local farey_list, n, d, x

   /k := 0				# default numerators

   k := integer(k) | fail

   farey_list := [farey(0.0, 0, 1)]

   every d := 1 to i do
      every n := 1 to d do {
         if gcd(n, d) = 1 then
            put(farey_list, farey(real(n) / d, n, d))
         }

   farey_list := sortf(farey_list, 1)

   case k of {
      0  :  every suspend (!farey_list).n	# numerator sequence
      1  :  every suspend (!farey_list).d	# denominator sequence
      }

end
      
procedure fareydseq(i)		#: Farey fraction denominator sequence
   local parity, j

   parity := 1

   every j := fareyseq(i) do {
      if parity < 0 then suspend j
      parity *:= -1
      }

end

procedure fareynseq(i)		#: Farey fraction numerator sequence
   local parity, j

   parity := 1

   every j := fareyseq(i) do {
      if parity > 0 then suspend j
      parity *:= -1
      }

end

procedure fareyn1seq(i)		#: Farey fraction numerator sequence, 1-based

   suspend fareynseq(i) + 1

end

procedure fibseq(i, j, k, m)	#: generalized Fibonacci sequence
   local n

   /i := 1
   /j := 1
   /k := 0

   if /m then {
      suspend i | j | |{
         n := i + j + k
         i := j
         j := n
         }
      }
   else {
      suspend i % m | j % m | |{
         n := (i + j + k) % m
         i := j
         j := n
         }
      }

end

#  Warning; if not all lines are generated from the input file, the
#  file is not closed until the next call of fileseq().

procedure fileseq(s, i)		#: sequence from file
   static input

   close(\input)

   input := dopen(s) | fail

   if /i then suspend !input
   else suspend !!input

   close(input)

   input := &null

end

procedure figurseq(i)		#: sequence of figurate numbers
   local j, k

   /i := 1

   suspend (j := 1, k := seq(i), j *:= k + 1, j /:= k + 1 - i)

end

procedure friendseq(k)		#: random friendly sequence
   local state

   state := ?k

   repeat {
      suspend state
      case state of {
         1        :  state +:= 1
         k        :  state -:= 1
         default  :  state +:= ?[1, -1]
         }
      }

end
      
procedure geomseq(i, j)		#: geometric sequence

   /i := 1
   /j := 1

   suspend seq(i) * j

end

procedure hailseq(i)		#: hailstone sequence

   /i := 1

   suspend |if i % 2 = 0 then i /:= 2 else i := 3 * i + 1

end

procedure irepl(i, j)		#: repeated sequence

   /i := 1
   /j := 1

   suspend |i \ j

end

procedure lindseq(f, i, p)			# generate symbols from L-system
   local input, gener

   /p := "lindsys"

   if \i then input := open(p || " -g " || i || " <" || f, "p")
   else input := open(p || " <" || f, "p")

   while gener := read(\input) do
      suspend !gener

   close(input)		# pipe will be left open if not all result are generated

   fail

end

procedure logmapseq(k, x)		# logistic map

   suspend x := k * x * (1 - |x)

end

procedure linrecseq(terms, coeffs)	#: synonym for lrrcseq
   linrecseq := lrrcseq

   suspend lrrcseq(terms, coeffs)

end

procedure lrrcseq(terms, coeffs)	#: linear recurrence sequence
   local i, term

   suspend !terms

   repeat {
      term := 0
      every i := 1 to *coeffs do
         term +:= terms[i] * coeffs[-i]
      suspend term
      get(terms)
      put(terms, term)
      }

end

procedure meanderseq(alpha, n)	#: generate meandering characters
   local sequence, trial, i, c

   i := *alpha

   sequence := repl(alpha[1], n - 1)			# base string

   while c := alpha[i] do {			# try a character
      trial := right(sequence, n - 1) || c
      if find(trial, sequence) then
         i -:= 1
      else {
         sequence ||:= c				# add it
         i := *alpha				# and start from end again
         suspend c
         }
      }

end

procedure mthueseq()		#: Morse-Thue sequence
   local s, t

   s := 0

   suspend s

   repeat {
      t := map(s, "01", "10")
      every suspend integer(!t)
      s ||:= t
      }

end

procedure mthuegseq(j)		#: generalized Morse-Thue sequence

   suspend adr(exbase10(seq(0), j)) % j		# only works through base 10
      
end
   
procedure multiseq(i, j, k)	#: sequence of repeated integers

   /i := 1
   /j := 1
   /k := 0

   suspend (i := seq(i), (|i \ (i * j + k)) & i)

end

procedure ngonalseq(i)		#: sequence of polygonal numbers
   local j, k

   /i := 2

   k := i - 2

   suspend ((j := 1) | (j +:= 1 + k * seq()))

end

procedure nibonacciseq(values[])	#: n-valued Fibonacci generalization
   local sum

   if *values = 0 then fail

   suspend !values

   repeat {
      sum := 0
      every sum +:= !values
      suspend sum
      get(values)
      put(values, sum)
      }

end

procedure partitseq(i, j, k)	#: sequence of integer partitions

   /i := 1
   /j := 1
   /k := i

   suspend !partit(i, j, k)

end

procedure pellseq(i, j, k)	#: generalized Pell sequence
   local m

   /i := 1
   /j := 2
   /k := 2

   suspend i | j | |{
      m := i + k * j
      i := j
      j := m
      }

end

procedure perrinseq()		#: perrin sequence
   local i, j, k, l

   suspend i := 0
   suspend j := 2
   suspend k := 3

   repeat {
      suspend l := i + j
      i := j
      j := k
      k := l
      }

end

procedure polyseq(coeff[])	#: sequence of polynomial evaluations
   local i, j, sum

   every i := seq() do {
      sum := 0
      every j := 1 to *coeff do
         sum +:= coeff[j] * i ^ (j - 1)
      suspend sum
      }

end

procedure primeseq()		#: sequence of prime numbers
   local i, k

   suspend 2 | ((i := seq(3, 2)) & (not(i = (k := (3 to sqrt(i) by 2)) *
      (i / k))) & i)

end

procedure powersofseq(i)		#: powers

   /i := 2

   suspend i ^ seq(i)

end

procedure powerseq(i)		#: powers sequence

   suspend seq() ^ i

end

procedure rabbitseq()		#: rabbit sequence
   local seq, i

   seq := [0]

   suspend 1

   repeat {
      i := get(seq)
      suspend i
      if i = 0 then put(seq, 1)
      else put(seq, 1, 0)
      }

end
   
procedure ratsseq(i, p)		#: reverse add and then sort sequence

   /p := 1

   repeat {
      i +:= reverse(i)
      i := integer(p(csort(i)))
      suspend i
      }

end

record entry(value, i, j)

procedure signaseq(r, n, m)	#: signature sequence
   local i, j, result

   /n := 100
   /m := 100

   result := []

   every j := 1 to n do 
      every i := 1 to m do
         put(result, entry(i + j * r, i, j))

   result := sortf(result, 1)

   suspend (!result)[2]

end

procedure spectseq(r)		#: spectral sequence

   /r := 1.0

   suspend integer(seq() * r)

end


procedure srpseq(n, m)		#: generate square-root palindrome
   local iter, count, okay, rat, j, pal

   if not (1 <= m <= 2 * n) then fail

   iter := 5

   repeat {
      pal := []
      count := 0
      okay := 1
      rat := Sqrt(n ^ 2 + m, iter)
      every j := cfseq(rat.numer, rat.denom) do {
         count +:= 1
         if count = 1 then next	# don't examine first term
         if j = 2 * n then {	# presumed end
            if not lequiv(pal, lreverse(pal)) then break
            okay := &null
            break
            }
         else if j > n then break			# too big; error
         else put(pal, j)
         }
      if \okay then {
         iter +:= 1			# back to repeat loop
         if iter > 12 then fail		# too many iterations required.
         next
         }
      break
      }

   suspend !pal

end

procedure versumseq(i, j)	#: generalized reversed-sum sequence

   /j := 0

   /i := 196

   repeat {
      i +:= reverse(i) + j
      suspend i
      }

end

procedure versumopseq(i, p, args[])	#: versum sequence with operator

   /i := 196

   /p := csort

   push(args, &null)			# make room for first argument

   repeat {
      i := reverse(i)
      args[1] := i		# make current i first argument
      i := integer(p ! args)
      suspend i
      }

end

procedure vishwanathseq(i, j)	#: random variation on Fibonacci sequence
   local m

   /i := 1
   /j := 1

   suspend i | j | |{
      m := case ?4 of {
         1  :  i + j
         2  :  i - j
         3  :  -i + j
         4  :  -i - j
         }
      i := j
      j := m
      }

end

procedure zebra(args[])		#: black and white bands
   local i, clr, clr_alt

   clr := 2			# light
   clr_alt := 1			# dark

   while i := get(args) do {
      suspend (1 to i) & clr
      clr :=: clr_alt
      }

end