c - Calculate multiplicative inverse of positive real number -
i wanted calculate multiplicative inverse of positive real number without using division. after reading through various pages on newton raphson method, implemented following code calculating inverse:
#define precision 0.1e-5 double inverse(double num) { double ans = num; while (ans * num > 1 + precision || ans * num < 1 - precision) { ans = ans * (2 - num * ans); } return ans; }
now, problem is, not inverse value. loop goes on infinitely.
so, 1st thing noticed was: value of ans
become negative, added statement if (ans < 0) ans *= -1;
so ans
remains positive always.
2nd point noticed: if initialization ans = 0.5
correct answer few values of num
= (1, 2, 3, 5)
.
so, assumption is, implementation isn't working because of inappropriate initialization of variable ans
.
so, questions:
1.can method used calculate inverse of positive real numbers?
2.if yes, there conditions required initial value assumed when using newton raphson method?
3.is there other better method calculate reciprocal without using division?
edit:
thanks answer. so, mentioned in answer, changed initial value of ans
precision
, works! also, initial value particular max limit on num
, ans
never becomes negative, no need negative checking condition had added initially.
so, working code (atleast works input have given.)
double inverse(double num) { // added make valid inputs. // large number under precision constraint, answer 0. if (num > 999999) return 0; double ans = precision; while (ans * num > 1 + precision || ans * num < 1 - precision) { ans = ans * (2 - num * ans); } return ans; }
you should pick initial approximation (0, 2/num). if pick other side of 2/num, method may not converge: approximation overshoot 0 , sequence tend -∞.
to arrive @ interval, let's see ans*(2-num*ans)
changes sign: becomes negative when 2-num*ans < 0, or when ans > 2/num. ans
should less 2/num.
to able pick initial approximation have know little how floating point numbers expressed. typically computers use x = s*m*2e, s sign, m ("mantissa") in (0.5, 1) , e ("exponent") integer. 1/x = s*1/m * 2-e, problem reduced calculating inverse of number in range (0.5, 1), , in range can use example 1 initial guess. apparently optimal initial guess in range 48/17 - 32/17*m.
one initial guess should work numbers s*m*2e s*2-e. in c can calculated as:
double ans = num; // initial guess: ans = 2**(-floor(log num)) long int *t = (long int*)(&ans); *t = (*t & (0x1l<<63)) | (~*t & (0x3ffl << 52)) - (0x1l << 52); /* sign */ /* negative exponent */
(beware: haven't checked edge conditions.)
Comments
Post a Comment