Это сложно, поэтому я собираюсь изложить подход, который я бы выбрал.Я не обещаю никаких ошибок, и вам останется много работы.
Я изменю переменные с того, что вы сказали, на exp(x, y, err)
на x^y
в пределах ошибки err
.If y
не находится в диапазоне 0 <= y < 1
, тогда мы можем легко умножить на соответствующее x^k
с k
целое число, чтобы сделать это так.Поэтому нам нужно беспокоиться только о дробном `y
. Если бы все числители и знаменатели были маленькими, было бы легко решить эту проблему, сначала взяв целочисленную степень, а затем получив корень, используя метод Ньютона.Но эта наивная идея развалится до боли, когда вы попытаетесь оценить что-то вроде (1000001/1000000)^(2000001/1000000)
.Поэтому задача состоит в том, чтобы не дать вам взорваться.
Я бы порекомендовал рассмотреть проблему вычисления x^y
как x^y = (x0^y0) * (x0^(y-y0)) * (x/x0)^y = (x0^y0) * e^((y-y0) * log(x0)) * e^(y * log(x/x0))
.И мы выберем x0
и y0
так, чтобы вычисления были проще и ошибки были ограничены.
Чтобы ограничить ошибки, мы сначала можем придумать наивную верхнюю границу b
для x0^y0
- что-то вроде «следующего наивысшего целого, чем x
, до степени следующего наивысшего целого, чем y
».Мы выберем x0
и y0
, чтобы быть достаточно близкими к x
и y
, чтобы последние термины были под 2
.И тогда нам просто нужно оценить три члена с точностью до err/12
, err/(6*b)
и err/(6*b)
.(Возможно, вы захотите сделать эти ошибки наполовину более жесткими, чтобы сделать окончательный ответ близким к рациональному.)
Теперь, когда мы выберем x0
и y0
, мы будем стремиться к "близкому рациональному с небольшим числителем /знаменатель".Для этого мы начнем вычислять продолжение доли .Это дает последовательность рациональных чисел, которая быстро сходится к целевому реальному.Если мы довольно скоро обрежем последовательность, мы сможем быстро найти рациональное число, которое находится на любом желаемом расстоянии от целевого действительного числа, сохранив при этом относительно небольшие числители и знаменатели.
Давайте работать от третьего слагаемого назад.
Мы хотим y * log(x/x0) < log(2)
.Но из серии Тейлора, если x/2 < x0 < 2x
, то log(x/x0) < x/x0 - 1
.Таким образом, мы можем найти в непрерывной дроби соответствующий x0
.
Как только мы его найдем, мы можем использовать ряд Тейлора для log(1+z)
, чтобы вычислить log(x/x0)
с точностью до err/(12*y*b)
.А затем ряд Тейлора для e^z
, чтобы вычислить член для нашей желаемой ошибки.
Второй член является более сложным.Нам нужно оценить log(x0)
.Мы находим подходящее целое число k
, такое что 1.1^k <= x0 < 1.1^(k+1)
.И тогда мы можем достаточно точно оценить и k * log(1.1)
, и log(x0 / 1.1^k)
.Найдите наивную верхнюю границу этого log
и используйте ее, чтобы найти достаточно близкое значение y0
, чтобы второе слагаемое находилось в пределах 2. А затем используйте ряд Тейлора, чтобы оценить e^((y-y0) * log(x0))
до желаемой точности.
Для первого члена мы используем наивный метод возведения x0
в целое число, а затем метод Ньютона, чтобы получить корень, чтобы дать x0^y0
с нашей желаемой точностью.
Затем умножить их вместе, иу нас есть ответ.(Если вы выбрали «более жесткие ошибки, более хороший ответ», то теперь вы будете продолжать дробную часть этого ответа, чтобы выбрать более рациональное возвращение.)