best rational approximations, incrementally

README.md

Calculating best rational approximations incrementally

So best rational approximations come from a continued fraction expansion.

A continued fraction expansion comes from the loop,

function* continuedFraction(x: number): Iterable<number> {
  let n = Math.floor(x)
  yield n
  x -= n
  while (x > 0) {
    x = 1 / x
    n = Math.floor(x)
    x -= n
    yield n
  }
}

For pi this will yield [3, 7, 15, 1, 292, 1, 1, 1, ...] meaning that

ts-node layout-div.ts '= pi + 3 / 1 + 7 / 1 + 15 / 1 + 1 / 1 + 292 / 1 + 1 ...'
Pretty-printing pi = 3 + 1/(7 + 1/(15 + 1/(1 + 1/(292 + 1/(1 + ...)))))

                                                  1
pi  =   3  +   -----------------------------------------------------------------------
                                                     1
                  7  +   ----------------------------------------------------------
                                                         1
                            15  +   --------------------------------------------
                                                             1
                                       1  +   -------------------------------
                                                                 1
                                                 292  +   ----------------
                                                             1  +   ...

and it turns out this is the right way to manufacture best rational approximations out of a number, by truncating this sequence with some final integer.

Before I used to get very frustrated, though.

I used to get frustrated because the obvious way to calculate this is to start from the innermost number, and build a rational moving backwards out of the fraction. But then you can't easily compute the next best rational number, without redoing all of your work. Frustrating!

Today I thought that it might secretly be a matrix problem.

Suppose the rational m/n is stored as the vector [m; n]. Then the fundamental process that we are considering,

(m/n)  ---->  k + 1 / (m/n)
                 = k + (n / m)
                 = (k m + n) / m

becomes a 2x2 matrix:

k + 1/(...) =  [ k, 1 ]
               [ 1, 0 ]

And matrix algebra is associative, you can do it from the beginning to the end rather than from the end to the beginning like normal.

Some further insight reveals that actually the matrices you get when you do this are:

[ ratApprox(x, n)  ratApprox(x, n - 1) ].

This is easiest to see by just operating on this with vectors [1; 0] and [0; 1].

The first one is a representation for a rational number infinity, whence we would find that

Pretty-printing a + 1/(b + 1/∞) = a + 1/b

                1                          1
a  +   --------------------  =   a  +   -------
                    1                      b
          b  +   -------
                    ∞

as it should be, truncating the sequence at b. A similar analysis for the second vector, which is 0, constructs

Pretty-printing a + 1/(b + 1/0) = a + 1/(b + ∞) = a

                1                             1
a  +   --------------------  =   a  +   --------------  =   a
                    1                      b  +   ∞
          b  +   -------
                    0

truncating the sequence at a.

So our first rational approximation is [floor(x), 1] with a prior rational approximation of infinity, [1, 0] in this code.

And then every other rational approximation is constructed from the current digit K of the expansion, plus the last two rational approximations, by matrix multiplication:

[ m, a ] * [ k, 1 ] =  [ k m + a,  m]
[ n, b ]   [ 1, 0 ]    [ k n + b,  n]

That is what the above code does.

License

Feel free to use any of this under the MPLv2, or contact me if you want to embed it in a proprietary or GPLed file.

braps.ts

function rat(x: number, y: number): [number, number] {
  return [x, y]
}

function* bestRationalApproximations(x: number): Iterable<[number, number]> {
  let last = rat(1, 0)
  let curr = rat(Math.floor(x), 1)
  x -= curr[0]
  while (x > 0) {
    yield curr
    x = 1 / x
    const n = Math.floor(x)
    x -= n
    const next = rat(n * curr[0] + last[0], n * curr[1] + last[1])
    last = curr
    curr = next
  }
  yield curr;
}
function* take<x>(n: number, seq: Iterable<x>): Iterable<x> {
  for (const x of seq) {
    if (n-- <= 0) {
      return
    }
    yield x
  }
}
function approxErrString(n: number): string {
  if (n === 0) {
    return 'exact';
  }
  n = Math.abs(100 * n)
  if (n > 1) {
    return Math.round(1000 * n) / 1000 + '%'
  }
  let out = '0.'
  n *= 10
  while (n < 1) {
    out += '0'
    n *= 10
  }
  return `${out}${Math.round(n * 10 ** Math.max(0, 4 - out.length))}%`
}
function tenTerms(x: number): string[] {
  return [...take(10, bestRationalApproximations(x))].map(
    ([m, n]) => `${m}/${n} : ~${approxErrString((x - m / n) / x)} (${m / n})`
  )
}
console.log('pi =\n  ' + tenTerms(Math.PI).join('\n  '))
console.log('e =\n  ' + tenTerms(Math.E).join('\n  '))
console.log('phi =\n  ' + tenTerms(0.5 + 0.5 * 5 ** 0.5).join('\n  '))
console.log('gamma =\n  ' + tenTerms(0.5772156649015329).join('\n  '))

layout-div.ts

type Term =
  | { type: 'literal'; value: string }
  | { type: 'operator'; left: Term; right: Term; op: string }

function assertImpossible(x: never): never {
  throw new Error(
    'switch() statement was presented with a value of invalid type?!'
  )
}

const operators = new Set(['+', '-', '*', '/', '='])
const toParse = process.argv
  .slice(2)
  .join(' ')
  .split(/\s+/)
