Module: Abst

Defined in:
lib/abst.rb,
lib/config.rb,
lib/include/ring.rb,
lib/include/cache.rb,
lib/include/graph.rb,
lib/include/group.rb,
lib/include/prime.rb,
lib/include/matrix.rb,
lib/include/vector.rb,
lib/include/integer.rb,
lib/include/residue.rb,
lib/include/sequence.rb,
lib/include/polynomial.rb,
lib/include/prime_mpqs.rb,
lib/include/fundamental.rb

Defined Under Namespace

Modules: Group, Ring Classes: Cache, IntegerIdeal, IntegerResidueField, IntegerResidueRing, MPQS, Matrix, MatrixSizeError, Polynomial, SquareMatrix, SystemCache, Vector, VectorSizeError

Constant Summary collapse

ABST_ROOT = 
File.dirname(__FILE__) + '/'
BASE_BYTE =

Integer block byte size

1.size
INFINITY = 
Float::INFINITY
I = 
Complex::I
PRIME_CACHE_LIMIT =

Precompute primes under the value

1_000_000
DATA_DIR =

System cache files will be created in this directory

ABST_ROOT + 'data/'
CACHE_DIR =

User cache files will be created in this directory

ABST_ROOT + 'cache/'
PRIMES_LIST =

Cache primes

DATA_DIR + 'primes_list'
FIXNUM_BIT_SIZE =

Bit size of Fixnum integer e s.t. (2 ** e) is not Fixnum and (2 ** e - 1) is Fixnum

30
THREAD_NUM = 
defined?(JRUBY_VERSION) ? 6 : 1
AUTO_INCLUDE = 
false
Z = 
Integer

Class Method Summary collapse

Instance Method Summary collapse

Class Method Details

.apr(n) ⇒ Object

Primality test

Param

positive integer n

Return

boolean whether n is prime or not

Raises:

  • (NotImplementedError) 


136
137
138
 
# File 'lib/include/prime.rb', line 136

def apr(n)
	raise NotImplementedError
end

.bhaskara_brouncker(n) ⇒ Object

Param

positive non-square integer n

Return

least positive integer a, b s.t. a**2 - n * b**2 == 1 or -1



824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
 
# File 'lib/include/fundamental.rb', line 824

def bhaskara_brouncker(n)
	t = sqrt = isqrt(n)
	u, u1 = 0, sqrt
	c, c1 = 1, n - sqrt ** 2
	a, a1 = 1, sqrt
	b, b1 = 0, 1

	until 1 == c1
		t = (sqrt + u1) / c1
		u1, u = t * c1 - u1, u1
		c1, c = c + t * (u - u1), c1
		a1, a = a + t * a1, a1
		b1, b = b + t * b1, b1
	end

	return a1, b1
end

.binary_gcd(a, b) ⇒ Object

Binary GCD

Param

integer a, b

Return

gcd of a and b



211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
 
# File 'lib/include/fundamental.rb', line 211

def binary_gcd(a, b)
	a = -a if a < 0
	b = -b if b < 0
	a, b = b, a if a < b
	return a if 0 == b

	# Reduce size once
	a, b = b, a % b
	return a if 0 == b

	# Compute powers of 2
	k = 0
	k += 1 while 0 == a[k] and 0 == b[k]
	if 0 < k
		a >>= k
		b >>= k
	end

	# Remove initial power of 2
	a >>= 1 while a.even?
	b >>= 1 while b.even?

	loop do
		# Subtract (Here a and b are both odd.)
		t = a - b
		return a << k if 0 == t

		count = 0
		count += 1 while 0 == t[count]
		t >>= count

		(0 < t) ? a = t : b = -t
	end
end

.chinese_remainder_theorem(list) ⇒ Object

Chinese Remainder Theorem using Algorithm 1.3.12 of [CCANT]

Param

array of pair integer s.t. [[x_1, m_1], [x_2, m_2], … , [x_k, m_k]] m_i are pairwise coprime

Return

integer x s.t. x = x_i (mod m_i) for all i and 0 <= x < m_1 * m_2 * … * m_k

Example: chinese_remainder_theorem() return 23



423
424
425
426
427
428
429
430
431
432
433
434
 
# File 'lib/include/fundamental.rb', line 423

def chinese_remainder_theorem(list)
	x, m = list.shift

	list.each do |xi, mi|
		u, v, = extended_gcd(m, mi)
		x = u * m * xi + v * mi * x
		m *= mi
		x %= m
	end

	return x
end

.continued_fraction(a, b, a_, b_) ⇒ Object

Param

integer a, b, a_, b_ s.t. (real number x s.t. a/b <= x <= a_/b_)

Return

continueed fraction expansion of x and lower and upper bounds for this next partial quotient. Format: [[continueed fraction expansion…], [lower, upper]]



440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
 
# File 'lib/include/fundamental.rb', line 440

def continued_fraction(a, b, a_, b_)
	# Initialize
	q = q_ = nil
	rslt = []

	until 0 == b or 0 == b_
		# Euclidean step
		q, r = a.divmod(b)

		r_ = a_ - b_ * q
		if r_ < 0 or b_ <= r_
			q_ = a_ / b_
			break
		end

		rslt.push(q)
		a , b  = b , r
		a_, b_ = b_, r_
	end

	if 0 == b
		return rslt, [] if 0 == b_
		return rslt, [a_ / b_, INFINITY]
	end

	return rslt, [a / b, INFINITY] if 0 == b_

	q, q_ = q_, q if q > q_
	return rslt, [q, q_]
end

.continued_fraction_of_sqrt(m) ⇒ Object

Param

positive integer m

Return

continued fraction of sqrt m

n, [a_1, a_2, …, a_i]

s.t. sqrt m = n + frac1+ frac{1+ …}

n is integer part and a_1…a_i are repeating part

Example::continued_fraction_of_sqrt(2) -> [1, [2]]

continued_fraction_of_sqrt(3) -> [1, [1, 2]]
continued_fraction_of_sqrt(4) -> [2, []]



792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
 
# File 'lib/include/fundamental.rb', line 792

