Module:Rational: Difference between revisions

← Older edit Newer edit →

@@ Line 658: / Line 658: @@
 	end
 	return max_factor
+end
+-- FJS: x = a * 2^n : x >= 1, x < 2
+local function red(a)
+	if type(a) == 'number' then
+		a = p.new(a)
+	end
+	local b = p.copy(a)
+	while p.lt(b, 1) do
+		b = p.mul(b, 2)
+	end
+	while p.geq(b, 2) do
+		b = p.div(b, 2)
+	end
+	return b
+end
+-- FJS: x = a * 2^n : x >= 1/sqrt(2), x < sqrt(2)
+local function reb(a)
+	local b = red(a)
+	if p.geq(p.mul(b, b), 2) then
+		b = p.div(b, 2)
+	end
+	return b
+end
+-- FJS: master algorithm
+local function FJS_master(prime)
+	prime = red(prime)
+	local tolerance = p.new(65, 63)
+	local fifth = p.new(3, 2)
+	local k = 0
+	while true do
+		local a = red(p.pow(fifth, k))
+		if math.abs(p.cents(p.div(prime, a))) < p.cents(tolerance) then
+			return k
+		end
+		if k == 0 then
+			k = 1
+		elseif k > 0 then
+			k = -k
+		else
+			k = -k + 1
+		end
+	end
+end
+-- FJS: formal comma
+local function formal_comma(prime)
+	local fifth_shift = FJS_master(prime)
+	return reb(p.div(prime, p.pow(3, fifth_shift)))
 end