mathness.rb

require 'madness' # https://gist.github.com/pumbur/2df08af2ad846ace75fa
need 'prime', 'bigdecimal/math', 'matrix'

{:prm=>Prime,:mth=>Math,:vec=>Vector,:mat=>Matrix}.each{|n,k|  Kernel.class_alias(n,k) }


Num.act :sqrt, ->{ Math.sqrt(self) },
        :cbrt, ->{ Math.cbrt(self) },
        :sin,  ->{ Math.sin(self) },
        :cos,  ->{ Math.cos(self) },
        :log,  ->(q=Math::E){ Math.log(self,q) },
        :sqrt?,->{ (Math.sqrt(self)**2)==self },
        :cbrt?,->{ (Math.cbrt(self)**3)==self }
Int.act :fac,   ->{ Math.fac(self) },
        :frc,   ->*a{ Math.frc(self,*a) },
        :flip,  ->b=2,n=0{ to_s(b).rjust(n,'0').flip.to_i(b) }, # digits position flip (todo)
        :move,  ->a,b{ n=self.abs % (2**b); a = a%b; (n >> a) | ((n % (2**a)) << (b-a)) }, # bit rotation
        :prefix,->i=0,n=0{ return 0 if self<1; i+=1 while (self[i]==n); i }, # ctz (todo)
        :digsum,->i=2{ (i==2) ? digsum2 : digits(i).sum }, # digit sum (popcount/hamming-weight for binary)
        :digsum2, ->{
          x = self & 0xffffffff
          t=0x55555555; x = (x & t) + (x >> 1 & t)
          t=0x33333333; x = (x & t) + (x >> 2 & t)
          t=0x0f0f0f0f; x = (x & t) + (x >> 4 & t)
          t=0x00ff00ff; x = (x & t) + (x >> 8 & t)
          t=0x0000ffff; x = (x & t) + (x >>16 & t)
          t=self>>32;   x + (t==0 ? 0 : t.digsum2)
        }





Prime.aka :div => :prime_division
# euler's totient, number of divisors, sum of divisors, sum of the proper divisors
Prime.slf :phi, ->q{ q.pos? && div(q).fold(q){|o,(q,_)| o*q.dec/q } || 0 },
          :tau, ->q{ q.pos? && div(q).fold(1){|o,(_,q)| o*q.inc } || 0 },
          :sig, ->q{ q.pos? && div(q).fold(1){|o,(a,b)| o*((a**b.inc).dec/a.dec) } || 0 },
          :alq, ->q{ q.pos? && (sig(q)-q) || 0 }
# coprimes, divisors
Prime.slf :cop, ->n{ (n<2) ? [1] : div(n).fold((1...n).to_a){|o,(w,_)| o-(n/w).map{|q| q*w } } },
          :dvs, ->n{ (n<2) ? [n] : all(n).fold([1]){|o,q| o+o.map{|w| w*q } }.sort.uniq }
# ...
Prime.slf :sqr, ->n{ (q=dvs(n))[q.size/2...-1].map{|q| (q-(n/q))/2 } },
          :cot, ->n{ (n=cop(n))[0,n.size.inc/2]  },
          :one, ->n{ div(n).map(&:head) },
          :all, ->n{ div(n).flat_map{|a,b| [a]*b } }


# factorial  =->n{ n.prod }
#             ->n{ Prime.each(n).prod{|q| q**((n-n.digsum(q))/(q-1)) } }
#             ->n{ Prime.each(n).prod{|q| a=b=n/q; b+=(a/=q) while a>0; q**b } }
Mth.slf :fac, ->n{
  a,b,i,s=1,1,n.bit_length-1,->(q,w){ (d=w-q) < 2 ? q : (d==2) ? q*w : s[q,k=(q+d/2-1)|1] * s[k+2,w] }
  b *= ((m=n>>i)<3 ? 1 : (a *= s[(m>>1).inc|1, (m-1)|1])) while (i-=1)>=0; b << (n-n.digsum)
}
# inverse factorial (todo)
Mth.slf :frc, ->(n,m=1,i=1){ m*=(i+=1) while n>=m; i-1 }
# subfactorial (aka num of derangements/full-permutations)
#             ->n{ (n<1) ? 1 : Perm.prt(n).sum{|w| Perm.pcl(w,n) } }
#             ->n{ (n<1) ? 1 : ((f=n.fac)/BigMath.E(f.bit_length)).round }
Mth.slf :der, ->n{ n.inc.prod{|q| q+~0**q } }


# binom, catalan
Mth.slf :bin, ->n,r{ (n>=r) && (n.fac * (r.fac**-1) * ((n-r).fac**-1)).to_i || 0 },
        :cat, ->n{ ((2*n).fac * (n.fac**-2) * (n.inc**-1)).to_i }

# fibonacci numbers (lucas sequences; U=n,[0,1],P,Q; V=n,[2,P],P,Q)
Mth.slf :fib, ->n,s=[0,1],a=1,b=-1{ (n.rng? ? n.max : n).fold(s){|o,_| o << (o[-1]*a)-(o[-2]*b) }[n]  }
# rencontres numbers (per row, starting from triangular number)
Mth.slf :ren, ->n{ (n<2) ? [] : (2...n).fold([(n**2-n)/2]){|o,w| o<<((o[-1]*(n-w)*(d=der(w+1)))/(d-d.par)) }}
# 1st and 2nd stirling numbers (per row, starting from triangular number)
Mth.slf :st1, ->n{ n.fold([1]){|o,n| o.mapi{|w,i| (i.zero? ? 0 : o[i-1]) + (w*n) } << 1 }[1..-2].flip },
        :st2, ->n{ n.fold([1]){|o,n| o.mapi{|w,i| (i.zero? ? 0 : o[i-1]) + (w*i) } << 1 }[1..-2].flip }
Mth.slf :psc, ->n{ (1..n).fold([1]){|o,k| o << (o[-1]*(n-k+1)) / k; o } }, # pascal triangle rows (bin(n, 0..n))
        :bel, ->n{ st2(n).sum.inc }, # bell numbers
        :plg, ->n,i{ n * (i*i.inc)/2 + i.inc }, # polygonal numbers (n-2)
        :plt, ->n,i{ Math.bin(i+n+1,n+1) } # polytopic numbers (n-2)

Mth::PHI = 5.sqrt.inc / 2 # golden ratio (todo)



Perm = Class.new
Perm.slf :[], ->a,*q{ Perm.method(a).curry[*q] }

### permutation <-> index
# lexicographic order
Perm.slf :pid, ->(r,s=r.max){ r.flip.eachi.fold([0,1]){|(q,f),(u,j)| [q+f*r.rest(s-j).count{|w|w<u},f*(j+1)] }[0] },
         :prd, ->(n,s,o=s.to_a){ s.map{|i| o.delete_at((_,n=n.divmod((s-1-i).fac))[0]) }}
# flipped lexicographic order (ignoring fixed points - same id among all symmetric groups)
Perm.slf :fid, ->(r,s=r.max){ (r=r.map{|q| s-1-q }).eachi.fold([0,1]){|(q,f),(u,j)| [q+f*r[0,j].count{|w|w<u},f*(j+1)] }[0] },
         :frd, ->(n,s=n.frc,o=s.inc.to_a){ o.clone.map{|i| s-o.delete_at((_,n=n.divmod((s-i).fac))[0]) }.flip }

### partitions
# num of partitions
Perm.slf :prn, ->n,k{ (k>0)&&(n>0) ? prn(n-k,k) + prn(n-1,k-1) : (k==0)&&(n==0) ? 1 : 0 }
# all partitions of size n
Perm.slf :prt, ->(n,m=2){ (n<m) ? [[]] : (m..(n/2)).fmap{|q| prt(n-q,q).qap(:push,q) } << [n] }
# num of permutations with spectific partition
Perm.slf :pcl, ->(b,n=b.sum){ n.fac / (b+=[1]*(n-b.sum)).uniq.prod{|q| (w=b.count(q)).fac*(q**w) } }
# num of cyclic subgroups of S_n with specific partition (sum of every $b in $n == A051625)
Perm.slf :pcg, ->(b,n=b.sum){ pcl(b,n)/Prime.phi(b.fold(1,&:lcm)) }

# orbits (cycles ids)
Perm.slf :orb, ->(r,f=1,o=[],l=-1){ r.eachi{|w,i| (s=->(h,*q){ o[h]||=( (w==i) ? (f&&(l+=1)) : q.inc?(h) ? (l+=1) : s[r[h],h,*q] ) })[i] }; o }
# cycles
Perm.slf :cyc, ->(r){ ((o=orb(r,nil)).compact.max||-1).inc.map{|i| t=[l=o.index(i)]; t << l while !t.inc?(l=r[l]); t } }
# inverse, reduce
Perm.slf :inv, ->(q){ q.ids(*q.size.to_a) }, :red, ->(r){ r[0].num? ? r.sort.ids(*r) : (r.map(&red)-[[0]]) }


# num of cyclic subgroups of S_n with specific partition (sum of every $b in $n == A051625)
Perm.slf :pcg, ->(b,n=b.sum){  pcl(b,n)/Prime.phi(b.fold(1,&:lcm))  }
# # partition of perm
Perm.slf :prr, ->(r,t=1){ Perm.cyc(r).map(&:size).lap{|r| t ? r.rsort : r }}

# # generate all permutations with specific partition
Perm.slf :pgl, ->(b,n=b.sum,o=nil){
  return [n.to_a] if (b.empty? || (n<2))
  a,b,o = Perm.pgl(b.rest,n-b.head), n.to_a.perms(b.head), (o||Prime.cop(b.head))
  a.flat_map{|a| o.flat_map{|o| b.map{|b| [b+(n.to_a-b), b.move(o)+(n.to_a-b).vals(*a)].turn.fold([]){|o,(q,w)| o[q]=w; o } }}}.uniq.sort
}
Perm.slf :pgg, ->(b,n=b.sum,o=nil){
  return [n.to_a] if (b.empty? || (n<2))
  a,b,o = Perm.pgl(b.rest,n-b.head), n.to_a.perms(b.head), (o||Prime.cop(b.head))
  a.flat_map{|a| o.flat_map{|o| b.map{|b| Perm.fid([b+(n.to_a-b), b.move(o)+(n.to_a-b).vals(*a)].turn.fold([]){|o,(q,w)| o[q]=w; o }) }}}.uniq.sort
}

# transposition-permutations products paritions statistics
Perm.slf :stp, ->(q,n=q.sum+2,t=2,c=q.count_by){ [].tap{|o|
  # ±2 elements, transposition deletion/addition (5 <- 5:2 -> 5:2:2)
  o << [(q-[2]+[2]*c[2].dec).rsort, 2] if c[2]
  o << [[*q,2].rsort, 2] if q.sum<=(n-2)
  # ±1 elements, cycle decrease/increase (4:2 <- 5:2 -> 6:2)
  (c-[2]).each{|k,v| o << [(q-[k]+[k]*v.dec+[k.dec]).rsort, 1] }
  (q.sum<=n.dec) && c.each{|k,v| o << [(q-[k]+[k]*v.dec+[k.inc]).rsort, 1] }
  # ±0 elements, cycle split/merge (3:2:2 <- 5:2 -> 7)
  c.each{|k,v| (k>3)&&(k/2).dec.each{|i| o << [[*(q-[k]+[k]*v.dec),2+i,k-2-i].rsort, 0] } }
  q.combs(2).uniq.map{|w| o << [[*c.fmap{|k,v| [k]*(v-w.count(k)) }, w.sum].rsort, 0] }
}.map{|k,v| [k, t<1 ? v : Perm.send("stp#{v}",k,q)*(t>1 ? Math.bin(n-2,[k.sum,q.sum].max-2) : 1) ] } }
Perm.slf :stp2, ->a,b{ (a.sum>b.sum)&&(a,b=b,a); Perm.pcl(a,a.sum) },
         :stp1, ->a,b{ (a.sum>b.sum)&&(a,b=b,a); 2 * a.sub(b).sum.dec.fac * Perm.pcl(a.and(b),a.sum-1) },
         :stp0, ->a,b{ (a.len>b.len)&&(a,b=b,a); ((z=b.sub(a)).same? ? 1 : 2) * (z=z.sum-2).fac * Perm.pcl(a=a.and(b),a.sum+z) }






