;============================================================
;  アルゴリズム 26〜30: 数論
;    26. Euclidean GCD + LCM
;    27. Extended Euclidean (逆元)
;    28. Modular Exponentiation
;    29. Sieve of Eratosthenes
;    30. Euler's Totient φ(n)
;============================================================

#include "hsp3cl_net_64.as"

; ----- 26. GCD / LCM -----
#module
#defcfunc gcd int _a, int _b, local _x, local _y, local _t
	_x = _a : _y = _b
	repeat
		if _y == 0 : break
		_t = _y : _y = _x \ _y : _x = _t
	loop
	return _x

#defcfunc lcm int _a, int _b
	return _a / gcd(_a, _b) * _b
#global

; ----- 27. Extended Euclidean -----
#module
#deffunc ext_gcd int _a, int _b, var _x, var _y, local _x1, local _y1
	if _b == 0 {
		_x = 1 : _y = 0 : return _a
	}
	ext_gcd _b, _a \ _b, _x1, _y1
	_x = _y1
	_y = _x1 - (_a / _b) * _y1
	return stat

#defcfunc mod_inv int _a, int _m, local _x, local _y, local _g
	ext_gcd _a, _m, _x, _y
	_g = stat
	if _g != 1 : return -1
	return (_x \ _m + _m) \ _m
#global

; ----- 28. Modular Exponentiation (a^b mod m) -----
#module
#defcfunc mod_pow int _a, int _b, int _m, local _r, local _x, local _y
	_r = 1 : _x = _a \ _m : _y = _b
	repeat
		if _y <= 0 : break
		if _y & 1 : _r = _r * _x \ _m
		_x = _x * _x \ _m
		_y >>= 1
	loop
	return _r
#global

; ----- 29. Sieve of Eratosthenes -----
#module
#deffunc sieve int _n, array _out, local _i, local _j
	dim _out, _n + 1
	repeat _n - 1
		_out(cnt + 2) = 1
	loop
	repeat _n - 1
		_i = cnt + 2
		if _out(_i) == 0 : continue
		_j = _i * _i
		repeat
			if _j > _n : break
			_out(_j) = 0
			_j += _i
		loop
	loop
	return
#global

; ----- 30. Euler's Totient φ(n) -----
#module
#defcfunc euler_phi int _n, local _nv, local _result, local _p
	_nv = _n : _result = _n
	_p = 2
	repeat
		if _p * _p > _nv : break
		if _nv \ _p == 0 {
			repeat
				if _nv \ _p != 0 : break
				_nv = _nv / _p
			loop
			_result -= _result / _p
		}
		_p++
	loop
	if _nv > 1 : _result -= _result / _nv
	return _result
#global

; ===== デモ =====
	mes "=== GCD / LCM ==="
	mes strf("gcd(48, 18) = %d", gcd(48, 18))
	mes strf("lcm(12, 18) = %d", lcm(12, 18))

	mes ""
	mes "=== Extended Euclidean ==="
	ext_gcd 30, 18, xx, yy
	mes strf("30*x + 18*y = gcd: x=%d y=%d gcd=%d", xx, yy, stat)
	mes strf("mod_inv(3, 11) = %d (期待: 4)", mod_inv(3, 11))

	mes ""
	mes "=== Modular Exponentiation ==="
	mes strf("mod_pow(2, 10, 1000) = %d", mod_pow(2, 10, 1000))
	mes strf("mod_pow(3, 200, 100) = %d", mod_pow(3, 200, 100))

	mes ""
	mes "=== Sieve 50 ==="
	sieve 50, primes
	sdim _line, 256
	_line = "primes: "
	repeat 51
		if primes(cnt) : _line += "" + cnt + " "
	loop
	mes _line

	mes ""
	mes "=== Euler's phi ==="
	mes strf("phi(1) = %d",  euler_phi(1))
	mes strf("phi(9) = %d (期待: 6)",  euler_phi(9))
	mes strf("phi(10) = %d (期待: 4)", euler_phi(10))
	mes strf("phi(36) = %d (期待: 12)", euler_phi(36))
	end 0