Numpy high precision

Question

I am using numpy and pyfits to manipulate spectra and I require high precision (something like 8-10 decimal places on a value which might go as high as 10^12). For that the data type "decimal" would be perfect (float64 is not good enough), but unfortunalely numpy.interp does not like it:

File ".../modules/manip_fits.py", line 47, in get_shift
pix_shift = np.interp(x, xp, fp)-fp
File "/usr/lib/python2.7/dist-packages/numpy/lib/function_base.py", line 1053, in interp
return compiled_interp(x, xp, fp, left, right)
TypeError: array cannot be safely cast to required type

A simplified version of the code I use:

fp = np.array(range(new_wave.shape[-1]),dtype=Decimal)
pix_shift = np.empty_like(wave,dtype=Decimal)
      x = wave
  xp = new_wave
 pix_shift = np.interp(x, xp, fp)-fp

where 'wave' and 'new_wave' are a one-dimension numpy array representing a 1D spectrum. This code is needed to shift my spectra along the x-axis (which is wavelenght)

My biggest issue is that further down the code I divide my spectra by a template spectrum constructed from the sum of all my spectra in order to analyse the differences and since I do not have enough decimal places I am getting rounding errors. Any ideas?

Thanks!

UPDATE:

Test Example:

import numpy as np
from decimal import *
getcontext().prec = 12

wave = np.array([Decimal(xx*np.pi) for xx in range(0,10)],dtype=np.dtype(Decimal))
new_wave = np.array([Decimal(xx*np.pi+0.5) for xx in range(0,10)],dtype=np.dtype(Decimal))

fp = np.array(range(new_wave.shape[-1]),dtype=Decimal)
pix_shift = np.empty_like(wave,dtype=Decimal)

x = wave
xp = new_wave
pix_shift = np.interp(x, xp, fp)-fp

The error is:

Traceback (most recent call last):
  File "untitled.py", line 16, in <module>
    pix_shift = np.interp(x, xp, fp)-fp
  File "/usr/lib/python2.7/dist-packages/numpy/lib/function_base.py", line 1053, in interp
    return compiled_interp(x, xp, fp, left, right)
TypeError: array cannot be safely cast to required type

this is the closest I can provide without using the real spectra in fits format.

UPDATE 2: Some typical values of my spectra, printed using Decimal:

  18786960689.118938446044921875
  18473926205.282184600830078125
  18325454516.792461395263671875
  18400241010.149127960205078125
2577901751996.03857421875
2571812230557.63330078125
2567431795280.80712890625

the problem I am getting is when I make operations between them, I get rounding up errors. For instance, I create a template for all spectra by summing all of them. Then I use this template to normalize every spectra. An example:

Spectra = np.array([Spectrum1, Spectrum2, ...])
Template = np.nansum(Spectra, axis= 0)

NormSpectra = Spectra/Template

This should return me only the noise on the spectra (assuming that the template is a good representation of the star). I tried normalizing each spectrum to its total flux

(Spectrum1 = Spectrum1/np.nansum(Spectrum1), ...) 

as well as the template, but would get even worse rounding up errors.

Using Decimal would work fine for me, but I need to "shift" my spectra so all spectral features/lines are aligned.

Hope this makes sense?