Skip to content

Instantly share code, notes, and snippets.

@cfsamson
Last active September 15, 2018 14:22
Show Gist options
  • Save cfsamson/6b176634b83e016932d0df51e186dd9f to your computer and use it in GitHub Desktop.
Save cfsamson/6b176634b83e016932d0df51e186dd9f to your computer and use it in GitHub Desktop.
crystal_gmp
# Fastest implementation in Crystal, but needs the following lines added to lib_gmp.cr
# and recompile the Crystal:
#
# fun submul_ui = __gmpz_submul_ui(rop : MPZ*, op1 : MPZ*, op2 : ULong)
# fun addmul_ui = __gmpz_addmul_ui(rop : MPZ*, op1 : MPZ*, op2 : ULong)
#
# ...or you could rebind to GMPlib and implement them in your own `lib` and use that
# in the #addmul and #submul methods instead if you don't have Crystal source code available.
# **SEE VERSION 5 for an implementation of this without changing the stdlib and confirm the perfomance**
#
require "big"
N = (ARGV[0]? || 100).to_i
struct BigInt
def add_with(other : BigInt) : self
LibGMP.add(self, self, other)
self
end
def add_with(other : Int) : self
if other < 0
sub(other.abs)
elsif other <= LibGMP::ULong::MAX
LibGMP.add_ui(self, self, other)
else
add_with(other.to_big_i)
end
self
end
def sub_with(other : Int) : self
if other < 0
add_with(other.abs)
elsif other <= LibGMP::ULong::MAX
LibGMP.sub_ui(self, self, other)
else
sub_with(other.to_big_i)
end
self
end
def sub_with(other : BigInt) : self
LibGMP.sub(self, self, other)
self
end
def mul_with(other : BigInt) : self
LibGMP.mul(self, self, other)
self
end
def mul_with(other : LibGMP::IntPrimitiveSigned) : self
LibGMP.mul_si(self, self, other)
self
end
def mul_with(other : LibGMP::IntPrimitiveUnsigned) : self
LibGMP.mul_ui(self, self, other)
self
end
# --------------------- NEW GMP METHODS -------------------
def add(op1 : BigInt, op2 : BigInt) : self
LibGMP.add(self, op1, op2)
self
end
def add(op1 : BigInt, op2 : Int) : self
if op2 < 0
sub(op1, op2.abs)
elsif op2 <= LibGMP::ULong::MAX
LibGMP.add_ui(self, op1, op2)
else
add(op1, op2.to_big_i)
end
self
end
def sub(op1 : BigInt, op2 : Int) : self
if op2 < 0
add(op1, op2.abs)
elsif op2 <= LibGMP::ULong::MAX
LibGMP.sub_ui(self, op1, op2)
else
sub(op1, op2.to_big_i)
end
self
end
def sub(op1 : BigInt, op2 : BigInt) : self
LibGMP.sub(self, op1, op2)
self
end
def mul(op1 : BigInt, op2 : BigInt) : self
LibGMP.mul(self, op1, op2)
self
end
def mul(op1 : BigInt, op2 : LibGMP::IntPrimitiveSigned) : self
LibGMP.mul_si(self, op1, op2)
self
end
def mul(op1 : BigInt, op2 : LibGMP::IntPrimitiveUnsigned) : self
LibGMP.mul_ui(self, op1, op2)
self
end
def tdiv_q(op1 : BigInt, op2 : BigInt)
check_division_by_zero op2
LibGMP.tdiv_q(self, op1, op2)
self
end
def submul(op1 : BigInt, op2 : UInt32)
LibGMP.submul_ui(self, op1, op2)
self
end
def addmul(op1 : BigInt, op2 : UInt32)
LibGMP.addmul_ui(self, op1, op2)
self
end
end
class Context
property k = 0
property tmp1 : BigInt = 0.to_big_i
property tmp2 : BigInt = 0.to_big_i
property acc : BigInt = 0.to_big_i
property den : BigInt = 1.to_big_i
property num : BigInt = 1.to_big_i
def extract_digit(nth : Int32) : UInt32
@tmp1.mul(@num, nth)
@tmp2.add(@tmp1, @acc)
@tmp1.tdiv_q(@tmp2, @den)
@tmp1.to_u32
end
def eliminate_digit(d)
@acc.submul(@den, d.to_u32)
@acc.mul_with(10)
@num.mul_with(10)
end
def next_term
@k += 1
k2 = @k * 2 + 1
@acc.addmul(@num, 2)
@acc.mul_with(k2)
@den.mul_with(k2)
@num.mul_with(@k)
end
end
buffer = String::Builder.new
context = Context.new
(1...(N + 1)).each do |i|
loop do
context.next_term
next if context.num > context.acc
d = context.extract_digit(3)
next if d != context.extract_digit(4)
context.eliminate_digit(d)
buffer << d
if i % 10 == 0
buffer << "\t:%d\n" % i
end
break
end
end
print buffer.to_s
# Does not require any additions to lib_gmp.cr but is more than 2x slower
# than version 1 since it's not using submul and addmul
require "big"
N = (ARGV[0]? || 100).to_i
struct BigInt
def add(other : BigInt) : self
LibGMP.add(self, self, other)
self
end
def add(other : Int) : self
if other < 0
sub(other.abs)
elsif other <= LibGMP::ULong::MAX
LibGMP.add_ui(self, self, other)
else
add(other.to_big_i)
end
self
end
def sub(other : Int) : self
if other < 0
add(other.abs)
elsif other <= LibGMP::ULong::MAX
LibGMP.sub_ui(self, self, other)
else
sub(other.to_big_i)
end
self
end
def sub(other : BigInt) : self
LibGMP.sub(self, self, other)
self
end
def mul(other : BigInt) : self
LibGMP.mul(self, self, other)
self
end
def mul(other : LibGMP::IntPrimitiveSigned) : self
LibGMP.mul_si(self, self, other)
self
end
def mul(other : LibGMP::IntPrimitiveUnsigned) : self
LibGMP.mul_ui(self, self, other)
self
end
# --------------------- NEW METHODS -----------------------------
def add(op1 : BigInt, op2 : BigInt) : self
LibGMP.add(self, op1, op2)
self
end
def add(op1 : BigInt, op2 : Int) : self
if op2 < 0
sub(op1, op2.abs)
elsif op2 <= LibGMP::ULong::MAX
LibGMP.add_ui(self, op1, op2)
else
add(op1, op2.to_big_i)
end
self
end
def sub(op1 : BigInt, op2 : Int) : self
if op2 < 0
add(op1, op2.abs)
elsif op2 <= LibGMP::ULong::MAX
LibGMP.sub_ui(self, op1, op2)
else
sub(op1, op2.to_big_i)
end
self
end
def sub(op1 : BigInt, op2 : BigInt) : self
LibGMP.sub(self, op1, op2)
self
end
def mul(op1 : BigInt, op2 : BigInt) : self
LibGMP.mul(self, op1, op2)
self
end
def mul(op1 : BigInt, op2 : LibGMP::IntPrimitiveSigned) : self
LibGMP.mul_si(self, op1, op2)
self
end
def mul(op1 : BigInt, op2 : LibGMP::IntPrimitiveUnsigned) : self
LibGMP.mul_ui(self, op1, op2)
self
end
def tdiv_q(op1 : BigInt, op2 : BigInt)
check_division_by_zero op2
LibGMP.tdiv_q(self, op1, op2)
self
end
end
class Context
property k = 0
property tmp1 : BigInt = 0.to_big_i
property tmp2 : BigInt = 0.to_big_i
property acc : BigInt = 0.to_big_i
property den : BigInt = 1.to_big_i
property num : BigInt = 1.to_big_i
def extract_digit(nth) : UInt32
@tmp1.mul(@num, nth)
@tmp2.add(@tmp1, @acc)
@tmp1.tdiv_q(@tmp2, @den)
@tmp1.to_u
end
def eliminate_digit(d : UInt32)
x = 0.to_big_i.mul(@den, d)
@acc.sub(x)
@acc.mul(10)
@num.mul(10)
end
def next_term
@k += 1
k2 = @k * 2 + 1
x = 0.to_big_i.mul(@num, 2)
@acc.add(x)
@acc.mul(k2)
@den.mul(k2)
@num.mul(@k)
end
end
# calculate size of buffer
buf_size = 0
order = 2
(0..N).each do |i|
if i % 10 == 0 && i > 0
buf_size += 13 + order
if (i + 10) % 10**order == 0
order += 1
end
end
end
buffer = String::Builder.new(buf_size)
context = Context.new
(1...(N + 1)).each do |i|
loop do
context.next_term
next if context.num > context.acc
d = context.extract_digit(3)
next if d != context.extract_digit(4)
context.eliminate_digit(d)
buffer << d
if i % 10 == 0
buffer << "\t:%d\n" % i
end
break
end
end
print buffer.to_s
# Imlemented using Python implementation as inspiration (that is similar to the PHP implementation)
# Faster than version 4, but still not matching the fastest implementations
require "big"
N = (ARGV[0]? || 100).to_i
struct BigInt
def add(other : BigInt) : self
LibGMP.add(self, self, other)
self
end
def add(other : Int) : self
if other < 0
sub(other.abs)
elsif other <= LibGMP::ULong::MAX
LibGMP.add_ui(self, self, other)
else
add(other.to_big_i)
end
self
end
def sub(other : Int) : self
if other < 0
add(other.abs)
elsif other <= LibGMP::ULong::MAX
LibGMP.sub_ui(self, self, other)
else
sub(other.to_big_i)
end
self
end
def sub(other : BigInt) : self
LibGMP.sub(self, self, other)
self
end
def mul(other : BigInt) : self
LibGMP.mul(self, self, other)
self
end
def mul(other : LibGMP::IntPrimitiveSigned) : self
LibGMP.mul_si(self, self, other)
self
end
def mul(other : LibGMP::IntPrimitiveUnsigned) : self
LibGMP.mul_ui(self, self, other)
self
end
# --------------------- NEW METHODS -----------------------------
def add(op1 : BigInt, op2 : BigInt) : self
LibGMP.add(self, op1, op2)
self
end
def add(op1 : BigInt, op2 : Int) : self
if op2 < 0
sub(op1, op2.abs)
elsif op2 <= LibGMP::ULong::MAX
LibGMP.add_ui(self, op1, op2)
else
add(op1, op2.to_big_i)
end
self
end
def sub(op1 : BigInt, op2 : Int) : self
if op2 < 0
add(op1, op2.abs)
elsif op2 <= LibGMP::ULong::MAX
LibGMP.sub_ui(self, op1, op2)
else
sub(op1. op2.to_big_i)
end
self
end
def sub(op1 : BigInt, op2 : BigInt) : self
LibGMP.sub(self, op1, op2)
self
end
def mul(op1 : BigInt, op2 : BigInt) : self
LibGMP.mul(self, op1, op2)
self
end
def mul(op1 : BigInt, op2 : LibGMP::IntPrimitiveSigned) : self
LibGMP.mul_si(self, op1, op2)
self
end
def mul(op1 : BigInt, op2 : LibGMP::IntPrimitiveUnsigned) : self
LibGMP.mul_ui(self, op1, op2)
self
end
end
#calculate size of buffer
buf_size = 0
order = 2
(0..N).each do |i|
if i % 10 == 0 && i > 0
buf_size += 13 + order
if (i + 10) % 10**order == 0
order += 1
end
end
end
buffer = String::Builder.new(buf_size)
w = 0.to_big_i
k = 1
n1 = 4.to_big_i
n2 = 3.to_big_i
d = 1.to_big_i
i = 0
produced = 0
extracted = 0
while true
u = n1.tdiv(d)
v = n2.tdiv(d)
if (u <=> v) == 0
buffer << u
i += 1
# echo
break if i == N
# extract
u = d * (-10 * u)
n1.mul(10)
n1.add(u)
n2.mul(10)
n2.add(u)
else
# produce
k2 = k << 1
u.mul(n1, (k2-1))
v.add(n2, n2)
w.mul(n1, (k -1))
n1.add(u, v)
u.mul(n2, (k + 2))
n2.add(w, u)
d.mul((k2 + 1))
k += 1
end
end
print buffer.to_s
# The original code with first suggestion of a mutating implementation of BigInt
require "big"
N = (ARGV[0]? || 100).to_i
struct BigInt
def add(other : BigInt) : self
LibGMP.add(self, self, other)
self
end
def add(other : Int) : self
if other < 0
sub(other.abs)
elsif other <= LibGMP::ULong::MAX
LibGMP.add_ui(self, self, other)
else
add(other.to_big_i)
end
self
end
def sub(other : Int) : self
if other < 0
add(other.abs)
elsif other <= LibGMP::ULong::MAX
LibGMP.sub_ui(self, self, other)
else
sub(other.to_big_i)
end
self
end
def sub(other : BigInt) : self
LibGMP.sub(self, self, other)
self
end
def mul(other : BigInt) : self
LibGMP.mul(self, self, other)
self
end
def mul(other : LibGMP::IntPrimitiveSigned) : self
LibGMP.mul_si(self, self, other)
self
end
def mul(other : LibGMP::IntPrimitiveUnsigned) : self
LibGMP.mul_ui(self, self, other)
self
end
end
#calculate size of buffer
buf_size = 0
order = 2
(0..N).each do |i|
if i % 10 == 0 && i > 0
buf_size += 13 + order
if (i + 10) % 10**order == 0
order += 1
end
end
end
buffer = String::Builder.new(buf_size)
i = 0
k = 0
ns = 0.to_big_i
a = 0.to_big_i
t = 0
u = 0.to_big_i
k1 = 1
n = 1.to_big_i
d = 1.to_big_i
while true
k += 1
t = n << 1
n.mul(k)
k1 += 2
a.add(t)
a.mul(k1)
d.mul(k1)
if a >= n
t, u = (n * 3 + a).divmod(d)
u.add(n)
if d > u
ns.mul(10)
ns.add(t)
i += 1
if i % 10 == 0
#print "%010d\t:%d\n" % {ns.to_u64, i}
buffer << "%010d\t:%d\n" % {ns.to_u64, i}
ns = 0.to_big_i
end
break if i >= N
a.sub(d*t)
a.mul(10)
n.mul(10)
end
end
end
print buffer.to_s
# Fastest implementation in a copy-paste format
@[Link("gmp")]
lib LibGMP2
alias ULong = LibC::ULong
alias BitcntT = ULong
fun submul_ui = __gmpz_submul_ui(rop : LibGMP::MPZ*, op1 : LibGMP::MPZ*, op2 : ULong)
fun addmul_ui = __gmpz_addmul_ui(rop : LibGMP::MPZ*, op1 : LibGMP::MPZ*, op2 : ULong)
end
require "big"
N = (ARGV[0]? || 100).to_i
struct BigInt
def add_with(other : BigInt) : self
LibGMP.add(self, self, other)
self
end
def add_with(other : Int) : self
if other < 0
sub_with(other.abs)
elsif other <= LibGMP::ULong::MAX
LibGMP.add_ui(self, self, other)
else
add_with(other.to_big_i)
end
self
end
def sub_with(other : Int) : self
if other < 0
add_with(other.abs)
elsif other <= LibGMP::ULong::MAX
LibGMP.sub_ui(self, self, other)
else
sub_with(other.to_big_i)
end
self
end
def sub_with(other : BigInt) : self
LibGMP.sub(self, self, other)
self
end
def mul_with(other : BigInt) : self
LibGMP.mul(self, self, other)
self
end
def mul_with(other : LibGMP::IntPrimitiveSigned) : self
LibGMP.mul_si(self, self, other)
self
end
def mul_with(other : LibGMP::IntPrimitiveUnsigned) : self
LibGMP.mul_ui(self, self, other)
self
end
# --------------------- NEW METHODS -----------------------------
def add(op1 : BigInt, op2 : BigInt) : self
LibGMP.add(self, op1, op2)
self
end
def add(op1 : BigInt, op2 : Int) : self
if op2 < 0
sub(op1, op2.abs)
elsif op2 <= LibGMP::ULong::MAX
LibGMP.add_ui(self, op1, op2)
else
add(op1, op2.to_big_i)
end
self
end
def sub(op1 : BigInt, op2 : Int) : self
if op2 < 0
add(op1, op2.abs)
elsif op2 <= LibGMP::ULong::MAX
LibGMP.sub_ui(self, op1, op2)
else
sub(op1, op2.to_big_i)
end
self
end
def sub(op1 : BigInt, op2 : BigInt) : self
LibGMP.sub(self, op1, op2)
self
end
def mul(op1 : BigInt, op2 : BigInt) : self
LibGMP.mul(self, op1, op2)
self
end
def mul(op1 : BigInt, op2 : LibGMP::IntPrimitiveSigned) : self
LibGMP.mul_si(self, op1, op2)
self
end
def mul(op1 : BigInt, op2 : LibGMP::IntPrimitiveUnsigned) : self
LibGMP.mul_ui(self, op1, op2)
self
end
def tdiv_q(op1 : BigInt, op2 : BigInt)
check_division_by_zero op2
LibGMP.tdiv_q(self, op1, op2)
self
end
def submul(op1 : BigInt, op2 : UInt32)
LibGMP2.submul_ui(self, op1, op2)
self
end
def addmul(op1 : BigInt, op2 : UInt32)
LibGMP2.addmul_ui(self, op1, op2)
self
end
end
class Context
property k = 0
property tmp1 : BigInt = 0.to_big_i
property tmp2 : BigInt = 0.to_big_i
property acc : BigInt = 0.to_big_i
property den : BigInt = 1.to_big_i
property num : BigInt = 1.to_big_i
def extract_digit(nth : UInt32)
@tmp1.mul(@num, nth)
@tmp2.add(@tmp1, @acc)
@tmp1.tdiv_q(@tmp2, @den)
@tmp1.to_i
end
def eliminate_digit(d)
@acc.submul(@den, d)
@acc.mul_with(10)
@num.mul_with(10)
end
def next_term
@k += 1
k2 = @k * 2 + 1
@acc.addmul(@num, 2)
@acc.mul_with(k2)
@den.mul_with(k2)
@num.mul_with(@k)
end
end
buffer = String::Builder.new
context = Context.new
(1...(N + 1)).each do |i|
loop do
context.next_term
next if context.num > context.acc
d = context.extract_digit(3)
next if d != context.extract_digit(4)
context.eliminate_digit(d.to_u)
buffer << d
if i % 10 == 0
buffer << "\t:%d\n" % i
end
break
end
end
print buffer.to_s
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment