http://mathling.com/core/sequences  library module

Variable: $ODIOUS as xs:integer*

First few odious numbers = numbers with odd numbers of 1s in binary
expansion
A000069

Variable: $KATYDID as xs:integer*

First few Katydid numbers
A060031

Variable: $GRASSHOPPER as xs:integer*

First few Grasshopper numbers

Function: fibonacci-sequence
declare function fibonacci-sequence($n as xs:integer) as xs:integer*

fibonacci-sequence()
All Fibonacci numbers to the nth (starting at 1)
A000045

Params

• n as xs:integer: up to how many to generate n >= 1

Returns

• xs:integer*: sequence

declare function this:fibonacci-sequence($n as xs:integer) as xs:integer*
{
  if ($n lt 1) then errors:error("ML-BADARGS", ("n", $n))
  else if ($n = 1) then 1
  else if ($n = 2) then (1,1)
  else (
    fold-left(
      1 to ($n - 2), (1,1),
      function($fib as xs:integer*, $i as xs:integer) {
        $fib, $fib[last()]+$fib[last()-1]
      }
    )
  )
}

Errors

ML-BADARGS if $n is less than 1

e.g.: $n=1 => (1)
e.g.: $n=7 => (1,1,2,3,5,8,13)

Function: fibonacci
declare function fibonacci($n as xs:integer) as xs:integer

fibonacci()
Return nth Fibonacci number

Params

• n as xs:integer: which item in the sequence to return

Returns

• xs:integer: integer

declare function this:fibonacci($n as xs:integer) as xs:integer
{
  if ($n lt 1) then errors:error("ML-BADARGS", ("n", $n))
  else if ($n = 1) then 1
  else if ($n = 2) then 1
  else this:fibonacci-sequence($n)[last()]
}

Errors

ML-BADARGS if $n is less than 1

Function: fibonacci-word
declare function fibonacci-word($n as xs:integer) as xs:integer*

fibonacci-word()
Return the first $n terms of the infinite Fibonacci word
e.g. fibonacci-word(10) = (0,1,0,0,1,0,1,0,0,1)
A003849

Params

• n as xs:integer: up to how many to generate n >= 1

Returns

• xs:integer*: sequence

declare function this:fibonacci-word($n as xs:integer) as xs:integer*
{
  if ($n lt 1) then errors:error("ML-BADARGS", ("n", $n))
  else if ($n = 1) then 0
  else if ($n = 2) then (0, 1)
  else (
    util:while(
      function($a as array(*), $b as array(*)) as xs:boolean {
        count(array:flatten($b)) <= $n
      },
      function($a as array(*), $b as array(*)) as array(*)* {
        $b,
        array { array:flatten($b), array:flatten($a) }
      },
      array {0}, array {0, 1}
    )=>tail()=>array:flatten()
  )[position() <= $n]
}

Errors

ML-BADARGS if $n is less than 1

Function: padovan
declare function padovan($n as xs:integer) as xs:integer*

padovan()
Return the first n terms of the Padovan sequence
a(n) = a(n-2) + a(n-3); a(0)=1, a(1)=a(2)=0
First few terms:
1,0,0,1,0,1,1,1,2,2,3,4,5,7,9,12,16,21,28,37,49,65,86,114,151,200,265,351

A000931

Params

• n as xs:integer: up to how many to generate n >= 1

Returns

• xs:integer*: sequence

declare function this:padovan($n as xs:integer) as xs:integer*
{
  if ($n lt 1) then errors:error("ML-BADARGS", ("n", $n))
  else if ($n = 1) then 1
  else if ($n = 2) then (1,0)
  else if ($n = 3) then (1,0,0)
  else (
    fold-left(
      1 to ($n - 3), (1,0,0),
      function($pad as xs:integer*, $i as xs:integer) {
        $pad, $pad[last()-1]+$pad[last()-2]
      }
    )
  )
}

Errors

ML-BADARGS if $n is less than 1

Function: fractal
declare function fractal($n as xs:integer) as xs:integer*

fractal()
Triangle read by rows
First few terms:
1,1,2,1,2,3,1,2,3,4,1,2,3,4,5,1,2,3,4,5,6,1,2,3,4,5,6,7,1,2,3,4,5,6,7,8

A002260

Params

• n as xs:integer: up to how many to generate

Returns

• xs:integer*: sequence

declare function this:fractal($n as xs:integer) as xs:integer*
{
  (: n = Σ[i=1:k]i = (1 + k) * k/2 => k = -1/2 + √(2n + 1/4) :)
  (: ceiling to make sure we get enough :)
  let $k := ceiling(math:sqrt(2 * $n + 0.25) - 0.25) cast as xs:integer
  return (
    for $i in 1 to $k
    for $j in 1 to $i
    return $j
  )[position() <= $n]
}

Function: golomb
declare function golomb($n as xs:integer) as xs:integer*

golomb()
a[n] is number of times n occurs
First few terms:
1,2,2,3,3,4,4,4,5,5,5,6,6,6,6,7,7,7,7,8,8,8,8,9,9,9,9,9,10,10,10,10,10,11

A001462

Params

• n as xs:integer: up to how many to generate

Returns

• xs:integer*: sequence

declare function this:golomb($n as xs:integer) as xs:integer*
{
  if ($n lt 1) then errors:error("ML-BADARGS", ("n", $n))
  else if ($n = 1) then 1
  else if ($n = 2) then (1, 2)
  else if ($n = 3) then (1, 2, 2)
  else (
    fold-left(
      3 to $n, (1, 2, 2),
      function($gol as xs:integer*, $i as xs:integer) {
        $gol, for $j in 1 to $gol[$i] return $i
      }
    )[position() <= $n]
  )
}

Errors

ML-BADARGS if $n is less than 1

Function: lazy-caterer
declare function lazy-caterer($n as xs:integer) as xs:integer*

lazy-caterer()
n(n+1)/2 + 1
First few terms:
1,2,4,7,11,16,22,29,37,46,56,67,79,92,106,121,137,154,172,191,211,232,254

