Is it possible to force exponent or significand of a float to match another float (Python)?

Question

This is an interesting question that I was trying to work through the other day. Is it possible to force the significand or exponent of one float to be the same as another float in Python?

The question arises because I was trying to rescale some data so that the min and max match another data set. However, my rescaled data was slightly off (after about 6 decimal places) and it was enough to cause problems down the line.

To give an idea, I have f1 and f2 (type(f1) == type(f2) == numpy.ndarray). I want np.max(f1) == np.max(f2) and np.min(f1) == np.min(f2). To achieve this, I do:

import numpy as np

f2 = (f2-np.min(f2))/(np.max(f2)-np.min(f2)) # f2 is now between 0.0 and 1.0
f2 = f2*(np.max(f1)-np.min(f1)) + np.min(f1)  # f2 is now between min(f1) and max(f1)

The result (just as an example) would be:

np.max(f1) # 5.0230593
np.max(f2) # 5.0230602 but I need 5.0230593 

My initial thought is that forcing the exponent of the float would be the correct solution. I couldn't find much on it, so I made a workaround for my need:

exp = 0
mm = np.max(f1)

# find where the decimal is
while int(10**exp*mm) == 0
  exp += 1

# add 4 digits of precision
exp += 4

scale = 10**exp

f2 = np.round(f2*scale)/scale
f1 = np.round(f1*scale)/scale

now np.max(f2) == np.max(f1)

However, is there a better way? Did I do something wrong? Is it possible to reshape a float to be similar to another float (exponent or other means)?

EDIT: as was suggested, I am now using:

scale = 10**(-np.floor(np.log10(np.max(f1))) + 4)

While my solution above will work (for my application), I'm interested to know if there's a solution that can somehow force the float to have the same exponent and/or significand so that the numbers will become identical.

Monophagous

Depending on the range you're working with, you might be able to finagle something using np.log10(x). With some dummy values I get ` np.log10(1.2312412) = 0.090343139501527295` and np.log10(12.312412) = 1.0903431395015273, but np.log10(.12312412) = -0.90965686049847261 – Corr

Or, so long as you're only scaling by powers of 10 10**(np.log10(1.2312412)) = 1.2312411999999999, 10**(np.log10(1.2312412)+2) = 123.12412000000005, and 10**(np.log10(1.2312412)-2) = 0.012312411999999998, meaning the error in the mantissa is pushed to the last digit. – Corr

@CeramicSheep Using log10 would have probably been a smarter idea than incrementing exp. That's a good point. Could modify it to scale = 10**(np.floor(np.log10(mm))+4). Definitely another workaround. Still curious if I could make np.max(f1) == np.max(f2) without doing any rounding and losing precision. – Monophagous

I don't understand how you're using the word "mantissa" in this context. Can you clarify? Do you mean that you want f1.max()/f2.max() to be a power of 10 (which is what the normal usage of mantissa would suggest). – Goodard

@MarkDickinson Floating point Wiki. I think I meant the exponent. Mantissa would be the significand (I thought it was the exponent). Either way, I'd be interested in know if you can "force" either part. – Monophagous

Ah, okay. FWIW, the floating-point use of "mantissa" tends to be frowned upon, since it conflicts with the older mathematical definition of mantissa as the fractional part of the (base 10) logarithm. Use "significand" instead. – Goodard

@MarkDickinson to further clarify, I want f1.max() == f2.max() without having to round to some number of significant digits. Mathematically speaking, they should be identical due to the equation. However, floating numbers gunk that up. – Monophagous

@MarkDickinson thanks. Did not know that! – Monophagous

re f2 *= (np.max(f1)-np.min(f1)) + np.min(f1). Should this not be f2 = f2 * (np.max(f1)-np.min(f1)) + np.min(f1). I.e. you want to rescale by the range of f1, and then add the minimum of f1. It looks to me that as it stands you are multiplying your 0 to 1 range by the whole of the term on the right of the *=. – Bondholder