def continued_fraction_of_sqrt(m)
	if r = m.square?
		return [r, []]
	end

	# Initialize
	rslt = [isqrt(m), []]

	a = 1
	_B = b = -rslt[0]
	c = 1

	loop do
		# inverse
		a, b, c = a * c, -b * c, a**2 * m - b**2

		# reduction
		t = gcd(gcd(a, b), c)
		a /= t
		b /= t
		c /= t

		t = (isqrt(m * a ** 2) + b) / c
		rslt[1].push(t)
		b -= c * t

		return rslt if 1 == a and _B == b and 1 == c
	end
end

.cornacchia(d, p) ⇒ Object

Param

positive integer d prime number p (d < p)

Return

integer solution (x, y) to the Diophantine equation x ** 2 + d * y ** 2 = p if exists else nil



616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
 
# File 'lib/include/fundamental.rb', line 616

def cornacchia(d, p)
	return nil if -1 == kronecker_symbol(-d, p)

	# Compute square root
	x0 = mod_sqrt(-d, p)
	x0 = p - x0 if x0 <= p >> 1
	a = p
	b = x0
	border = isqrt(p)

	# Euclidean algorithm
	while border < b
		a, b = b, a % b
	end

	# Test solution
	c, r = (p - b ** 2).divmod(d)
	return nil if 0 != r

	q = c.square?
	return nil unless q

	return [b, q]
end

.create_matrix(coef_class, height, width = nil) ⇒ Object Also known as: Matrix 



121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
 
# File 'lib/include/matrix.rb', line 121

def create_matrix(coef_class, height, width = nil)
	if height.kind_of?(Array)
		elems = height
		width = height[0].size
		height = height.size
	end

	super_class = height == width ? SquareMatrix : Matrix
	matrix = Class.new(super_class) do
		@coef_class = coef_class
		@height = height
		@width = width
		@coef_vector = Abst::Vector(coef_class, width)
	end

	return matrix.new(elems) if elems
	return matrix
end

.create_polynomial(coef_class, coef = nil) ⇒ Object Also known as: Polynomial

Param

immutable class coef_class



118
119
120
121
122
123
124
125
126
127
 
# File 'lib/include/polynomial.rb', line 118

def create_polynomial(coef_class, coef = nil)
	poly = Class.new(Polynomial) do
		@coef_class = coef_class
		@zero = self.new([coef_class.zero]).freeze
		@one = self.new([coef_class.one]).freeze
	end

	return poly.new(coef) if coef
	return poly
end

.create_vector_space(coef_class, size) ⇒ Object Also known as: Vector 



72
73
74
75
76
77
78
79
80
81
82
83
84
85
 
# File 'lib/include/vector.rb', line 72

def create_vector_space(coef_class, size)
	if size.kind_of?(Array)
		elems = size
		size = size.size
	end

	vector = Class.new(Vector) do
		@coef_class = coef_class
		@size = size
	end

	return vector.new(elems) if elems
	return vector
end

.differences_of_squares(n, limit = 1_000_000) ⇒ Object

Param

positive integer n

Return

factor a and b (a <= b) if found else false



234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
 
# File 'lib/include/prime.rb', line 234

def differences_of_squares(n, limit = 1_000_000)
	sqrt = isqrt(n)
	dx = (sqrt << 1) + 1
	dy = 1
	r = sqrt ** 2 - n
	terms = 0

	# find x and y s.t. x**2 - y**2 == n
	loop do
		while 0 < r
			r -= dy
			dy += 2
		end

		break if 0 == r

		r += dx
		dx += 2

		terms += 1
		return false if limit == terms
	end

	return (dx - dy) >> 1, (dx + dy - 2) >> 1
end

.dlog_icm(base, m, p) ⇒ Object

Param

Integer base (2 <= base < p) Integer m (2 <= m < p) prime p

Return

Integer e s.t. base ** e == m (mod p)

Raises:

  • (NotImplementedError) 


890
891
892
893
894
895
896
897
898
899
900
901
902
 
# File 'lib/include/fundamental.rb', line 890

def dlog_icm(base, m, p)
	# Select factor base

	# Select smooth elements on factor base

	# Solve DL of factor base for base

	# Find smooth element h s.t. h == m * base ** f where f is integer

	# Solve DL of h for factor base

	raise NotImplementedError
end

.dlog_rho(base, m, p) ⇒ Object

Param

Integer base (2 <= base < p) Integer m (2 <= m < p) prime p

Return

Integer e s.t. base ** e == m (mod p)

Raises:

  • (NotImplementedError) 


882
883
884
 
# File 'lib/include/fundamental.rb', line 882

def dlog_rho(base, m, p)
	raise NotImplementedError
end

.element_order(g, group_order) ⇒ Object 



18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
 
# File 'lib/include/group.rb', line 18

def element_order(g, group_order)
	one = g.class.one
	order = group_order
	factor = factorize(group_order)

	factor.each do |f, e|
		order /= f ** e

		t = g ** order
		until one == t
			t **= f
			order *= f
		end
	end

	return order
end

.eratosthenes_sieve(n) {|2| ... } ⇒ Object

generation

Yields:

  • (2) 


455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
 
# File 'lib/include/prime.rb', line 455

def eratosthenes_sieve(n)
	return Enumerator.new(self, :eratosthenes_sieve, n) unless block_given?

	return if n < 2

	yield 2
	return if 2 == n

	yield 3

	# make list for sieve
	sieve_len_max = (n + 1) >> 1
	sieve = [true, false, true]
	sieve_len = 3
	k = 5
	i = 2
	while sieve_len < sieve_len_max
		if sieve[i]
			yield k
			sieve_len *= k
			if sieve_len_max < sieve_len
				sieve_len /= k
				# adjust sieve list length
				sieve *= sieve_len_max / sieve_len
				sieve += sieve[0...(sieve_len_max - sieve.size)]
				sieve_len = sieve_len_max
			else
				sieve *= k
			end

			i.step(sieve_len - 1, k) do |j|
				sieve[j] = false
			end
		end

		k += 2
		i += 1
	end

	# sieve
	limit = (isqrt(n) - 1) >> 1
	while i <= limit
		if sieve[i]
			yield k = (i << 1) + 1
			j = (k ** 2) >> 1
			while j < sieve_len_max
				sieve[j] = false
				j += k
			end
		end

		i += 1
	end

	# output result
	limit = (n - 1) >> 1
	while i <= limit
		yield((i << 1) + 1) if sieve[i]
		i += 1
	end