A000124

Params

• n as xs:integer: up to how many to generate

Returns

• xs:integer*: sequence

declare function this:lazy-caterer($n as xs:integer) as xs:integer*
{
  if ($n lt 1) then errors:error("ML-BADARGS", ("n", $n))
  else for $i in 0 to $n - 1 return (floor($i * ($i + 1) div 2) + 1) cast as xs:integer
}

Errors

ML-BADARGS if $n is less than 1

Function: jacobsthal
declare function jacobsthal($n as xs:integer) as xs:integer*

jacobsthal()
a(n)=a(n-1) + 2*a(n-2) starting 0, 1, 1
First few terms:
0,1,1,3,5,11,21,43,85,171,341,683,1365,2731,5461,10923,21845,43691,87381

A001045

Params

• n as xs:integer: up to how many to generate

Returns

• xs:integer*: sequence

declare function this:jacobsthal($n as xs:integer) as xs:integer*
{
  if ($n lt 1) then errors:error("ML-BADARGS", ("n", $n))
  else if ($n = 1) then 0
  else if ($n = 2) then (0, 1)
  else if ($n = 3) then (0, 1, 1)
  else (
    fold-left(
      4 to $n, (0, 1, 1),
      function($jac as xs:integer*, $i as xs:integer) {
        $jac, $jac[last()] + 2 * $jac[last()-1]
      }
    )
  )
}

Errors

ML-BADARGS if $n is less than 1

Function: jacobsthal-lucas
declare function jacobsthal-lucas($n as xs:integer) as xs:integer*

jacobsthal-lucas()
a(n+1) = 2*a(n) - (-1)^n*3 = 2^n + (-1)^n
First few terms:
2,1,5,7,17,31,65,127,257,511,1025,2047,4097,8191,16385,32767,65537,131071

A014551

Params

• n as xs:integer: up to how many to generate

Returns

• xs:integer*: sequence

declare function this:jacobsthal-lucas($n as xs:integer) as xs:integer*
{
  if ($n lt 1) then errors:error("ML-BADARGS", ("n", $n))
  else for $i in 0 to $n - 1 return (math:pow(2, $i) + math:pow(-1, $i)) cast as xs:integer
}

Errors

ML-BADARGS if $n is less than 1

Function: recaman
declare function recaman($n as xs:integer) as xs:integer*

recaman()
Recamán's sequence. a(0)=0; a(n)=a(n)-n if >=0 and not already in sequence;
a(n)=a(n-1)+n otherwise
First few terms:
0,1,3,6,2,7,13,20,12,21,11,22,10,23,9,24,8,25,43,62,42,63,41,18,42,17,43,16

A005132

Params

• n as xs:integer: up to how many to generate

Returns

• xs:integer*: sequence

declare function this:recaman($n as xs:integer) as xs:integer*
{
  if ($n lt 1) then errors:error("ML-BADARGS", ("n", $n))
  else if ($n = 1) then 0
  else (
    fold-left(
       1 to $n - 1, 0,
       function($seq as xs:integer*, $i as xs:integer) {
         $seq,
         let $an := $seq[last()] - $i
         return (
           if ($an >= 0 and not($an=$seq)) then $an
           else $seq[last()] + $i
         )
       }
    )
  )
}

Errors

ML-BADARGS if $n is less than 1

Function: gerrymander2
declare function gerrymander2($n as xs:integer) as xs:integer*

gerrymander2()
2-party gerrymander: a(n) is min number of total votes for one party
to win when there are n² voters divided into equal districts.
First few terms:
1,3,4,8,9,14,16,24,25,33,36,45,49,60,64,80,81,95,100,117,121,138,144,165
A341578

Params

• n as xs:integer: up to how many to generate

Returns

• xs:integer*: sequence

declare function this:gerrymander2($n as xs:integer) as xs:integer*
{
  if ($n lt 1) then errors:error("ML-BADARGS", ("n", $n))
  else (
    for $i in 1 to $n return (
      let $i2 := $i*$i
      let $divisors := this:factors($i2)
      return (
        min((
          for $d in $divisors return (
            (floor($d div 2)+1)*(floor($i2 div (2*$d))+1)
          ),
          if ($i mod 2 = 0) then (
            (($i div 2) + 1)*(($i div 2) + 1) - 1
          ) else ()
        )) cast as xs:integer
      )
    )
  )
}

Errors

ML-BADARGS if $n is less than 1

Function: inventory
declare function inventory($n as xs:integer) as xs:integer*

inventory()

First few terms:
0,1,1,0,2,2,2,0,3,2,4,1,1,0,4,4,4,1,4,0,5,5,4,1,6,2,1,0,6,7,5,1,6,3,3,1
A342585

Params

• n as xs:integer: up to how many to generate

Returns

• xs:integer*: sequence

declare function this:inventory($n as xs:integer) as xs:integer*
{
  if ($n lt 1) then errors:error("ML-BADARGS", ("n", $n)) 
  else if ($n = 1) then 0
  else (
    fold-left(2 to $n, (0, 0),
      function($seq as xs:integer*, $i as xs:integer) {
        let $last := $seq[last()]
        let $looking-for := head($seq)
        return (
          if ($last = 0) then (
            let $an := count(tail($seq)[. = 0])
            return (
              1,
              tail($seq),
              $an
            )
          ) else (
            let $an := count(tail($seq)[. = $looking-for])
            return (
              $looking-for + 1,
              tail($seq),
              $an
            )
          )
        )
      }
    )=>tail()
  )
}

Errors

ML-BADARGS if $n is less than 1

Function: pascal-row
declare function pascal-row($n as xs:integer) as xs:integer*

pascal-row()
The nth row of Pascal's triangle

Params

• n as xs:integer: which row

Returns

• xs:integer*: sequence

declare function this:pascal-row($n as xs:integer) as xs:integer*
{
  for $i in 0 to $n - 1 return util:binomial($n - 1, $i)
}

Function: pascal
declare function pascal($n as xs:integer) as xs:integer*

