MIPS How can I calculate numbers greater than 1.111 x2^1023?

Question

I'm working with MIPS and MARS 4.5. Is there any way to calculate numbers that extend the range of 1.111.. * 2^1023?

For example (factorial func):

.data
     dval:    .double   171.0  # its working up to 170!
     one:     .double   1.0  

.text
     l.d      $f2, dval
     l.d      $f4, one
     l.d      $f12, one
lo:  c.eq.d   $f2, $f4          # calc factorial of 171 = Infinity 
     bc1t     ex
     mul.d    $f12, $f12, $f2
     sub.d    $f2, $f2, $f4
     j        lo
ex:  li       $v0, 3
     syscall

How can I calculate and print the factorial of 171?

Implement arbitrary precision integer arithmetic for example. — Jester, Jan 11 '19 at 18:24
"is there a way?" yes, you can calculate numbers almost up to "your size of memory in bits"-th power of two (minus memory for code and you may need multiple values of large precision so then the spare memory divided by amount of arbitrary-precision numbers), or even more if you swap to disk. "how?" - quite some extra code, as native integer types are 32b max, and native float types are 64b max (53b for digits), so the larger amount of bits requires all the arithmetic to be done manually (the amount of bits is decisive in "how much information can I encode"). — Ped7g, Jan 12 '19 at 06:15
And as noted below in answer, your code does use only about 53 bits to store precise value, so you are already losing digits (at right end) much sooner, at 171! you just run also out of exponent, but the value is probably skewed a lot already by the missing bottom part. (so your comment "working up to 170!" is "working" as in "does not trap/crash", not as "working correctly") — Ped7g, Jan 12 '19 at 06:17

score 6 · Accepted Answer · answered Jan 11 '19 at 19:51

It is surprising to compute a factorial in FP. Precision of FP number is limited by their mantissa and actually above ~20!, you need more that 54 bits to store the (integer) value and the FP result is just an approximation and hence is incorrect.

But, if you only need an imprecise value in double, I would suggest:

1/ using a simpler formula, such as Stirling formula. https://en.wikipedia.org/wiki/Stirling%27s_approximation#Speed_of_convergence_and_error_estimates If you keep 4 or 5 numbers in the Stirling serie, the loss of approximation due to using this formula, will not be larger than the one due to the limited mantissa. (and it will be much faster).

2/ whether you use the factorial or Stirling, you can use a trick to extend the exponent in the double. The idea just to keep the exponent below 2^512 and after every mult, you do the following check

 int extra_exponent=0;
 double theshold = 2^512;
 double threshold_inv=2^-512;
 ....
 # check exponent on n
 if (n>threshold) {
   n *= threshold_inv ; # reduce exponent
   extra_exponent++ ; # but keep track of it
 }

You will never overflow a long as you multiply by a number < 2^512. At the end, to get the "real" value you have to multiply but 2^(extra*512) (which may not be so simple). And your precision is limited by the the mantissa size.

But the best proposal is to use infinite precision arithmetic. And this is the only way to get an exact result.

MIPS How can I calculate numbers greater than 1.111 x2^1023?

1 Answers1