end

.extended_binary_gcd(a, b) ⇒ Object

Param

non-negative integer a, b

Return

(u, v, d) s.t. a*u + b*v = gcd(a, b) = d



327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
 
# File 'lib/include/fundamental.rb', line 327

def extended_binary_gcd(a, b)
	if a < b
		a, b = b, a
		exchange_flag_1 = true
	end

	if 0 == b
		return 0, 1, a if exchange_flag_1
		return 1, 0, a
	end

	# Reduce size once
	_Q, r = a.divmod(b)
	if 0 == r
		return 1, 0, b if exchange_flag_1
		return 0, 1, b
	end
	a, b = b, r

	# Compute power of 2
	_K = 0
	_K += 1 while 0 == a[_K] and 0 == b[_K]
	if 0 < _K
		a >>= _K
		b >>= _K
	end

	if b.even?
		a, b = b, a
		exchange_flag_2 = true
	end

	# Initialize
	u = 1
	d = a	# d == a * u + b * v, (v = 0)
	u_ = 0
	d_ = b	# d_ == a * u_ + b * v_, (v_ = 1)

	# Remove intial power of 2
	while d.even?
		d >>= 1
		u += b if u.odd?
		u >>= 1
	end

	loop do
		# Substract
		next_u = u - u_
		next_d = d - d_		# next_d == a * next_u + b * next_v
		next_u += b if next_u < 0

		break if 0 == next_d

		# Remove powers of 2
		while next_d.even?
			next_d >>= 1
			next_u += b if next_u.odd?
			next_u >>= 1
		end

		if 0 < next_d
			u = next_u
			d = next_d
		else
			u_ = b - next_u
			d_ = -next_d
		end
	end

	v = (d - a * u) / b

	u, v = v, u if exchange_flag_2
	d <<= _K
	u, v = v, u - v * _Q
	u, v = v, u if exchange_flag_1

	return u, v, d
end

.extended_gcd(a, b) ⇒ Object

Param

a and b are member of a Euclidean domain

Return

(u, v, d) s.t. a*u + b*v = gcd(a, b) = d



248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
 
# File 'lib/include/fundamental.rb', line 248

def extended_gcd(a, b)
	u0 = a.class.one
	u1 = a.class.zero

	return u0, u1, a if b.zero?

	d0 = a	# d0 = a * u0 + b * v0
	d1 = b	# d1 = a * u1 + b * v1

	loop do
		q, r = d0.divmod(d1)

		return u1, (d1 - a * u1) / b, d1 if r.zero?

		d0, d1 = d1, r
		u0, u1 = u1, u0 - q * u1
	end
end

.extended_lehmer_gcd(a, b) ⇒ Object

Param

non-negative integer a, b

Return

(u, v, d) s.t. a*u + b*v = gcd(a, b) = d



269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
 
# File 'lib/include/fundamental.rb', line 269

def extended_lehmer_gcd(a, b)
	d0 = a
	u0 = 1	# d0 = a * u0 + b * v0
	d1 = b
	u1 = 0	# d1 = a * u1 + b * v1

	loop do
		if d1.instance_of?(Fixnum)
			_u, _v, d = extended_gcd(d0, d1)

			# here
			# d == _u * d0 + _v * d1
			# d0 == u0 * a + v0 * b
			# d1 == u1 * a + v1 * b

			u = _u * u0 + _v * u1
			v = (d - u * a) / b

			return u, v, d
		end

		# Get most significant digits of d0 and d1
		shift_size = (d0 < d1 ? d1 : d0).bit_size - FIXNUM_BIT_SIZE
		a_ = d0 >> shift_size
		b_ = d1 >> shift_size

		# Initialize (Here a_ and b_ are next value of d0, d1)
		_A = 1
		_B = 0	# a_ == msd(d0) * _A + msd(d1) * _B
		_C = 0
		_D = 1	# b_ == msd(d0) * _C + msd(d1) * _D

		# Test Quotient
		until 0 == b_ + _C or 0 == b_ + _D
			q1 = (a_ + _B) / (b_ + _D)
			q2 = (a_ + _A) / (b_ + _C)
			break if q1 != q2

			# Euclidean step
			_A, _C = _C, _A - q1 * _C
			_B, _D = _D, _B - q1 * _D
			a_, b_ = b_, a_ - q1 * b_
		end

		# Multi-precision step
		if 0 == _B
			q, r = d0.divmod(d1)
			d0, d1  = d1, r
			u0, u1 = u1, u0 - q * u1
		else
			d0, d1  = d0 * _A + d1 * _B, d0 * _C + d1 * _D
			u0, u1 =u0 *  _A + u1 * _B, u0 * _C + u1 * _D
		end
	end
end

.factorize(n, return_hash = false) ⇒ Object

Param

integer n boolean return_hash

Return

prime factorization of n s.t. [[p_1, e_1], [p_2, e_2], …] n == p_1**e_1 * p_2**e_2 * … (p_1 < p_2 < …) if |n| <= 1 then return [[n, 1]]. if n < 0 then [-1, 1] is added as a factor.



368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
 
# File 'lib/include/prime.rb', line 368

def factorize(n, return_hash = false)
	unless return_hash
		return factorize(n, true).to_a.sort
	end

	factor = Hash.new(0)
	if n <= 1
		return {n => 1} if -1 <= n
		n = -n
		factor[-1] = 1
	end

	found_factor, n = trial_division(n, td_lim = 10_000)
	found_factor.each {|k, v| factor[k] += v}
	td_lim_square = td_lim ** 2

	check_finish = lambda do
		if n <= td_lim_square or prime?(n)
			factor[n] += 1 unless 1 == n
			return true
		end
		return false
	end

	divide = lambda do |f|
		f.size.times do |i|
			d = f[i][0]
			loop do
				q, r = n.divmod(d)
				break unless 0 == r
				n = q
				f[i][1] += 1
			end
		end

		return f
	end

	return factor if check_finish.call

	# pollard_rho
	loop do
		c = nil
		loop do
			c = rand(n - 3) + 1
			break unless c.square?
		end
		s = rand(n)
		f = pollard_rho(n, c, s, 50_000)
		break unless f

		# f is prime?
		n /= f
		if f <= td_lim_square or prime?(f)
			f = divide.call([[f, 1]])
		else
			f = divide.call(factorize(f))
		end

		f.each {|k, v| factor[k] += v}
		return factor if check_finish.call
	end

	# MPQS
	loop do
		f = mpqs(n)
		break unless f

		# f is prime?
		n /= f
		if f <= td_lim_square or prime?(f)
			f = divide.call([[f, 1]])
		else
			f = divide.call(factorize(f))
		end

		f.each {|k, v| factor[k] += v}
		return factor if check_finish.call
	end

	raise [factor, n].to_s
end

.fibonacci(n) ⇒ Object

Param

non-negative integer n

Return

the n-th Fibonacci number

effect $fibonacci = fibonacci(n)



9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
 
# File 'lib/include/sequence.rb', line 9

def fibonacci(n)
	return $fibonacci[n] if $fibonacci.include?(n)

	m = n >> 1
	if n.even?
		f1 = fibonacci(m - 1)
		f2 = fibonacci(m)
		$fibonacci[n] = (f1 + f1 + f2) * f2
	else
		f1 = fibonacci(m)
		f2 = fibonacci(m + 1)
		$fibonacci[n] = f1 ** 2 + f2 ** 2
	end

	return $fibonacci[n]
end

.gcd(a, b) ⇒ Object

GCD

Param

a and b are member of a Euclidean domain

Return

gcd of a and b



146
147
148
149
150
151
152
 
# File 'lib/include/fundamental.rb', line 146

def gcd(a, b)
	until b.zero?
		a, b = b, a % b
	end

	return a
end

.heptagonal(n) ⇒ Object 



47
48
49
 
# File 'lib/include/sequence.rb', line 47

def heptagonal(n)
	return n * (5 * n - 3) >> 1
end

.hexagonal(n) ⇒ Object

Hexagonal numbers are generated by the formula, H_n = n * (2 * n - 1) The first ten hexagonal numbers are:

1, 6, 15, 28, 45, 66, 91, 120, 153, 190, ...



43
44
45
 
# File 'lib/include/sequence.rb', line 43

def hexagonal(n)
	return n * ((n << 1) - 1)
end

.ilog2(n) ⇒ Object

Param

positive integer n

Return

integer e s.t. 2 ** e <= n < 2 ** (e + 1)



767
768
769
770
 
# File 'lib/include/fundamental.rb', line 767

def ilog2(n)
	bits = (n.size - BASE_BYTE) << 3
	return bits + Bisect.bisect_right(powers_of_2, n >> bits) - 1
end

.include?(n) ⇒ Boolean

Returns:

  • (Boolean) 


33
34
35
 
# File 'lib/include/prime.rb', line 33

def $primes.include?(n)
	return Bisect.index(self, n)
end

.inverse(n, mod) ⇒ Object

Param

integer n positive integer mod which is reratively prime with n

Return

inverse of n modulo mod. (1 <= inverse < mod)

Raises:

  • (ArgumentError) 


409
410
411
412
413
414
415
 
# File 'lib/include/fundamental.rb', line 409

def inverse(n, mod)
	u, _, d = extended_gcd(n, mod)

	raise ArgumentError, 'n and mod must be reratively prime.' unless 1 == d

	return u % mod
end

.iroot(n, pow, return_power = false) ⇒ Object

Integer Power Root

Param

positive integer n positive integer pow

Return

integer part of the power root of n i.e. the number m s.t. m ** pow <= n < (m + 1) ** pow if return_power is true then return n ** pow



695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
 
# File 'lib/include/fundamental.rb', line 695

def iroot(n, pow, return_power = false)
	# get integer e s.t. (2 ** (e - 1)) ** pow <= n < (2 ** e) ** pow
	e = ilog2(n) / pow + 1		# == Rational(ilog2(n) + 1, pow).ceil

	x = 1 << e					# == 2 ** e
	z = nil
	q = n >> (e * (pow - 1))	# == n / (x ** (pow -  1))

	loop do
		# Newtonian step
		x += (q - x) / pow
		z = x ** (pow - 1)
		q = n / z

		break if x <= q
	end

	return x, x * z if return_power
	return x
end

.isqrt(n) ⇒ Object

Integer Square Root

Param

positive integer n

Return

integer part of the square root of n

i.e. the number m s.t. m ** 2 <= n < (m + 1) ** 2



677
678
679
680
681
682
683
684
685
686
687
 
# File 'lib/include/fundamental.rb', line 677

def isqrt(n)
	x = 1 << ((ilog2(n) >> 1) + 1)

	loop do
		# Newtonian step
		next_x = (x + n / x) >> 1
		return x if x <= next_x

		x = next_x
	end
end

.kronecker_symbol(n, m) ⇒ Object Also known as: legendre_symbol, jacobi_symbol

Param

integer n, m

Return

kronecker symbol (n|m)



499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
 
# File 'lib/include/fundamental.rb', line 499

def kronecker_symbol(n, m)
	if 0 == m
		return (-1 == n or 1 == n) ? 1 : 0
	end

	return 0 if n.even? and m.even?

	m8 = [0, 1, 0, -1, 0, -1, 0, 1]

	# Remove 2's from m
	count = 0
	count += 1 while 0 == m[count]
	m >>= count
	rslt = count.even? ? 1 : m8[n & 7]

	if m < 0
		m = -m
		rslt = -rslt if n < 0
	end
	n %= m

	until 0 == n
		count = 0
		count += 1 while 0 == n[count]
		n >>= count
		rslt *= m8[m & 7] if count.odd?

		# Apply reciprocity
		rslt = -rslt if 1 == n[1] and 1 == m[1]
		n, m = m % n, n
	end

	return (1 == m) ? rslt : 0
