Skip to content

Instantly share code, notes, and snippets.

@etale

etale/alg.rb

Created November 24, 2009 09:28
Show Gist options
  • Save etale/241759 to your computer and use it in GitHub Desktop.
Save etale/241759 to your computer and use it in GitHub Desktop.
Download ZIP
module Algebraic
Raw
alg.rb
$radix ||= 10 # 2 <= $radix <= 36
$adic ||= 10**4 # 2 <= $adic
$delim ||= ',' # or '|', etc.
$type ||= :arch # or :adic
$prec ||= 44 # default precision
$term_order ||= :gr # or :lex
def scale; $adic ** $prec; end # default scale
def prec_2; (scale - 1).size 2; end # default precision in bits
def prec_calc; ($prec / $adic.log).ceil; end
module Algebraic
include Comparable
def algebraic?; true; end
def numeric?; false; end
def integer?; false; end
def rational?; false; end
def residual?; false; end
def adic?; false; end
def arch?; false; end
def polynomial?; false; end
Unity = Object.new
def Unity.=== a; 1.equal? a; end
@@pool = {}
def self.pool; @@pool; end
def identity; @@pool[self] ||= self; end
def finalize; self; end
def new *args, &block; self.class.new *args, &block; end
def get *args, &block; self.class.get *args, &block; end
def hash; self.class.hash; end
def eql? a; self.class.eql? a.class; end
def coerce a
_ = self
case a
when NilClass
_ = a
when Fraction
_ = Fraction.get _, 1
when Quotient
_ = Quotient.get _, a.r
else
raise TypeError.new "#{self}.coerce #{a}"
end
[a, _]
end
def cast a; a, = coerce a; a; end
def <=> a
case a
when Algebraic
case
when (equal? a)
0
when numeric? && ! a.numeric?
-1
when ! numeric? && a.numeric?
1
else
a, _ = coerce a
_.compare a
end
else
_, a = a.coerce self
_ <=> a
end
end
def zero; 0; end
def unity; 1; end
def sgn; self <=> zero; end # root of unity
def abs; zero? ? self : self / sgn; end
def unit; zero? ? unity : sgn; end
def puni; zero? ? unity : unit / sgn; end # principal unit
def body; self / unit; end
def zero?; self == zero; end
def unity?; self == unity; end
def sgn?; self == sgn; end
def abs?; self == abs; end
def unit?; self == unit; end
def puni?; self == puni; end
def body?; self == body; end
def +@; self; end
def + a
case a
when Algebraic
a, _ = coerce a
(_.add a).finalize
else
_, a = a.coerce self
_ + a
end
end
def -@; self * -1; end
def - a; self + -a; end
def * a
case a
when Algebraic
a, _ = coerce a
(_.multiply a).finalize
else
_, a = a.coerce self
_ * a
end
end
def inverse
Fraction.get 1, self
end
def / a; self * a.inverse; end
def % a
self +
case a
when Numeric
Residual.get 0, a.body
when Symbol, Term, Poly
_ = Ideal.get a => 0
Quotient.get 0, _
when Ideal
Quotient.get 0, a
else
raise NotImplementedError
end
end
def ** a
case a
when Integer
return inverse ** -a if a < 0
_, r = self, unity
until a.zero?
r *= _ if a.odd?
_ *= _
a >>= 1
end
r
else
(log * a).exp
end
end
def binomial a
case a
when Integer
a.inject {|_, i| _ * (self - i) / (a - i) } || unity
else
raise TypeError.new "#{self}.binomial #{a}"
end
end
def divmod a = $adic
case a
when Algebraic
a, _ = coerce a
_._divmod a
else
_, a = a.coerce self
_.divmod a
end
end
def orduni a = $adic # necessary?
case a
when Algebraic
a, _ = coerce a
_._orduni a
else
_, a = a.coerce self
_.orduni a
end
end
def div a = $adic; q, r = divmod a; q; end
def mod a = $adic; q, r = divmod a; r; end
def ord a = $adic; e, u = orduni a; e; end
def uni a = $adic; e, u = orduni a; u; end
def gcd a
case a
when Algebraic
a, _ = coerce a
_._gcd a
else
_, a = a.coerce self
_.gcd a
end
end
def _gcd a
_ = self
until a.zero?
_, a = a, (_.mod a)
end
_
end
def lcm a
(self * a).div gcd a
end
def inv a
_, x, z = self, 1, 0
until a.zero?
q, r = _.divmod a
_, a, x, z = a, r, z, x - q * z
end
x
end
def euclid a
_, x, y, z, w = self, 1, 0, 0, 1
until a.zero?
q, r = _.divmod a
_, a, x, y, z, w = a, r, z, w, x - q * z, y - q * w
end
[_, x, y]
end
def to_adic prec = ($prec / (Real.get $adic<<32, 32).log).ceil, adic = $adic
u = 1 % adic ** prec
self * (Adic.get 0, u, prec, adic)
end
def to_real prec = ($prec / (Real.get 2<<32, 32).log).ceil
self + (Real.get 0, prec)
end
def to_arch prec = ($prec / (Real.get 2<<32, 32).log).ceil
ord = Real.get 0, prec
arg = 0 % 2 ** prec
self * (Arch.get ord, arg)
end
end
class Numeric
include Algebraic
def numeric?; true; end
def polynomial?; true; end
def -@; super; end
def <=> a; super; end
def abs; super; end
def ceil; super; end
def floor; super; end
def round; super; end
def truncate; super; end
def div *a; super; end
def divmod *a; super; end
def modulo a; super; end
def quo a; super; end
def reminder a; super; end
# def eql? a; super; end
def integer?; super; end
def nonzero?; super; end
def zero?; super; end
# def step a, b; super; end
def to_int; super; end
def method_missing method, *args, &block
_ = Poly.get 1 => self
_.send method, *args, &block
end
def coerce a
_ = self
case a
when Numeric
_, a = a.coerce _
when Symbol, Term
a = Poly.get a => 1
_ = Poly.get 1 => _
when Poly
_ = Poly.get 1 => _
else
return super
end
[a, _]
end
end
class Integer
include Enumerable
def ord *a; super; end # Integer#ord is already defined
def integer?; true; end
def rational?; true; end
def residual?; true; end
def adic?; true; end
def arch?; true; end
def r; self; end
def s; 1; end
def n; 0; end
def method_missing method, *args, &block
# if unity?
# Term.get
# else
# Poly.get 1 => self
# end.send method, *args, &block
(Rational.get self, 1).send method, *args, &block
end
def each
(- self).times {|i| yield self + i } if self < 0
times {|i| yield i }
end
def coerce a
_ = self
case a
when Integer
when Numeric
_, a = a.coerce _
else
return super unless unity?
case a
when Symbol
a = Term.get a => 1
_ = Term.get
when Term
_ = Term.get
when Poly
_ = Poly.get 1 => 1
else
return super
end
end
[a, _]
end
def inverse
Rational.get 1, self
end
def _divmod a
return [0, self] if a.zero?
__divmod a
end
def to_a a = $adic
_, arr = abs, []
until _.zero?
_, i = _.divmod a
arr << i
end
arr
end
def size a = $adic
(to_a a).size
end
def orduni a = $adic
return [:infty, 1] if zero?
e, _ = 0, self
while true
q, r = _.divmod a
break unless r.zero?
e, _ = e + 1, q
end
[e, _]
end
def factor a = 2
return if unit?
bound = 2 << ((size 2).div 2)
while a < bound
return a if (mod a).zero?
a += 1
end
self
end
def factorize
return self if zero?
f, _ = unit, body
p = 2
until _.unity?
p = _.factor p
e, _ = _.orduni p
f *= Term.get p => e
end
f
end
# 3203431780337.fermat? => false ???
def fermat? a = 2
((a % self) ** (self - 1)).unity?
end
def rho
t = 2 % self
h = t**2 + 1
while true
f = (t - h).rn
return f unless f == 1
t, h = t**2 + 1, (h**2 + 1)**2 + 1
end
end
def ordfloat
return [-:infty, 1.0] if zero?
e = size 2
_ = abs / (1 << e)
_ = ((_.to_a 64, 2, :arch).to_num 0, 2) * 2.0 ** -64
_ = - _ if self < 0
[e, _]
end
def ordreal prec = prec_2
return [-:infty, (Real.get 1 << prec, prec)] if zero?
e = size 2
[e, (Real.get self << prec - e, prec)]
end
def factorial
map {|i| i + 1 }.product
end
end
class Fixnum
alias __div div
alias __divmod divmod
alias __to_s to_s
def ** a; super; end
def % a; super; end
def / a; super; end
def [] a; super; end
def div *a; super; end
def divmod *a; super; end
def ord *a; super; end
def size *a; super; end
def to_s *a; super; end
end
class Bignum
alias __div div
alias __divmod divmod
alias __to_s to_s
def ** a; super; end
def % a; super; end
def / a; super; end
def [] a; super; end
def div *a; super; end
def divmod *a; super; end
def ord *a; super; end
def size *a; super; end
def to_s *a; super; end
end
# inductive system of Z[1/s]
# inductive limit is Q, i.e. rational numbers
# \frac{r}{s}
# r, s \in \Z
class Rational < Numeric
def rational?; true; end
def adic?; true; end
def arch?; true; end
def initialize r = 0, s = 1
@r, @s = r * s.sgn, s.abs
end
def finalize
return if _s.zero?
return _r if _s.unity?
(get _r, _s).identity
end
attr_reader :r, :s
def n; s.zero? ? 1 : 0; end
def floor; r.div s; end
def fract; r.mod s; end
def rs; r.gcd s; end
def _r; r.div rs; end
def _s; s.div rs; end
def to_rat; self; end
def to_res; Residual.get r, s; end
def hash; super ^ r.hash ^ s.hash; end
def eq? a; s.eql? a.s; end
def eql? a; super and r.eql? a.r and eq? a; end
def method_missing method, *args, &block
case $type
when :arch
to_real
when :adic
to_adic
else
raise ''
end.send method, *args, &block
end
def coerce a
_ = self
case a
when _.class
return a, _ if eq? a
_s = s.lcm a.s
_ = get r * (_s.div s), _s unless _s.eql? s
a = get a.r * (_s.div a.s), _s unless _s.eql? a.s
when Integer
a = get a * s, s
when Residual
_, a = a.coerce _
else
return super
end
[a, _]
end
def compare a; r <=> a.r; end
def zero; get 0, s; end
def unity; get s, s; end
def sgn; r.sgn * s.sgn; end
def abs; get r * sgn, s; end
def unit; zero? ? unity : (get rs * sgn, s); end
def body; _r * sgn; end
def add a; get r + a.r, s; end
def -@; get -r, s; end
def multiply a; get r * a.r, s * s; end
def inverse; get s, r; end
def _divmod a
q, _ = r.divmod a.r
[q, (get _, s)]
end
def _gcd a
_r = r.gcd a.r
get _r, s
end
def orduni a = $adic
return [:infty, unity] if zero?
e_r, _r = r.orduni a
e_s, _s = s.orduni a
[e_r - e_s, _r / _s]
end
def factor; r.factor or s.factor; end
def factorize; r.factorize / s.factorize; end
def to_a prec = $prec, a = $adic, type = $type
scale = a ** prec
case type
when :arch
(self * scale).floor
when :adic
pole_ord = - [0, (ord a)].min
pole_scale = a ** pole_ord
((self * pole_scale) % (scale * pole_scale)).r
else
raise ArgumentError.new "#{type} is unknown type. use :adic or :arch"
end.to_a a
end
def ordfloat
return 0.ordfloat if zero?
e_r, _r = r.ordfloat
e_s, _s = s.ordfloat
[e_r - e_s, _r / _s]
end
def ordreal prec = prec_2
return 0.ordreal prec if zero?
e_r, _r = r.ordreal prec
e_s, _s = s.ordreal prec
[e_r - e_s, _r / _s]
end
def teich
to_adic.sgn
end
end
# projective system of Z/n
# projective limit is \hat{Z}, i.e. finite adeles
# r \pmod n
# r, n \in \Z
class Residual < Numeric
def residual?; true; end
def adic?; true; end
def initialize r = 0, n = 0
@r, @n = r, n
end
def finalize
_n = n.abs
_r = r.mod _n
return _r if _n.zero?
return if _n.unity?
(get _r, _n).identity
end
attr_reader :r, :n
def s; 1; end
def rn; r.gcd n; end
def _r; r.div rn; end
def _n; n.div rn; end
def to_rat; Rational.get r, n; end
def to_res; self; end
def hash; super ^ r.hash ^ n.hash; end
def eq? a; n.eql? a.n; end
def eql? a; super and r.eql? a.r and eq? a; end
def method_missing method, *args, &block
to_adic.send method, *args, &block
end
def coerce a
_ = self
case a
when _.class
return a, _ if eq? a
_n = n.gcd a.n
_ = get r, _n unless _n.eql? n
a = get a.r, _n unless _n.eql? a.n
when Integer
a = get a, n
when Rational
a = (cast a.r) / (cast a.s)
return coerce a
else
return super
end
[a, _]
end
def compare a; r <=> a.r; end
def zero; get 0, n; end
def unity; get 1, n; end
def sgn; get _r, n; end
def abs; get rn, n; end
def add a; get r + a.r, n; end
def -@; get -r, n; end
def multiply a; get r * a.r, n; end
def inverse
get (_r.inv _n), _n
end
def _divmod a
q, _ = r.divmod a.r
[q, (get _, n)]
end
def _gcd a
_r = r.gcd a.r
get _r, n
end
def orduni a = $adic
e, _ = r.orduni a
[e, (get _, n)]
end
def factor
rn.factor
end
def factorize
_ = (Term.get _r => 1) * rn.factorize
get _, n.factorize
end
def generator
g = unity
f = (n - 1).factorize
g += 1 while f.any? {|l, | (g ** ((n - 1)/l)).unity? }
g
end
end
# projective system of Z/p^n[1/p]
# projective limit is Qp, i.e. p-adic numbers
# p^\infty 1_\Z/p^prec
# p^ord u ord \in \Z, u \in \Gm(\Z/p^prec)
class Adic < Numeric
def adic?; true; end
def initialize ord = :infty, u = 1 % scale, prec = $prec, adic = $adic
@ord, @u, @prec, @adic = ord, u, prec, adic
end
attr_reader :ord, :u, :prec, :adic
def finalize; identity; end
def r; self * s; end
def s
return adic ** (- ord) if ord < 0
1
end
def n; 0; end
def scale; adic ** prec; end
def to_rat
adic ** ord * u.r
end
def to_res
if ord < 0
nil
else
adic ** ord * u
end
end
def hash; super ^ ord.hash ^ u.hash; end
def eq? a; u.eq? a.u; end
def eql? a; super and ord.eql? a.ord and u.eql? a.u; end
def coerce a
_ = self
case a
when _.class
return a, _ if eq? a
a_u, _u = u.coerce a.u
return nil, nil if _u.n.unity?
_ = get ord, _u, a.prec, adic unless _u.eq? u
a = get a.ord, a_u, prec, adic unless _u.eq? a_u
when Integer, Rational
a_ord, a_u = a.orduni adic
a_u = u.cast a_u
a = get a_ord, a_u, prec, adic
when Residual
_ = zero? ? a.zero : u.r % a.n
# a, _u = u.coerce a
# return nil, nil if _u.n.unity?
# a_prec = a.n.ord adic
# a_ord, a_u = a.orduni adic
# a = get a_ord, a_u, a_prec, adic
# _ = get ord, _u, a_prec, adic if a_prec < prec
when Real
# through continued fraction according to $type
return coerce a.to_rat
when Arch
return nil, nil
else
return super
end
[a, _]
end
def compare a
[a.ord, u] <=> [ord, a.u]
end
def zero; get :infty, u.unity, prec, adic; end
def unity; get 0, u.unity, prec, adic; end
def unit; get 0, u, prec, adic; end
def sgn
return self if zero?
_ = unit
a = _.frob
_, a = a, a.frob until _.eql? a
_
end
def zero?; ord.equal? :infty; end
def add a
return a if zero?
return self if a.zero?
return zero if eql? -a
return a.add self unless a.ord >= ord
_ = u + a.u * adic ** (a.ord - ord)
get ord + (_.rn.ord adic), _.unit, prec, adic
end
def -@
return self if zero?
get ord, - u, prec, adic
end
def multiply a
return self if zero?
return a if a.zero?
get ord + a.ord, u * a.u, prec, adic
end
def inverse
return if zero?
get -ord, u.inverse, prec, adic
end
def _divmod a
return [zero, self] if ord < a.ord
[(get ord - a.ord, u / a.u, prec, adic), zero]
end
def _orduni a
raise
end
def factor
return if unit?
adic
end
def factorize
return self if zero?
unit * (Term.get adic => ord)
end
def frob; self ** adic; end
def exp; 1 + exp0; end
def log
raise RangeError.new 'log' if zero?
- (1 - puni).log0
end
def exp_1
raise RangeError.new 'exp' unless zero? or (adic - 1) * ord > 1
return unity if zero?
rep = u.r
l = adic ** ord
r = scale / l
_, t, k = 1, rep, 1
until t.zero?
# puts (t * l).decimally prec
_ += t * l
k += 1
e, _k = k.orduni adic
s = adic ** (ord - e)
r /= s
l *= s
t = (t % r * rep % r * (_k % r).inverse).r
end
get 0, _ % scale, prec, adic
end
def exp_0
raise RangeError.new 'exp' unless zero? or (adic - 1) * ord > 1
return unity if zero?
_, t, k = unity, self, 1
while t.ord < prec + ord
# puts (t.u.r * adic ** t.ord).decimally prec
_ += t
k += 1
t *= self / k
end
_
end
def exp0
raise RangeError.new 'exp' unless zero? or (adic - 1) * ord > 1
r = self
unless zero?
n = 2
t = r * r / n
while t.ord < prec + ord
# puts r
r += t
n += 1
t *= self / n
end
end
r
end
def log0
r = self
unless zero?
n = 2
t = r * r / n
while t.ord < prec + ord
r += t
t *= n
n += 1
t *= self / n
end
end
r
end
def factorial_series
r = zero
n = zero
t = unity
while t.ord < prec
r += t
n += 1
t *= n
end
r
end
end
def list_fs
[2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31].map {|p|
$adic = p
puts 0.to_adic.factorial_series
}
end
class Frml < Numeric
# def adic?; true; end
def initialize ord = :infty, u = [1 % $adic]
@ord, @u = ord, u
end
attr_reader :ord, :u
def prec; u.size; end
def char; u[0].n; end
def adic
_u = Array.new prec, u[0].zero
_u[0] = u[0].unity
get 1, _u
end
def finalize; identity; end
def r; self * s; end
def s
_ord = - [ord, 0].min
_u = Array.new prec, u[0].zero
_u[0] = u[0].unity
get _ord, _u
end
def n; char; end
def scale; adic ** prec; end
def to_rat
adic ** ord * u.r
end
def hash; super ^ ord.hash ^ u.hash; end
def eq? a; char.eql? a.char and prec.eql? a.prec; end
def eql? a; super and ord.eql? a.ord and u.eql? a.u; end
def coerce a
_ = self
case a
when _.class
return a, _ if eq? a
return nil, nil unless char.eql? a.char
_prec = [prec, a.prec].min
range = 0..._prec
if _prec < prec
_ = get ord, u[range]
else
a = get a.ord, a.u[range]
end
when Integer, Rational, Residual
a_u = unity.u.map {|c| c * a }
a = if a_u[0].zero?
zero
else
get 0, a_u
end
when Real
return nil, nil
when Arch
return nil, nil
else
return super
end
[a, _]
end
def compare a
[a.ord, u] <=> [ord, a.u]
end
def zero
_u = Array.new prec, u[0].zero
_u[0] = u[0].unity
get :infty, _u
end
def unity
_u = Array.new prec, u[0].zero
_u[0] = u[0].unity
get 0, _u
end
def unit
return unity if zero?
get 0, u
end
def sgn
return self if zero?
_u = Array.new prec, u[0].zero
_u[0] = u[0]
get 0, _u
end
def zero?; ord.equal? :infty; end
def add a
return a if zero?
return self if a.zero?
return zero if eql? -a
return a.add self unless a.ord >= ord
d = a.ord - ord
_u = u[0...d] + (d...prec).map {|i| u[i] + a.u[i - d] }
i = 0
i += 1 while _u[i].zero?
_ord = ord + i
_u = _u[i...prec] + (Array.new i, u[0].zero)
get _ord, _u
end
def -@
return self if zero?
_u = u.map {|c| -c }
get ord, _u
end
def multiply a
return self if zero?
return a if a.zero?
get ord + a.ord, prec.map {|i| (i + 1).map {|j| u[j] * a.u[i - j] }.sum }
end
def inverse
return if zero?
s = u[0].inverse
_u = unit.u.map {|c| c * s }
x = 1 - (get 0, _u)
prec.map {|i| x**i }.sum
end
def _divmod a
return [zero, self] if ord < a.ord
[(get ord - a.ord, u / a.u, prec, adic), zero]
end
def _orduni a
raise
end
def factor
return if unit?
adic
end
def factorize
return self if zero?
unit * (Term.get adic => ord)
end
end
class Float
def method_missing method, *args, &block
to_rat.send method, *args, &block
end
def rep
return [-:infty, 0] if zero? # ???
e, _ = 0, abs
e, _ = e + 1, _ / 2 until _ < 1.0
e, _ = e - 1, _ * 2 while _ < 0.5
_ = (_ * 2.0 ** 64).floor
_ = - _ if self < 0
[e, _]
end
def to_exact
return Rational.get if zero?
e, _ = rep
2 ** (e - 64) * _
end
def to_real prec = prec_2 # use rep
(Real.get 0, prec) + to_rat
end
def to_arch prec = prec_2
ord = (to_real prec).log
arg = 0 % 2 ** prec
arg += arg.n.div 2 if self < 0
Arch.get ord, arg
end
def to_rat
_, a = to_exact, []
exact = _
error_scale = 1 << 32
while true
k = _.floor
a << k
r = a.to_rat
break if (exact - r).abs * error_scale < 1
_ = 1 / (_ - k)
end
r
end
end
# projective system of (1/2^prec)\Z
# projective limit is Ga(R), i.e. real numbers
# \frac{m}{2^prec}
class Real < Numeric
def arch?; true; end
def initialize m = 0, prec = prec_2
@m, @prec = m, prec
end
def finalize; identity; end
attr_reader :m, :prec
def r; self; end
def s; 1; end
def n; 0; end
def epsilon; get 1, prec; end
def scale; 1 << prec; end
def floor; m >> prec; end
def ceil; -(-self).floor; end
def to_exact; Rational.get m, scale; end
def to_res; Residual.get m, scale; end
def hash; super ^ m.hash ^ prec.hash; end
def eq? a; prec.eql? a.prec; end
def eql? a; super and m.eql? a.m and eq? a; end
def coerce a
_ = self
case a
when _.class
return a, _ if eq? a
_prec = [prec, a.prec].min
if prec > a.prec
_m = m >> prec - a.prec
_ = get _m, _prec
else
a_m = a.m >> a.prec - prec
a = get a_m, _prec
end
when Integer
a = get a << prec, prec
when Rational
a = (a * scale).floor
a = get a, prec
when Residual, Adic
# through continued fraction according to $type
return nil, nil # case of 2?
when Arch
return (m / scale).coerce a
else
return super
end
[a, _]
end
def compare a
m <=> a.m
end
def zero; get 0, prec; end
def unity; get scale, prec; end
def zero?; m.zero?; end
def add a
get m + a.m, prec
end
def -@
get -m, prec
end
def multiply a
_m = m * a.m >> prec
get _m, prec
end
def inverse
return if zero?
_m = (1 << prec * 2).div m
get _m, prec
end
def orduni_old
return [-:infty, unity] if zero?
e = (m.size 2) - prec
[e, (get m >> e, prec)]
end
def orduni
return [-:infty, unity] if zero?
e = (m.size 2) - prec
[e, (get m, prec + e)]
end
def _log
return zero if unity?
raise RangeError.new 'log' if zero?
e, _ = abs.orduni
_.__log + e * (get 0, prec + (e.size 2)).log2
end
def log2
_, r, k = 0, scale, 1
until (t = r.div k * 3 ** k).zero?
# puts new t, prec
_ += t
k += 2
end
get _ * 2, prec
end
def __log # change series?
_m = scale - m
_, r, k = 0, _m, 1
until (t = r.div k).zero?
# puts (new t, prec).abs
_ += t
r = r * _m >> prec
k += 1
end
get -_, prec
end
def _atan
return zero if zero?
return tau / 4 - inverse.atan if abs > 1
_, r, k = 0, m, 1
_m = - m * m >> prec
until (t = r.div k).zero?
# puts (new t, prec).abs
_ += t
r = r * _m >> prec
k += 2
end
get _, prec
end
def _tau
(4 * (unity / 5)._atan - (unity / 239)._atan) * 8
end
def _exp
return unity if zero?
_, r, s, k = scale, m, 1, 1
until (t = r.div s).zero?
# puts new t, prec
_ += t
r = r * m >> prec
k += 1
s *= k
end
get _, prec
end
def _cos
return unity if zero?
_, r, s, k = 0, m, 1, 0
_m = m * m >> prec
until (t = r.div s).zero?
# puts (new t, prec).abs
_ += t
r = - r * _m >> prec
k += 2
s *= k * (k - 1)
end
get _, prec
end
def _sin
return zero if zero?
_, r, s, k = 0, m, 1, 0
_m = m * m >> prec
until (t = r.div s).zero?
# puts (new t, prec).abs
_ += t
r = - r * _m >> prec
k += 2
s *= k * (k + 1)
end
get _, prec
end
def _tan
_sin / _cos
end
def ordfloat
return [-:infty, 1.0] if zero?
e_r, _r = m .ordfloat
e_s, _s = scale.ordfloat
[e_r - e_s, _r / _s]
end
def to_f
return 0.0 if zero?
e, _ = ordfloat
2.0 ** e * _
end
def use_math method
(Math.send method, to_f).to_real prec
end
def refine
get m * 2**32, prec + 32
end
def exp; cast refine._exp; end
def log; cast refine._log; end
def atan; cast refine._atan; end
def cos; cast refine._cos; end
def sin; cast refine._sin; end
def tan; cast refine._tan; end
def tau; cast refine._tau; end
def to_rat
return - (- self).to_rat if self < 0
exact, a = to_exact, []
_ = exact
while true
k = _.floor
a << k
r = a.to_rat
break if (exact - r).abs < 1 / scale
_ = 1 / (_ - k)
end
r
end
end
# projective system of \upsilon^\Z \sigma^\Z
# projective limit is \Gm(\C) and zero, i.e. complex numbers
# (e^{1/2^prec})^m (e^{\tau/2^prec})^r
class Arch < Numeric
def arch?; true; end
def initialize ord = (Real.new 0, prec_2), arg = 0 % 2 ** prec_2
@ord, @arg = ord, arg
end
attr_reader :ord, :arg
def finalize; identity; end
def r; self; end
def s; 1; end
def n; 0; end
def prec; arg.n.ord 2; end
def i
new rarg.zero, arg.zero + (arg.zero.n.div 4)
end
def tau
(new rarg.tau.log, arg.zero) * i
end
def rarg
Real.new arg.r, prec
end
def hash; super ^ ord.hash ^ arg.hash; end
def eq? a; arg.eq? a.arg; end
def eql? a
super and ord.eql? a.ord and arg.eql? a.arg
end
def coerce a
_ = self
case a
when _.class
a_arg, _arg = arg.coerce a.arg
if zero? or a.zero?
_ord = a_ord = -:infty
_ord = rarg.cast ord unless zero?
a_ord = rarg.cast a.ord unless a.zero?
else
a_ord, _ord = ord.coerce a.ord
end
a = new a_ord, a_arg unless a_arg.eq? a.arg
_ = new _ord, _arg unless _arg.eq? arg
when Integer
a_arg = arg.zero
a_arg += arg.n.div 2 if a < 0
a_ord = if a.zero?
-:infty
else
(rarg.cast a.abs).log
end
a = new a_ord, a_arg
when Rational
a = (cast a.r) / (cast a.s)
when Residual, Adic
return nil, nil
when Real
s = a.sgn
a = new a.log, arg.zero
a *= s
return _.coerce a
else
return super
end
[a, _]
end
def compare a
[ord, arg] <=> [a.ord, a.arg]
end
def zero
return self if zero?
get -:infty, arg.zero
end
def unity
return self if unity?
new rarg.zero, arg.zero
end
def zero?; ord.eql? -:infty and arg.zero?; end
def unity?; ord.zero? and arg.zero?; end
def repsilon; new rarg.epsilon, arg.zero; end
def iepsilon; new rarg.zero, arg.zero + 1; end
def add a
return a if zero?
return self if a.zero?
return zero if eql? -a
self * (a / self).inc
end
def -@
return self if zero?
_arg = arg + (arg.n.div 2)
get ord, _arg
end
def multiply a
return self if zero?
return a if a.zero?
new ord + a.ord, arg + a.arg
end
def inverse
return if zero?
new -ord, -arg
end
def inc
_ = exp
(new _.ord + 1, _.arg).log
end
def exp
return unity if zero?
t = ord.tau
targ = t * rarg
eord = ord.exp
_ord = eord * targ.cos
_arg = (eord * targ.sin / t).to_res
new _ord, _arg
end
def log
raise if zero?
t = ord.tau
targ = t * rarg
_ord = (ord ** 2 + targ ** 2).log / 2
_arg = ord.zero? ? i.arg : ((targ / ord).atan / t).to_res
new _ord, _arg
end
def reim
return [rarg.zero, rarg.zero] if zero?
t = ord.tau
targ = t * rarg
eord = ord.exp
[eord * targ.cos, eord * targ.sin]
end
# test
def it
self - 1 + (- self).exp
end
end
class NilClass
include Algebraic
def method_missing method, *args, &block
(Residual.new 0, 1).send method, *args, &block
end
def coerce a
case a
when Algebraic
[self, self]
else
raise TypeError.new "nil.coerce #{a}"
end
end
def compare a; 0; end
def zero; self; end
def unity; self; end
def sgn; self; end
end
class Symbol
include Algebraic
def polynomial?; true; end
def method_missing method, *args, &block
_ = Term.new self => 1
_.send method, *args, &block
end
def coerce a
_ = self
case a
when _.class
when Unity
a = Term.new
_ = Term.new _ => 1
when Term
_ = Term.new _ => 1
when Poly
_ = Poly.new _ => 1
when Numeric
a = Poly.new 1 => a
_ = Poly.new _ => 1
else
return super
end
[a, _]
end
def compare a; to_s <=> a.to_s; end
def add a
_ = Poly.new self => 1
a = Poly.new a => 1
_.add a
end
def multiply a
_ = Term.new self => 1
a = Term.new a => 1
_.multiply a
end
def inverse
Term.new self => -1
end
end
# inductive system of multiplicative free abelian groups
# inductive limit is free abelian group of countable generators,
# i.e. terms \prod_v v^{e_v}
class Term < Hash
include Algebraic
def polynomial?; true; end
def initialize a = {}
a.each {|v, e| e.zero? or self[v] = e }
self.default = 0
end
def finalize
return 1 if unity?
return lv if size.unity? and le.unity? and Symbol === lv
identity
end
alias vars keys
def lv; vars.max; end; def le; self[lv]; end
def ev; vars.min; end; def ee; self[ev]; end
def flatten; map {|v, e| v ** e }.product; end
def expand; map {|v, e| v.flatten ** e.flatten}.product; end
def hash; super ^ map {|v, e| v.hash ^ e.hash }.xor; end
def eql? a
super and size.eql? a.size and all? {|v, e| e.eql? a[v] }
end
def method_missing method, *args, &block
_ = Poly.new self => 1
_.send method, *args, &block
end
def coerce a
_ = self
case a
when _.class
when Unity
a = new
when Symbol
a = new a => 1
when Numeric
a = Poly.new 1 => a
_ = Poly.new _ => 1
when Poly
_ = Poly.new _ => 1
else
return super
end
[a, _]
end
alias unity? empty?
def compare a; (self / a).weight; end
def add a
_ = Poly.new self => 1
a = Poly.new a => 1
_.add a
end
def multiply a; new gen(a) {|e, a_e| e + a_e }; end
def inverse; new gen {|e| - e }; end
def ** a; get gen {|e| e * a }; end
def _lcm a; get gen(a) {|e, a_e| [e, a_e].max }; end
def _gcd a; get gen(a) {|e, a_e| [e, a_e].min }; end
def effective?; all? {|v, e| e > 0 }; end
def deg; map {|v, e| e }.sum; end
def weight term_order = $term_order
return 0 if unity?
case term_order
when :gr
d = deg
return d <=> 0 unless d.zero?
0 <=> ee
when :lex
le <=> 0
else
raise ArgumentError.new "#{term_order} is unknown order. use :gr or :lex"
end
end
def dif
f = {}
each {|v, e|
f[('d'+v.to_s).intern] = der v
}
Poly.get f
end
def der a
if include? a
_, c = {}, 0
each {|v, e|
if v == a
_[v], c = e - 1, e
else
_[v] = e
end
}
t = get _
Poly.get t => c
else
0
end
end
end
# inductive system of free modules over ground algebra
# inductive limit is free algebra of locally finite type,
# i.e. polynomials \sum_t c_t t
class Poly < Hash
include Algebraic
def polynomial?; true; end
def initialize a = {}
_ = Hash.new 0
a.each {|t, c| _[t.finalize] += c }
a = _
_ = 0
a.each {|t, c| _ = c.cast _ }
a.each {|t, c| (c = _.cast c).zero? or self[t] = c }
self.default = _
end
def finalize
return zero if zero?
if size.unity?
return lc if lt.unity?
return lt if Unity === lc
end
identity
end
alias terms keys
def vars; terms.map {|t| t.vars }.flatten.uniq; end
def deg; terms.map {|t| t.deg }.max; end
def lt; terms.max; end; def lc; self[lt]; end; def lm; lc * lt; end
def et; terms.min; end; def ec; self[et]; end; def em; ec * et; end
def flatten; map {|t, c| c * t }.sum; end
def expand; map {|t, c| c.flatten * t.flatten }.sum; end
def hash
super ^ map {|t, c| t.hash ^ c.hash }.xor
end
def eq? a
zero.eql? a.zero
end
def eql? a
super and eq? a and size.eql? a.size and all? {|t, c| c.eql? a[t] }
end
def coerce a
_ = self
case a
when _.class
return a, _ if eq? a
_zero = zero.cast a.zero
_ = change_base _zero unless _zero.eql? zero
a = a.change_base _zero unless _zero.eql? a.zero
when Symbol, Term
a = new a => 1
when Numeric
a = new 1 => a
else
return super
end
[a, _]
end
def change_base a
r = {}
each {|t, c| (c = a.cast c).zero? or r[t] = c }
new r
end
alias zero? empty?
def zero; default; end
def unity; default.unity; end
def compare a
return 0 if eql? a
ts, a_ts = [self, a].map {|_| _.terms.sort.reverse }
_ = ts <=> a_ts
return _ unless _.zero?
ts.each {|t|
_ = self[t] <=> a[t]
return _ unless _.zero?
}
raise "fail: #{self}.compare #{a}"
end
def add a; new gen(a) {|c, a_c| c + a_c }; end
def -@; get gen {|c| - c }; end
def multiply a
_ = Hash.new 0
each {|t, c|
a.each {|a_t, a_c|
_[t * a_t] += c * a_c
}
}
new _
end
def divmod a
_, q, r = self, new, 0
until _.zero?
f, = a.find {|f, | (_.lt / f.lt).effective? }
if f
c = (_.lc / f.lc) * (_.lt / f.lt)
_ -= c * f
q += new f => c
else
g = _.lm
_ -= g
r += g
end
end
[q, r]
end
def divmod_old a
_, q, r = self, new, 0
until _.zero?
f, = a.find {|f, | (_.lt / f.lt).effective? }
if f
q += new f => (_.lc / f.lc) * (_.lt / f.lt)
else
r += _.lm
end
_ -= _.lm
end
[q, r]
end
def syz a
t = lt.lcm a.lt
f = (lc/a.lc) * (t/lt) * self - (a.lc/lc) * (t/a.lt) * a
(1/f.lc) * f
end
def dif
f = {}
vars.each {|v|
f[('d'+v.to_s).intern] = der v
}
get f
end
def der a
_ = 0
each {|t, c|
_ += c * (t.der a)
}
_
end
end
class Ideal < Hash
include Algebraic
def initialize a = {}
self.default = 0
a.each {|l, r|
f = l - r
self[f] = 0
self.default += f.zero
}
end
def finalize
return zero if zero?
if size.unity?
return lc if lt.unity?
return lt if Unity === lc
end
identity
end
alias polys keys
def vars; terms.map {|t| t.vars }.flatten.uniq; end
def deg; terms.map {|t| t.deg }.max; end
def lt; terms.max; end; def lc; self[lt]; end; def lm; lc * lt; end
def et; terms.min; end; def ec; self[et]; end; def em; ec * et; end
def flatten; map {|t, c| c * t }.sum; end
def expand; map {|t, c| c.flatten * t.flatten }.sum; end
def hash
super ^ map {|f, | f.hash }.xor
end
def eq? a
zero.eql? a.zero
end
def eql? a
super and eq? a and size.eql? a.size and all? {|f, | a.include? f }
end
def coerce a
_ = self
case a
when _.class
return a, _ if eq? a
_zero = zero.cast a.zero
_ = change_base _zero unless _zero.eql? zero
a = a.change_base _zero unless _zero.eql? a.zero
when Symbol, Term
a = new a => 0
when Numeric
a = new a => 0
else
return super
end
[a, _]
end
def change_base a
r = {}
each {|t, c| (c = a.cast c).zero? or r[t] = c }
new r
end
alias zero? empty?
def zero; default; end
def unity; default.unity; end
def compare a
return 0 if eql? a
ts, a_ts = [self, a].map {|_| _.terms.sort.reverse }
_ = ts <=> a_ts
return _ unless _.zero?
ts.each {|t|
_ = self[t] <=> a[t]
return _ unless _.zero?
}
raise "fail: #{self}.compare #{a}"
end
def add a; get gen(a) {|f, a_f| 0 }; end
def -@; get gen {|c| 0 }; end
def multiply a
_ = Hash.new 0
each {|f, |
a.each {|a_f, |
_[f * a_f] = 0
}
}
get _
end
end
class Fraction
include Algebraic
def initialize r = 0, s = 1
@r, @s = r, s
end
attr_reader :r, :s
def finalize
d = r.gcd s
_r, _s = (r.div d), (s.div d)
return nil if _s.zero?
return _r if _s.unity?
(new _r, _s).identity
end
def add a
new r * a.s + s * a.r, s * a.s
end
def multiply a
new r * a.r, s * a.s
end
def inverse
new s, r
end
end
class Quotient
include Algebraic
def initialize r = 0, n = {}
@r, @n = r, n
end
attr_reader :r, :n
def finalize
_n = n
_r = r.mod _n
return _r if _n.zero?
return nil if _n.unity?
(new _r, _n).identity
end
def add a
get r + a.r, n
end
def -@
get -r, n
end
def multiply a
get r * a.r, n
end
def inverse
get r.inv n, n
end
end
class Class
def get *args, &block
(new *args, &block).identity
end
end
class Array
def sum; reduce 0, :+; end
def xor; reduce 0, :^; end
def product; reduce 1, :*; end
def to_num prec = 0, a = $adic
(reverse.inject {|s, i| s * a + i } || 0) / a ** prec
end
def pair_to_s adic = $adic, radix = $radix, delim = $delim
delim = adic > radix ? delim : ''
map {|a|
a.map {|i|
i.__to_s radix
}.join delim
}.join '.'
end
def to_rat
reverse.inject {|_, a| 1 / _ + a } || 0
end
end
class Hash
def gen *args
r = {}
case args.size
when 0
each {|k, v| r[k] = yield v }
when 1
a, = *args
(keys | a.keys).each {|k|
r[k] = yield self[k], a[k]
}
else # can enhance spec
raise ArgumentError, "#{self}.gen #{args}"
end
r
end
end
class String
def to_num adic = $adic, type = $type, radix = $radix, delim = $delim
case
when self[0..0] == '-'
- (self[1..-1].to_num adic, type, radix, delim)
when (include? '/')
r, s = split '/'
r = r.to_num adic, type, radix, delim
s = s.to_num adic, type, radix, delim
r / s
when (include? '%')
r, n = split '%'
r = r.to_num adic, type, radix, delim
n = n.to_num adic, type, radix, delim
r % n
when (include? '...')
_, = split '...'
l, r = _.split '.'
prec = r.size + (type == :adic ? 1 : 0)
_ = _.to_num adic, type, radix, delim
case type
when :arch
prec_2 = (adic ** prec - 1).size 2
_.to_real prec_2
when :adic
_.to_adic prec
else
raise
end
else
delim = adic > radix ? delim : ''
l, r = split '.'
l, r = [l, r].map {|s|
s.nil? ? [] : (s.split delim).map {|r| r.to_i radix }
}
a = l + r
case type
when :arch
a.reverse.to_num r.size, adic
when :adic
a.to_num l.size - 1, adic
else
raise
end
end
end
def to_adic
raise NotImplementedError
end
def to_arch
raise NotImplementedError
end
end
module Algebraic; def inspect; to_s; end; end
class Integer
def to_s adic = $adic, type = $type, radix = $radix, delim = $delim
return '-' + ((- self).to_s adic, type, radix, delim) if self < 0
s = (to_a adic).map {|i| i.__to_s radix }
s << '0' if zero?
delim = adic > radix ? delim : ''
case type
when :arch
s.reverse * delim
when :adic
s[0] + '.' + s[1..-1] * delim
end
end
end
class Rational
def to_s; "#@r/#@s"; end
def decimally prec = $prec, adic = $adic, type = $type, radix = $radix, delim = $delim
s = to_a prec, adic, type
prec.map {|i| s[i] ||= 0 }
delim = adic > radix ? delim : ''
case type
when :arch
s[prec] ||= 0
(sgn < 0 ? '-' : '') +
([s[prec..-1].reverse, s[0..prec-1].reverse].pair_to_s adic, radix, delim)
when :adic
e = - [0, (ord adic)].min
[s[0..e], s[(e+1)..-1]].pair_to_s adic, radix, delim
else
raise ArgumentError.new "#{type} is unknown type. use :adic or :arch"
end
end
end
class Residual;
def to_s; "#@r%#@n"; end;
def decimally prec = $prec, adic = $adic, type = $type, radix = $radix, delim = $delim
s = r.to_a adic
prec.map {|i| s[i] ||= 0 }
delim = adic > radix ? delim : ''
case type
when :arch
s[prec] ||= 0
([s[0..prec-1], s[prec..-1]].pair_to_s adic, radix, delim).reverse
when :adic
e = - [0, (ord adic)].min
[s[0..e], s[(e+1)..-1]].pair_to_s adic, radix, delim
else
raise ArgumentError.new "#{type} is unknown type. use :adic or :arch"
end
end
end
class Adic
def to_s radix = $radix, delim = $delim
delim = adic > radix ? delim : ''
s = zero? ? 0 : u.r
"#{s.decimally prec, adic, :adic, radix, delim}..." +
if zero? or unit?
''
else
' * ' +
if ord > 0
'0.1'
else
"1#{delim}0"
end +
if ord.unit?
''
else
"^#{ord.abs}"
end
end
end
def decimally prec = $prec, adic = $adic, type = :adic, radix = $radix, delim = $delim
to_rat.decimally prec, adic, type, radix, delim
end
end
class Frml
def to_s
return '0.' + '0' * (prec - 1) if zero?
u[0].r.to_s + '.' + u[1..-1].map {|c| c.r.to_s}.join + '...' +
if zero? or unit?
''
else
" * " +
if ord > 0
'0.1'
else
'10'
end +
if ord.unit?
''
else
"^#{ord.abs}"
end
end
end
end
class Real
def to_s adic = $adic
prec_adic = (scale.size adic) - 1
"#{to_exact.decimally prec_adic, adic, :arch}..."
end
def decimally prec = $prec, adic = $adic, type = :arch, radix = $radix, delim = $delim
to_exact.decimally prec, adic, type, radix, delim
end
end
class Arch
def to_s adic = $adic
return '0' if zero?
prec_adic = (arg.n.size adic) - 1
sarg = "#{arg.to_rat.decimally prec_adic, adic, :arch}"
sord = "#{ord.to_exact.decimally prec_adic, adic, :arch}"
sord + sarg[1..-1] + '... X'
end
=begin
def to_s adic = $adic
return '0' if zero?
prec_adic = (arg.n.size adic) - 1
sarg = "#{arg.to_rat.decimally prec_adic, adic, :arch}..."
case
when ord.zero?
return '1' if unity?
"etau#{sarg}"
when arg.zero?
"e#{ord.to_s adic}"
else
"e(#{ord.to_s adic}+tau#{sarg})"
end
end
=end
end
#=begin
class Symbol
alias inspect to_s
end
class Term
def to_s
return "1:#{self.class}" if unity?
vars.sort.map {|v|
e = self[v]
case v
when Symbol, Numeric
"#{v}"
else
"(#{v})"
end +
(e.unity? ? '' : "^#{e}")
}.join(' ')
end
end
class Poly
def r; self; end
def n; 0; end
def s; 1; end
=begin
(3.14 * :s).to_s => 157s
314/100 = 157/50
=end
def to_s
r = {}
each {|t, c| r[t] = c.r }
s = (new r)._to_s
return "(#{s})%#{zero.n}" unless zero.n.zero?
return "(#{s})/#{zero.s}" unless zero.s.unity?
s
end
def _to_s
return "0:#{self.class}" if zero?
terms.sort.reverse.map {|t|
c = self[t]
if c == 1 and t != 1
if t == lt
''
else
' + '
end
elsif c < 0
c = -c
' - ' +
if c == 1 and t != 1
''
else
c.to_s
end
else
if t == lt
''
else
' + '
end +
case c
when Numeric
c.to_s
else
"(#{c})"
end
end +
if t == 1
''
else
case t
when Symbol, Term
t.to_s
else
"(#{t})"
end
end
}.join
end
end
class Ideal
def to_s
_ = polys.sort.map {|f| f.to_s }.join ', '
"(#{_})"
end
end
class Fraction
def to_s
"(#{r})/(#{s})"
end
end
class Quotient
def to_s
"#{r}%#{n}"
end
end
#=end
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment