Overflow/Underflow problem

[fanta2]

I was adding the sum of the square of numbers as shown below:

do i=1, num
myresult = x(i) * x(i)
end do

However some values of x(i) when squared they become underflow/overflow how can I scale these numbers to get the correct result? Thanks in advance for any help.

 

[gummibaer]

Check max and min limits of numbers in your environment.

in 32 bit environment these would be

integer: +/- 2,147,483,647
real: something of between 1.0E-37 and 1.0 E+39 both positive and negative.

so you should limit your x to the squareroot of these values.

Norbert

 

[xwb]

By default myresult is an integer unless declared otherwise

 

[fanta2]

The variable myresult is declared as double precision. The problem is when I check the numbers there are values in the order of 1.0E+/-22. When this numbers are squared they overflow/underflow....so am asking how I can scale them?

 

[gummibaer]

????

Trusting that you work on a 32 bit system double precision would allow results in the magnitude of 2E+/-300 which should be large enough to hold the result of any of your values.

Tip #1:
Make sure that all the values involved in your arithmetics are declared real*2 (or double precision). Maybe it is not myresult that is overflowing but some other variable on the way.

Tip #2:
Just do your own arithmetics: save mantissa and exponent in two different variables and do your operations for these two variables seperately, for instance:

program square_big_number

  real figure,mant,result_mant
  integer exp,result_exp
  figure = 1.234567e+25    ! just a big number

  mant = fraction(figure)  !standard intrinsic to get matissa of real number
  exp = exponent(figure)   !standard intrinsic to get exponent of real number

! evaluate figure * figure

  result_mant = mant * mant
  result_exp = exp + exp

! print result

  write(5,10) result_mant,result_exp
10 format(f10.7,'E',i4)

end

Multiplication of two figures:
multiply both mantissas and add both exponents

Division of two figures
divide the mantissas and substract the exponents

Addition and subtraction however is somewhat tricky as you have to observe the relative positions of the two mantissas given by the exponents.

Squaring is as above, computing squareroot you should make sure that your exponent is an even number (just dividing your mantissa by 10 and adding 1 to your exponent when it is an odd number).


New idea:
Depending of what you want to do with the result afterwards maybe it would be wise to use logarithms intead of the numbers itself.

Norbert

 

[xwb]

I was thinking of logs too but Nobert beat me to it. Note that ln/exp is faster than log/10**.

 

[ArkM]

Floating point data is scaled by definition. No need to invent another scaling mechanics.
Common source for under/overflow: uninitialized data (garbage). Check x array contents in your code.

 

[fanta2]

Thanks gummibaer and xwb as usual! I tried to identify where it overflows and underflows but couldn't solve the problem.

Without enabling the floating point exception (visual fortran) ... it calculates without error and gives me a wrong result. But when I enable the floating point exception it terminates in the middle it says <65> floating invalid.

In this line I have a division - (verysmall number)/(Very small number) - the ratio is around 1.0.

The denominator is sum of a product of two columns. Sum(x*y) --- x and y can vary from large number to very small numbers in the order of 1E-200. The product of these two small numbers become zero in the fortran when it underflows --- eventhough they are not zero.

I cannot transform them using ln or log as they have also negative values.

I tried using separating mantissa and exponent but I couldn't get it right as follows.

Real(8)::x(lastcol), y(lastcol)
Real(8):: mant1, mant2, mant_res, finalsum
Integer:: exp1, exp2, exp_res

finalsum = 0.0D0
do i=1, lastCol
mant1 = fraction(x(i))
exp1 = exponent(x(i))
mant2 = fraction(y(i))
exp2 = exponent(y(i))

mant_res = mant1*mant2
exp_res = exp1 + exp2

finalsum = finalsum + { How can I put the exp and mant together}

End do

I appreciate any kind of help!

 

[fanta2]

I use this function set_exponent(mant, exp) from the manual but still it doesn't give me the correct one. I think it is not possible to put the result in finalsum variable.

 

[ArkM]

If x(i) is real (32-bit) then x(i)*x(i) treated as 32-bit multiplication. Try to cast one (or both) operand to double precision before multiplication. You may try DBLE(x) intrinsic function of DPROD(x,y) intrinsic function (get the double precision value of multiplying two reals) or simply rewrite (if myresult is double precision variable):

do i=1, num
myresult = x(i)
myresult = myresult * myresult
end do

No matter what are left hand attributes when right hand expression calculated. So overflows/underflows and other troubles...

 

[xwb]

Note that the definition of DP is different from that of 16-bit minis. In those days, it really was double precision i.e. bigger mantissa.

Nowadays, it is both double the precision (bigger mantissa) and roughly 6 times the range of the exponent as well.

 

[gummibaer]

finalsum = finalsum + { How can I put the exp and mant together}

finalsum = mant_res*10.0**exp_res


How do you know about your results being valid or not ? From what you write down here it seems you have a very severe problem of what we called a 'poorly conditioned problem' (do not know if this is the proper English term). This means any operation of division and subtraction - or a combination thereof - where numbers of very different magnitudes are involved.

Look at the following for illustration:

10,002 - 10.000 = 2

If any of the values on the lefthandside change by only a very minute fraction your result could go astray completely.

10,003 - 10,000 = 3

A change of less than a percent of a percent on the left brought a change of 50 % on the right of this evaluation.

So the difference you encounter checking on the results of your PC's evaluation might come fromt the different number of significant digits that you use in the calculation.

Just an idea. But if this is the nature of your problem, then you should contact some folks with in depth knowledge of computer maths (which unfortunately I do not have)

Norbert