pascal()
Pascal's triangle read by rows
First few terms:
1,1,1,1,2,1,1,3,3,1,1,4,6,4,1,1,5,10,10,5,1,1,6,15,20,15,6,1,1,7,21,35,35,21

A007318

Params

• n as xs:integer: up to how many to generate

Returns

• xs:integer*: sequence

declare function this:pascal($n as xs:integer) as xs:integer*
{
  if ($n lt 1) then errors:error("ML-BADARGS", ("n", $n)) else
  (: n = Σ[i=1:k]i = (1 + k) * k/2 => k = -1/2 + √(2n + 1/4) :)
  (: ceiling to make sure we get enough :)
  let $k := ceiling(math:sqrt(2 * $n + 0.25) - 0.25) cast as xs:integer
  return (
    for $i in 1 to $k
    return this:pascal-row($i)
  )[position() <= $n]
}

Errors

ML-BADARGS if $n is less than 1

Function: katydid
declare function katydid($n as xs:integer) as xs:integer*

katydid()
Katydid sequence; closed under n=> 2n+2 and 7n+7
Only have the first 53 terms
A060031

Params

• n as xs:integer: up to how many to generate

Returns

• xs:integer*: sequence

declare function this:katydid($n as xs:integer) as xs:integer*
{
  if ($n lt 1) then errors:error("ML-BADARGS", ("n", $n))
  else $this:KATYDID[position() <= $n]
}

Errors

ML-BADARGS if $n is less than 1

Function: grasshopper
declare function grasshopper($n as xs:integer) as xs:integer*

grasshopper()
Grasshopper sequence; closed under n=> 2n+2 and 6n+6
Only have the first 100 terms

Params

• n as xs:integer: up to how many to generate

Returns

• xs:integer*: sequence

declare function this:grasshopper($n as xs:integer) as xs:integer*
{
  if ($n lt 1) then errors:error("ML-BADARGS", ("n", $n))
  else $this:GRASSHOPPER[position() <= $n]
}

Errors

ML-BADARGS if $n is less than 1

Function: thue-morse
declare function thue-morse($n as xs:integer) as xs:integer*

thue-morse()
Return the first n terms of the Thue-Morse sequence.
First few terms:
0,1,1,0,1,0,0,1,1,0,0,1,0,1,1,0,1,0,0,1,0,1,1,0,0,1,1,0,1,0,0,1,1,0,0,1,0,1

A010060

Params

• n as xs:integer: up to how many to generate

Returns

• xs:integer*: sequence

declare function this:thue-morse($n as xs:integer) as xs:integer*
{
  if ($n lt 1) then errors:error("ML-BADARGS", ("n", $n)) else
  let $generations := util:count-bits($n)
  return (
    fold-left(1 to $generations, 0,
      function($seq as xs:integer*, $i as xs:integer) {
        for $bit in $seq
        return
          if ($bit=0) then (0, 1) else (1, 0)
      }
    )[position() <= $n]
  )
}

Errors

ML-BADARGS if $n is less than 1

Function: thue-morse
declare function thue-morse($n as xs:integer, $k as xs:integer) as xs:integer*

thue-morse()
Return the first n terms of the Thue-Morse sequence over alphabet of size k
See A010060, A053838, etc.

First few terms of k=3
0,1,2,1,2,0,2,0,1,1,2,0,2,0,1,0,1,2,2,0,1,0,1,2,1,2,0,1,2,0,2,0,1,0,1,2,2

Params

• n as xs:integer: up to how many to generate
• k as xs:integer

Returns

• xs:integer*: sequence

declare function this:thue-morse($n as xs:integer, $k as xs:integer) as xs:integer*
{
  if ($n lt 1) then errors:error("ML-BADARGS", ("n", $n)) else
  let $generations := util:count-digits($n, $k)
  let $mappings := (
    array { 0 to $k - 1 },
    for $i in 1 to $k - 1 return array { this:rotate(0 to $k - 1, $i + 1) }
  )
  return (
    fold-left(1 to $generations, 0,
      function($seq as xs:integer*, $i as xs:integer) {
        for $digit in $seq
        return array:flatten($mappings[$digit + 1])
      }
    )[position() <= $n]
  )
}

Errors

ML-BADARGS if $n is less than 1

Function: feigenbaum1
declare function feigenbaum1($n as xs:integer) as xs:integer*

feigenbaum1()
Return first n terms in the Feigenbaum symbolic sequence.
First few terms:
1,0,1,1,1,0,1,0,1,0,1,1,1,0,1,1,1,0,1,1,1,0,1,0,1,0,1,1,1,0,1,0,1,0,1,1,1

A035263

Params

• n as xs:integer: up to how many to generate

Returns

• xs:integer*: sequence

declare function this:feigenbaum1($n as xs:integer) as xs:integer*
{
  if ($n lt 1) then errors:error("ML-BADARGS", ("n", $n))
  else if ($n = 1) then 1
  else if ($n = 2) then (1, 0)
  else (
    let $bits := util:count-bits($n - 2)
    return (
      fold-left(
        1 to $bits, (1, 0),
        function($fg as xs:integer*, $i as xs:integer) {
          $fg, $fg[position()<last()],
          if ($fg[last()]=1) then 0 else 1
        }
      )[position() <= $n]
    )
  )
}

Errors

ML-BADARGS if $n is less than 1

Function: ruler-sequence
declare function ruler-sequence($n as xs:integer) as xs:integer*

ruler-sequence()
First n terms of the ruler function.
nth term = exponent of highest power of 2 dividing n
e.g. 4th term = 2 because 2^2 mod 4 = 0
First few terms:
0,1,0,2,0,1,0,3,0,1,0,2,0,1,0,4,0,1,0,2,0,1,0,3,0,1,0,2,0,1,0,5,0,1,0,2,0,1

A007814

Params

• n as xs:integer: up to how many to generate

Returns

• xs:integer*: sequence

