title: Numpy high precision
body: <p>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:</p>
  
  <pre><code>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
  </code></pre>
  
  <p>A simplified version of the code I use:</p>
  
  <pre><code>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
  </code></pre>
  
  <p>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)</p>
  
  <p>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?</p>
  
  <p>Thanks!</p>
  
  <p><strong>UPDATE:</strong></p>
  
  <p>Test Example:</p>
  
  <pre><code>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
  </code></pre>
  
  <p>The error is:</p>
  
  <pre><code>Traceback (most recent call last):
    File "untitled.py", line 16, in &lt;module&gt;
      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
  </code></pre>
  
  <p>this is the closest I can provide without using the real spectra in fits format.</p>
  
  <p>UPDATE 2:
  Some typical values of my spectra, printed using Decimal:</p>
  
  <pre><code>  18786960689.118938446044921875
    18473926205.282184600830078125
    18325454516.792461395263671875
    18400241010.149127960205078125
  2577901751996.03857421875
  2571812230557.63330078125
  2567431795280.80712890625
  </code></pre>
  
  <p>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:</p>
  
  <pre><code>Spectra = np.array([Spectrum1, Spectrum2, ...])
  Template = np.nansum(Spectra, axis= 0)
  
  NormSpectra = Spectra/Template
  </code></pre>
  
  <p>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 </p>
  
  <pre><code>(Spectrum1 = Spectrum1/np.nansum(Spectrum1), ...) 
  </code></pre>
  
  <p>as well as the template, but would get even worse rounding up errors. </p>
  
  <p>Using Decimal would work fine for me, but I need to "shift" my spectra so all spectral features/lines are aligned.</p>
  
  <p>Hope this makes sense?</p>