Реализация функции вычисления корня

Question

Реализация математических функций для различных вещей достаточно проста. int mul(int,int);, int pow(int,int);, даже double div(float,float); легко сделать, и их можно реализовать с помощью циклов или рекурсии. (Это те же методы, которые используются для выполнения этих функций вручную или в голове.) Чтобы умножить, просто несколько раз добавьте число. Чтобы разделить, неоднократно вычитайте это. Чтобы получить власть, многократно умножьте. И так далее.

Одна математическая функция, о которой я всегда задумывался, это корни. Например, как бы вы написали функцию для вычисления квадратного (или куба и т. Д.) Корня числа (т. Е. double root(float num, float root);)? Я попытался осмотреться и не смог найти алгоритм или метод для этого.

Когда я пытаюсь вычислить корень вручную, я обычно использую метод предположения (начните с приблизительного числа, добавьте дробь, умножьте, посмотрите, как далеко это расстояние, добавьте меньшую дробь, умножьте, проверьте еще раз и повторить пока не выполнится). Я полагаю, что это может сработать, но, безусловно, есть лучший и более быстрый метод (независимо от того, насколько быстрее компьютер может сделать это, чем вручную).

Очевидно, что LUT не являются релевантными, поскольку они должны быть достаточно универсальными, чтобы брать любые операнды (если вы не пишете игру с конечным набором данных). В статье Википедии упоминается метод угадывания и перечисляются некоторые древние (задолго до того, как были изобретены компьютеры), а также некоторые методы чистой математики и даже исчисления (в том числе те, которые имеют «бесконечность» в качестве компонента). Единственные, кто, кажется, имеет какое-либо отношение к электронике, используют трюки или логирифмы. (И это только для квадратных корней , не говоря уже о кубических корнях и подобных.)

Нет ли простого способа вычисления корня? Как калькуляторы это делают? Как компьютеры делают это? (Нет, просто выполнение double pow(a,0.5); не сработает, потому что тогда как будет реализовано double pow(float,float)?)

Я просто неправильно группирую корневые функции с более простыми функциями? Они сложнее, чем кажутся?

Answer 1

Есть несколько возможностей. Есть несколько разных итеративных подходов, таких как бисекция или ньютоновский подход. Что касается использования pow, то на некоторых компьютерах (например, x86) есть инструкция, чтобы сделать (по крайней мере, часть) возведение числа в степень, так что это просто вопрос написания некоторой структуры вокруг этого.

Вот реализация на языке ассемблера метода Ньютона для квадратного корня, в данном случае работающая только с 16-разрядными целыми числами, но та же основная идея применима и к другим типам. Я написал это около 20 лет назад, так что это было для 16-битных процессоров без аппаратного обеспечения с плавающей запятой.

isqrt proc uses di, number:word
;
; uses bx, cx, dx
;
    mov di,number
    mov ax,255
start_loop:
    mov bx,ax
    xor dx,dx
    mov ax,di
    div bx
    add ax,bx
    shr ax,1
    mov cx,ax
    sub cx,bx
    cmp cx,2
    ja  start_loop
    ret
isqrt endp

Вот некоторый код для x87 для вычисления произвольных степеней:

pow proc
    ; x^y = 2^(log2(x) * y)
    fyl2x    
    fld st(0)
    frndint  
    fld1     
    fscale   
    fxch st(2)
    fsubrp   
    f2xm1    
    fld1     
    faddp    
    fmulp    
    ret
endp

Обратите внимание, однако, что вы обычно не хотите реализовывать умножение путем простого многократного сложения или деления путем повторного вычитания. Скорее, вы хотите сдвинуть и сложить / вычесть последовательные степени двух, чтобы получить результат намного быстрее.

Вот некоторый код, который показывает общую идею:

mult proc
; multiplies eax by ebx and places result in edx:ecx
    xor ecx, ecx
    xor edx, edx
mul1:
    test ebx, 1
    jz  mul2
    add ecx, eax
    adc edx, 0
mul2:
    shr ebx, 1
    shl eax, 1
    test ebx, ebx
    jnz  mul1
done:
    ret
mult endp

Это довольно бессмысленно для x86, потому что в него встроены инструкции умножения, но на старых процессорах (PDP-11, 8080, 6502 и т. Д.) Такой код был вполне распространенным.

Answer 2

Это зависит от того, насколько вы хотите быть в целом. Если вы хотите вычислить, например, (-4.2) ^0.23 , вам понадобится сложная арифметика. Как указывает Мэт, существуют быстрые алгоритмы для вычисления x ^{1 / n} для целого числа n и положительного значения x. Если вы хотите x ^y для положительного x и любого y, то журналы и экспоненты будут работать.

Answer 3

Существует простой для реализации N-й корневой алгоритм , связанный непосредственно со статьей, которую вы заявляете.Он получен из метода Ньютона .