end

.lcm(a, b) ⇒ Object 



154
155
156
 
# File 'lib/include/fundamental.rb', line 154

def lcm(a, b)
	return a * b / gcd(a, b)
end

.left_right_base2k_power(g, n, mod = nil, k = nil) ⇒ Object

Left-Right Base 2**k Power

Param

group element g integer n Euclidean domain element mod base bit size k

Return

g ** n % mod



74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
 
# File 'lib/include/fundamental.rb', line 74

def left_right_base2k_power(g, n, mod = nil, k = nil)
	rslt = g.class.one
	return rslt if 0 == n

	# Initialize
	if n < 0
		n = -n
		g **= (-1)
	end

	g %= mod if mod

	unless k
		e = ilog2(n)
		optim = [0, 8, 24, 69, 196, 538, 1433, 3714]

		if e <= optim.last
			k = Bisect.bisect_left(optim, e)
		else
			k = optim.size
			k += 1 until e <= k * (k + 1) * (1 << (k << 1)) / ((1 << (k + 1)) - k - 2)
		end
	end

	# convert n into base 2**k
	digits = []
	mask = (1 << k) - 1
	until 0 == n
		digits.unshift(n & mask)
		n >>= k
	end

	# Precomputations
	z_powers = [nil, g]
	g_square = g * g
	g_square %= mod if mod
	3.step((1 << k) - 1, 2) do |i|
		z_powers[i] = z_powers[i - 2] * g_square
		z_powers[i] %= mod if mod
	end

	digits.each do |a|
		# Multiply
		if 0 == a
			k.times do
				rslt = rslt ** 2
				rslt %= mod if mod
			end
		else
			t = 0
			t += 1 while 0 == a[t]
			a >>= t if 0 < t

			(k - t).times do
				rslt = rslt ** 2
				rslt %= mod if mod
			end
			rslt *= z_powers[a]
			rslt %= mod if mod
			t.times do
				rslt = rslt ** 2
				rslt %= mod if mod
			end
		end
	end

	return rslt
end

.left_right_power(g, n, mod = nil) ⇒ Object Also known as: power

Left-Right Binary Power

Param

group element g integer n Euclidean domain element mod

Return

g ** n % mod



39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
 
# File 'lib/include/fundamental.rb', line 39

def left_right_power(g, n, mod = nil)
	return g.class.one if 0 == n

	if n < 0
		n = -n
		g **= (-1)
	end

	g %= mod if mod
	e = ilog2(n)

	rslt = g
	while 0 != e
		e -= 1

		rslt *= rslt
		rslt %= mod if mod

		if 1 == n[e]
			rslt *= g
			rslt %= mod if mod
		end
	end

	return rslt
end

.lehmer_gcd(a, b) ⇒ Object

Param

integer a, b

Return

gcd of a and b



160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
 
# File 'lib/include/fundamental.rb', line 160

def lehmer_gcd(a, b)
	a = -a if a < 0
	b = -b if b < 0
	a, b = b, a if a < b

	until 0 == b
		return gcd(a, b) if b.instance_of?(Fixnum)

		# Get most significant digits of a and b
		shift_size = (a < b ? b : a).bit_size - FIXNUM_BIT_SIZE
		a_ = a >> shift_size
		b_ = b >> shift_size

		_A = 1
		_B = 0	# a_ == msd(a) * _A + msd(b) * _B
		_C = 0
		_D = 1	# b_ == msd(a) * _C + msd(b) * _D

		# Always
		# a_ + _B <= msd(a * _A + b * _B) < a_ + _A AND
		# b_ + _C <= msd(a * _C + b * _D) < a_ + _D
		# OR
		# a_ + _B > msd(a * _A + b * _B) >= a_ + _A AND
		# b_ + _C > msd(a * _C + b * _D) >= a_ + _D

		# Test quotient
		until 0 == b_ + _C or 0 == b_ + _D
			q1 = (a_ + _A) / (b_ + _C)
			q2 = (a_ + _B) / (b_ + _D)
			break if q1 != q2

			# Euclidean step
			_A, _C = _C, _A - q1 * _C
			_B, _D = _D, _B - q1 * _D
			a_, b_ = b_, a_ - q1 * b_
		end

		# Multi-precision step
		if 0 == _B
			a, b = b, a % b
		else
			a, b = a * _A + b * _B, a * _C + b * _D
		end
	end

	return a
end

.load_precomputed_primes ⇒ Object 



18
19
20
 
# File 'lib/include/prime.rb', line 18

def load_precomputed_primes
	open(PRIMES_LIST) {|io| return io.read.split("\n").map(&:to_i)}
end

.lucas_lehmer(p) ⇒ Object

Lucas-Lehmer primality test

Param

odd prime p

Return

boolean whether M(p) == 2**p - 1 is prime or not there M(p) means p-th Mersenne number



170
171
172
173
174
175
 
# File 'lib/include/prime.rb', line 170

def lucas_lehmer(p)
	s = 4
	m = (1 << p) - 1
	(p - 2).times {s = (s * s - 2) % m}
	return 0 == s
end

.miller_rabin(n, trials = 20, return_witness = false) ⇒ Object

Miller-Rabin pseudo-primality test

Param

odd integer n >= 3

Return

boolean whether n passes trials times random base strong pseudoprime test or not integer witness (or nil) if return_witness



85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
 
# File 'lib/include/prime.rb', line 85

def miller_rabin(n, trials = 20, return_witness = false)
	# Precomputation
	n_minus_1 = n - 1
	e = 0
	e += 1 while 0 == n_minus_1[e]
	k = n_minus_1 >> e

	trials.times do
		base = rand(n - 2) + 2
		unless strong_pseudoprime_test(n, base, e, k)
			return return_witness ? [false, base] : false
		end
	end

	return return_witness ? [true, nil] : true
end

.mobius(n) ⇒ Object Also known as: moebius

Param

positive integer n

Return

the value of moebius function for n



871
872
873
874
875
 
# File 'lib/include/fundamental.rb', line 871