@Bondholder good catch! I've edited my question. – Monophagous

@MarkDickinson I've edited the question to remove "mantissa" references. I hope that clears it up. – Monophagous

re the edit to fix the *=, does the # 5.0230602 but I need 5.0230593 still hold? – Bondholder

@TooTone, yes. I am actually doing the math in one line in my code. I made it 2 lines for SO and messed up. However, keep in mind the numbers I used my example are not the numbers I was working with. I am just using them as an explanation. – Monophagous

the relative error you are seeing is −0.000000179. That is not consistent with double precision arithmetic. Maybe there is some error in your code, or the data has some peculiar properties. Please can you post code with a sample of actual data run through that code that shows the problem? You don't need to post the whole data set, for obvious reasons.. Just the highest and lowest number in each data set, plus one or two random numbers inbetween. Otherwise it is very difficult to see what is going on. In general terms I also have some issues with using rounding, but need to see data as well. – Bondholder

@Bondholder I can. It may take me a little time to parse through and make a nice clean example with real numbers. I probably won't get to it today. Stay tuned. – Monophagous

@TooTone: It could well be single-precision arithmetic (i.e., NumPy arrays of dtype float32). – Goodard

@Mark thanks, I've been sadly out of the Python numpy world for a while. Would still be good to see the data set but point taken. – Bondholder

@Bondholder sadly (and sort of embarrassingly) I cannot replicate my original problem. I don't think the data has changed... I must have been doing something wrong. I will leave my question posed as a fundamental question. – Monophagous

Answer 1

TL;DR

Use

f2 = f2*np.max(f1)-np.min(f1)*(f2-1)  # f2 is now between min(f1) and max(f1)

and make sure you're using double precision, compare floating point numbers by looking at absolute or relative differences, avoid rounding for adjusting (or comparing) floating point numbers, and don't set the underlying components of floating point numbers manually.

Details

This isn't a very easy error to reproduce, as you have discovered. However, working with floating numbers is subject to error. E.g., adding together 1 000 000 000 + 0 . 000 000 000 1 gives 1 000 000 000 . 000 000 000 1, but this is too many significant figures even for double precision (which supports around 15 significant figures), so the trailing decimal is dropped. Moreover, some "short" numbers can't be represented exactly, as noted in @Kevin's answer. See, e.g., here, for more. (Search for something like "floating point truncation roundoff error" for even more.)

Here's an example which does demonstrate a problem:

import numpy as np

numpy.set_printoptions(precision=16)

dtype=np.float32                     
f1 = np.linspace(-1000, 0.001, 3, dtype=dtype)
f2 = np.linspace(0, 1, 3, dtype=dtype)

f2 = (f2-np.min(f2))/(np.max(f2)-np.min(f2)) # f2 is now between 0.0 and 1.0
f2 = f2*(np.max(f1)-np.min(f1)) + np.min(f1)  # f2 is now between min(f1) and max(f1)

print (f1)
print (f2)

output

[ -1.0000000000000000e+03  -4.9999951171875000e+02   1.0000000474974513e-03]
[ -1.0000000000000000e+03  -4.9999951171875000e+02   9.7656250000000000e-04]

Following @Mark Dickinson's comment, I have used 32 bit floating point. This is consistent with the error you reported, a relative error of around 10^-7, around the 7th significant figure

In: (5.0230602 - 5.0230593) / 5.0230593
Out: 1.791736760621852e-07

Going to dtype=np.float64 makes things better but it still isn't perfect. The program above then gives

[ -1.0000000000000000e+03  -4.9999950000000001e+02   1.0000000000000000e-03]
[ -1.0000000000000000e+03  -4.9999950000000001e+02   9.9999999997635314e-04]

This isn't perfect, but is generally close enough. When comparing floating point numbers you almost never want to use strict equality because of the possibility of small errors as noted above. Instead subtract one number from the other and check the absolute difference is less than some tolerance, and/or look at the relative error. See, e.g., numpy.isclose.