declare function this:ruler-sequence($n as xs:integer) as xs:integer*
{
  if ($n lt 1) then errors:error("ML-BADARGS", ("n", $n))
  else if ($n = 1) then 0
  else if ($n = 2) then (0, 1)
  else (
    let $bits := util:count-bits($n - 2)
    return (
      fold-left(
        1 to $bits, (0, 1),
        function($fg as xs:integer*, $i as xs:integer) {
          $fg, $fg[position()<last()],
          $fg[last()] + 1
        }
      )[position() <= $n]
    )
  )
}

Errors

ML-BADARGS if $n is less than 1

Function: ruler
declare function ruler($n as xs:integer) as xs:integer

ruler()
nth term of the ruler function = exponent of highest power of 2 dividing n
e.g. ruler(4) = 2 because 2^2 mod 4 = 0
A007814

Params

• n as xs:integer: which term to return

Returns

• xs:integer: integer

declare function this:ruler($n as xs:integer) as xs:integer
{
  if ($n lt 1) then errors:error("ML-BADARGS", ("n", $n))
  else if ($n = 1) then 0
  else if ($n = 2) then 1
  else this:ruler-sequence($n)[last()]
}

Errors

ML-BADARGS if $n is less than 1

Function: debruijn
declare function debruijn($n as xs:integer) as xs:integer*

debruijn()
Compute a binary deBruijn sequence of order n
Cyclic sequence where every possible sequence of length n occurs exactly
once as a contiguous subsequence in the cycle.

Params

• n as xs:integer the order

Returns

• xs:integer*: sequence

declare function this:debruijn(
  $n as xs:integer
) as xs:integer*
{
  if ($n lt 1) then errors:error("ML-BADARGS", ("n", $n))
  else this:lempel($n, ())
}

Errors

ML-BADARGS if $n is less than 1

Function: debruijn
declare function debruijn($alphabet as item()*,
  $n as xs:integer)

debruijn()
Construct a deBruijn sequence of order n using the given alphabet
Cyclic sequence where every possible sequence of length n occurs exactly
once as a contiguous subsequence in the cycle.

Example: debruijn(("a","b"), 3) = ("a","a","a","b","a","b","b","b")
Example: debruijn(2, 3) = (0,0,0,1,0,1,1,1)

Params

• alphabet as item()*: either the sequence of alphabet values to use or an integer indicating the numerical range to use, e.g. 2 => (0,1)
• n as xs:integer the order

declare function this:debruijn(
  $alphabet as item()*,
  $n as xs:integer
)
{
  if ($n lt 1) then errors:error("ML-BADARGS", ("n", $n)) else
  let $alphabet :=
    if ($alphabet instance of xs:integer) then 0 to $alphabet - 1
    else $alphabet
  let $k := count($alphabet)
  let $a := sparse:array()=>sparse:put(1,0)
  let $a := this:ruskey($a, $k, $n, 1, 1, 0)
  let $res := $a=>map:get("seq")
  for $ix in $res return $alphabet[$ix + 1]
}

Errors

ML-BADARGS if $n is less than 1

Function: prime-sequence
declare function prime-sequence($n as xs:integer,
  $k as xs:integer,
  $q as xs:integer,
  $r as xs:integer) as xs:integer*

prime-sequence()
A sequence of digits in range [0, k-1] from the division of two
mutually prime numbers (use actual primes or rand:mutual-prime-pair())
If the numbers aren't mutually prime, the sequence gets boring.
You want q > r, in general.

Params

• n as xs:integer: length of sequence
• k as xs:integer: number range (e.g. 2 = binary 0 and 1)
• q as xs:integer: first number
• r as xs:integer: second number

Returns

• xs:integer*: sequence

declare function this:prime-sequence(
  $n as xs:integer,
  $k as xs:integer,
  $q as xs:integer,
  $r as xs:integer
) as xs:integer*
{
  if ($n lt 1) then errors:error("ML-BADARGS", ("n", $n)) else
  for $i in 1 to $n
  let $dx := (($k*$q - $r) * $i) mod ($k*$q)
  return
    (
      for $i in 1 to $k - 1
      return if ($dx < $i*$q) then $i - 1 else (),
      $k
    )[1]
}

Errors

ML-BADARGS if $n is less than 1

Function: permutations
declare function permutations($input as item()*) as array(*)*

permutations()
All permutations of the input. Very naive algorithm: Don't use
on large sequences. Doesn't eliminate duplicates if there are
duplicate values within the input.

Params

• input as item()*: the input sequence of values to permute

Returns

• array(*)*: sequence of arrays, one per permutation

declare function this:permutations($input as item()*) as array(*)*
{
  let $N := count($input)
  return (
    switch ($N)
    case 0 return $input
    case 1 return array {$input}
    case 2 return (
      array {$input},
      array {reverse($input)}
    )
    case 3 return (
      array {$input},
      array {($input[1], $input[3], $input[2])},
      array {($input[2], $input[1], $input[3])},
      array {($input[2], $input[3], $input[1])},
      array {($input[3], $input[1], $input[2])},
      array {reverse($input)}
    )
    default return (
      for $i in 1 to $N return (
        for $c in this:permutations($input[position()!=$i])
        return
          array {(
            $input[$i],
            util:array-values($c)
         )}
      )
    )
  )
}

Function: foreach-permutation
declare function foreach-permutation($input as item()*,
  $f as function(item()*) as item()*) as item()*

foreach-permutation()
Execute a function over all permutations of the input. Very naive algorithm:
Don't use on large sequences. Doesn't eliminate duplicates if there are
duplicate values within the input.

Params

• input as item()*: the input sequence of values to permute
• f as function(item()*)asitem()*: the function to call for each permutation

Returns

• item()*: values of function for each permutation

