regression/harmonic.sml

   1 fun die str = (
   2        print (str ^ "\n");
   3        raise Overflow
   4        )
   5
   6 local fun loop (big: IntInf.int, small: IntInf.int): IntInf.int =
   7              if small = 0
   8                 then big
   9                 else loop (small,
  10                            IntInf.rem (big, small))
  11 in fun gcd (x: IntInf.int, y: IntInf.int): IntInf.int =
  12           let val x = IntInf.abs x
  13               val y = IntInf.abs y
  14               val (x, y) = if x >= y
  15                               then (x, y)
  16                               else (y, x)
  17           in loop (x, y)
  18           end
  19 end
  20
  21 fun reduce (num: IntInf.int, den: IntInf.int) : IntInf.int * IntInf.int =
  22        let val g = gcd (num, den)
  23            val gs = if den >= 0
  24                        then g
  25                        else ~ g
  26        in if gs = 1
  27              then (num, den)
  28           else let val rnum = IntInf.quot (num, gs)
  29                    val badn = IntInf.rem (num, gs)
  30                    val rden = IntInf.quot (den, gs)
  31                    val badd = IntInf.rem (den, gs)
  32                in if badn <> 0
  33                   orelse num <> rnum * gs
  34                   orelse badd <> 0
  35                   orelse den <> rden * gs
  36                      then die ("Bad: num " ^ (IntInf.toString num)
  37                                ^ ", den " ^ (IntInf.toString den)
  38                                ^ ", gcds " ^ (IntInf.toString gs)
  39                                ^ ", rnum " ^ (IntInf.toString rnum)
  40                                ^ ", rden " ^ (IntInf.toString rden)
  41                                ^ ", badn " ^ (IntInf.toString badn)
  42                                ^ ", badd " ^ (IntInf.toString badd))
  43                   else ();
  44                        (rnum, rden)
  45                end
  46        end
  47
  48 fun addrecip (xxx: int, (num: IntInf.int, den: IntInf.int))
  49              : IntInf.int * IntInf.int =
  50        let val xxx = IntInf.fromInt xxx
  51            val xnum = xxx * num + den
  52            val xden = xxx * den
  53        in reduce (xnum, xden)
  54        end
  55
  56 fun printRat (num: IntInf.int, den: IntInf.int): unit =
  57         print (IntInf.toString num ^ "/" ^ IntInf.toString den ^ "\n")
  58
  59 fun spin (limit: int): IntInf.int * IntInf.int =
  60        let fun loop (n: int, res: IntInf.int * IntInf.int)
  61                     : IntInf.int * IntInf.int =
  62                   if n = limit
  63                      then res
  64                   else loop (n + 1,
  65                              addrecip (n, res))
  66        in if limit <= 0
  67              then die "Bad limit"
  68              else loop (1, (0, 1))
  69        end
  70
  71 val (n, d) = spin 3000
  72 val _ = printRat (n, d)
  73 val _ = printRat (reduce (n * d * n, d * d * n))
  74 val _ = printRat (reduce (n + 1, d + 1))
  75 val _ = printRat (reduce (n * (d + 1) + d * (n + 1), d * (d + 1)))