Returning to your problem, it seems like it should be possible to do better. After all, f2 has the range 0 to 1, so you should be able to replicate the maximum in f1. The problem comes in the line

f2 = f2*(np.max(f1)-np.min(f1)) + np.min(f1)  # f2 is now between min(f1) and max(f1)

because when an element of f2 is 1 you're doing a lot more to it than just multiplying 1 by the max of f1, leading to the possibility of floating point arithmetic errors occurring. Notice that you can multiply out the brackets f2*(np.max(f1)-np.min(f1)) to f2*np.max(f1) - f2*np.min(f1), and then factorize the resulting - f2*np.min(f1) + np.min(f1) to np.min(f1)*(f2-1) giving

f2 = f2*np.max(f1)-np.min(f1)*(f2-1)  # f2 is now between min(f1) and max(f1)

So when an element of f2 is 1, we have 1*np.max(f1) - np.min(f1)*0. Conversely when an element of f2 is 0, we have 0*np.max(f1) - np.min(f1)*1. The numbers 1 and 0 can be exactly represented so there should be no errors.

The modified program outputs

[ -1.0000000000000000e+03  -4.9999950000000001e+02   1.0000000000000000e-03]
[ -1.0000000000000000e+03  -4.9999950000000001e+02   1.0000000000000000e-03]

i.e. as desired.

Nevertheless I would still strongly recommend only using inexact floating point comparison (with tight bounds if you need) unless you have a very good reason not to do so. There are all sorts of subtle errors that can occur in floating point arithmetic and the easiest way to avoid them is never to use exact comparison.

An alternative approach to that given above, that might be preferable, would be to rescale both arrays to between 0 and 1. This might be the most suitable form to use within the program. (And both arrays could be multiplied by a scaling factor such the original range of f1, if necessary.)

Re using rounding to solve your problem, I would not recommend this. The problem with rounding -- apart from the fact that it unnecessary reduces the accuracy of your data -- is that numbers that are very close can round in different directions. E.g.

f1 = np.array([1.000049])
f2 = np.array([1.000051])
print (f1)
print (f2)
scale = 10**(-np.floor(np.log10(np.max(f1))) + 4)
f2 = np.round(f2*scale)/scale
f1 = np.round(f1*scale)/scale
print (f1)
print (f2)

Output

[ 1.000049]
[ 1.000051]
[ 1.]
[ 1.0001]

This is related to the fact that although it's common to discuss numbers matching to so many significant figures, people don't actually compare them this way in the computer. You calculate the difference and then divide by the correct number (for a relative error).

Re mantissas and exponents, see math.frexp and math.ldexp, documented here. I would not recommend setting these yourself however (consider two numbers that are very close but have different exponents, for example -- do you really want to set the mantissa). Much better to just directly set the maximum of f2 explicitly to the maximum of f1, if you want to ensure the numbers are exactly the same (and similarly for the minimum).

Bondholder

To add context, this data is used for interpolation (of something through time). I don't necessarily use equality operations, but when I try to map an f1 time position to an f2 time position, wonky things start to happen. Then things like np.unique return more than 2 values when time is infinite. Does that make sense? This is a good discussion. – Monophagous

@Jason. That context helps. Do you mean that f1[i] is something like the time that particle 1 reaches position i, or that f1[t] is the distance that particular 1 has reached at time t? (Similarly for particle 2.) – Bondholder

The two data pieces can be thought of like an acceleration or de-acceleration curves. However, I will jump between the curves at different points in time. As time approached infinity, the extrema should always be one of two values. In my original discussion, that wouldn't happen. Jumping curves was bad because one extrema didn't match the other. – Monophagous

Answer 2

It depends what you mean by "mantissa."

Internally, floats are stored using scientific notation in base 2. So if you mean the base 2 mantissa, it is actually very easy: Just multiply or divide by powers of two (not powers of 10), and the mantissa will stay the same (provided the exponent doesn't get out of range; if it does, you'll get clamped to infinity or zero, or possibly go into denormal numbers depending on architectural details). It's important to understand that the decimal expansions will not match up when you rescale on powers of two. It's the binary expansion that's preserved with this method.

But if you mean the base 10 mantissa, no, it's not possible with floats, because the rescaled value may not be exactly representable. For example, 1.1 cannot be represented exactly in base 2 (with a finite number of digits) in much the same way as 1/3 cannot be represented in base 10 (with a finite number of digits). So rescaling 11 down by 1/10 cannot be done perfectly accurately:

>>> print("%1.29f" % (11 * 0.1))
1.10000000000000008881784197001

You can, however, do the latter with decimals. Decimals work in base 10, and will behave as expected in terms of base 10 rescaling. They also provide a fairly large amount of specialized functionality to detect and handle various kinds of loss of precision. But decimals don't benefit from NumPy speedups, so if you have a very large volume of data to work with, they may not be efficient enough for your use case. Since NumPy depends on hardware support for floating point, and most (all?) modern architectures provide no hardware support for base 10, this is not easily remedied.

Gipon

IBM System z machines (from the z10 onwards) have hardware support for decimal arithmetic. But there's no facility that I'm aware of for either NumPy or Python to make use of that hardware support. – Goodard

Answer 3

Try replacing the second line by

f2 = f2*np.max(f1) + (1.0-f2)*np.min(f1)

Explanation: There are 2 places where the difference could creep in:

Step 1) f2 = (f2-np.min(f2))/(np.max(f2)-np.min(f2))

When you inspect np.min(f2) and np.max(f2), do you get exactly 0 and 1 or something like 1.0000003?

Step 2) f2 = f2*(np.max(f1)-np.min(f1)) + np.min(f1)

The expression like (a-b)+b doesn't always produce exactly a, due to rounding error. The suggested expression is slightly more stable.

For a very detailed explanation, please see What Every Computer Scientist Should Know About Floating-Point Arithmetic by David Goldberg.

Subtract

note that + (1.0-f2)*np.min(f1) is equivalent to -np.min(f1)*(f2-1) in my answer – Bondholder

Answer 4

def rescale(val, in_min, in_max, out_min, out_max):
    return out_min + (val - in_min) * ((out_max - out_min) / (in_max - in_min))

value_to_rescale = 5
current_scale_min = 0
current_scale_max = 10
target_scale_min = 100
target_scale_max = 200

new_value = rescale(value_to_rescale, current_scale_min, current_scale_max, target_scale_min, target_scale_max)
print(new_value)

new_value = rescale(10, 0, 10, 0, 100)
print(new_value)

answer:

150 100

Joellenjoelly

Answer 5

Here is one with decimals

from decimal import Decimal, ROUND_05UP
num1 = Decimal('{:.5f}'.format(5.0230593))  ## Decimal('5.02306')
num2 = Decimal('{}'.format(5.0230602))  ## Decimal('5.0230602')
print num2.quantize(num1, rounding=ROUND_05UP) ## 5.02306

EDIT** i am slightly confused of why I get so much negative feedback, so here is yet another solution not using decimals:

a = 5.0230593
b = 5.0230602
if abs(a - b) < 1e-6:
    b = a

Columbite

I think you are are getting negative feedback because this a workaround and inexact. In my case, this would cause discontinuities in the data. If a[1] - b[1] < 1e-6 and a[2] - b[2] > 1e-6 there could be a discontinuity of up to 1e-6 from b[1] to b[2]. Also Decimal would be incredibly slow in my case. – Monophagous

well maybe you haven't understood the problem but we have to make discontinuities because we want the final floats to be identical. the question has not defined how the rounding will be dealt with, just that it would round upto 5 or 6 decimals. you can certainly make a chained if statement that does 5 decimals for some difference, and 6 decimals for some other difference using just the same logic. and decimal is the fastest way to access the mantisa bits and operate on them without having to write machine code. – Columbite

The other solutions would not give discontinuities. Decimal is incredibly slow if you are performing many operations using. – Monophagous