declare function this:foreach-permutation(
  $input as item()*,
  $f as function(item()*) as item()*
) as item()*
{
  let $N := count($input)
  return (
    switch ($N)
    case 0 return $f($input)
    case 1 return $f($input)
    case 2 return (
      $f($input),
      $f(reverse($input))
    )
    case 3 return (
      $f($input),
      $f(($input[1], $input[3], $input[2])),
      $f(($input[2], $input[1], $input[3])),
      $f(($input[2], $input[3], $input[1])),
      $f(($input[3], $input[1], $input[2])),
      $f(reverse($input))
    )
    default return (
      for $i in 1 to $N return (
        for $c in this:permutations($input[position()!=$i])
        return
          $f((
            $input[$i],
            util:array-values($c)
         ))
      )
    )
  )
}

Function: powerset
declare function powerset($input as item()*) as array(*)*

powerset()
Get all combinations of the input.

Params

• input as item()*: the input sequence of values to get combinations from

Returns

• array(*)*: sequence of arrays, one per combination

declare function this:powerset($input as item()*) as array(*)*
{
  let $N := count($input)
  for $i in 0 to (math:pow(2, count($input)) cast as xs:integer) - 1
  return
    array {
      for $j in 1 to $N
      where $i idiv math:pow(2, $j - 1) mod 2 = 1
      return $input[$j]
    }
}

Function: foreach-combination
declare function foreach-combination($input as item()*,
  $f as function(item()*) as item()*) as item()*

foreach-combination()
Execute a function over all combinations of the input.

Params

• input as item()*: the input sequence of values to get combinations from
• f as function(item()*)asitem()*: the function to call for each combination

Returns

• item()*: values of function over each combinatuion

declare function this:foreach-combination(
  $input as item()*,
  $f as function(item()*) as item()*
) as item()*
{
  let $N := count($input)
  for $i in 0 to (math:pow(2, count($input)) cast as xs:integer) - 1
  return (
    $f(
      for $j in 1 to $N
      where $i idiv math:pow(2, $j - 1) mod 2 = 1
      return $input[$j]
    )
  )
}

Function: k-combinations
declare function k-combinations($input as item()*,
  $k as xs:integer) as array(*)*

k-combinations()
Members of the powerset of $items that have cardinality $k

Params

• input as item()*: the input sequence of values to get combinations from
• k as xs:integer: cardinality constraint

Returns

• array(*)*: sequence of arrays, one per combination of cardinality $k

declare function this:k-combinations(
  $input as item()*,
  $k as xs:integer
) as array(*)*
{
  let $N := count($input)
  for $i in 0 to (math:pow(2, count($input)) cast as xs:integer) - 1
  let $members := (
    for $j in 1 to $N
    where $i idiv math:pow(2, $j - 1) mod 2 = 1
    return $input[$j]
  )
  where count($members)=$k
  return array { $members }
}

Function: foreach-k-combination
declare function foreach-k-combination($input as item()*,
  $k as xs:integer,
  $f as function(item()*) as item()*) as item()*

foreach-k-combination()
Execute a function over all combinations of the input of cardinality k.

Params

• input as item()*: the input sequence of values to get combinations from
• k as xs:integer: cardinality constraint
• f as function(item()*)asitem()*: the function to call for each combination

Returns

• item()*: values of function over each combinary

declare function this:foreach-k-combination(
  $input as item()*,
  $k as xs:integer,
  $f as function(item()*) as item()*
) as item()*
{
  let $N := count($input)
  for $i in 0 to (math:pow(2, count($input)) cast as xs:integer) - 1
  let $members := (
    for $j in 1 to $N
    where $i idiv math:pow(2, $j - 1) mod 2 = 1
    return $input[$j]
  )
  where count($members)=$k
  return $f($members)
}

Function: rotate
declare function rotate($seq as item()*,
  $k as xs:integer) as item()*

rotate()
Rotate a sequence by k slots
Example this:rotate((1,0,1,1,1,0,0,0), 2) = (1,1,1,0,0,0,1,0)

Params

• seq as item()*: input sequence
• k as xs:integer: position in input sequence to start new sequence

Returns

• item()*: rotated sequence

declare function this:rotate(
  $seq as item()*,
  $k as xs:integer
) as item()*
{
  $seq[position()>=$k], $seq[position() < $k]
}

Function: lodumo
declare function lodumo($seq as xs:integer*,
  $m as xs:integer) as xs:integer*

lodumo()
Modify the input sequence s(n) where a(n) = smallest non-negative integer
not yet in list such that a(n) = s(n) mod m

See: https://oeis.org/wiki/Lodumo_transform

Params

• seq as xs:integer*: input sequence
• m as xs:integer: modulo

Returns

• xs:integer*: lodumo of input sequence

declare function this:lodumo(
  $seq as xs:integer*,
  $m as xs:integer
) as xs:integer*
{
  if ($m lt 1) then errors:error("ML-BADARGS", ("m", $m))
  else if (some $a in $seq satisfies $a < 0) then errors:error("ML-BADARGS", ("seq", $seq))
  else (
    fold-left(
       $seq, (),
       function($output as xs:integer*, $a as xs:integer) {
         $output,
         util:while(
           function($an as xs:integer) { $an = $output },
           function($an as xs:integer) { $an + $m },
           $a mod $m
         )
       }
    )
  )
}

Function: primes
declare function primes($n as xs:integer) as xs:integer*

primes()
Set of primes less than or equal $n

Params

• n as xs:integer: maximum prime

Returns

• xs:integer*: sequence of primes

declare function this:primes($n as xs:integer) as xs:integer*
{
  if ($n <= $util:CANVAS-PRIMES[last()]) then $util:CANVAS-PRIMES[. <= $n]
  else (
    $util:CANVAS-PRIMES,
    (: Brute force, may take a while :)
    for $i in $util:CANVAS-PRIMES[last()] + 1 to $n
    where util:is-prime($i)
    return $i
  )
}

Function: prime-factors
declare function prime-factors($n as xs:integer) as xs:integer*

prime-factors()
Prime factors of the number, excluding 1, in ascending order
So factor(1) = (); factor(12) = (2, 2, 3)
n > 0

Params

• n as xs:integer: number to factor

Returns

• xs:integer*: sequence of prime factors

