diff --git a/src/lisp/s-expr.pla b/src/lisp/s-expr.pla
index ef23131..9cf911b 100644
--- a/src/lisp/s-expr.pla
+++ b/src/lisp/s-expr.pla
@@ -1,5 +1,7 @@
 include "inc/cmdsys.plh"
 include "inc/int32.plh"
+include "inc/fpstr.plh"
+include "inc/fpu.plh"
 
 const TYPE_MASK  = $70
 const NIL        = $00
@@ -10,6 +12,7 @@ const SYM_TYPE   = $20
 const SYM_LEN    = $0F
 const NUM_TYPE   = $30
 const NUM_INT    = $31
+const NUM_FLOAT  = $32
 const MARK_BIT   = $80
 const MARK_MASK  = $7F
 
@@ -32,17 +35,24 @@ struc t_numint
   res[t_elem]
   word[2] intval
 end
+struc t_numfloat
+  res[t_elem]
+  res[t_fpureg] floatval
+end
 
 predef eval_expr(expr)
 
 var sym_quote, sym_lambda, sym_cond
-res[t_elem] pred_true  = 0, 0, BOOL_TRUE
-res[t_elem] pred_false = 0, 0, BOOL_FALSE
+res[t_elem]   pred_true  = 0, 0, BOOL_TRUE
+res[t_elem]   pred_false = 0, 0, BOOL_FALSE
+res[t_numint] nan = 0, 0, NUM_INT, 0, 0, 0, 128 // NaN
 
 var cons_list  = NULL
 var cons_free  = NULL
 var int_list   = NULL
 var int_free   = NULL
+var float_list = NULL
+var float_free = NULL
 var sym_list   = NULL
 var assoc_list = NULL // SYM->value association list
 
@@ -60,6 +70,7 @@ end
 def mark_elems#0
   mark_list(cons_list)
   mark_list(int_list)
+  mark_list(float_list)
 end
 
 def sweep_expr(expr)#0
@@ -114,8 +125,9 @@ def collect_list(listhead, freehead)#2
 end
 
 def collect_unused#0
-  cons_list, cons_free = collect_list(cons_list, cons_free)
-  int_list,  int_free  = collect_list(int_list,  int_free)
+  cons_list,  cons_free  = collect_list(cons_list,  cons_free)
+  int_list,   int_free   = collect_list(int_list,   int_free)
+  float_list, float_free = collect_list(float_list, float_free)
 end
 
 export def gc#0
@@ -177,6 +189,43 @@ export def new_int(intlo, inthi)#1
   return intptr
 end
 
+def match_float(extptr)
+  var floatptr
+  byte i
+
+  floatptr = float_list
+  while floatptr
+    for i = 0 to 4
+      if floatptr=>floatval[i] <> extptr=>[i]
+        break
+      fin
+    next
+    if i > 4
+      return floatptr
+    fin
+    floatptr = floatptr=>link
+  loop
+  return NULL
+end
+
+export def new_float(extptr)#1
+  var floatptr
+
+  floatptr = match_float(extptr)
+  if floatptr; return floatptr; fin
+  if float_free
+    floatptr   = float_free
+    float_free = float_free=>link
+  else
+    floatptr = heapalloc(t_numfloat)
+  fin
+  floatptr=>link = float_list
+  float_list     = floatptr
+  floatptr->type = NUM_FLOAT
+  memcpy(floatptr + floatval, extptr, t_fpureg)
+  return floatptr
+end
+
 def match_sym(symstr)
   var symptr
   byte len, typelen, i
@@ -293,6 +342,9 @@ def print_atom(atom)#0
           is NUM_INT
             puti32(atom + intval)
             break
+          is NUM_FLOAT
+            puts(ext2str(atom + floatval, @prstr, 6, 4, FPSTR_FIXED|FPSTR_STRIP|FPSTR_FLEX))
+            break
         wend
         break
       is SYM_TYPE
@@ -336,26 +388,51 @@ end
 // Parse textual representation of S-expression
 //
 
-def is_int(c); return c >= '0' and c <= '9'; end
+def is_num(c); return c >= '0' and c <= '9'; end
 
 def is_alphasym(c)
   return (c >= '*' and c <= 'z') and (c <> '.') and (c <> ',')
 end
 
-def parse_int(evalptr)#2 // return evalptr, intptr
-  var int[2]
+def parse_num(evalptr)#2 // return evalptr, intptr
+  var startptr
+  var int[2], ext[5]
   byte sign
 
-  zero32
   sign = FALSE
   if ^evalptr == '-'
     sign = TRUE
     evalptr++
   fin
+  startptr = evalptr
   while ^evalptr >= '0' and ^evalptr <= '9'
-    muli16(10); addi16(^evalptr - '0')
     evalptr++
   loop
+  if (evalptr - startptr > 10) or ^evalptr == '.' or toupper(^evalptr) == 'E'
+    if ^evalptr == '.'
+      evalptr++
+      while ^evalptr >= '0' and ^evalptr <= '9'
+        evalptr++
+      loop
+    fin
+    if toupper(^evalptr) == 'E'
+      ^evalptr = 'E'
+      evalptr++
+      if ^evalptr == '-' or ^evalptr == '+'; evalptr++; fin
+      while ^evalptr >= '0' and ^evalptr <= '9'
+        evalptr++
+      loop
+    fin
+    if sign; startptr--; fin
+    ^(startptr - 1) = evalptr - startptr
+    str2ext(startptr - 1, @ext)
+    return evalptr, new_float(@ext)
+  fin
+  zero32
+  while startptr <> evalptr
+    muli16(10); addi16(^startptr - '0')
+    startptr++
+  loop
   if sign; neg32; fin
   store32(@int)
   return evalptr, new_int(int[0], int[1])
@@ -435,8 +512,8 @@ export def parse_expr(evalptr, level, refill)#2 // return evalptr, exprptr
         elemptr = NULL
         break
       otherwise
-        if is_int(^evalptr) or (^evalptr == '-' and is_int(^(evalptr+1)))
-          evalptr, elemptr = parse_int(evalptr)
+        if is_num(^evalptr) or (^evalptr == '-' and is_num(^(evalptr+1)))
+          evalptr, elemptr = parse_num(evalptr)
         elsif is_alphasym(^evalptr)
           evalptr, elemptr = parse_sym(evalptr)
         else
@@ -677,101 +754,228 @@ def natv_setq(expr)
   return valptr
 end
 
-export def eval_num(expr)#2
+export def eval_num(expr)
   var result
 
   result = eval_expr(expr=>car)
-  if result->type == NUM_INT
-    return result=>intval[0], result=>intval[1]
+  if result and (result->type & TYPE_MASK == NUM_TYPE)
+    return result
   fin
   puts("Not an number\n")
-  return 0, 0
+  return @nan
+end
+
+def push_int32(intptr)#0
+  var[2] int
+  byte isneg
+
+  isneg  = FALSE
+  if intptr=>[1] < 0
+    load32(intptr)
+    isneg = TRUE
+    neg32
+    store32(@int)
+  else
+    int[0] = intptr=>[0]
+    int[1] = intptr=>[1]
+  fin
+  fpu:pushInt(@int[1])
+  fpu:scalebXInt(16)
+  fpu:pushInt(@int[0])
+  fpu:addXY()
+  if isneg
+    fpu:negX()
+  fin
+end
+
+def push_num(numptr)#0
+  var int
+
+  if numptr->type == NUM_FLOAT
+    fpu:pushExt(numptr + floatval)
+  elsif numptr->type == NUM_INT
+    push_int32(numptr + intval)
+  else
+    puts("Pushing non number!\n")
+    int = 0
+    fpu:pushInt(@int)
+  fin
 end
 
 def natv_add(expr)
-  var[2] sum, num
+  var num
+  var[2] intsum
+  var[5] extsum
 
-  sum[0] = 0
-  sum[1] = 0
-  while expr
-    num[0], num[1] = eval_num(expr)
-    load32(@num)
-    add32(@sum)
-    store32(@sum)
-    expr = expr=>cdr
-  loop
-  return new_int(sum[0], sum[1])
+  intsum[0] = 0
+  intsum[1] = 0
+  num       = eval_num(expr)
+  expr      = expr=>cdr
+  if num->type == NUM_INT
+    //
+    // Sum as integers unless a float is encountered
+    //
+    intsum[0] = num=>intval[0]
+    intsum[1] = num=>intval[1]
+    while expr
+      num  = eval_num(expr)
+      expr = expr=>cdr
+      if num->type == NUM_FLOAT
+        break
+      fin
+      load32(@intsum)
+      add32(num + intval)
+      store32(@intsum)
+    loop
+  fin
+  if num->type == NUM_FLOAT
+    //
+    // Sum as floating point numbers
+    //
+    push_int32(@intsum)
+    push_num(num)
+    fpu:addXY()
+    while expr
+      num = eval_num(expr)
+      push_num(num)
+      fpu:addXY()
+      expr = expr=>cdr
+    loop
+    fpu:pullExt(@extsum)
+    return new_float(@extsum)
+  fin
+  return new_int(intsum[0], intsum[1])
 end
 
 def natv_sub(expr)
-  var[2] dif, num
+  var num1, num2
+  var[2] dif
+  var[5] ext
 
-  num[0], num[1] = eval_num(expr)
-  dif[0], dif[1] = eval_num(expr=>cdr)
-  load32(@num)
-  sub32(@dif)
-  store32(@dif)
-  return new_int(dif[0], dif[1])
+  num1 = eval_num(expr)
+  num2 = eval_num(expr=>cdr)
+  if num1->type == NUM_INT and num2->type == NUM_INT
+    load32(num1 + intval)
+    sub32(num2 + intval)
+    store32(@dif)
+    return new_int(dif[0], dif[1])
+  fin
+  push_num(num1)
+  push_num(num2)
+  fpu:subXY()
+  fpu:pullExt(@ext)
+  return new_float(@ext)
 end
 
 def natv_mul(expr)
-  var[2] mul, num
+  var num1, num2
+  var[2] mul
+  var[5] ext
 
-  num[0], num[1] = eval_num(expr)
-  mul[0], mul[1] = eval_num(expr=>cdr)
-  load32(@num)
-  mul32(@mul)
-  store32(@mul)
-  return new_int(mul[0], mul[1])
+  num1 = eval_num(expr)
+  num2 = eval_num(expr=>cdr)
+  if num1->type == NUM_INT and num2->type == NUM_INT
+    load32(num1 + intval)
+    mul32(num2 + intval)
+    store32(@mul)
+    return new_int(mul[0], mul[1])
+  fin
+  push_num(num1)
+  push_num(num2)
+  fpu:mulXY()
+  fpu:pullExt(@ext)
+  return new_float(@ext)
 end
 
 def natv_div(expr)
-  var[2] num, div
+  var num1, num2
+  var[2] div
+  var[5] ext
 
-  num[0], num[1] = eval_num(expr)
-  div[0], div[1] = eval_num(expr=>cdr)
-  load32(@num)
-  div32(@div)
-  store32(@div)
-  return new_int(div[0], div[1])
+  num1 = eval_num(expr)
+  num2 = eval_num(expr=>cdr)
+  if num1->type == NUM_INT and num2->type == NUM_INT
+    load32(num1 + intval)
+    div32(num2 + intval)
+    store32(@div)
+    return new_int(div[0], div[1])
+  fin
+  push_num(num1)
+  push_num(num2)
+  fpu:divXY()
+  fpu:pullExt(@ext)
+  return new_float(@ext)
 end
 
 def natv_rem(expr)
-  var[2] num, div
+  var num1, num2
+  var[2] rem, div
+  var[5] ext
 
-  num[0], num[1] = eval_num(expr)
-  div[0], div[1] = eval_num(expr=>cdr)
-  load32(@num)
-  num[1], num[0] = div32(@div)
-  return new_int(num[0], num[1])
+  num1 = eval_num(expr)
+  num2 = eval_num(expr=>cdr)
+  if num1->type == NUM_INT and num2->type == NUM_INT
+    load32(num1 + intval)
+    rem[1], rem[0] = div32(num2 + intval)
+    return new_int(rem[0], rem[1])
+  fin
+  push_num(num1)
+  push_num(num2)
+  fpu:remXY()
+  fpu:pullExt(@ext)
+  return new_float(@ext)
 end
 
 def natv_neg(expr)
-  var num[2]
+  var num
+  var[2] neg
+  var[5] ext
 
-  num[0], num[1] = eval_num(expr)
-  load32(@num)
-  neg32
-  store32(@num)
-  return new_int(num[0], num[1])
+  num = eval_num(expr)
+  if num=>type == NUM_INT
+    load32(num + intval)
+    neg32
+    store32(@neg)
+    return new_int(neg[0], neg[1])
+  fin
+  push_num(num)
+  fpu:negX()
+  fpu:pullExt(@ext)
+  return new_float(@ext)
 end
 
 def natv_gt(expr)
-  var[2] num1, num2
+  var num1, num2
+  var[5] ext
 
-  num1[0], num1[1] = eval_num(expr)
-  num2[0], num2[1] = eval_num(expr=>cdr)
-  load32(@num1)
-  return bool_pred(isgt32(@num2))
+  num1 = eval_num(expr)
+  num2 = eval_num(expr=>cdr)
+  if num1->type == NUM_INT and num2->type == NUM_INT
+    load32(num1 + intval)
+    return bool_pred(isgt32(num2 + intval))
+  fin
+  push_num(num2)
+  push_num(num1)
+  fpu:subXY()
+  fpu:pullExt(@ext)
+  return bool_pred(ext[4] < 0)
 end
 
 def natv_lt(expr)
-  var[2] num1, NUM2
+  var num1, num2
+  var[5] ext
 
-  num1[0], num1[1] = eval_num(expr)
-  num2[0], num2[1] = eval_num(expr=>cdr)
-  load32(@num1)
-  return bool_pred(islt32(@num2))
+  num1 = eval_num(expr)
+  num2 = eval_num(expr=>cdr)
+  if num1->type == NUM_INT and num2->type == NUM_INT
+    load32(num1 + intval)
+    return bool_pred(islt32(num2 + intval))
+  fin
+  push_num(num1)
+  push_num(num2)
+  fpu:subXY()
+  fpu:pullExt(@ext)
+  return bool_pred(ext[4] < 0)
 end
 
 def natv_print(expr)
@@ -814,5 +1018,6 @@ new_sym("NEG")=>natv    = @natv_neg
 new_sym(">")=>natv      = @natv_gt
 new_sym("<")=>natv      = @natv_lt
 new_sym("PRINT")=>natv  = @natv_print
+fpu:reset()
 return modkeep | modinitkeep
 done