def mobius(n)
	return 1 if 1 == n
	return 0 unless n.squarefree?
	return factorize(n).size.odd? ? -1 : 1
end

.mod_sqrt(n, p, exp = 1, return_list = false) ⇒ Object

Param

integer n odd prime p positive integer exp (if 1 < exp then n must be relatively prime with p)

Return

the square root of n mod (p ** exp) if exists else nil



542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
 
# File 'lib/include/fundamental.rb', line 542

def mod_sqrt(n, p, exp = 1, return_list = false)
	if 1 < exp or return_list
		x = mod_sqrt(n, p)
		return x unless x
		return [x] if 1 == exp
		raise ArgumentError, "if 1 < exp then n must be relatively prime with p" if 0 == x

		rslt = [x] if return_list
		p_power = p
		z = extended_lehmer_gcd(x << 1, p)[0]
		(exp - 1).times do
			x += (n - x ** 2) / p_power * z % p * p_power
			p_power *= p
			rslt.push(x) if return_list
		end

		return return_list ? rslt : x
	end

	unless (k = kronecker_symbol(n, p)) == 1
		return nil if -1 == k
		return 0
	end

	if 0 < p & 6
		return power(n, (p >> 2) + 1, p) if p[1] == 1
		n %= p
		x = power(n, (p >> 3) + 1, p)
		return x if x ** 2 % p == n
		return x * power(2, p >> 2, p) % p
	end

	# get q and e s.t. p - 1 == 2**e * q with q odd
	e = 0
	q = p - 1
	e += 1 while 0 == q[e]
	q >>= e

	# Find generator
	g = 2
	g += 1 until -1 == kronecker_symbol(g, p)
	z = power(g, q, p)	# |<z>| == 2 ** e

	# Initialize
	temp = power(n, q >> 1, p)
	x = n * temp % p	# n ** ((q + 1) / 2) mod p
	b = x * temp % p	# n ** q mod p

	# always
	# n * b == x ** 2
	until 1 == b
		# Find exponent f s.t. b ** (2 ** f) == 1 (mod p)
		f = 0
		b_ = b
		until 1 == b_
			b_ = b_ ** 2 % p
			f += 1
		end

		# Reduce exponent
		(e - f - 1).times { z = z ** 2 % p }
		e = f
		x = x * z % p
		z = z ** 2 % p
		b = b * z % p
	end

	return x
end

.modified_cornacchia(d, p) ⇒ Object

Param

d is a negative integer s.t. d = 0 or 1 mod 4 prime number p (|d| < 4p)

Return

integer solution (x, y) to the Diophantine equation x ** 2 + |d| * y ** 2 = 4p if exists else nil



645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
 
# File 'lib/include/fundamental.rb', line 645

def modified_cornacchia(d, p)
	# Case p == 2
	if 2 == p
		return nil unless q = (d + 8).square?
		return q, 1
	end

	# Test if residue
	return nil if -1 == kronecker_symbol(d, p)

	# Compute square root
	x0 = mod_sqrt(d, p)
	x0 = p - x0 if 0 != x0 - d % 2

	a = p << 1
	b = x0
	l = isqrt(p << 2)

	# Euclidean algorithm
	a, b = b, a % b while b < l

	# Test solution
	c, r = (4 * p - b**2).divmod(-d)
	return nil if 0 != r
	return nil unless q = c.square?
	return b, q
end

.mpqs(n, thread_num = Abst::THREAD_NUM) ⇒ Object 



479
480
481
482
 
# File 'lib/include/prime_mpqs.rb', line 479

def mpqs(n, thread_num = Abst::THREAD_NUM)
	mpqs = MPQS.new(n, thread_num)
	return mpqs.find_factor
end

.n_minus_1(n, factor = nil) ⇒ Object

Primality test

Param

positive integer n > 2 factor is factorization of n - 1

Return

boolean whether n is prime or not



106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
 
# File 'lib/include/prime.rb', line 106

def n_minus_1(n, factor = nil)
	factor = factorize(n - 1) unless factor
	factor.shift if 2 == factor[0][0]

	n_1 = n - 1
	half_n_1 = n_1 >> 1

	primes = primes_list.each
	find_base = proc do
		b = primes.next
		until (t = power(b, half_n_1, n)) == n_1
			return false unless t == 1
			b = primes.next
		end
		b
	end

	base = find_base.call
	factor.each do |prime, e|
		while power(base, half_n_1 / prime, n) == n_1
			base = find_base.call
		end
	end

	return true
end

.next_prime(n) ⇒ Object

Param

integer n

Return

The least prime greater than n



519
520
521
522
523
524
525
526
 
# File 'lib/include/prime.rb', line 519

def next_prime(n)
	return Bisect.find_gt(primes_list, n) if n < primes_list.last

	n += (n.even? ? 1 : 2)
	n += 2 until prime?(n)

	return n
end

.octagonal(n) ⇒ Object 



51
52
53
 
# File 'lib/include/sequence.rb', line 51

def octagonal(n)
	return n * (3 * n - 2)
end

.p_minus_1(n, bound = 10_000, m = 2) ⇒ Object

Param

positive integer n positive integer bound <= PRIME_CACHE_LIMIT positive integer m (2 <= m < n)

Return

a factor f (1 < f < n) if found else nil



301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
 
# File 'lib/include/prime.rb', line 301

def p_minus_1(n, bound = 10_000, m = 2)
	plist = primes_list

	p = nil
	old_m = m
	old_i = i = -1

	loop do
		i += 1
		p = plist[i]
		break if nil == p or bound < p

		# Compute power
		p_pow = p
		lim = bound / p
		p_pow *= p while p_pow <= lim
		m = power(m, p_pow, n)

		next unless 15 == i & 15

		# Compute GCD
		g = lehmer_gcd(m - 1, n)
		if 1 == g
			old_m = m
			old_i = i
			next
		end
		return g unless n == g
		break
	end

	if nil == p or bound < p
		return nil if 0 == i & 15

		g = lehmer_gcd(m - 1, n)
		return nil if 1 == g
		return g unless n == g
	end

	# Backtrack
	i = old_i
	m = old_m

	loop do
		i += 1
		p_pow = p = plist[i]

		loop do
			m = power(m, p, n)
			g = lehmer_gcd(m - 1, n)
			if 1 == g
				p_pow *= p
				break if bound < p_pow
				next
			end
			return nil if n == g
			return g unless 1 == g
		end
	end
