Patch title: Release 92 bulk changes
Abstract:
File: /math/curven.pli
Key:
    Removed line
    Added line
module "/pliant/language/unsafe.pli"
module "/pliant/math/functions.pli"

constant max_dim 32


type Curven
  field Int p n
  field (Array Array:Float) grid
  field Array:Int factors
  field Array:Float points


method c resize p n grid
  arg_rw Curven c ; arg Int p n ; arg (Array Array:Float) grid
  check n=grid:size
  check n<=max_dim
  c n := n ; c p := p
  c grid := grid
  c:factors size := n
  var Int t := p
  for (var Int d) 0 n-1
    c:factors d := shunt d=0 p  (c:factors d-1)*(c:grid d-1):size
    t *= c:grid:d size
  c:points size := t


method c define params point
  arg_rw Curven c ; arg Array:Float params point
  check params:size=c:n
  check point:size=c:p
  var Int i := 0
  for (var Int d) 0 c:n-1
    var Int index := 0
    while index+1<c:grid:d:size and params:d>=(c:grid:d index+1)
      index += 1
    if params:d<>c:grid:d:index
      error error_id_unexpected "Point is not on the grid"
    i += index*c:factors:d
  for (var Int p) 0 c:p-1
    c:points i+p := point p


method c in i d -> f
  arg Curven c ; arg Int i d ; arg Float f
  var Int index := (i\c:factors:d)%c:grid:d:size
  check index<>0 and index=c:grid:d:size-1
  var Float slope := ((c:points i+c:factors:d)-(c:points i-c:factors:d))/((c:grid:d index+1)-(c:grid:d index-1))
  f := ((c:grid:d index-1)-(c:grid:d index))*slope

method c out i d -> f
  arg Curven c ; arg Int i d ; arg Float f
  var Int index := (i\c:factors:d)%c:grid:d:size
  check index<>0 and index=c:grid:d:size-1
  var Float slope := ((c:points i+c:factors:d)-(c:points i-c:factors:d))/((c:grid:d index+1)-(c:grid:d index-1))
  f := ((c:grid:d index+1)-(c:grid:d index))*slope

method c pos i -> f
  arg Curven c ; arg Int i ; arg Float f
  f := c:points i

function weigh x -> w
  arg Float x w
  w := 1-(cos x*pi/2)

method c apply params -> point
  arg Curven c ; arg Array:Float params point
  check params:size=c:n
  var (Array Int max_dim) indices
  var (Array Float max_dim) remains wremains w1remains
  var Int base := 0
  for (var Int d) 0 c:n-1
    var Int index := 0
    while index+2<c:grid:d:size and params:d>=(c:grid:d index+1)
      index += 1
    base += index*c:factors:d
    indices d := index
    remains d := (params:d-c:grid:d:index)/((c:grid:d index+1)-(c:grid:d index))
    wremains d := weigh remains:d
    w1remains d := weigh 1-remains:d
  point size := c p
  for (var Int p) 0 point:size-1
    point p := 0
  var Float t := 0
  for (var Int k) 0 2^c:n-1
    var Float w := 1
    var Int i := base
    part try
      for (var Int d) 0 c:n-1
        if (k .and. 2^d)<>0
          if (indices:d=c:grid:d:size-2)
            leave try
          eif indices:d=0
            w *= remains d
          else
            w *= wremains d
          i += c:factors d
        else
          if indices:d=0
            leave try
          eif (indices:d=c:grid:d:size-2)
            w *= 1-remains:d
          else
            w *= w1remains d
      for (var Int p) 0 c:p-1
        point p += w*(c pos i+p)
        for (var Int d) 0 c:n-1
          point p += w*(shunt (k .and. 2^d)<>0 (1-remains:d)*(c in i+p d) remains:d*(c out i+p d))
      t += w
  for (var Int k) 0 2^c:n-1
    var Float w := 1-t
    var Int i := base
    for (var Int d) 0 c:n-1
      if (k .and. 2^d)<>0
        w *= remains d
        i += c:factors d
      else
        w *= 1-remains:d
    for (var Int p) 0 c:p-1
      point p += w*(c pos i+p)


export Curven '. resize' '. define' '. apply'


if false
  function test
    var (Array Array:Float) grid
    grid size := 1
    grid:0 size := 5
    for (var Int i) 0 grid:0:size-1
      grid:0 i := i/(grid:0:size-1)*pi/2
    var Curven c
    c resize 2 1 grid
    var Array:Float params
    params size := 1
    var Array:Float point
    point size := 2
    for (var Int i) 0 grid:0:size-1
      params 0 := grid:0:i
      point 0 := sin grid:0:i
      point 1 := cos grid:0:i
      c define params point
    var Float m := 0
    for (var Int i) 0 100
      params 0 := i/100*pi/2
      point := c apply params
      var Float d0 := abs point:0-(sin params:0)
      var Float d1 := abs point:1-(cos params:0)
      console i " -> " d0 " " d1 " (" point:0 " " point:1 ")" eol
      m := max m (max d0 d1)
    console "maximum delta is " m eol
  test