declare function this:prime-factors($n as xs:integer) as xs:integer*
{
  if ($n lt 1) then errors:error("ML-BADARGS", ("n", $n))
  else if ($n = 1) then ()
  else if ($n = $util:CANVAS-PRIMES) then $n
  else (
    let $factors :=
      let $primes := reverse(this:primes(math:sqrt($n) cast as xs:integer))
      return (
        fold-left($primes, $n,
          function($factors as xs:integer*, $prime as xs:integer) as xs:integer* {
            util:while(
              function($residue as xs:integer*) { head($residue) mod $prime = 0 },
              function($residue as xs:integer*) { head($residue) idiv $prime, $prime, tail($residue) },
              head($factors)
            ),
            tail($factors)
          }
        )
      )
    for $f in $factors order by $f ascending where $f!=1 return $f
  )
}

Errors

ML-BADARGS if $n is less than 1

Function: factors
declare function factors($n as xs:integer) as xs:integer*

factors()
All distinct factors of the number, including 1
n >= 0

Params

• n as xs:integer: number to factor

Returns

• xs:integer*: set of distinct factors

declare function this:factors($n as xs:integer) as xs:integer*
{
  if ($n < 0) then errors:error("ML-BADARGS", ("n", $n))
  else if ($n = 0) then 0
  else (
    for $f in 
    distinct-values(
      for $p in this:powerset(this:prime-factors($n))
      return fold-left(array:flatten($p), 1, function($prod as xs:integer, $f as xs:integer) {$prod*$f})
    )
    order by $f ascending
    return $f
  )
}

Errors

ML-BADARGS if $n is less than zero

Function: continued-fraction
declare function continued-fraction($n as xs:integer,
  $x as xs:double) as xs:integer*

continued-fraction()
Generated the continued fraction sequence for the given number.
$x=π => 3, 7, 15, 1, 292, 1, 1, 1, 2, 1, 3,...
$x=e => 2, 1, 2, 1, 1, 4, 1, 1, 6, 1, 1, 8, ...
If $x is rational this is finite but this function will return 0s
to fill out the sequence (technically wrong: should be INF)

Params

• n as xs:integer: now many terms to generate
• x as xs:double: the number to use

Returns

• xs:integer*

declare function this:continued-fraction(
  $n as xs:integer,
  $x as xs:double
) as xs:integer*
{
  fold-left(1 to $n, ($x, ()),
    function($data as xs:numeric*, $i as xs:integer) {
      let $fraction := head($data) - floor(head($data))
      return (
        if ($fraction = 0) then 0
        else 1 div $fraction
      ),
      tail($data),
      floor(head($data)) cast as xs:integer
    }
  )=>tail()
}

Function: convergents
declare function convergents($seq as xs:integer*) as array(xs:integer)*

convergents()
Get the numerator/denominator pairs of the convergents of the sequence,
taken as a continued fraction. Each pair $pair(1)/$pair(2) represents
an approximation of the value represented by the sequence.

Params

• seq as xs:integer*: sequence to interpret as continued fraction sequence

Returns

• array(xs:integer)*: sequence of the successive rational approximations

declare function this:convergents($seq as xs:integer*) as array(xs:integer)*
{
  let $h :=
    fold-left(1 to count($seq), (0, 1),
      function($h as xs:integer*, $i as xs:integer) {
        $h,
        $seq[$i]*$h[2 + $i - 1] + $h[2 + $i - 2]
      }
    )=>tail()=>tail()
  let $k :=
    fold-left(1 to count($seq), (1, 0),
      function($k as xs:integer*, $i as xs:integer) {
        $k,
        $seq[$i]*$k[2 + $i - 1] + $k[2 + $i - 2]
      }
    )=>tail()=>tail()
  for $n in 1 to count($seq) return (
    array { $h[$n], $k[$n] }
  )
}

Function: numerators
declare function numerators($seq as xs:integer*) as xs:integer*

numerators()
Get the numerators of the convergents of the sequence,
taken as a continued fraction.

Params

• seq as xs:integer*: sequence to interpret as continued fraction sequence

Returns

• xs:integer*: sequence of the numerators of successive rational approximations

declare function this:numerators($seq as xs:integer*) as xs:integer*
{
  fold-left(1 to count($seq), (0, 1),
    function($h as xs:integer*, $i as xs:integer) {
      $h,
      $seq[$i]*$h[2 + $i - 1] + $h[2 + $i - 2]
    }
  )=>tail()=>tail()
}

Function: denominators
declare function denominators($seq as xs:integer*) as xs:integer*

denominators()
Get the denominators of the convergents of the sequence,
taken as a continued fraction.

Params

• seq as xs:integer*: sequence to interpret as continued fraction sequence

Returns

• xs:integer*: sequence of the denominators of successive rational approximations

declare function this:denominators($seq as xs:integer*) as xs:integer*
{
  fold-left(1 to count($seq), (1, 0),
    function($k as xs:integer*, $i as xs:integer) {
      $k,
      $seq[$i]*$k[2 + $i - 1] + $k[2 + $i - 2]
    }
  )=>tail()=>tail()
}

Function: μ
declare function μ($n as xs:integer) as xs:integer

μ()
Möbius μ function:
0: if n has 1 or more repeated prime factors
1: if n = 1
(-1)^k: if n is product of k distinct primes

Params

• n as xs:integer: number to calculate Möbius μ function for

Returns

• xs:integer: μ(n)

declare function this:μ($n as xs:integer) as xs:integer
{
  if ($n = 1) then 1 else (
    let $factors := this:prime-factors($n)
    let $k := count($factors)
    let $d := count(distinct-values($factors))
    return (
      if ($k = $d) then (
        if ($k mod 2 = 0) then 1 else -1
      ) else 0
    )
  )
}

Function: mobius
declare function mobius($seq as xs:integer*) as xs:integer*

mobius()
Möbius transform b[n] = Σ (μ(n/d)a[d])
d|n d divides into n

Params

• seq as xs:integer*: sequence to apply transform to

Returns

• xs:integer*: Möbius transform of the sequence