end

.pentagonal(n) ⇒ Object

Pentagonal numbers are generated by the formula, P_n = n * (3 * n - 1) / 2. The first ten pentagonal numbers are:

1, 5, 12, 22, 35, 51, 70, 92, 117, 145, ...



36
37
38
 
# File 'lib/include/sequence.rb', line 36

def pentagonal(n)
	return n * (3 * n - 1) >> 1
end

.phi(n) ⇒ Object 



528
529
530
531
532
533
 
# File 'lib/include/prime.rb', line 528

def phi(n)
	return 1 if 1 == n
	return n - 1 if prime?(n)

	return factorize(n).inject(1) {|r, i| r * i[0] ** (i[1] - 1) * (i[0] - 1)}
end

.pollard_rho(n, c = 1, s = 2, max = 10_000) ⇒ Object

Param

positive integer n integer c which is used recurrence formula x ** 2 + c (mod n) integer s starting value of the recurrence formula integer max number of trials

Return

a factor f (1 < f < n) if found else nil



265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
 
# File 'lib/include/prime.rb', line 265

def pollard_rho(n, c = 1, s = 2, max = 10_000)
	u = s
	v = (s ** 2 + c) % n
	range = 1
	product = 1
	terms = 0

	loop do
		range.times do
			v = (v ** 2 + c) % n
			temp = product * (u - v) % n
			if 0 == temp
				g = lehmer_gcd(n, product)
				return (1 < g) ? g : nil
			end
			product = temp
			terms += 1

			if terms & 1023 == 0
				g = lehmer_gcd(n, product)
				return g if 1 < g
				return nil if max <= terms
			end
		end

		# reset
		u = v
		range <<= 1
		range.times { v = (v ** 2 + c) % n }
	end
end

.power_detection(n, largest_exp = true) ⇒ Object

Param

positive integer n boolean largest_exp

Return

integer x, k s.t. n == x ** k if n == 1 then x, k == 1, 1 (2 <= k if exist else x, k == n, 1) if largest_exp is true then return largest k



744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
 
# File 'lib/include/fundamental.rb', line 744

def power_detection(n, largest_exp = true)
	x, k = n, 1

	limit = ilog2(n)
	(2..limit).each_prime do |exp|
		break if limit < exp

		root, pow = iroot(n, exp, true)
		if pow == n
			return root, exp unless largest_exp

			n = x = root
			k *= exp
			limit = ilog2(n)
			redo
		end
	end

	return x, k
end

.powers_of_2 ⇒ Object

Return

array power of 2 s.t. [1, 2, 4, 8, 16, 32, …]



774
775
776
777
778
779
780
781
782
783
 
# File 'lib/include/fundamental.rb', line 774

def powers_of_2
	unless $powers_of_2
		$powers_of_2 = [1]
		((BASE_BYTE << 3) - 1).times do |i|
			$powers_of_2[i + 1] = $powers_of_2[i] << 1
		end
	end

	return $powers_of_2
end

.precompute_primes ⇒ Object

precompute primes by eratosthenes



9
10
11
12
13
14
15
16
 
# File 'lib/include/prime.rb', line 9

def precompute_primes
	primes = eratosthenes_sieve(PRIME_CACHE_LIMIT).to_a

	Dir::mkdir(DATA_DIR) unless FileTest.exist?(DATA_DIR)
	open(PRIMES_LIST, "w") {|io| io.write(primes.map(&:to_s).join("\n"))}

	return primes
end

.prime?(n) ⇒ Boolean

Param

positive integer n

Return

boolean whether n is prime or not

Returns:

  • (Boolean) 


142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
 
# File 'lib/include/prime.rb', line 142

def prime?(n)
	if n <= 3
		return false if n <= 1
		return true
	end

	if n <= primes_list.last
		return Bisect.index(primes_list, n) ? true : false
	end

	factor = trial_division(n, 257)[0]
	return factor[0][0] == n unless factor.empty?

	if n < 341_550_071_728_321
		[2, 3, 5, 7, 11, 13, 17].each do |i|
			return false unless strong_pseudoprime_test(n, i)
		end
		return true
	end

	return false unless miller_rabin(n)
	return n_minus_1(n)
end

.prime_power?(n) ⇒ Boolean

Param

positive integer n > 1

Return

p if n is of the form p^k with p prime else false

Returns:

  • (Boolean) 


718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
 
# File 'lib/include/fundamental.rb', line 718

def prime_power?(n)
	if n.even?
		return n.power_of?(2) ? 2 : false
	end

	p = n
	loop do
		rslt, witness = miller_rabin(p, 10, true)
		if rslt
			# Final test
			redo unless prime?(p)
			return n.power_of?(p) ? p : false
		end

		d = lehmer_gcd(power(witness, p, p) - witness, p)
		return false if 1 == d or d == p
		p = d
	end
end

.primes_list ⇒ Object 



23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
 
# File 'lib/include/prime.rb', line 23

def primes_list
	return $primes if $primes

	# precomputed?
	if FileTest.exist?(PRIMES_LIST)
		$primes = load_precomputed_primes
	else
		$primes = precompute_primes
	end

	def $primes.include?(n)
		return Bisect.index(self, n)
	end

	$primes.freeze

	return $primes
end

.primitive_root(p) ⇒ Object

Param

odd prime p

Return

primitive root modulo p



473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
 
# File 'lib/include/fundamental.rb', line 473

def primitive_root(p)
	p_m1 = p - 1
	check_order = factorize(p_m1)
	check_order.shift
	half_p_m1 = p_m1 >> 1
	check_order.map!{|f| half_p_m1 / f.first}

	a = 1
	loop do
		a += 1

		next unless -1 == kronecker_symbol(a, p)
		check = true
		check_order.each do |e|
			if p_m1 == power(a, e, p)
				check = false
				break
			end
		end

		return a if check
	end