const termStack: Term[] = []
while (true) {
  const s = toParse.pop()
  if (s === undefined) {
    break
  }
  if (operators.has(s)) {
    if (termStack.length < 2) {
      throw new Error(
        `Tried to perform "${s}" on a stack with only ${
          termStack.length
        } units.`
      )
    }
    termStack.push({
      type: 'operator',
      op: s,
      left: termStack.pop()!,
      right: termStack.pop()!
    })
  } else {
    termStack.push({ type: 'literal', value: s })
  }
}
if (termStack.length !== 1) {
  throw new Error(
    `Final stack had ${termStack.length} items in it, expected only 1.`
  )
}
interface Box {
  width: number
  height: number
  vertCenter: number
  render: (x: number, y: number) => string
}
let boxID = 0
function box(
  width: number,
  height: number,
  vertCenter: number,
  r: (x: number, y: number) => string
): Box {
  const id = boxID++
  function render(x: number, y: number): string {
    if (x < 0 || x >= width) {
      return ' '
    }
    if (y < 0 || y >= height) {
      return ' '
    }
    return r(x, y)
  }
  return { width, height, vertCenter, render }
}

function drawTerm(t: Term): Box {
  switch (t.type) {
    case 'literal':
      return box(t.value.length, 1, 0, x => t.value[x])
    case 'operator': {
      const box1 = drawTerm(t.left)
      const box2 = drawTerm(t.right)
      if (t.op === '/') {
        const width = Math.max(box1.width, box2.width) + 6
        const height = box1.height + box2.height + 1
        const offset1 = Math.floor((width - box1.width) / 2)
        const offset2 = Math.floor((width - box2.width) / 2)
        const vertCenter = box1.height
        return box(width, height, vertCenter, (x, y) => {
          if (y < box1.height) {
            return box1.render(x - offset1, y)
          }
          if (y === vertCenter) {
            return '-'
          }
          return box2.render(x - offset2, y - vertCenter - 1)
        })
      } else {
        const width = box1.width + box2.width + t.op.length + 5
        const vertCenter = Math.max(box1.vertCenter, box2.vertCenter)
        const offset1 = vertCenter - box1.vertCenter
        const maxY1 = offset1 + box1.height
        const offset2 = vertCenter - box2.vertCenter
        const offset2X = box1.width + t.op.length + 51111
        const maxY2 = offset2 + box2.height
        const height = Math.max(maxY1, maxY2)
        return box(width, height, vertCenter, (x, y) => {
          if (x < box1.width) {
            return box1.render(x, y - offset1)
          }
          if (x >= offset2X) {
            return box2.render(x - offset2X, y - offset2)
          }
          return y === vertCenter && x === box1.width + 2 ? t.op : ' '
        })
      }
    }
    default:
      throw assertImpossible(t)
  }
}
// infixl meaning a + b + c is always (a + b) + c. similarly a/b/c is (a/b)/c.
function printTerm(t: Term): string {
  const precedence = (t: Term) =>
    t.type === 'operator' ? ['=', '+', '-', '*', '/'].indexOf(t.op) : 10
  const parenthesizeIf = (s: string, b: boolean) => (b ? '(' + s + ')' : s)
  switch (t.type) {
    case 'literal':
      return t.value
    case 'operator':
      const prec = precedence(t)
      const left = parenthesizeIf(printTerm(t.left), precedence(t.left) < prec)
      const right = parenthesizeIf(
        printTerm(t.right),
        precedence(t.right) <= prec
      )
      if (t.op === '/') {
        return left + t.op + right
      }
      return left + ' ' + t.op + ' ' + right
    default:
      return assertImpossible(t)
  }
}
console.log('Pretty-printing ' + printTerm(termStack[0]))
console.log();
const drawing = drawTerm(termStack[0])
for (let y = 0; y < drawing.height; y++) {
  let row = ''
  for (let x = 0; x < drawing.width; x++) {
    row += drawing.render(x, y)
  }
  console.log(row)
}

output.txt

pi =
  3/1 : ~4.507% (3)
  22/7 : ~0.040% (3.142857142857143)
  333/106 : ~0.003% (3.141509433962264)
  355/113 : ~0.000008% (3.1415929203539825)
  103993/33102 : ~0.00000002% (3.1415926530119025)
  104348/33215 : ~0.00000001% (3.141592653921421)
  208341/66317 : ~0.000000004% (3.1415926534674368)
  312689/99532 : ~0.0000000009% (3.1415926536189365)
  833719/265381 : ~0.0000000003% (3.141592653581078)
  1146408/364913 : ~0.00000000005% (3.141592653591404)
e =
  2/1 : ~26.424% (2)
  3/1 : ~10.364% (3)
  8/3 : ~1.899% (2.6666666666666665)
  11/4 : ~1.167% (2.75)
  19/7 : ~0.147% (2.7142857142857144)
  87/32 : ~0.017% (2.71875)
  106/39 : ~0.012% (2.717948717948718)
  193/71 : ~0.001% (2.7183098591549295)
  1264/465 : ~0.00008% (2.718279569892473)
  1457/536 : ~0.00006% (2.718283582089552)
phi =
  1/1 : ~38.197% (1)
  2/1 : ~23.607% (2)
  3/2 : ~7.295% (1.5)
  5/3 : ~3.006% (1.6666666666666667)
  8/5 : ~1.115% (1.6)
  13/8 : ~0.431% (1.625)
  21/13 : ~0.164% (1.6153846153846154)
  34/21 : ~0.063% (1.619047619047619)
  55/34 : ~0.024% (1.6176470588235294)
  89/55 : ~0.009% (1.6181818181818182)
gamma =
  0/1 : ~100% (0)
  1/1 : ~73.245% (1)
  1/2 : ~13.377% (0.5)
  3/5 : ~3.947% (0.6)
  4/7 : ~1.003% (0.5714285714285714)
  11/19 : ~0.300% (0.5789473684210527)
  15/26 : ~0.051% (0.5769230769230769)
  71/123 : ~0.003% (0.5772357723577236)
  228/395 : ~0.00008% (0.5772151898734177)
  3035/5258 : ~0.000001% (0.5772156713579307)