declare function this:mobius($seq as xs:integer*) as xs:integer*
{
  for $n in 1 to count($seq)
  let $factors := distinct-values(
    for $p in this:powerset(this:prime-factors($n))
    return fold-left(array:flatten($p), 1, function($prod as xs:integer*, $f as xs:integer) {$prod*$f})
  )
  return (
    if ($seq[$n] = 0) then 0
    else sum(for $d in $factors return this:μ($n idiv $d) * $seq[$d])
  )
}

Function: binomial
declare function binomial($seq as xs:integer*) as xs:integer*

binomial()
n
Binomial transform of a sequence: b[n] = Σ (n choose k) a[k]
k=0

Params

• seq as xs:integer*: sequence to apply transform to

Returns

• xs:integer*: binomial transform of the sequence

declare function this:binomial($seq as xs:integer*) as xs:integer*
{
  for $n in 0 to count($seq) - 1 return (
    sum(
      for $k in 0 to $n return util:binomial($n, $k)*$seq[$k + 1]
    ) cast as xs:integer
  )  
}

Function: ordinal
declare function ordinal($seq as xs:integer*) as xs:integer*

ordinal()

Ordinal transformation: b[n] = count of number of times a[n] appears in a[i]
sequence up to and include a[n]

Params

• seq as xs:integer*: sequence to apply transform to

Returns

• xs:integer*: ordinal transform of the sequence

declare function this:ordinal($seq as xs:integer*) as xs:integer*
{
  for $a at $i in $seq return count($seq[position() <= $i and .=$a])
}

Function: mirror
declare function mirror($seq as xs:integer*) as xs:integer*

mirror()

Mirror transformation: reversed negation of first half of sequence
followed by first half of sequence

Params

• seq as xs:integer*: sequence to apply transform to

Returns

• xs:integer*: mirror transform of the sequence

declare function this:mirror($seq as xs:integer*) as xs:integer*
{
  let $n := count($seq)
  return (
    reverse($seq[position() <= $n div 2])!(-.),
    $seq[position() <= $n div 2]
  )
}

Function: alternate
declare function alternate($seq as xs:integer*) as xs:integer*

alternate()

Alternation transformation: b[n] = a[n] is n mod 2 = 1 else -a[n]

Params

• seq as xs:integer*: sequence to apply transform to

Returns

• xs:integer*: alternate transform of the sequence

declare function this:alternate($seq as xs:integer*) as xs:integer*
{
  for $val at $i in $seq return (
    if ($i mod 2 = 1) then abs($val) else -abs($val)
  )
}

Function: shift
declare function shift($seq as xs:integer*) as xs:integer*

shift()

Shift transformation: b[n] = a[n]+1
A lot of the sequences produce 0 values which we may not want.

Params

• seq as xs:integer*: sequence to apply transform to

Returns

• xs:integer*: shift transform of the sequence

declare function this:shift($seq as xs:integer*) as xs:integer*
{
  $seq!(. + 1)
}

Function: modix
declare function modix($seq as xs:integer*,
  $m as xs:integer) as xs:integer*

modix()
Modify the input sequence s(n) where a(n) = util:modix(s(n), m)
i.e. s(n) mod m shifted to 1 to m range (maps negative values as well)

Params

• seq as xs:integer*: input sequence
• m as xs:integer: modulo

Returns

• xs:integer*: modix of input sequence

declare function this:modix(
  $seq as xs:integer*,
  $m as xs:integer
) as xs:integer*
{
  if ($m lt 1) then errors:error("ML-BADARGS", ("m", $m))
  else (
    $seq!util:modix(., $m)
  )
}

Function: modp
declare function modp($seq as xs:integer*,
  $m as xs:integer) as xs:integer*

modp()
Modify the input sequence s(n) where a(n) = s(n) if s(n) mod m != 0, 0
otherwise.

Params

• seq as xs:integer*: input sequence
• m as xs:integer: modulo

Returns

• xs:integer*: input sequence when value is positive mod m, 0 when not

declare function this:modp(
  $seq as xs:integer*,
  $m as xs:integer
) as xs:integer*
{
  if ($m lt 1) then errors:error("ML-BADARGS", ("m", $m))
  else (
    for $a in $seq return if ($a mod $m = 0) then 0 else $a
  )
}

Function: modz
declare function modz($seq as xs:integer*,
  $m as xs:integer) as xs:integer*

modz()
Modify the input sequence s(n) where a(n) = s(n) if s(n) mod m = 0, 0
otherwise.

Params

• seq as xs:integer*: input sequence
• m as xs:integer: modulo

Returns

• xs:integer*: input sequence when value is zero mod m, 0 when not

declare function this:modz(
  $seq as xs:integer*,
  $m as xs:integer
) as xs:integer*
{
  if ($m lt 1) then errors:error("ML-BADARGS", ("m", $m))
  else (
    for $a in $seq return if ($a mod $m = 0) then $a else 0
  )
}

Function: unmodp
declare function unmodp($seq as xs:integer*,
  $m as xs:integer) as xs:integer*

unmodp()
Modify the input sequence s(n) where a(n) = $n if s(n) mod m != 0, 0
otherwise.

Params

• seq as xs:integer*: input sequence
• m as xs:integer: modulo

Returns

• xs:integer*: index when value is positive mod m, 0 when not

declare function this:unmodp(
  $seq as xs:integer*,
  $m as xs:integer
) as xs:integer*
{
  if ($m lt 1) then errors:error("ML-BADARGS", ("m", $m))
  else (
    for $i in 1 to count($seq) return if ($seq[$i] mod $m = 0) then 0 else $i
  )
}

Function: unmodz
declare function unmodz($seq as xs:integer*,
  $m as xs:integer) as xs:integer*

unmodz()
Modify the input sequence s(n) where a(n) = $n if s(n) mod m = 0, 0
otherwise.

Params

• seq as xs:integer*: input sequence
• m as xs:integer: modulo

Returns

• xs:integer*: index when value is zero mod m, 0 when not