end

.pseudoprime_test(n, base) ⇒ Object

Param

positive integer n positive integer base

Return

boolean whether n passes a pseudoprime test for the base or not When n and base are not relatively prime, this algorithm may judge a prime number n to be a composite number



51
52
53
 
# File 'lib/include/prime.rb', line 51

def pseudoprime_test(n, base)
	return power(base, n - 1, n) == 1
end

.pythagorean(max_c) ⇒ Object

primitive Pythagorean number

Param

integer max_c

Return

array of primitive Pythagorean numbers s.t. [[a, b, c], …] a**2 + b**2 == c**2 and a < b < c <= max_c



846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
 
# File 'lib/include/fundamental.rb', line 846

def pythagorean(max_c)
	return Enumerator.new(self, :pythagorean, max_c) unless block_given?
	return [] if max_c <= 4

	1.upto(isqrt(max_c - 1)) do |m|
		mm = m ** 2
		s = m.even? ? 1 : 2
		s.step(m, 2) do |n|
			next unless gcd(m, n) == 1

			nn = n ** 2
			c = mm + nn
			break if max_c < c

			a = mm - nn
			b = (m * n) << 1
			a, b = b, a if b < a

			yield [a, b, c]
		end
	end
end

.residue_class(ideal) ⇒ Object 



4
5
6
7
8
9
10
 
# File 'lib/include/residue.rb', line 4

def residue_class(ideal)
	if prime?(ideal.n)
		return residue_class_field(ideal)
	else
		return residue_class_ring(ideal)
	end
end

.residue_class_field(maximal_ideal) ⇒ Object 



20
21
22
23
24
25
26
 
# File 'lib/include/residue.rb', line 20

def residue_class_field(maximal_ideal)
	residue_class = Class.new(IntegerResidueField) do
		@mod = maximal_ideal.n
	end

	return residue_class
end

.residue_class_ring(ideal) ⇒ Object 



12
13
14
15
16
17
18
 
# File 'lib/include/residue.rb', line 12

def residue_class_ring(ideal)
	residue_class = Class.new(IntegerResidueRing) do
		@mod = ideal.n
	end

	return residue_class
end

.right_left_power(g, n, mod = nil) ⇒ Object

Right-Left Binary Power

Param

group element g integer n Euclidean domain element mod

Return

g ** n % mod



9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
 
# File 'lib/include/fundamental.rb', line 9

def right_left_power(g, n, mod = nil)
	rslt = g.class.one
	return rslt if 0 == n

	if n < 0
		n = -n
		g **= (-1)
	end

	g %= mod if mod

	loop do
		rslt *= g if n.odd?
		rslt %= mod if mod

		n >>= 1
		break if 0 == n

		g = g ** 2
		g %= mod if mod
	end

	return rslt
end

.root_mod_p ⇒ Object

Param
Return

Raises:

  • (NotImplementedError) 


134
135
136
 
# File 'lib/include/polynomial.rb', line 134

def root_mod_p
	raise NotImplementedError
end

.strong_pseudoprime_test(n, b, e = nil, k = nil) ⇒ Object

Param

positive odd integer n positive integer b integer e and k s.t. n = 2 ** e * k and k is odd.

Return

boolean whether n passes a strong pseudoprime test for the base b or not When n and b are not relatively prime, this algorithm may judge a prime number n to be a composite number



61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
 
# File 'lib/include/prime.rb', line 61

def strong_pseudoprime_test(n, b, e = nil, k = nil)
	n_minus_1 = n - 1

	unless e
		e = 0
		e += 1 while 0 == n_minus_1[e]
		k = n_minus_1 >> e
	end

	z = power(b, k, n)

	return true if 1 == z or n_minus_1 == z
	(e - 1).times do
		z = z ** 2 % n
		return true if z == n_minus_1
	end

	return false
end

.trial_division(n, limit = INFINITY) ⇒ Object

Param

positive integer n >= 2 positive integer limit

Return

factorization up to limit and remainder of n i.e.

[[a, b], [c, d], …], r

s.t.

n == a**b * c**d * … * r, (r have no factor less than or equal to limit ** 2)



186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
 
# File 'lib/include/prime.rb', line 186

def trial_division(n, limit = INFINITY)
	factor = []
	lim = [limit, isqrt(n)].min

	divide = lambda do |d|
		n /= d
		div_count = 1
		loop do
			q, r = n.divmod(d)
			break unless 0 == r

			n = q
			div_count += 1
		end

		factor.push([d, div_count])
		lim = [lim, isqrt(n)].min

		lim < d
	end

	(plist = primes_list).each do |d|
		break if lim < d
		break if 0 == n % d and divide.call(d)
	end

	if plist.last < lim
		d = plist.last - plist.last % 30 - 1
		while d <= lim
			# [2, 6, 4, 2, 4, 2, 4, 6] are difference of [1, 7, 11, 13, 17, 19, 23, 29]
			# which are prime bellow 30 == 2 * 3 * 5 except 2, 3, 5
			[2, 6, 4, 2, 4, 2, 4, 6].each do |diff|
				d += diff
				break if 0 == n % d and divide.call(d)
			end
		end
	end

	# n is prime?
	if 1 < n and isqrt(n) <= lim
		factor.push([n, 1])
		n = 1
	end
	return factor, n
end

.triangle(n) ⇒ Object

Triangle numbers are generated by the formula, T_n = n * (n + 1) / 2. The first ten triangle numbers are:

1, 3, 6, 10, 15, 21, 28, 36, 45, 55, ...



29
30
31
 
# File 'lib/include/sequence.rb', line 29

def triangle(n)
	return n * (n + 1) >> 1
end

Instance Method Details

#show_gif_image(path) ⇒ Object 



2
3
4
5
6
7
8
9
10
 
# File 'lib/include/graph.rb', line 2

def show_gif_image(path)
	require 'tk'
	TkButton.new do
		image TkPhotoImage.new(file: path)
		command lambda {TkRoot.destroy }
		pack
	end
	Tk.mainloop
end