--* From bmt@spadserv.watson.ibm.com  Sat Mar 20 14:14:08 1993
--* Received: from spadserv.watson.ibm.com by radical.watson.ibm.com (AIX 3.2/UCB 5.64/900524)
--*           id AA26266; Sat, 20 Mar 1993 14:14:08 -0500
--* Received: by spadserv.watson.ibm.com (AIX 3.2/UCB 5.64/900524)
--*           id AA23925; Sat, 20 Mar 1993 14:04:51 -0500
--* Date: Sat, 20 Mar 1993 14:04:51 -0500
--* From: bmt@spadserv.watson.ibm.com
--* X-External-Networks: yes
--* Message-Id: <9303201904.AA23925@spadserv.watson.ibm.com>
--* To: axc-bug@radical.watson.ibm.com, sutor@watson.vnet.ibm.com
--* Subject: asharp -M no-warnings says: Bug: Bad case 14 (line 680 in .. scobind.c) [xxfloat.as][v27.3 (current)]

--@ Fixed  by: RSS Wed Oct 13 14:59:00 1993
--@ Tested by: <name of new or existing file in test directory>
--@ Summary:   <One line description of real problem and the fix>


I ==> Integer
PI ==> Integer
NonNegativeInteger ==> Integer

Float:Type ==
  BASE ==> 2
  BITS:Reference PI := ref 68
  LENGTH ==> INTEGER_-LENGTH$Lisp
  ISQRT ==> approxSqrt$IntegerRoots(I)
  Rep := Record(mantissa:I,exponent:I)
  StoredConstant ==> Record(precision:PI,value:%)
  UCA ==> Record(unit:%,coef:%,associate:%)
  inc ==> increasePrecision
  dec ==> decreasePrecision
  times:(%,%) -> %
  itimes:(I,%) -> %
  chop:(%,PI) -> %
  dvide:(%,%) -> %
  square:(%,I) -> %
  power:(%,I) -> %
  plus:(%,%) -> %
  sub:(%,%) -> %
  negate:% -> %
  ceillog10base2:PI -> PI
  floorln2:PI -> PI
  atanSeries:% -> %
  atanInverse:I -> %
  expInverse:I -> %
  expSeries:% -> %
  logSeries:% -> %
  sinSeries:% -> %
  cosSeries:% -> %
  piRamanujan:() -> %
  asin (x:%):% ==
    zero? x => 0
    negative? x => - asin (- x)
    one? x => pi() / 2
    1 < x => error "asin: argument > 1 in magnitude"
    increasePrecision 5
    r := atan (x / sqrt sub(1,times(x,x)))
    decreasePrecision 5
    normalize r
  acos (x:%):% ==
    zero? x => pi() / 2
    negative? x =>
      increasePrecision 3
      r := pi() - acos (- x)
      decreasePrecision 3
      normalize r
    one? x => 0
    1 < x => error "acos: argument > 1 in magnitude"
    increasePrecision 5
    r := atan (sqrt sub(1,times(x,x)) / x)
    decreasePrecision 5
    normalize r
  atan(x:%,y:%):% ==
    x = 0 =>
      0 < y => pi() / 2
      y < 0 => - pi() / 2
      0
    theta := atan abs (y / x)
    if x < 0 then theta := pi() - theta
    if y < 0 then theta := - theta
    theta
  atan (x:%):% ==
    zero? x => 0
    negative? x => - atan (- x)
    if 1 < x then
      increasePrecision 4
      r :=
        zero? fractionPart x and x < per [bits()::Integer,0] =>
          atanInverse wholePart x
        atan (1 / x)
      r := pi / 2 - r
      decreasePrecision 4
      return normalize r
    k := approxSqrt(bits()::Integer - 100)$IntegerRoots(Integer)::Integer quo 5
    k := max(0,2 + k + order x)
    increasePrecision (2 * k)
    for i in 1..k repeat (x := x / (1 + sqrt (1 + x * x)))
    t := atanSeries x
    decreasePrecision (2 * k)
    t := shift(t,k)
    normalize t
  atanSeries (x:%):% ==
    p := bits() + INTEGER_-LENGTH(bits())$Lisp + 2
    s:Integer := d:Integer := shift(1,p::Integer)
    y := times(x,x)
    t := m := - shift((rep y).mantissa,(rep y).exponent + p::Integer)
    for i in 3.. by 2 while t ^= 0 repeat
      s := s + t quo i::Integer
      t := m * t quo d
    x * per [s,- p::Integer]
  atanInverse (n:Integer):% ==
    n2 := - n * n
    e:Integer :=
      (bits() + INTEGER_-LENGTH(bits())$Lisp + INTEGER_-LENGTH(n)$Lisp + 1)::
        Integer
    s:Integer := shift(1,e) quo n
    t:Integer := s quo n2
    for k in 3.. by 2 while t ^= 0 repeat
      s := s + t quo k::Integer
      t := t quo n2
    normalize per [s,- e]
  sin (x:%):% ==
    s := sign x
    x := abs x
    p := bits()
    increasePrecision 4
    if per [6,0] < x then
      increasePrecision (p::Integer)
      x := 2 * pi * fractionPart (x / pi / 2)
      bits p
    if per [3,0] < x then
      increasePrecision (p::Integer)
      s := - s
      x := x - pi
      bits p
    if per [3,- 1] < x then
      increasePrecision (p::Integer)
      x := pi - x
      decreasePrecision (p::Integer)
    k := approxSqrt(bits()::Integer - 100)$IntegerRoots(Integer)::Integer quo 4
    k := max(0,2 + k + order x)
    if 0 < k then
      increasePrecision k
      x := x / (3 ** k::NonNegativeInteger)::Integer
    r := sinSeries x
    for i in 1..k repeat (r := itimes(3,r) - shift(r ** 3,2))
    bits p
    s * r
  sinSeries (x:%):% ==
    p := bits() + INTEGER_-LENGTH(bits())$Lisp + 2
    y := times(x,x)
    s:Integer := d:Integer := shift(1,p::Integer)
    m:Integer := - shift((rep y).mantissa,(rep y).exponent + p::Integer)
    t:Integer := m quo 6
    for i in 4.. by 2 while t ^= 0 repeat
      s := s + t
      t := m * t quo (i * (i + 1))::Integer
      t := t quo d
    x * per [s,- p::Integer]
  cos (x:%):% ==
    s:Integer := 1
    x := abs x
    p := bits()
    increasePrecision 4
    if per [6,0] < x then
      increasePrecision (p::Integer)
      x := 2 * pi * fractionPart (x / pi / 2)
      decreasePrecision (p::Integer)
    if per [3,0] < x then
      increasePrecision (p::Integer)
      s := - s
      x := x - pi
      decreasePrecision (p::Integer)
    if per [1,0] < x then
      increasePrecision (p::Integer)
      x := pi / 2 - x
      bits p
      x := normalize x
      return s * sin x
    k := approxSqrt(bits()::Integer - 100)$IntegerRoots(Integer)::Integer quo 3
    k := max(0,2 + k + order x)
    if 0 < k then (increasePrecision k; x := shift(x,- k))
    r := cosSeries x
    for i in 1..k repeat (r := shift(r * r,1) - 1)
    bits p
    s * r
  cosSeries (x:%):% ==
    p := bits() + INTEGER_-LENGTH(bits())$Lisp + 1
    y := times(x,x)
    s:Integer := d:Integer := shift(1,p::Integer)
    m:Integer := - shift((rep y).mantissa,(rep y).exponent + p::Integer)
    t:Integer := m quo 2
    for i in 3.. by 2 while t ^= 0 repeat
      s := s + t
      t := m * t quo (i * (i + 1))::Integer
      t := t quo d
    normalize per [s,- p::Integer]
  tan (x:%):% ==
    s := sign x
    x := abs x
    p := bits()
    increasePrecision 6
    if per [3,0] < x then
      increasePrecision (p::Integer)
      x := pi() * fractionPart (x / pi())
      decreasePrecision (p::Integer)
    if per [3,- 1] < x then
      increasePrecision (p::Integer)
      x := pi() - x
      s := - s
      decreasePrecision (p::Integer)
    if 1 < x then (c := cos x; t := sqrt (1 - c * c) / c)
    else (c := sin x; t := c / sqrt (1 - c * c))
    bits p
    s * t
  P:StoredConstant := [1,[1,2]]
  pi():% ==
    free P
    P precision < bits() => (P := [bits(),piRamanujan()]).value
    normalize P value
  piRamanujan():% ==
    n := bits() + INTEGER_-LENGTH(bits())$Lisp + 11
    t:Integer := shift(1,n::Integer) quo 882
    d:Integer := (4 * 882 ** 2)::Integer
    s:Integer := 0
    for i in 2.. by 2 for j in 1123.. by 21460 while t ^= 0 repeat
      s := s + j::Integer * t
      m := - (i::Integer - 1) * ((2 * i)::Integer - 1) * ((2 * i)::Integer - 3)
      t := m * t quo (d * (i ** 3)::Integer)
    1 / per [s,- n::Integer - 2]
  sinh (x:%):% ==
    zero? x => 0
    lost:Integer := max(- order x,0)
    bits()::Integer < 2 * lost => x
    increasePrecision (5 + lost)
    e := exp x
    s := (e - 1 / e) / 2
    decreasePrecision (5 + lost)
    normalize s
  cosh (x:%):% ==
    increasePrecision 5
    e := exp x
    c := (e + 1 / e) / 2
    decreasePrecision 5
    normalize c
  tanh (x:%):% ==
    zero? x => 0
    lost:Integer := max(- order x,0)
    bits()::Integer < 2 * lost => x
    increasePrecision (6 + lost)
    e := exp x
    e := e * e
    t := (e - 1) / (e + 1)
    decreasePrecision (6 + lost)
    normalize t
  asinh (x:%):% ==
    p := min(0,order x)
    if zero? x or 2 * p < - bits()::Integer then return x
    increasePrecision (5 - p)
    r := log (x + sqrt (1 + x * x))
    decreasePrecision (5 - p)
    normalize r
  acosh (x:%):% ==
    if x < 1 then error "invalid argument to acosh"
    increasePrecision 5
    r := log (x + sqrt sub(times(x,x),1))
    decreasePrecision 5
    normalize r
  atanh (x:%):% ==
    if 1 < x or x < - 1 then error "invalid argument to atanh"
    p := min(0,order x)
    if zero? x or 2 * p < - bits()::Integer then return x
    increasePrecision (5 - p)
    r := log ((x + 1) / (1 - x)) / 2
    decreasePrecision (5 - p)
    normalize r
  log (x:%):% ==
    negative? x => error "negative log"
    zero? x => error "log 0 generated"
    p := bits()
    increasePrecision 5
    if (n := order x) < 0 then n := n + 1
    l :=
      n = 0 => 0
      x := shift(x,- n)
      n * log2
    k := approxSqrt(p::Integer - 100)$IntegerRoots(Integer)::Integer quo 3
    if 1 < k then
      k := max(1,k + order (x - 1))
      increasePrecision k
      ek := expInverse ((2 ** k::NonNegativeInteger)::Integer)
      decreasePrecision ((p::NonNegativeInteger quo 2)::Integer)
      m := order square(x,k)
      increasePrecision ((p::NonNegativeInteger quo 2)::Integer)
      m := 6847196937 * m quo 9878417065
      x := x * ek ** (- m)
      l := l + per [m,- k]
    l := l + logSeries x
    bits p
    normalize l
  logSeries (x:%):% ==
    p := bits() + (g := INTEGER_-LENGTH(bits())$Lisp + 3)
    increasePrecision (g::Integer)
    y := (x - 1) / (x + 1)
    decreasePrecision (g::Integer)
    s:Integer := d:Integer := shift(1,p::Integer)
    z := times(y,y)
    t := m := shift((rep z).mantissa,(rep z).exponent + p::Integer)
    for i in 3.. by 2 while t ^= 0 repeat
      s := s + t quo i::Integer
      t := m * t quo d
    y * per [s,1 - p::Integer]
  L2:StoredConstant := [1,1]
  log2():% ==
    free L2
    n := bits()::NonNegativeInteger::NonNegativeInteger
    L2 precision::NonNegativeInteger < n =>
      n := n + INTEGER_-LENGTH(n)$Lisp + 3
      s:Integer := shift(1,(n + 1)::Integer) quo 3
      t:Integer := s quo 9
      for k in 3.. by 2 while t ^= 0 repeat
        s := s + t quo k::Integer
        t := t quo 9
      L2 := [bits(),per [s,- n::Integer]]
      normalize L2 value
    normalize L2 value
  L10:StoredConstant := [1,[1,1]]
  log10():% ==
    n := bits()::NonNegativeInteger::NonNegativeInteger
    free L10
    L10 precision::NonNegativeInteger < n =>
      n := n + INTEGER_-LENGTH(n)$Lisp + 5
      s:Integer := shift(1,(n + 1)::Integer) quo 9
      t:Integer := s quo 81
      for k in 3.. by 2 while t ^= 0 repeat
        s := s + t quo k::Integer
        t := t quo 81
      increasePrecision 2
      L10 := [bits(),per [s,- n::Integer] + 3 * log2]
      decreasePrecision 2
      normalize L10 value
    normalize L10 value
  log2 (x:%):% ==
    increasePrecision 2
    r := log x / log2
    decreasePrecision 2
    normalize r
  log10 (x:%):% ==
    increasePrecision 2
    r := log x / log10
    decreasePrecision 2
    normalize r
  exp (x:%):% ==
    p := bits()
    increasePrecision 5
    e1:% := 1
    if not ((n := wholePart x) = 0) then
      increasePrecision INTEGER_-LENGTH(n)$Lisp
      e1 := exp1 ** n
      decreasePrecision INTEGER_-LENGTH(n)$Lisp
      x := fractionPart x
    if zero? x then (bits p; return normalize e1)
    k := approxSqrt(p::Integer - 100)$IntegerRoots(Integer)::Integer quo 3
    k := max(0,2 + k + order x)
    if 0 < k then (increasePrecision k; x := shift(x,- k))
    e := expSeries x
    if 0 < k then e := square(e,k)
    bits p
    e * e1
  expSeries (x:%):% ==
    p := bits() + INTEGER_-LENGTH(bits())$Lisp + 1
    s:Integer := d:Integer := shift(1,p::Integer)
    t:Integer :=
      n:Integer := shift((rep x).mantissa,(rep x).exponent + p::Integer)
    for i in 2.. while t ^= 0 repeat
      s := s + t
      t := n * t quo i::Integer
      t := t quo d
    normalize per [s,- p::Integer]
  expInverse (k:Integer):% ==
    p0:Integer := 2 * k + 1
    p1:Integer := 6 * k * p0 + 1
    q0:Integer := 2 * k - 1
    q1:Integer := 6 * k * q0 + 1
    for i in 10 * k.. by 4 * k while 2 * INTEGER_-LENGTH(p0)$Lisp < bits()
      repeat
        (p0,p1) := (p1,i * p1 + p0)
        (q0,q1) := (q1,i * q1 + q0)
    dvide(per [p1,0],per [q1,0])
  E:StoredConstant := [1,[1,1]]
  exp1():% ==
    free E
    if E precision < bits() then E := [bits(),expInverse 1]
    normalize E value
  sqrt (x:%):% ==
    negative? x => error "negative sqrt"
    m := (rep x).mantissa
    e := (rep x).exponent
    l := INTEGER_-LENGTH(m)$Lisp
    p := (2 * bits())::Integer - l + 2
    if odd? (e - l) then p := p - 1
    i := shift((rep x).mantissa,p)
    i := approxSqrt(i)$IntegerRoots(Integer)
    normalize per [i,(e - p) quo 2]
  bits():PositiveInteger == BITS()
  bits (n:PositiveInteger):PositiveInteger ==
    t := bits()
    BITS() := n
    t
  precision():PositiveInteger == bits()
  precision (n:PositiveInteger):PositiveInteger == bits n
  increasePrecision (n:Integer):PositiveInteger ==
    b := bits()
    bits ((b::Integer + n)::PositiveInteger)
    b
  decreasePrecision (n:Integer):PositiveInteger ==
    b := bits()
    bits ((b::Integer - n)::PositiveInteger)
    b
  ceillog10base2 (n:PositiveInteger):PositiveInteger ==
    ((13301 * n + 4003)::NonNegativeInteger quo 4004)::PositiveInteger
  digits():PositiveInteger ==
    max(1,4004 * (bits()::Integer - 1) quo 13301)::PositiveInteger
  digits (n:PositiveInteger):PositiveInteger ==
    t := digits()
    bits (1 + ceillog10base2 n)
    t
  order (a:%):Integer ==
    INTEGER_-LENGTH((rep(a)).mantissa)$Lisp + (rep a).exponent - 1
  relerror(a:%,b:%):Integer == order ((a - b) / b)
  0:% == per [0,0]
  1:% == per [1,0]
  base():PositiveInteger == 2
  mantissa (x:%):Integer == (rep x).mantissa
  exponent (x:%):Integer == (rep x).exponent
  one? (a:%):Boolean == a = 1
  zero? (a:%):Boolean == zero? (rep a).mantissa
  negative? (a:%):Boolean == negative? (rep a).mantissa
  positive? (a:%):Boolean == positive? (rep a).mantissa
  chop(x:%,p:PositiveInteger):% ==
    e:Integer := INTEGER_-LENGTH((rep(x)).mantissa)$Lisp - p::Integer
    if 0 < e then x := per [shift((rep x).mantissa,- e),(rep x).exponent + e]
    x
  float(m:Integer,e:Integer):% == normalize per [m,e]
  float(m:Integer,e:Integer,b:PositiveInteger):% ==
    m = 0 => 0
    increasePrecision 4
    r := m * per [b::Integer,0] ** e
    decreasePrecision 4
    normalize r
  normalize (x:%):% ==
    m := (rep x).mantissa
    m = 0 => 0
    e:Integer := INTEGER_-LENGTH(m)$Lisp - bits()::Integer
    if 0 < e then
      y := shift(m,1 - e)
      if odd? y
         then
           y := (0 < y => y + 1; y - 1) quo 2
           bits() < INTEGER_-LENGTH(y)$Lisp =>
             y := y quo 2
             e := e + 1
         else y := y quo 2
      x := per [y,(rep x).exponent + e]
    x
  shift(x:%,n:Integer):% == per [(rep x).mantissa,(rep x).exponent + n]
  (x:% = y:%):Boolean == order x = order y and sign x = sign y and zero? (x - y)
  (x:% < y:%):Boolean ==
    (rep y).mantissa = 0 => (rep x).mantissa < 0
    (rep x).mantissa = 0 => 0 < (rep y).mantissa
    negative? x and positive? y => true
    negative? y and positive? x => false
    order x < order y => positive? x
    order y < order x => negative? x
    negative? (x - y)
  abs (x:%):% ==
    negative? x => - x
    normalize x
  ceiling (x:%):% ==
    if negative? x then return - floor (- x)
    zero? fractionPart x => x
    truncate x + 1
  wholePart (x:%):Integer == shift((rep x).mantissa,(rep x).exponent)
  floor (x:%):% ==
    negative? x => - ceiling (- x)
    truncate x
  round (x:%):% ==
    half := per [sign x,- 1]
    truncate (x + half)
  sign (x:%):Integer ==
    (rep x).mantissa < 0 => - 1
    1
  truncate (x:%):% ==
    if not ((rep x).exponent < 0) then return x
    normalize per [shift((rep x).mantissa,(rep x).exponent),0]
  recip (x:%):Union(%,"failed") ==
    x = 0 => "failed"
    (1 / x)::Union(%,"failed")
  differentiate (x:%):% == 0
  (- x:%):% == normalize negate x
  negate (x:%):% == per [- (rep x).mantissa,(rep x).exponent]
  (x:% + y:%):% == normalize plus(x,y)
  (x:% - y:%):% == normalize plus(x,negate y)
  sub(x:%,y:%):% == plus(x,negate y)
  plus(x:%,y:%):% ==
    mx := (rep x).mantissa
    my := (rep y).mantissa
    mx = 0 => y
    my = 0 => x
    ex := (rep x).exponent
    ey := (rep y).exponent
    ex = ey => per [mx + my,ex]
    de := ex + INTEGER_-LENGTH(mx)$Lisp - ey - INTEGER_-LENGTH(my)$Lisp
    (bits() + 1)::Integer < de => x
    de < - (bits() + 1)::Integer => y
    if ex < ey then (mx,my,ex,ey) := (my,mx,ey,ex)
    mw := my + shift(mx,ex - ey)
    per [mw,ey]
  (x:% * y:%):% == normalize times(x,y)
  (x:Integer * y:%):% ==
    bits() < INTEGER_-LENGTH(x)$Lisp => normalize per [x,0] * y
    normalize per [x * (rep y).mantissa,(rep y).exponent]
  (x:% / y:%):% == normalize dvide(x,y)
  (x:% / y:Integer):% ==
    bits() < INTEGER_-LENGTH(y)$Lisp => x / normalize per [y,0]
    x / per [y,0]
  inv (x:%):% == 1 / x
  times(x:%,y:%):% ==
    per [(rep x).mantissa * (rep y).mantissa,(rep x).exponent + (rep y).exponent
      ]
  itimes(n:Integer,y:%):% == per [n * (rep y).mantissa,(rep y).exponent]
  dvide(x:%,y:%):% ==
    ew :=
      INTEGER_-LENGTH((rep(y)).mantissa)$Lisp -
        INTEGER_-LENGTH((rep(x)).mantissa)$Lisp + bits()::Integer + 1
    mw := shift((rep x).mantissa,ew) quo (rep y).mantissa
    ew := (rep x).exponent - (rep y).exponent - ew
    per [mw,ew]
  square(x:%,n:Integer):% ==
    ma := (rep x).mantissa
    ex := (rep x).exponent
    for k in 1..n repeat
      ma := ma * ma
      ex := ex + ex
      l:Integer := bits()::Integer::Integer - INTEGER_-LENGTH(ma)$Lisp
      ma := shift(ma,l)
      ex := ex - l
    per [ma,ex]
  power(x:%,n:Integer):% ==
    y:% := 1
    z:% := x
    repeat
      if odd? n then y := chop(times(y,z),bits())
      if (n := n quo 2) = 0 then return y
      z := chop(times(z,z),bits())
    never
  (x:% ** y:%):% ==
    x = 0 =>
      y = 0 => error "0**0 is undefined"
      y < 0 => error "division by 0"
      0 < y => 0
      never
    y = 0 => 1
    y = 1 => x
    x = 1 => 1
    p := abs order y + 5
    increasePrecision p
    r := exp (y * log x)
    decreasePrecision p
    normalize r
  (x:% ** r:Fraction Integer):% ==
    x = 0 =>
      r = 0 => error "0**0 is undefined"
      r < 0 => error "division by 0"
      0 < r => 0
      never
    r = 0 => 1
    r = 1 => x
    x = 1 => 1
    n := numer r
    d := denom r
    negative? x =>
      odd? d =>
        odd? n => return - (- x) ** r
        return (- x) ** r
      error "negative root"
    if d = 2 then
      increasePrecision INTEGER_-LENGTH(n)$Lisp
      y := sqrt x
      y := y ** n
      decreasePrecision INTEGER_-LENGTH(n)$Lisp
      return normalize y
    y := per [n,0] / per [d,0]
    x ** y
  (x:% ** n:Integer):% ==
    x = 0 =>
      n = 0 => error "0**0 is undefined"
      n < 0 => error "division by 0"
      0 < n => 0
      never
    n = 0 => 1
    n = 1 => x
    x = 1 => 1
    p := bits()
    bits (p + INTEGER_-LENGTH(n)$Lisp + 2)
    y := power(x,abs n)
    if n < 0 then y := dvide(1,y)
    bits p
    normalize y
  ceilLength10:I -> I
  chop10:(%,I) -> %
  convert10:(%,I) -> %
  floorLength10:I -> I
  length10:I -> I
  normalize10:(%,I) -> %
  quotient10:(%,%,I) -> %
  power10:(%,I,I) -> %
  times10:(%,%,I) -> %
  convert10(x:%,d:Integer):% ==
    m := (rep x).mantissa
    e := (rep x).exponent
    b := bits()
    (q,r) := divide(abs e,b::Integer)
    b := 2 ** b::NonNegativeInteger::NonNegativeInteger
    r := (2 ** r::NonNegativeInteger)::Integer
    h := power10(per [b::Integer,0],q,d + 5)
    h := chop10(per [r * (rep h).mantissa,(rep h).exponent],d + 5)
    e < 0 => h := quotient10(per [m,0],h,d)
    times10(per [m,0],h,d)
  ceilLength10 (n:Integer):Integer ==
    ((146 * INTEGER_-LENGTH(n)$Lisp)::NonNegativeInteger quo 485 + 1)::Integer
  floorLength10 (n:Integer):Integer ==
    ((643 * INTEGER_-LENGTH(n)$Lisp)::NonNegativeInteger quo 2136)::Integer
  length10 (n:Integer):Integer ==
    ln := INTEGER_-LENGTH(n := abs(n))$Lisp
    upper := (76573 * ln)::NonNegativeInteger quo 254370
    lower := 21306 * (ln - 1) quo 70777
    upper::Integer = lower => (upper + 1)::Integer
    n := n quo (10 ** lower::NonNegativeInteger)::Integer
    while n >= 10 repeat
      n := n quo 10
      lower := lower + 1
    lower + 1
  chop10(x:%,p:Integer):% ==
    e:Integer := floorLength10 (rep x).mantissa - p
    if 0 < e then x :=
      per [(rep x).mantissa quo (10 ** e::NonNegativeInteger)::Integer,
        (rep x).exponent + e]
    x
  normalize10(x:%,p:Integer):% ==
    ma := (rep x).mantissa
    ex := (rep x).exponent
    e:Integer := length10 ma - p
    if 0 < e then
      ma := ma quo (10 ** (e - 1)::NonNegativeInteger)::Integer
      ex := ex + e
      (ma,r) := divide(ma,10)
      4 < r =>
        ma := ma + 1
        ma = (10 ** p::NonNegativeInteger)::Integer =>
          ma := 1
          ex := ex + p
    per [ma,ex]
  times10(x:%,y:%,p:Integer):% == normalize10(times(x,y),p)
  quotient10(x:%,y:%,p:Integer):% ==
    ew := floorLength10 (rep y).mantissa - ceilLength10 (rep x).mantissa + p + 2
    if ew < 0 then ew := 0
    mw :=
      (rep x).mantissa * (10 ** ew::NonNegativeInteger)::Integer quo
        (rep y).mantissa
    ew := (rep x).exponent - (rep y).exponent - ew
    normalize10(per [mw,ew],p)
  power10(x:%,n:Integer,d:Integer):% ==
    x = 0 => 0
    n = 0 => 1
    n = 1 => x
    x = 1 => 1
    p:Integer := d + INTEGER_-LENGTH(n)$Lisp + 1
    e:Integer := n
    y:% := 1
    z:% := x
    repeat
      if odd? e then y := chop10(times(y,z),p)
      if (e := e quo 2) = 0 then return y
      z := chop10(times(z,z),p)
    never
  zero() ==> char "0"
  separator() ==> space()$Character
  SPACING:Reference N := ref 10
  OUTMODE:Reference S := ref "general"
  OUTPREC:Reference I := ref (- 1)
  fixed:% -> S
  floating:% -> S
  general:% -> S
  padFromLeft (s:String):String ==
    zero? SPACING() => s
    n:Integer := #s::Integer - 1
    t :=
      new((n + 1 + n quo SPACING()::Integer)::NonNegativeInteger,space()$
        Character)
    for i in 0..n for j in minIndex t.. repeat
      t j := s (i::Integer + minIndex s)
      (i + 1) rem SPACING() = 0 => j := j + 1
    t
  padFromRight (s:String):String ==
    SPACING() = 0 => s
    n:Integer := #s::Integer - 1
    t :=
      new((n + 1 + n quo SPACING()::Integer)::NonNegativeInteger,space()$
        Character)
    for i in n..0 by - 1 for j in maxIndex t.. by - 1 repeat
      t j := s (i + minIndex s)
      (n - i + 1) rem SPACING()::Integer = 0 => j := j - 1
    t
  fixed (f:%):String ==
    zero? f => "0.0"
    zero? exponent f => padFromRight concat(convert mantissa f @ String,".0")
    negative? f => concat("-",fixed abs f)
    d :=
      OUTPREC() = - 1 => digits::Integer::Integer
      OUTPREC()
    g := convert10(abs f,d)
    m := (rep g).mantissa
    e := (rep g).exponent
    if not (OUTPREC() = - 1) then
      l := length10 m
      OUTPREC() < - e and - e < 2 * digits::Integer::Integer =>
        g := normalize10(g,l + e + OUTPREC())
        m := (rep g).mantissa
        e := (rep g).exponent
    s := convert m @ String
    n := #s
    o := e + n::Integer
    p :=
      OUTPREC() = - 1 => n::Integer::Integer
      OUTPREC()
    t:String
    if e < 0
       then
         0 < o =>
           t := s (o + minIndex s..n::Integer + minIndex s - 1)
           s := s (minIndex s..o + minIndex s - 1)
         t := concat(new((- o)::NonNegativeInteger,char "0"),s)
         s := "0"
       else (s := concat(s,new(e::NonNegativeInteger,char "0")); t := "")
    n := #t
    if OUTPREC() = - 1 then (t := rightTrim(t,char "0"); t = "" => t := "0")
    else
      p < n::Integer => t := t (minIndex t..p + minIndex t - 1)
      t := concat(t,new((p - n::Integer)::NonNegativeInteger,char "0"))
    concat(padFromRight s,concat(".",padFromLeft t))
  floating (f:%):String ==
    zero? f => "0.0"
    negative? f => concat("-",floating abs f)
    t:String :=
      zero? SPACING() => "E"
      " E "
    zero? exponent f =>
      s := convert mantissa f @ String
      concat ["0.",padFromLeft s,t,convert (#s::Integer) @ String]
    d :=
      OUTPREC() = - 1 => digits::Integer::Integer
      OUTPREC()
    g := convert10(f,d)
    m := (rep g).mantissa
    e := (rep g).exponent
    s := convert m @ String
    n := #s
    o := e + n::Integer
    s := padFromLeft s
    concat ["0.",s,t,convert o @ String]
  general (f:%):String ==
    zero? f => "0.0"
    negative? f => concat("-",general abs f)
    d :=
      OUTPREC() = - 1 => digits::Integer::Integer
      OUTPREC()
    zero? exponent f =>
      d := d + 1
      s := convert mantissa f @ String
      d < (e := #s::Integer) =>
        t:String :=
          zero? SPACING() => "E"
          " E "
        concat ["0.",padFromLeft s,t,convert (e::Integer) @ String]
      padFromRight concat(s,".0")
    g := convert10(f,d)
    m := (rep g).mantissa
    e := (rep g).exponent
    s := convert m @ String
    n := #s
    o := n::Integer + e
    0 < o and max(n::Integer,d) >= o =>
      if n::Integer < o then s :=
        concat(s,new((o - n::Integer)::NonNegativeInteger,char "0"))
      t := rightTrim(s (o + minIndex s..n::Integer + minIndex s - 1),char "0")
      if t = "" then t := "0" else t := padFromLeft t
      s := padFromRight s (minIndex s..o + minIndex s - 1)
      concat(s,concat(".",t))
    0 >= o and o >= - 5 =>
      concat("0.",
        padFromLeft
          concat(new((- o)::NonNegativeInteger,char "0"),rightTrim(s,char "0")))
    t := padFromLeft rightTrim(s,char "0")
    s :=
      zero? SPACING() => "E"
      " E "
    concat ["0.",t,s,convert (e + n::Integer) @ String]
  outputSpacing (n:NonNegativeInteger):Void == SPACING() := n
  outputFixed():Void ==
    OUTMODE() := "fixed"
    OUTPREC() := - 1
  outputFixed (n:NonNegativeInteger):Void ==
    OUTMODE() := "fixed"
    OUTPREC() := n::Integer::Integer
  outputGeneral():Void ==
    OUTMODE() := "general"
    OUTPREC() := - 1
  outputGeneral (n:NonNegativeInteger):Void ==
    OUTMODE() := "general"
    OUTPREC() := n::Integer::Integer
  outputFloating():Void ==
    OUTMODE() := "floating"
    OUTPREC() := - 1
  outputFloating (n:NonNegativeInteger):Void ==
    OUTMODE() := "floating"
    OUTPREC() := n::Integer::Integer
  convert (f:%):String ==
    b:Integer :=
      OUTPREC() = - 1 and not (zero? f) =>
        bits (length abs mantissa f::PositiveInteger)::Integer
      0
    s :=
      OUTMODE() = "fixed" => fixed f
      OUTMODE() = "floating" => floating f
      OUTMODE() = "general" => general f
      empty()$String
    if 0 < b then bits (b::PositiveInteger)
    s = empty()$String => error "bad output mode"
    s
  coerce (f:%):OutputForm ==
    f < 0 => - coerce (- f) @ OutputForm
    message (convert f @ String)
  convert (f:%):InputForm ==
    convert [convert("float"::Symbol),convert(mantissa(f)),convert(exponent(f)),
      convert(base()::Integer)]$List(InputForm)
  convert (x:%):Float == x pretend Float
  convert (x:%):SmallFloat == makeSF((rep(x)).mantissa,(rep(x)).exponent)$Lisp
  coerce (x:%):SmallFloat == convert x @ SmallFloat
  convert (sf:SmallFloat):% == float(mantissa sf,exponent sf,base()$SmallFloat)
  retract (f:%):Fraction Integer ==
    rationalApproximation(f,(bits()::NonNegativeInteger - 1)::NonNegativeInteger
      ,2)
  retractIfCan (f:%):Union(Fraction Integer,"failed") ==
    rationalApproximation(f,(bits()::NonNegativeInteger - 1)::NonNegativeInteger
      ,2)::Union(Fraction Integer,"failed")
  retract (f:%):Integer ==
    f = (n := wholePart f)::% => n
    error "Not an integer"
  retractIfCan (f:%):Union(Integer,"failed") ==
    f = (n::Union(Integer,"failed") := wholePart f)::% => n
    "failed"
  rationalApproximation(f:%,d:NonNegativeInteger):Fraction Integer ==
    rationalApproximation(f,d,10)
  rationalApproximation(f:%,d:NonNegativeInteger,b:NonNegativeInteger):Fraction Integer ==
    t:Integer
    nu := (rep f).mantissa
    ex := (rep f).exponent
    if not (ex < 0) then return nu * (2 ** ex::NonNegativeInteger)::Integer / 1
    de := 2 ** (- ex)::NonNegativeInteger
    if b < 2 then error "base must be > 1"
    tol := b ** d
    s := nu
    t := de::Integer
    (p0,p1,q0,q1):Integer
    p0 := 0
    p1 := 1
    q0 := 1
    q1 := 0
    repeat
      (q,r) := divide(s,t)
      p2 := q * p1 + p0
      q2 := q * q1 + q0
      if r = 0 or tol::Integer * abs (nu * q2 - de::Integer * p2) < de::Integer
        * abs p2 then return p2 / q2
      (p0,p1) := (p1,p2)
      (q0,q1) := (q1,q2)
      (s,t) := (t,r)
    never