Mth.slf :parabola, ->d{ # y = a*(x**2) + b*x + c
  xs,ys,xy,x2,x3,x4,x2y = d.map{|y,x| [x,y,x*y,x**2,x**3,x**4,(x**2)*y] }.turn.map(&:sum)
  rx,ry = Mat[[x2,xs,d.size],[x3,x2,xs],[x4,x3,x2]], Mat[[ys,xy,x2y]].transpose
  ((rx.transpose * rx).inverse * (rx.transpose * ry)).to_a.flat
}
Math.slf :paraboloid, ->d{ # z = a*(x**2) + b*(x*y) + c*(y**2) + d*x + e*y + f
  dd = d.turn.map{|d| 4.fold([d]){|o,_| o << d.mapi{|q,i| q*o[-1][i] } } }.flip
  ww = %w'201 111 021 101 011 001  400 310 220 300 210 200  310 220 130 210 120 110  220 130 040 120 030 020  300 210 120 200 110 100  210 120 030 110 020 010  200 110 020 100 010'
  dd = ww.fold({nil=>d.size}){|o,q| o[q] ||= q.split('').mapi{|q,i| (q=q.to_i).zero? ? [1]*d.size : dd[i][q-1] }.turn.sum(&:prod); o }
  a,b = Mat[*(dd.vals(*ww[6..-1]) << dd[nil]).slices(6)].transpose, Mat[dd.vals(*ww[0,6])].transpose
  ((a*a.transpose).inverse*(a*b)).to_a.flat
}
Mth.slf :bezier_cubic, ->d{
  a = [[2,-6,6,-2,-6,18,-18,6,6,-18,18,-6,-2,6,-6,2],[0,1,2,3,1,2,3,4,2,3,4,5,3,4,5,6],[6,5,4,3,5,4,3,2,4,3,2,1,3,2,1,0]].turn
  b = [[2,-6,6,-2],[0,1,2,3],[3,2,1,0]].turn
  a = Mat[*a.map{|q| d.sum{|x,i| q[0]*i**q[1]*(i-1)**q[2] } }.slices(4)]
  b = Mat[*b.map{|q|-d.sum{|x,i| q[0]*i**q[1]*(i-1)**q[2]*x } }.slices(1)]
  (a.inv * b).column(0).to_a
}
Mth.slf :bezier_cubic_m, ->d{
  a = [[18,-18,-18,18],[2,3,3,4],[4,3,3,2]].turn
  b = [[-6,6],[1,2],[2,1]].turn
  a = Mat[*a.map{|q| d.sum{|x,i| q[0]*i**q[1]*(i-1)**q[2] } }.slices(2)]
  b = Mat[*b.map{|q|-d.sum{|x,i| q[0]*i**q[1]*(i-1)**q[2]*(x-i**3) } }.slices(1)]
  (a.inv * b).column(0).to_a
}
Mth.slf :bezier_quartic, ->d{
  a = [[2,-8,12,-8,2,-8,32,-48,32,-8,12,-48,72,-48,12,-8,32,-48,32,-8,2,-8,12,-8,2],
      [0,1,2,3,4,1,2,3,4,5,2,3,4,5,6,3,4,5,6,7,4,5,6,7,8], [8,7,6,5,4,7,6,5,4,3,6,5,4,3,2,5,4,3,2,1,4,3,2,1,0]].turn
  b = [[2,-8,12,-8,2],[0,1,2,3,4],[4,3,2,1,0]].turn
  a = Mat[*a.map{|q| d.sum{|x,i| q[0]*i**q[1]*(i-1)**q[2]          }}.slices(5)]
  b = Mat[*b.map{|q| d.sum{|x,i| q[0]*i**q[1]*(i-1)**q[2]*(x-i**4) }}.slices(1)]
  (a.inv * b).column(0).to_a
}
Mth.slf :bezier_quartic_m, ->d{
  a = [[32,-48,32,-48,72,-48,32,-48,32],[2,3,4,3,4,5,4,5,6],[6,5,4,5,4,3,4,3,2]].turn
  b = [[-8,12,-8],[1,2,3],[3,2,1]].turn
  a = Mat[*a.map{|q| d.sum{|x,i| q[0]*i**q[1]*(i-1)**q[2]          }}.slices(3)]
  b = Mat[*b.map{|q| d.sum{|x,i| q[0]*i**q[1]*(i-1)**q[2]*(x-i**4) }}.slices(1)]
  (a.inv * b).column(0).to_a
}



Ary.act :der, ->i=1,f=nil{ i.fold(f ? [*[0]*i,*self] : self){|o,_| o.cons(2).map{|a,b| b-a }} } # derivative (finite differences, deltas) (todo)
Ary.act :int, ->i=1,d=[0]{ i.fold(self){|o,_| o.fold(d.clone){|o,q| o << (o.last+q) }}.rest(d.size) } # integral (prefix sum, scan) (todo)
Ary.act :gauss, ->r=1{ r.fold(self){|o,_| [(o[0]*3+o[1])/4.0,*o[1..-2].mapi{|w,i| (o[i]+o[i+2]+w+w)/4.0 },(o[-2]+o[-1]*3)/4.0] } } # gaussian(binomial) blur
#       :gauss, ->r=1{ r.fold(self){|o,_| [o[0],*o,o[-1]].cons(3).map{|a,b,c| (a+b+b+c)/4.0 } } }