ratfuncfield.myr

   1 use std
   2
   3 use t
   4
   5 use "traits"
   6 use "groebner"
   7 use "polyring"
   8 use "Q"
   9 use "Z"
  10
  11 /*
  12    Note that you probably want to use reduceratfunc liberally. As is,
  13    this won't even simplify "x^2 / (x - 1)" + "-1 / (x - 1)" for you.
  14    This is because our groebner implementation is so damn slow.
  15  */
  16
  17 pkg yakmo =
  18         type ratfunc = struct
  19                 num : polynomial#
  20                 den : polynomial#
  21                 reduced : bool /* For saving work */
  22         ;;
  23         impl disposable ratfunc#
  24         impl division_struct ratfunc#
  25         impl group_struct ratfunc#
  26         impl field_struct ratfunc#
  27         impl ring_struct ratfunc#
  28         impl std.equatable ratfunc#
  29         impl t.dupable ratfunc#
  30
  31         const ratfuncfrompoly : (p : polynomial# -> ratfunc#)
  32         const ratfuncfromQ : (Q : Q# -> ratfunc#)
  33         const ratfuncfromS : (n : byte[:], d : byte[:] -> std.result(ratfunc#, byte[:]))
  34         const reduceratfunc : (rf : ratfunc# -> void)
  35
  36         const polynomialfromratfunc : (rf : ratfunc# -> std.option(polynomial#))
  37 ;;
  38
  39 const __init__ = {
  40         var rf : ratfunc#
  41         rf = rf
  42         std.fmtinstall(std.typeof(rf), ratfuncfmt)
  43 }
  44
  45 const ratfuncfmt = {sb, ap, opts
  46         var rf : ratfunc# = std.vanext(ap)
  47         var one : Q# = rid()
  48         auto one
  49
  50         var den_is_one = rf.den.terms.len == 1 && rf.den.terms[0].vars.len == 0 && std.eq(rf.den.terms[0].coeff, one)
  51
  52         if rf.num.terms.len == 0
  53                 std.sbputs(sb, "0")
  54                 -> void
  55         elif rf.num.terms.len == 1 || den_is_one
  56                 std.sbfmt(sb, "{}", rf.num)
  57         else
  58                 std.sbfmt(sb, "({})", rf.num)
  59         ;;
  60
  61         if rf.den.terms.len == 1
  62                 if rf.den.terms[0].vars.len == 0
  63                         if std.eq(rf.den.terms[0].coeff, one)
  64                                 -> void
  65                         ;;
  66                         if std.bigeqi(rf.den.terms[0].coeff.q, 1)
  67                                 std.sbfmt(sb, " / {}", rf.den.terms[0].coeff.p)
  68                                 -> void
  69                         ;;
  70                 ;;
  71         ;;
  72         std.sbfmt(sb, " / ({})", rf.den)
  73 }
  74
  75 impl disposable ratfunc# =
  76         __dispose__ = {x : ratfunc#
  77                 __dispose__(x.num)
  78                 __dispose__(x.den)
  79                 std.bytefree((x : byte#), sizeof(ratfunc))
  80         }
  81 ;;
  82
  83 impl division_struct ratfunc# =
  84         div_maybe = {a, b
  85                 if b.num.terms.len == 0
  86                         -> `std.None
  87                 ;;
  88                 var fake_bi : ratfunc = [ .num = b.den, .den = b.num ]
  89
  90                 -> `std.Some rmul(a, &fake_bi)
  91         }
  92 ;;
  93
  94 impl group_struct ratfunc# =
  95         gid = {
  96                 var num : polynomial# = gid()
  97                 var den : polynomial# = rid()
  98                 var base : ratfunc = [ .num = num, .den = den, .reduced = true ]
  99                 -> std.mk(base)
 100         }
 101
 102         gadd_ip = {a, b
 103                 var summand : polynomial# = rmul(b.num, a.den)
 104                 auto summand
 105                 rmul_ip(a.den, b.den)
 106                 rmul_ip(a.num, b.den)
 107                 gadd_ip(a.num, summand)
 108                 a.reduced = false
 109         }
 110
 111         gadd = {a, b
 112                 var da = t.dup(a)
 113                 gadd_ip(da, b)
 114                 -> da
 115         }
 116
 117         gneg_ip = {a
 118                 gneg_ip(a.num)
 119         }
 120
 121         gneg = {a
 122                 var da = t.dup(a)
 123                 gneg_ip(da)
 124                 -> da
 125         }
 126
 127         eq_gid = {a
 128                 -> eq_gid(a.num)
 129         }
 130 ;;
 131
 132 impl field_struct ratfunc# =
 133         finv = {a
 134                 if a.num.terms.len == 0
 135                         -> `std.None
 136                 ;;
 137
 138                 var base = [ .num = t.dup(a.den), .den = t.dup(a.num), .reduced = a.reduced ]
 139                 reduce_numeric_fractions(&base)
 140                 -> `std.Some std.mk(base)
 141         }
 142 ;;
 143
 144 impl ring_struct ratfunc# =
 145         rid = {
 146                 var num : polynomial# = rid()
 147                 var den : polynomial# = rid()
 148                 var base : ratfunc = [ .num = num, .den = den ]
 149                 -> std.mk(base)
 150         }
 151
 152         rmul_ip = {a, b
 153                 rmul_ip(a.num, b.num)
 154                 rmul_ip(a.den, b.den)
 155                 a.reduced = false
 156         }
 157
 158         rmul = {a, b
 159                 var da = t.dup(a)
 160                 rmul_ip(da, b)
 161                 -> da
 162         }
 163 ;;
 164
 165 impl std.equatable ratfunc# =
 166         eq = {a, b
 167                 var lhs : polynomial# = rmul(a.num, b.den)
 168                 auto lhs
 169                 var rhs : polynomial# = rmul(b.num, a.den)
 170                 auto rhs
 171
 172                 -> std.eq(lhs, rhs)
 173         }
 174 ;;
 175
 176 impl t.dupable ratfunc# =
 177         dup = {a
 178                 var dnum = t.dup(a.num)
 179                 var dden = t.dup(a.den)
 180                 var base = [ .num = dnum, .den = dden, .reduced = a.reduced ]
 181                 -> std.mk(base)
 182         }
 183 ;;
 184
 185 const ratfuncfrompoly = {p : polynomial#
 186         var num = t.dup(p)
 187         var den = rid()
 188         var base = [ .num = num, .den = den, .reduced = true ]
 189         reduce_numeric_fractions(&base)
 190         reduce_variable_fractions(&base)
 191         -> std.mk(base)
 192 }
 193
 194 const ratfuncfromQ = {q : Q#
 195         var num : polynomial# = polynomialfromQ(q)
 196         var den : polynomial# = rid()
 197         var base : ratfunc = [ .num = num, .den = den ]
 198         -> std.mk(base)
 199 }
 200
 201 const ratfuncfromS = {n : byte[:], d : byte[:]
 202         match polynomialfromS(n)
 203         | `std.Err e:
 204                 auto (e : t.doomed_str)
 205                 -> `std.Err std.fmt("cannot parse numerator: {}", e)
 206         | `std.Ok num:
 207                 match polynomialfromS(d)
 208                 | `std.Err e:
 209                         auto (e : t.doomed_str)
 210                         -> `std.Err std.fmt("cannot parse denominator: {}", e)
 211                 | `std.Ok den :
 212                         if den.terms.len == 0
 213                                 -> `std.Err std.fmt("division by zero")
 214                         ;;
 215                         var rfbase = [ .num = num, .den = den, .reduced = false ]
 216                         reduceratfunc(&rfbase)
 217                         -> `std.Ok std.mk(rfbase)
 218                 ;;
 219         ;;
 220 }
 221
 222 const reduceratfunc = {x
 223         /*
 224            LCM takes a really fucking long time, so be a little clever
 225            about not doing work.
 226          */
 227         if x.reduced
 228                 -> void
 229         ;;
 230
 231         if x.num.terms.len == 0
 232                 __dispose__(x.den)
 233                 x.den = rid()
 234                 x.reduced = true
 235                 -> void
 236         ;;
 237         reduce_numeric_fractions(x)
 238         reduce_variable_fractions(x)
 239
 240         /*
 241            This is polynomial lcm, so it always returns something monic.
 242            This may screw up the coefficients again.
 243          */
 244         var d = lcm(x.num, x.den)
 245         var leading_coeff = x.num.terms[x.num.terms.len - 1].coeff
 246         for var j = 0; j < d.terms.len; ++j
 247                 rmul_ip(d.terms[j].coeff, leading_coeff)
 248         ;;
 249
 250         match div_maybe(d, x.num)
 251         | `std.None:
 252         | `std.Some new_den:
 253                 match div_maybe(d, x.den)
 254                 | `std.None:
 255                         __dispose__(new_den)
 256                 | `std.Some new_num:
 257                         __dispose__(x.den)
 258                         __dispose__(x.num)
 259                         x.num = new_num
 260                         x.den = new_den
 261                 ;;
 262         ;;
 263         __dispose__(d)
 264
 265         x.reduced = true
 266 }
 267
 268 const reduce_numeric_fractions = {x
 269         /*
 270            Collect (integer) denominators of numerator and denominator,
 271            multiply through, hope to cancel some things.
 272          */
 273         var one = std.mkbigint(1)
 274         for var j = 0; j < x.num.terms.len; ++j
 275                 if !std.bigeqi(x.num.terms[j].coeff.q, 1)
 276                         var p = std.bigdup(x.num.terms[j].coeff.q)
 277                         var factor : Q = [ .p = p, .q = one ]
 278                         multiply_by_oneQ(x, &factor)
 279                         std.bigfree(p)
 280                 ;;
 281         ;;
 282         for var j = 0; j < x.den.terms.len; ++j
 283                 if !std.bigeqi(x.den.terms[j].coeff.q, 1)
 284                         var p = std.bigdup(x.den.terms[j].coeff.q)
 285                         var factor : Q = [ .p = p, .q = one ]
 286                         multiply_by_oneQ(x, &factor)
 287                         std.bigfree(p)
 288                 ;;
 289         ;;
 290
 291         /* Maybe we've got (-ax + b) / (-cx + d) */
 292         if x.num.terms.len > 0 && x.num.terms[0].coeff.p.sign < 0
 293                         var mone = std.mkbigint(-1)
 294                         var factor : Q = [ .p = mone, .q = one ]
 295                         multiply_by_oneQ(x, &factor)
 296                         std.bigfree(mone)
 297         ;;
 298
 299         /* Maybe we've got (ax + b) / (-1) */
 300         if x.den.terms.len == 1 && x.den.terms[0].vars.len == 0
 301                         var factor : Q = [ .p = x.den.terms[0].coeff.q, .q = x.den.terms[0].coeff.p ]
 302                         multiply_by_oneQ(x, &factor)
 303         ;;
 304         std.bigfree(one)
 305 }
 306
 307 const reduce_variable_fractions = {x
 308         /* As above, but now with variables */
 309
 310         /* Bail on 0 */
 311         if x.num.terms.len < 1
 312                 -> void
 313         ;;
 314
 315         var again = true
 316         while again
 317                 again = false
 318                 for var j = 0; j < x.num.terms.len; ++j
 319                         for var k = 0; k < x.num.terms[j].vars.len; ++k
 320                                 if x.num.terms[j].vars[k].power.sign < 0
 321                                         again = true
 322                                         var p = t.dup(x.num.terms[j].vars[k].power)
 323                                         abs_ip(p)
 324                                         var v = std.sldup(x.num.terms[j].vars[k].variable)
 325                                         var vp = [ .variable = v, .power = p ]
 326                                         var mon = [ .coeff = rid(), .vars = std.sldup([ vp ][:]) ]
 327                                         multiply_by_onemono(x, &mon)
 328                                         __dispose__(mon)
 329                                 ;;
 330                         ;;
 331                 ;;
 332         ;;
 333
 334         /*
 335            Perhaps the opposite has happened: both num and den are
 336            multiples of x^2 y or something.
 337          */
 338         if x.num.terms[0].vars.len < 1
 339                 -> void
 340         ;;
 341
 342         var oneQ : Q# = rid()
 343         auto oneQ
 344         for var j = 0; j < x.num.terms[0].vars.len; ++j
 345                 var vp : var_power# = &x.num.terms[0].vars[j]
 346                 var min_power = t.dup(vp.power)
 347                 var this_var_divides_all = true
 348                 for var k = 1; k < x.num.terms.len; ++k
 349                         var cur_mon : monomial# = &x.num.terms[k]
 350                         var seen_this_var = false
 351                         for var l = 0; l < cur_mon.vars.len; ++l
 352                                 var cur_vp : var_power# = &cur_mon.vars[l]
 353                                 match t.cmp(vp.variable, cur_vp.variable)
 354                                 | `std.Before:
 355                                 | `std.After: break
 356                                 | `std.Equal:
 357                                         seen_this_var = true
 358                                         match t.cmp(min_power, cur_vp.power)
 359                                         | `std.After:
 360                                                 __dispose__(min_power)
 361                                                 min_power = t.dup(cur_vp.power)
 362                                         | _:
 363                                         ;;
 364                                 ;;
 365                         ;;
 366                         if !seen_this_var
 367                                 this_var_divides_all = false
 368                                 break
 369                         ;;
 370                 ;;
 371
 372                 if !this_var_divides_all
 373                         __dispose__(min_power)
 374                         continue
 375                 ;;
 376
 377                 for var k = 0; k < x.den.terms.len; ++k
 378                         var cur_mon : monomial# = &x.den.terms[k]
 379                         var seen_this_var = false
 380                         for var l = 0; l < cur_mon.vars.len; ++l
 381                                 var cur_vp : var_power# = &cur_mon.vars[l]
 382                                 match t.cmp(vp.variable, cur_vp.variable)
 383                                 | `std.Before:
 384                                 | `std.After: break
 385                                 | `std.Equal:
 386                                         seen_this_var = true
 387                                         match t.cmp(min_power, cur_vp.power)
 388                                         | `std.After:
 389                                                 __dispose__(min_power)
 390                                                 min_power = t.dup(cur_vp.power)
 391                                         | _:
 392                                         ;;
 393                                 ;;
 394                         ;;
 395                         if !seen_this_var
 396                                 this_var_divides_all = false
 397                                 break
 398                         ;;
 399                 ;;
 400
 401                 if !this_var_divides_all
 402                         __dispose__(min_power)
 403                         continue
 404                 ;;
 405
 406                 /* We've found that x^min_power is contained in everything */
 407                 gneg_ip(min_power)
 408                 var dup_name = std.sldup(vp.variable)
 409                 var killer_vp = [
 410                         .variable = dup_name,
 411                         .power = min_power,
 412                 ]
 413                 var killer_monomial = [
 414                         .coeff = oneQ,
 415                         .vars = [ killer_vp ][:],
 416                 ]
 417                 mul_poly_mon_ip(x.num, &killer_monomial)
 418                 mul_poly_mon_ip(x.den, &killer_monomial)
 419                 __dispose__(min_power)
 420                 std.slfree(dup_name)
 421                 j--
 422         ;;
 423 }
 424
 425 const multiply_by_oneQ = {rf : ratfunc#, factor : Q#
 426         for var j = 0; j < rf.num.terms.len; ++j
 427                 rmul_ip(rf.num.terms[j].coeff, factor)
 428         ;;
 429         for var j = 0; j < rf.den.terms.len; ++j
 430                 rmul_ip(rf.den.terms[j].coeff, factor)
 431         ;;
 432 }
 433
 434 const multiply_by_onemono = {rf : ratfunc#, factor : monomial#
 435         mul_poly_mon_ip(rf.num, factor)
 436         mul_poly_mon_ip(rf.den, factor)
 437         rf.reduced = false
 438 }
 439
 440 const polynomialfromratfunc = {rf
 441         if rf.den.terms.len > 1
 442                 /* 1 / (x + 1) is not a polynomial - it's not finite */
 443                 -> `std.None
 444         ;;
 445
 446         var m = &rf.den.terms[0]
 447
 448         match finv(m.coeff)
 449         | `std.None: -> `std.None
 450         | `std.Some yi:
 451                 var p = t.dup(rf.num)
 452                 for var j = 0; j < p.terms.len; ++j
 453                         rmul_ip(p.terms[j].coeff, yi)
 454
 455                         for var k = 0; k < m.vars.len; ++k
 456                                 std.slpush(&p.terms[j].vars, [
 457                                         .variable = std.sldup(m.vars[k].variable),
 458                                         .power = gneg(m.vars[k].power)
 459                                 ])
 460                         ;;
 461                         reducemonomial(&p.terms[j])
 462                 ;;
 463
 464                 reducepoly(p)
 465                 -> `std.Some p
 466         ;;
 467 }