declare function this:unmodz(
  $seq as xs:integer*,
  $m as xs:integer
) as xs:integer*
{
  if ($m lt 1) then errors:error("ML-BADARGS", ("m", $m))
  else (
    for $i in 1 to count($seq) return if ($seq[$i] mod $m = 0) then $i else 0
  )
}

Function: compose
declare function compose($seq1 as xs:integer*,
  $seq2 as xs:integer*,
  $op as function(xs:integer,xs:integer) as xs:integer) as xs:integer*

compose()
Compose two integer sequences using the particular operation.
If the sequences are different lengths, loop back to start of shorter seq.
e.g. compose((1,2,3,4),(-1,1),function($a,$b) {$a * $b}) = (-1,2,-3,4)

Params

• seq1 as xs:integer*: one sequence
• seq2 as xs:integer*: other sequence
• op as function(xs:integer,xs:integer)asxs:integer composition function

Returns

• xs:integer*

declare function this:compose(
  $seq1 as xs:integer*,
  $seq2 as xs:integer*,
  $op as function(xs:integer,xs:integer) as xs:integer
) as xs:integer*
{
  let $n1 := count($seq1)
  let $n2 := count($seq2)
  for $i in 1 to (max(($n1,$n2)) cast as xs:integer)
  let $a := $seq1[util:modix($i, $n1)]
  let $b := $seq2[util:modix($i, $n2)]
  return $op($a, $b)
}

Function: sequence
declare function sequence($name as xs:string,
  $n as xs:integer) as xs:integer*

sequence()
Get an integer sequence by name. The name is either the simple name of
the sequence, or the sequence modified.

Bare names:
fibonacci (Fibonacci sequence), fibonacci-word, padovan, fractal, golomb,
lazy-caterer, jacobsthal, jacobsthal-lucas, recaman, gerrymander2,
inventory, pascal, katydid, grasshopper, thue-morse, feigenbaum1,
ruler (Ruler sequence), primes, counter

Parameterized names:
thue-morse@k = thue-morse(n,k)
prime-sequence@k@q@r = prime-sequence(n, k, q, r)
continued-fraction@d = continued-fraction(n, d)

Operators: apply the operator to the sequence so far
mobius = seq=>mobius()
binomial = seq=>binomial()
ordinal = seq=>ordinal()
numerators = seq=>numerators()
denominators = seq=>denominators()
mirror = seq=>mirror()
alternate = seq=>alternate()
shift = seq=>shift()

Parameterized operators:
rotate@k = seq=>rotate(k)
lodumo@m = seq=>lodumo(m)
modix@m = seq=>modix(m)
modp@m = seq=>modp(m)
modz@m = seq=>modz(m)
unmodp@m = seq=>unmodp(m)
unmodz@m = seq=>unmodz(m)

seq1+operator = apply operator to seq1
seq1+seq1 = seq1 followed by seq2 (so sequence length will be 2*$n)

Params

• name as xs:string: name of sequence
• n as xs:integerame: name of sequence: length of each component sequence

Returns

• xs:integer*: sequence

declare function this:sequence(
  $name as xs:string,
  $n as xs:integer
) as xs:integer*
{
  fold-left(tokenize($name, "[+]")[. ne ""], (),
    function($seq as xs:integer*, $part as xs:string) {
      switch($part)
      case "none" return $seq
      case "counter" return ($seq, 1 to $n)
      case "fibonacci" return ($seq, this:fibonacci-sequence($n))
      case "fibonacci-word" return ($seq, this:fibonacci-word($n))
      case "padovan" return ($seq, this:padovan($n))
      case "fractal" return ($seq, this:fractal($n))
      case "golomb" return ($seq, this:golomb($n))
      case "lazy-caterer" return ($seq, this:lazy-caterer($n))
      case "jacobsthal" return ($seq, this:jacobsthal($n))
      case "jacobsthal-lucas" return ($seq, this:jacobsthal-lucas($n))
      case "recaman" return ($seq, this:recaman($n))
      case "gerrymander2" return ($seq, this:gerrymander2($n))
      case "inventory" return ($seq, this:inventory($n))
      case "pascal" return ($seq, this:pascal($n))
      case "katydid" return ($seq, this:katydid($n))
      case "grasshopper" return ($seq, this:grasshopper($n))
      case "thue-morse" return ($seq, this:thue-morse($n))
      case "feigenbaum1" return ($seq, this:feigenbaum1($n))
      case "ruler" return ($seq, this:ruler-sequence($n))
      case "primes" return ($seq, this:primes($n))
      case "mobius" return $seq=>this:mobius()
      case "binomial" return $seq=>this:binomial()
      case "ordinal" return $seq=>this:ordinal()
      case "numerators" return $seq=>this:numerators()
      case "denominators" return $seq=>this:denominators()
      case "mirror" return $seq=>this:mirror()
      case "alternate" return $seq=>this:alternate()
      case "shift" return $seq=>this:shift()
      default return ( (: More parsing required :)
        let $pieces := tokenize($part,"@")[. ne ""]
        return switch ($pieces[1])
        case "thue-morse" return ($seq, this:thue-morse($n, xs:integer($pieces[2])))
        case "prime-sequence" return ($seq, this:prime-sequence($n, xs:integer($pieces[2]), xs:integer($pieces[3]), xs:integer($pieces[4])))
        case "rotate" return $seq=>this:rotate(xs:integer($pieces[2]))
        case "lodumo" return $seq=>this:lodumo(xs:integer($pieces[2]))
        case "modix" return $seq=>this:modix(xs:integer($pieces[2]))
        case "continued-fraction" return ($seq, this:continued-fraction($n, this:double($pieces[2])))
        case "modp" return $seq=>this:modp(xs:integer($pieces[2]))
        case "modz" return $seq=>this:modz(xs:integer($pieces[2]))
        case "unmodp" return $seq=>this:unmodp(xs:integer($pieces[2]))
        case "unmodz" return $seq=>this:unmodz(xs:integer($pieces[2]))
        default return $seq
      )
    }
  )
}