source: thirdparty/SZ/sz/src/szd_double_pwr.c @ 9ee2ce3

Revision 9ee2ce3, 41.8 KB checked in by Hal Finkel <hfinkel@…>, 6 years ago (diff)

importing new SZ files
  • Property mode set to 100644
RevLine 
[2c47b73]1/**
2 *  @file szd_double_pwr.c
3 *  @author Sheng Di
4 *  @date May, 2016
5 *  @brief
6 *  (C) 2016 by Mathematics and Computer Science (MCS), Argonne National Laboratory.
7 *      See COPYRIGHT in top-level directory.
8 */
9
10#include <stdlib.h>
11#include <stdio.h>
12#include <string.h>
13#include "TightDataPointStorageD.h"
14#include "sz.h"
15#include "Huffman.h"
[9ee2ce3]16#include "utility.h"
[2c47b73]17//#include "rw.h"
18
19#pragma GCC diagnostic push
20#pragma GCC diagnostic ignored "-Wchar-subscripts"
21
22void decompressDataSeries_double_1D_pwr(double** data, size_t dataSeriesLength, TightDataPointStorageD* tdps) 
23{
24        updateQuantizationInfo(tdps->intervals);
25        unsigned char tmpPrecBytes[8] = {0}; //used when needing to convert bytes to double values
26        unsigned char* bp = tdps->pwrErrBoundBytes;
27        size_t i, j, k = 0, p = 0, l = 0; // k is to track the location of residual_bit
28                                                                // in resiMidBits, p is to track the
29                                                                // byte_index of resiMidBits, l is for
30                                                                // leadNum
31
32        unsigned char* leadNum;
33        double interval = 0;// = (double)tdps->realPrecision*2;
34       
35        convertByteArray2IntArray_fast_2b(tdps->exactDataNum, tdps->leadNumArray, tdps->leadNumArray_size, &leadNum);
36        *data = (double*)malloc(sizeof(double)*dataSeriesLength);
37
38        int* type = (int*)malloc(dataSeriesLength*sizeof(int));
39
40        HuffmanTree* huffmanTree = createHuffmanTree(tdps->stateNum);
41        decode_withTree(huffmanTree, tdps->typeArray, dataSeriesLength, type);
42        SZ_ReleaseHuffman(huffmanTree); 
43       
44        unsigned char preBytes[8];
45        unsigned char curBytes[8];
46       
47        memset(preBytes, 0, 8);
48
49        size_t curByteIndex = 0;
50        int reqLength = 0, reqBytesLength = 0, resiBitsLength = 0, resiBits = 0; 
51        unsigned char leadingNum;       
52        double medianValue, exactData, predValue = 0, realPrecision = 0;
53       
54        medianValue = tdps->medianValue;
55       
56        int type_, updateReqLength = 0;
57        for (i = 0; i < dataSeriesLength; i++) 
58        {
59                if(i%tdps->segment_size==0)
60                {
61                        tmpPrecBytes[0] = *(bp++);
62                        tmpPrecBytes[1] = *(bp++);
63                        memset(&tmpPrecBytes[2], 0, 6*sizeof(unsigned char));
64
65                        realPrecision = bytesToDouble(tmpPrecBytes);
66                        interval = realPrecision*2;
67                        updateReqLength = 0;
68                }
69                type_ = type[i];
70                switch (type_) {
71                case 0:
72                        // compute resiBits
73                        if(updateReqLength==0)
74                        {
75                                computeReqLength_double(realPrecision, tdps->radExpo, &reqLength, &medianValue);
76                                reqBytesLength = reqLength/8;
77                                resiBitsLength = reqLength%8;   
78                                updateReqLength = 1;   
79                        }
80                        resiBits = 0;
81                        if (resiBitsLength != 0) {
82                                int kMod8 = k % 8;
83                                int rightMovSteps = getRightMovingSteps(kMod8, resiBitsLength);
84                                if (rightMovSteps > 0) {
85                                        int code = getRightMovingCode(kMod8, resiBitsLength);
86                                        resiBits = (tdps->residualMidBits[p] & code) >> rightMovSteps;
87                                } else if (rightMovSteps < 0) {
88                                        int code1 = getLeftMovingCode(kMod8);
89                                        int code2 = getRightMovingCode(kMod8, resiBitsLength);
90                                        int leftMovSteps = -rightMovSteps;
91                                        rightMovSteps = 8 - leftMovSteps;
92                                        resiBits = (tdps->residualMidBits[p] & code1) << leftMovSteps;
93                                        p++;
94                                        resiBits = resiBits
95                                                        | ((tdps->residualMidBits[p] & code2) >> rightMovSteps);
96                                } else // rightMovSteps == 0
97                                {
98                                        int code = getRightMovingCode(kMod8, resiBitsLength);
99                                        resiBits = (tdps->residualMidBits[p] & code);
100                                        p++;
101                                }
102                                k += resiBitsLength;
103                        }
104
105                        // recover the exact data
106                        memset(curBytes, 0, 8);
107                        leadingNum = leadNum[l++];
108                        memcpy(curBytes, preBytes, leadingNum);
109                        for (j = leadingNum; j < reqBytesLength; j++)
110                                curBytes[j] = tdps->exactMidBytes[curByteIndex++];
111                        if (resiBitsLength != 0) {
112                                unsigned char resiByte = (unsigned char) (resiBits << (8 - resiBitsLength));
113                                curBytes[reqBytesLength] = resiByte;
114                        }
115                       
116                        exactData = bytesToDouble(curBytes);
117                        (*data)[i] = exactData + medianValue;
118                        memcpy(preBytes,curBytes,8);
119                        break;
120                default:
121                        //predValue = 2 * (*data)[i-1] - (*data)[i-2];
122                        predValue = (*data)[i-1];
123                        (*data)[i] = predValue + (type_-exe_params->intvRadius)*interval;
124                        break;
125                }
126                //printf("%.30G\n",(*data)[i]);
127        }
128        free(leadNum);
129        free(type);
130        return;
131}
132
133double* extractRealPrecision_2D_double(size_t R1, size_t R2, int blockSize, TightDataPointStorageD* tdps)
134{
135        size_t i,j,k=0, I;
136        unsigned char* bytes = tdps->pwrErrBoundBytes;
137        unsigned char tmpBytes[8] = {0};
138        double* result = (double*)malloc(sizeof(double)*R1*R2);
139        for(i=0;i<R1;i++)
140        {
141                I = i*R2;
142                for(j=0;j<R2;j++)
143                {
144                        tmpBytes[0] = bytes[k++];
145                        tmpBytes[1] = bytes[k++];
146                        result[I+j]=bytesToDouble(tmpBytes);
147                }
148        }
149        return result;
150}
151
152void decompressDataSeries_double_2D_pwr(double** data, size_t r1, size_t r2, TightDataPointStorageD* tdps) 
153{
154        updateQuantizationInfo(tdps->intervals);
155        //printf("tdps->intervals=%d, exe_params->intvRadius=%d\n", tdps->intervals, exe_params->intvRadius);
156       
157        size_t j, k = 0, p = 0, l = 0; // k is to track the location of residual_bit
158        // in resiMidBits, p is to track the
159        // byte_index of resiMidBits, l is for
160        // leadNum
161        size_t dataSeriesLength = r1*r2;
162        //      printf ("%d %d\n", r1, r2);
163
164        unsigned char* leadNum;
165
166        convertByteArray2IntArray_fast_2b(tdps->exactDataNum, tdps->leadNumArray, tdps->leadNumArray_size, &leadNum);
167
168        *data = (double*)malloc(sizeof(double)*dataSeriesLength);
169
170        int* type = (int*)malloc(dataSeriesLength*sizeof(int));
171
172        HuffmanTree* huffmanTree = createHuffmanTree(tdps->stateNum);
173        decode_withTree(huffmanTree, tdps->typeArray, dataSeriesLength, type);
174        SZ_ReleaseHuffman(huffmanTree); 
175
176        unsigned char preBytes[8];
177        unsigned char curBytes[8];
178
179        memset(preBytes, 0, 8);
180
181        size_t curByteIndex = 0;
182        int reqLength, reqBytesLength, resiBitsLength, resiBits; 
183        unsigned char leadingNum;       
184        double medianValue, exactData, realPrecision;
185        int type_;
186        double pred1D, pred2D;
187        size_t ii, jj, II = 0, JJ = 0, updateReqLength = 1;
188
189        int blockSize = computeBlockEdgeSize_2D(tdps->segment_size);
190        size_t R1 = 1+(r1-1)/blockSize;
191        size_t R2 = 1+(r2-1)/blockSize;         
192        double* pwrErrBound = extractRealPrecision_2D_double(R1, R2, blockSize, tdps);
193
194        realPrecision = pwrErrBound[0]; 
195        computeReqLength_double(realPrecision, tdps->radExpo, &reqLength, &medianValue);
196        reqBytesLength = reqLength/8;
197        resiBitsLength = reqLength%8;
198
199        /* Process Row-0, data 0 */
200
201        // compute resiBits
202        resiBits = 0;
203        if (resiBitsLength != 0) {
204                int kMod8 = k % 8;
205                int rightMovSteps = getRightMovingSteps(kMod8, resiBitsLength);
206                if (rightMovSteps > 0) {
207                        int code = getRightMovingCode(kMod8, resiBitsLength);
208                        resiBits = (tdps->residualMidBits[p] & code) >> rightMovSteps;
209                } else if (rightMovSteps < 0) {
210                        int code1 = getLeftMovingCode(kMod8);
211                        int code2 = getRightMovingCode(kMod8, resiBitsLength);
212                        int leftMovSteps = -rightMovSteps;
213                        rightMovSteps = 8 - leftMovSteps;
214                        resiBits = (tdps->residualMidBits[p] & code1) << leftMovSteps;
215                        p++;
216                        resiBits = resiBits
217                                        | ((tdps->residualMidBits[p] & code2) >> rightMovSteps);
218                } else // rightMovSteps == 0
219                {
220                        int code = getRightMovingCode(kMod8, resiBitsLength);
221                        resiBits = (tdps->residualMidBits[p] & code);
222                        p++;
223                }
224                k += resiBitsLength;
225        }
226
227        // recover the exact data
228        memset(curBytes, 0, 8);
229        leadingNum = leadNum[l++];
230        memcpy(curBytes, preBytes, leadingNum);
231        for (j = leadingNum; j < reqBytesLength; j++)
232                curBytes[j] = tdps->exactMidBytes[curByteIndex++];
233        if (resiBitsLength != 0) {
234                unsigned char resiByte = (unsigned char) (resiBits << (8 - resiBitsLength));
235                curBytes[reqBytesLength] = resiByte;
236        }
237
238        exactData = bytesToDouble(curBytes);
239        (*data)[0] = exactData + medianValue;
240        memcpy(preBytes,curBytes,8);
241
242        /* Process Row-0, data 1 */
243        type_ = type[1]; 
244        if (type_ != 0)
245        {
246                pred1D = (*data)[0];           
247                (*data)[1] = pred1D + 2 * (type_ - exe_params->intvRadius) * realPrecision;
248        }
249        else
250        {
251                // compute resiBits
252                resiBits = 0;
253                if (resiBitsLength != 0) {
254                        int kMod8 = k % 8;
255                        int rightMovSteps = getRightMovingSteps(kMod8, resiBitsLength);
256                        if (rightMovSteps > 0) {
257                                int code = getRightMovingCode(kMod8, resiBitsLength);
258                                resiBits = (tdps->residualMidBits[p] & code) >> rightMovSteps;
259                        } else if (rightMovSteps < 0) {
260                                int code1 = getLeftMovingCode(kMod8);
261                                int code2 = getRightMovingCode(kMod8, resiBitsLength);
262                                int leftMovSteps = -rightMovSteps;
263                                rightMovSteps = 8 - leftMovSteps;
264                                resiBits = (tdps->residualMidBits[p] & code1) << leftMovSteps;
265                                p++;
266                                resiBits = resiBits
267                                                | ((tdps->residualMidBits[p] & code2) >> rightMovSteps);
268                        } else // rightMovSteps == 0
269                        {
270                                int code = getRightMovingCode(kMod8, resiBitsLength);
271                                resiBits = (tdps->residualMidBits[p] & code);
272                                p++;
273                        }
274                        k += resiBitsLength;
275                }
276
277                // recover the exact data
278                memset(curBytes, 0, 8);
279                leadingNum = leadNum[l++];
280                memcpy(curBytes, preBytes, leadingNum);
281                for (j = leadingNum; j < reqBytesLength; j++)
282                        curBytes[j] = tdps->exactMidBytes[curByteIndex++];
283                if (resiBitsLength != 0) {
284                        unsigned char resiByte = (unsigned char) (resiBits << (8 - resiBitsLength));
285                        curBytes[reqBytesLength] = resiByte;
286                }
287               
288                exactData = bytesToDouble(curBytes);
289                (*data)[1] = exactData + medianValue;
290                memcpy(preBytes,curBytes,8);
291        }
292
293        /* Process Row-0, data 2 --> data r2-1 */
294        for (jj = 2; jj < r2; jj++)
295        {
296                if(jj%blockSize==0)
297                {
298                        II = 0;
299                        JJ++;
300                        realPrecision = pwrErrBound[JJ];
301                        updateReqLength = 0;                   
302                }               
303               
304                type_ = type[jj];
305                if (type_ != 0)
306                {                       
307                        pred1D = 2*(*data)[jj-1] - (*data)[jj-2];
308                        (*data)[jj] = pred1D + 2 * (type_ - exe_params->intvRadius) * realPrecision;
309                }
310                else
311                {
312                        if(updateReqLength==0)
313                        {
314                                computeReqLength_double(realPrecision, tdps->radExpo, &reqLength, &medianValue);
315                                reqBytesLength = reqLength/8;
316                                resiBitsLength = reqLength%8;                           
317                                updateReqLength = 1;
318                        }                       
319                        // compute resiBits
320                        resiBits = 0;
321                        if (resiBitsLength != 0) {
322                                int kMod8 = k % 8;
323                                int rightMovSteps = getRightMovingSteps(kMod8, resiBitsLength);
324                                if (rightMovSteps > 0) {
325                                        int code = getRightMovingCode(kMod8, resiBitsLength);
326                                        resiBits = (tdps->residualMidBits[p] & code) >> rightMovSteps;
327                                } else if (rightMovSteps < 0) {
328                                        int code1 = getLeftMovingCode(kMod8);
329                                        int code2 = getRightMovingCode(kMod8, resiBitsLength);
330                                        int leftMovSteps = -rightMovSteps;
331                                        rightMovSteps = 8 - leftMovSteps;
332                                        resiBits = (tdps->residualMidBits[p] & code1) << leftMovSteps;
333                                        p++;
334                                        resiBits = resiBits
335                                                        | ((tdps->residualMidBits[p] & code2) >> rightMovSteps);
336                                } else // rightMovSteps == 0
337                                {
338                                        int code = getRightMovingCode(kMod8, resiBitsLength);
339                                        resiBits = (tdps->residualMidBits[p] & code);
340                                        p++;
341                                }
342                                k += resiBitsLength;
343                        }
344
345                        // recover the exact data
346                        memset(curBytes, 0, 8);
347                        leadingNum = leadNum[l++];
348                        memcpy(curBytes, preBytes, leadingNum);
349                        for (j = leadingNum; j < reqBytesLength; j++)
350                                curBytes[j] = tdps->exactMidBytes[curByteIndex++];
351                        if (resiBitsLength != 0) {
352                                unsigned char resiByte = (unsigned char) (resiBits << (8 - resiBitsLength));
353                                curBytes[reqBytesLength] = resiByte;
354                        }
355
356                        exactData = bytesToDouble(curBytes);
357                        (*data)[jj] = exactData + medianValue;
358                        memcpy(preBytes,curBytes,8);
359                }
360        }
361
362        size_t index;
363        /* Process Row-1 --> Row-r1-1 */
364        for (ii = 1; ii < r1; ii++)
365        {
366                /* Process row-ii data 0 */
367                if(ii%blockSize==0)
368                        II++;
369                JJ = 0;
370                realPrecision = pwrErrBound[II*R2+JJ];                         
371                updateReqLength = 0;
372
373                index = ii*r2;
374
375                type_ = type[index];
376                if (type_ != 0)
377                {
378                        pred1D = (*data)[index-r2];
379                        (*data)[index] = pred1D + 2 * (type_ - exe_params->intvRadius) * realPrecision;
380                }
381                else
382                {
383                        // compute resiBits
384                        if(updateReqLength==0)
385                        {
386                                computeReqLength_double(realPrecision, tdps->radExpo, &reqLength, &medianValue);
387                                reqBytesLength = reqLength/8;
388                                resiBitsLength = reqLength%8;                           
389                                updateReqLength = 1;
390                        }
391                       
392                        resiBits = 0;
393                        if (resiBitsLength != 0) {
394                                int kMod8 = k % 8;
395                                int rightMovSteps = getRightMovingSteps(kMod8, resiBitsLength);
396                                if (rightMovSteps > 0) {
397                                        int code = getRightMovingCode(kMod8, resiBitsLength);
398                                        resiBits = (tdps->residualMidBits[p] & code) >> rightMovSteps;
399                                } else if (rightMovSteps < 0) {
400                                        int code1 = getLeftMovingCode(kMod8);
401                                        int code2 = getRightMovingCode(kMod8, resiBitsLength);
402                                        int leftMovSteps = -rightMovSteps;
403                                        rightMovSteps = 8 - leftMovSteps;
404                                        resiBits = (tdps->residualMidBits[p] & code1) << leftMovSteps;
405                                        p++;
406                                        resiBits = resiBits
407                                                        | ((tdps->residualMidBits[p] & code2) >> rightMovSteps);
408                                } else // rightMovSteps == 0
409                                {
410                                        int code = getRightMovingCode(kMod8, resiBitsLength);
411                                        resiBits = (tdps->residualMidBits[p] & code);
412                                        p++;
413                                }
414                                k += resiBitsLength;
415                        }
416
417                        // recover the exact data
418                        memset(curBytes, 0, 8);
419                        leadingNum = leadNum[l++];
420                        memcpy(curBytes, preBytes, leadingNum);
421                        for (j = leadingNum; j < reqBytesLength; j++)
422                                curBytes[j] = tdps->exactMidBytes[curByteIndex++];
423                        if (resiBitsLength != 0) {
424                                unsigned char resiByte = (unsigned char) (resiBits << (8 - resiBitsLength));
425                                curBytes[reqBytesLength] = resiByte;
426                        }
427
428                        exactData = bytesToDouble(curBytes);
429                        (*data)[index] = exactData + medianValue;
430                        memcpy(preBytes,curBytes,8);
431                }
432
433                /* Process row-ii data 1 --> r2-1*/
434                for (jj = 1; jj < r2; jj++)
435                {
436                        index = ii*r2+jj;
437
438                        if(jj%blockSize==0)
439                                JJ++;
440                        realPrecision = pwrErrBound[II*R2+JJ];                 
441                        updateReqLength = 0;
442
443                        type_ = type[index];
444                        if (type_ != 0)
445                        {
446                                pred2D = (*data)[index-1] + (*data)[index-r2] - (*data)[index-r2-1];
447                                (*data)[index] = pred2D + 2 * (type_ - exe_params->intvRadius) * realPrecision;
448                        }
449                        else
450                        {
451                                // compute resiBits
452                                if(updateReqLength==0)
453                                {
454                                        computeReqLength_double(realPrecision, tdps->radExpo, &reqLength, &medianValue);
455                                        reqBytesLength = reqLength/8;
456                                        resiBitsLength = reqLength%8;                           
457                                        updateReqLength = 1;
458                                }                                               
459                               
460                                resiBits = 0;
461                                if (resiBitsLength != 0) {
462                                        int kMod8 = k % 8;
463                                        int rightMovSteps = getRightMovingSteps(kMod8, resiBitsLength);
464                                        if (rightMovSteps > 0) {
465                                                int code = getRightMovingCode(kMod8, resiBitsLength);
466                                                resiBits = (tdps->residualMidBits[p] & code) >> rightMovSteps;
467                                        } else if (rightMovSteps < 0) {
468                                                int code1 = getLeftMovingCode(kMod8);
469                                                int code2 = getRightMovingCode(kMod8, resiBitsLength);
470                                                int leftMovSteps = -rightMovSteps;
471                                                rightMovSteps = 8 - leftMovSteps;
472                                                resiBits = (tdps->residualMidBits[p] & code1) << leftMovSteps;
473                                                p++;
474                                                resiBits = resiBits
475                                                                | ((tdps->residualMidBits[p] & code2) >> rightMovSteps);
476                                        } else // rightMovSteps == 0
477                                        {
478                                                int code = getRightMovingCode(kMod8, resiBitsLength);
479                                                resiBits = (tdps->residualMidBits[p] & code);
480                                                p++;
481                                        }
482                                        k += resiBitsLength;
483                                }
484
485                                // recover the exact data
486                                memset(curBytes, 0, 8);
487                                leadingNum = leadNum[l++];
488                                memcpy(curBytes, preBytes, leadingNum);
489                                for (j = leadingNum; j < reqBytesLength; j++)
490                                        curBytes[j] = tdps->exactMidBytes[curByteIndex++];
491                                if (resiBitsLength != 0) {
492                                        unsigned char resiByte = (unsigned char) (resiBits << (8 - resiBitsLength));
493                                        curBytes[reqBytesLength] = resiByte;
494                                }
495
496                                exactData = bytesToDouble(curBytes);
497                                (*data)[index] = exactData + medianValue;
498                                memcpy(preBytes,curBytes,8);
499                        }
500                }
501        }
502
503        free(pwrErrBound);
504        free(leadNum);
505        free(type);
506        return;
507}
508
509double* extractRealPrecision_3D_double(size_t R1, size_t R2, size_t R3, int blockSize, TightDataPointStorageD* tdps)
510{
511        size_t i,j,k=0, IR, JR, p = 0;
512        size_t R23 = R2*R3;
513        unsigned char* bytes = tdps->pwrErrBoundBytes;
514        unsigned char tmpBytes[4] = {0};
515        double* result = (double*)malloc(sizeof(double)*R1*R2*R3);
516        for(i=0;i<R1;i++)
517        {
518                IR = i*R23;
519                for(j=0;j<R2;j++)
520                {
521                        JR = j*R3;
522                        for(k=0;k<R3;k++)
523                        {
524                                tmpBytes[0] = bytes[p++];
525                                tmpBytes[1] = bytes[p++];
526                                result[IR+JR+k]=bytesToDouble(tmpBytes);                               
527                        }
528                }
529        }
530        return result;
531}
532
533void decompressDataSeries_double_3D_pwr(double** data, size_t r1, size_t r2, size_t r3, TightDataPointStorageD* tdps) 
534{
535        updateQuantizationInfo(tdps->intervals);
536        size_t j, k = 0, p = 0, l = 0; // k is to track the location of residual_bit
537        // in resiMidBits, p is to track the
538        // byte_index of resiMidBits, l is for
539        // leadNum
540        size_t dataSeriesLength = r1*r2*r3;
541        size_t r23 = r2*r3;
542//      printf ("%d %d %d\n", r1, r2, r3);
543
544        unsigned char* leadNum;
545
546        convertByteArray2IntArray_fast_2b(tdps->exactDataNum, tdps->leadNumArray, tdps->leadNumArray_size, &leadNum);
547
548        *data = (double*)malloc(sizeof(double)*dataSeriesLength);
549
550        int* type = (int*)malloc(dataSeriesLength*sizeof(int));
551
552        HuffmanTree* huffmanTree = createHuffmanTree(tdps->stateNum);
553        decode_withTree(huffmanTree, tdps->typeArray, dataSeriesLength, type);
554        SZ_ReleaseHuffman(huffmanTree); 
555
556        unsigned char preBytes[8];
557        unsigned char curBytes[8];
558
559        memset(preBytes, 0, 8);
560
561        size_t curByteIndex = 0;
562        int reqLength, reqBytesLength, resiBitsLength, resiBits; 
563        unsigned char leadingNum;
564        double medianValue, exactData, realPrecision;
565        int type_;
566        double pred1D, pred2D, pred3D;
567        size_t ii, jj, kk, II = 0, JJ = 0, KK = 0, updateReqLength = 1;
568
569        int blockSize = computeBlockEdgeSize_3D(tdps->segment_size);
570        size_t R1 = 1+(r1-1)/blockSize;
571        size_t R2 = 1+(r2-1)/blockSize;         
572        size_t R3 = 1+(r3-1)/blockSize;
573        size_t R23 = R2*R3;
574        double* pwrErrBound = extractRealPrecision_3D_double(R1, R2, R3, blockSize, tdps);
575
576        realPrecision = pwrErrBound[0]; 
577        computeReqLength_double(realPrecision, tdps->radExpo, &reqLength, &medianValue);
578        reqBytesLength = reqLength/8;
579        resiBitsLength = reqLength%8;
580
581        ///////////////////////////     Process layer-0 ///////////////////////////
582        /* Process Row-0 data 0*/
583        // compute resiBits
584        resiBits = 0;
585        if (resiBitsLength != 0) {
586                int kMod8 = k % 8;
587                int rightMovSteps = getRightMovingSteps(kMod8, resiBitsLength);
588                if (rightMovSteps > 0) {
589                        int code = getRightMovingCode(kMod8, resiBitsLength);
590                        resiBits = (tdps->residualMidBits[p] & code) >> rightMovSteps;
591                } else if (rightMovSteps < 0) {
592                        int code1 = getLeftMovingCode(kMod8);
593                        int code2 = getRightMovingCode(kMod8, resiBitsLength);
594                        int leftMovSteps = -rightMovSteps;
595                        rightMovSteps = 8 - leftMovSteps;
596                        resiBits = (tdps->residualMidBits[p] & code1) << leftMovSteps;
597                        p++;
598                        resiBits = resiBits
599                                        | ((tdps->residualMidBits[p] & code2) >> rightMovSteps);
600                } else // rightMovSteps == 0
601                {
602                        int code = getRightMovingCode(kMod8, resiBitsLength);
603                        resiBits = (tdps->residualMidBits[p] & code);
604                        p++;
605                }
606                k += resiBitsLength;
607        }
608
609        // recover the exact data
610        memset(curBytes, 0, 8);
611        leadingNum = leadNum[l++];
612        memcpy(curBytes, preBytes, leadingNum);
613        for (j = leadingNum; j < reqBytesLength; j++)
614                curBytes[j] = tdps->exactMidBytes[curByteIndex++];
615        if (resiBitsLength != 0) {
616                unsigned char resiByte = (unsigned char) (resiBits << (8 - resiBitsLength));
617                curBytes[reqBytesLength] = resiByte;
618        }
619
620        exactData = bytesToDouble(curBytes);
621        (*data)[0] = exactData + medianValue;
622        memcpy(preBytes,curBytes,8);
623
624        /* Process Row-0, data 1 */
625        pred1D = (*data)[0];
626
627        type_ = type[1];
628        if (type_ != 0)
629        {
630                (*data)[1] = pred1D + 2 * (type_ - exe_params->intvRadius) * realPrecision;
631        }
632        else
633        {
634                // compute resiBits
635                resiBits = 0;
636                if (resiBitsLength != 0) {
637                        int kMod8 = k % 8;
638                        int rightMovSteps = getRightMovingSteps(kMod8, resiBitsLength);
639                        if (rightMovSteps > 0) {
640                                int code = getRightMovingCode(kMod8, resiBitsLength);
641                                resiBits = (tdps->residualMidBits[p] & code) >> rightMovSteps;
642                        } else if (rightMovSteps < 0) {
643                                int code1 = getLeftMovingCode(kMod8);
644                                int code2 = getRightMovingCode(kMod8, resiBitsLength);
645                                int leftMovSteps = -rightMovSteps;
646                                rightMovSteps = 8 - leftMovSteps;
647                                resiBits = (tdps->residualMidBits[p] & code1) << leftMovSteps;
648                                p++;
649                                resiBits = resiBits
650                                                | ((tdps->residualMidBits[p] & code2) >> rightMovSteps);
651                        } else // rightMovSteps == 0
652                        {
653                                int code = getRightMovingCode(kMod8, resiBitsLength);
654                                resiBits = (tdps->residualMidBits[p] & code);
655                                p++;
656                        }
657                        k += resiBitsLength;
658                }
659
660                // recover the exact data
661                memset(curBytes, 0, 8);
662                leadingNum = leadNum[l++];
663                memcpy(curBytes, preBytes, leadingNum);
664                for (j = leadingNum; j < reqBytesLength; j++)
665                        curBytes[j] = tdps->exactMidBytes[curByteIndex++];
666                if (resiBitsLength != 0) {
667                        unsigned char resiByte = (unsigned char) (resiBits << (8 - resiBitsLength));
668                        curBytes[reqBytesLength] = resiByte;
669                }
670
671                exactData = bytesToDouble(curBytes);
672                (*data)[1] = exactData + medianValue;
673                memcpy(preBytes,curBytes,8);
674        }
675
676        /* Process Row-0, data 2 --> data r3-1 */
677        for (jj = 2; jj < r3; jj++)
678        {
679                if(jj%blockSize==0)
680                {
681                        KK = 0;//dimension 1 (top)
682                        II = 0;//dimension 2 (mid)
683                        JJ++;
684                        realPrecision = pwrErrBound[JJ];
685                        updateReqLength = 0;                   
686                }               
687                type_ = type[jj];
688                if (type_ != 0)
689                {
690                        pred1D = 2*(*data)[jj-1] - (*data)[jj-2];
691                        (*data)[jj] = pred1D + 2 * (type_ - exe_params->intvRadius) * realPrecision;
692                }
693                else
694                {
695                        // compute resiBits
696                        if(updateReqLength==0)
697                        {
698                                computeReqLength_double(realPrecision, tdps->radExpo, &reqLength, &medianValue);
699                                reqBytesLength = reqLength/8;
700                                resiBitsLength = reqLength%8;                           
701                                updateReqLength = 1;
702                        }
703
704                        resiBits = 0;
705                        if (resiBitsLength != 0) {
706                                int kMod8 = k % 8;
707                                int rightMovSteps = getRightMovingSteps(kMod8, resiBitsLength);
708                                if (rightMovSteps > 0) {
709                                        int code = getRightMovingCode(kMod8, resiBitsLength);
710                                        resiBits = (tdps->residualMidBits[p] & code) >> rightMovSteps;
711                                } else if (rightMovSteps < 0) {
712                                        int code1 = getLeftMovingCode(kMod8);
713                                        int code2 = getRightMovingCode(kMod8, resiBitsLength);
714                                        int leftMovSteps = -rightMovSteps;
715                                        rightMovSteps = 8 - leftMovSteps;
716                                        resiBits = (tdps->residualMidBits[p] & code1) << leftMovSteps;
717                                        p++;
718                                        resiBits = resiBits
719                                                        | ((tdps->residualMidBits[p] & code2) >> rightMovSteps);
720                                } else // rightMovSteps == 0
721                                {
722                                        int code = getRightMovingCode(kMod8, resiBitsLength);
723                                        resiBits = (tdps->residualMidBits[p] & code);
724                                        p++;
725                                }
726                                k += resiBitsLength;
727                        }
728
729                        // recover the exact data
730                        memset(curBytes, 0, 8);
731                        leadingNum = leadNum[l++];
732                        memcpy(curBytes, preBytes, leadingNum);
733                        for (j = leadingNum; j < reqBytesLength; j++)
734                                curBytes[j] = tdps->exactMidBytes[curByteIndex++];
735                        if (resiBitsLength != 0) {
736                                unsigned char resiByte = (unsigned char) (resiBits << (8 - resiBitsLength));
737                                curBytes[reqBytesLength] = resiByte;
738                        }
739
740                        exactData = bytesToDouble(curBytes);
741                        (*data)[jj] = exactData + medianValue;
742                        memcpy(preBytes,curBytes,8);
743                }
744        }
745
746        size_t index;
747        /* Process Row-1 --> Row-r2-1 */
748        for (ii = 1; ii < r2; ii++)
749        {
750                /* Process row-ii data 0 */             
751                if(ii%blockSize==0)
752                        II++;           
753                JJ = 0;
754                realPrecision = pwrErrBound[II*R3+JJ];
755                updateReqLength = 0;           
756
757                index = ii*r3;
758               
759                type_ = type[index];
760                if (type_ != 0)
761                {
762                        pred1D = (*data)[index-r3];                     
763                        (*data)[index] = pred1D + 2 * (type_ - exe_params->intvRadius) * realPrecision;
764                }
765                else
766                {
767                        // compute resiBits
768                        if(updateReqLength==0)
769                        {
770                                computeReqLength_double(realPrecision, tdps->radExpo, &reqLength, &medianValue);
771                                reqBytesLength = reqLength/8;
772                                resiBitsLength = reqLength%8;                           
773                                updateReqLength = 1;
774                        }
775                       
776                        resiBits = 0;
777                        if (resiBitsLength != 0) {
778                                int kMod8 = k % 8;
779                                int rightMovSteps = getRightMovingSteps(kMod8, resiBitsLength);
780                                if (rightMovSteps > 0) {
781                                        int code = getRightMovingCode(kMod8, resiBitsLength);
782                                        resiBits = (tdps->residualMidBits[p] & code) >> rightMovSteps;
783                                } else if (rightMovSteps < 0) {
784                                        int code1 = getLeftMovingCode(kMod8);
785                                        int code2 = getRightMovingCode(kMod8, resiBitsLength);
786                                        int leftMovSteps = -rightMovSteps;
787                                        rightMovSteps = 8 - leftMovSteps;
788                                        resiBits = (tdps->residualMidBits[p] & code1) << leftMovSteps;
789                                        p++;
790                                        resiBits = resiBits
791                                                        | ((tdps->residualMidBits[p] & code2) >> rightMovSteps);
792                                } else // rightMovSteps == 0
793                                {
794                                        int code = getRightMovingCode(kMod8, resiBitsLength);
795                                        resiBits = (tdps->residualMidBits[p] & code);
796                                        p++;
797                                }
798                                k += resiBitsLength;
799                        }
800
801                        // recover the exact data
802                        memset(curBytes, 0, 8);
803                        leadingNum = leadNum[l++];
804                        memcpy(curBytes, preBytes, leadingNum);
805                        for (j = leadingNum; j < reqBytesLength; j++)
806                                curBytes[j] = tdps->exactMidBytes[curByteIndex++];
807                        if (resiBitsLength != 0) {
808                                unsigned char resiByte = (unsigned char) (resiBits << (8 - resiBitsLength));
809                                curBytes[reqBytesLength] = resiByte;
810                        }
811
812                        exactData = bytesToDouble(curBytes);
813                        (*data)[index] = exactData + medianValue;
814                        memcpy(preBytes,curBytes,8);
815                }
816
817                /* Process row-ii data 1 --> r3-1*/
818                for (jj = 1; jj < r3; jj++)
819                {
820                        index = ii*r3+jj;
821
822                        if(jj%blockSize==0)
823                                JJ++;
824                        realPrecision = pwrErrBound[II*R3+JJ];                 
825                        updateReqLength = 0;                   
826                       
827                        type_ = type[index];
828                        if (type_ != 0)
829                        {
830                                pred2D = (*data)[index-1] + (*data)[index-r3] - (*data)[index-r3-1];                           
831                                (*data)[index] = pred2D + 2 * (type_ - exe_params->intvRadius) * realPrecision;
832                        }
833                        else
834                        {
835                                // compute resiBits
836                                if(updateReqLength==0)
837                                {
838                                        computeReqLength_double(realPrecision, tdps->radExpo, &reqLength, &medianValue);
839                                        reqBytesLength = reqLength/8;
840                                        resiBitsLength = reqLength%8;                           
841                                        updateReqLength = 1;
842                                }
843
844                                resiBits = 0;
845                                if (resiBitsLength != 0) {
846                                        int kMod8 = k % 8;
847                                        int rightMovSteps = getRightMovingSteps(kMod8, resiBitsLength);
848                                        if (rightMovSteps > 0) {
849                                                int code = getRightMovingCode(kMod8, resiBitsLength);
850                                                resiBits = (tdps->residualMidBits[p] & code) >> rightMovSteps;
851                                        } else if (rightMovSteps < 0) {
852                                                int code1 = getLeftMovingCode(kMod8);
853                                                int code2 = getRightMovingCode(kMod8, resiBitsLength);
854                                                int leftMovSteps = -rightMovSteps;
855                                                rightMovSteps = 8 - leftMovSteps;
856                                                resiBits = (tdps->residualMidBits[p] & code1) << leftMovSteps;
857                                                p++;
858                                                resiBits = resiBits
859                                                                | ((tdps->residualMidBits[p] & code2) >> rightMovSteps);
860                                        } else // rightMovSteps == 0
861                                        {
862                                                int code = getRightMovingCode(kMod8, resiBitsLength);
863                                                resiBits = (tdps->residualMidBits[p] & code);
864                                                p++;
865                                        }
866                                        k += resiBitsLength;
867                                }
868
869                                // recover the exact data
870                                memset(curBytes, 0, 8);
871                                leadingNum = leadNum[l++];
872                                memcpy(curBytes, preBytes, leadingNum);
873                                for (j = leadingNum; j < reqBytesLength; j++)
874                                        curBytes[j] = tdps->exactMidBytes[curByteIndex++];
875                                if (resiBitsLength != 0) {
876                                        unsigned char resiByte = (unsigned char) (resiBits << (8 - resiBitsLength));
877                                        curBytes[reqBytesLength] = resiByte;
878                                }
879
880                                exactData = bytesToDouble(curBytes);
881                                (*data)[index] = exactData + medianValue;
882                                memcpy(preBytes,curBytes,8);
883                        }
884                }
885        }
886
887        ///////////////////////////     Process layer-1 --> layer-r1-1 ///////////////////////////
888
889        for (kk = 1; kk < r1; kk++)
890        {
891                /* Process Row-0 data 0*/
892                index = kk*r23;         
893                if(kk%blockSize==0)
894                        KK++;
895                II = 0;
896                JJ = 0;
897
898                realPrecision = pwrErrBound[KK*R23];                   
899                updateReqLength = 0;                   
900
901                type_ = type[index];
902                if (type_ != 0)
903                {
904                        pred1D = (*data)[index-r23];                   
905                        (*data)[index] = pred1D + 2 * (type_ - exe_params->intvRadius) * realPrecision;
906                }
907                else
908                {
909                        // compute resiBits
910                        if(updateReqLength==0)
911                        {
912                                computeReqLength_double(realPrecision, tdps->radExpo, &reqLength, &medianValue);
913                                reqBytesLength = reqLength/8;
914                                resiBitsLength = reqLength%8;                           
915                                updateReqLength = 1;
916                        }
917
918                        resiBits = 0;
919                        if (resiBitsLength != 0) {
920                                int kMod8 = k % 8;
921                                int rightMovSteps = getRightMovingSteps(kMod8, resiBitsLength);
922                                if (rightMovSteps > 0) {
923                                        int code = getRightMovingCode(kMod8, resiBitsLength);
924                                        resiBits = (tdps->residualMidBits[p] & code) >> rightMovSteps;
925                                } else if (rightMovSteps < 0) {
926                                        int code1 = getLeftMovingCode(kMod8);
927                                        int code2 = getRightMovingCode(kMod8, resiBitsLength);
928                                        int leftMovSteps = -rightMovSteps;
929                                        rightMovSteps = 8 - leftMovSteps;
930                                        resiBits = (tdps->residualMidBits[p] & code1) << leftMovSteps;
931                                        p++;
932                                        resiBits = resiBits
933                                                        | ((tdps->residualMidBits[p] & code2) >> rightMovSteps);
934                                } else // rightMovSteps == 0
935                                {
936                                        int code = getRightMovingCode(kMod8, resiBitsLength);
937                                        resiBits = (tdps->residualMidBits[p] & code);
938                                        p++;
939                                }
940                                k += resiBitsLength;
941                        }
942
943                        // recover the exact data
944                        memset(curBytes, 0, 8);
945                        leadingNum = leadNum[l++];
946                        memcpy(curBytes, preBytes, leadingNum);
947                        for (j = leadingNum; j < reqBytesLength; j++)
948                                curBytes[j] = tdps->exactMidBytes[curByteIndex++];
949                        if (resiBitsLength != 0) {
950                                unsigned char resiByte = (unsigned char) (resiBits << (8 - resiBitsLength));
951                                curBytes[reqBytesLength] = resiByte;
952                        }
953
954                        exactData = bytesToDouble(curBytes);
955                        (*data)[index] = exactData + medianValue;
956                        memcpy(preBytes,curBytes,8);
957                }
958
959                /* Process Row-0 data 1 --> data r3-1 */
960                for (jj = 1; jj < r3; jj++)
961                {
962                        index = kk*r23+jj;
963
964                        if(jj%blockSize==0)
965                                JJ++;
966
967                        realPrecision = pwrErrBound[KK*R23+JJ];                 
968                        updateReqLength = 0;                   
969                       
970                        type_ = type[index];
971                        if (type_ != 0)
972                        {
973                                pred2D = (*data)[index-1] + (*data)[index-r23] - (*data)[index-r23-1];                 
974                                (*data)[index] = pred2D + 2 * (type_ - exe_params->intvRadius) * realPrecision;
975                        }
976                        else
977                        {
978                                // compute resiBits
979                                if(updateReqLength==0)
980                                {
981                                        computeReqLength_double(realPrecision, tdps->radExpo, &reqLength, &medianValue);
982                                        reqBytesLength = reqLength/8;
983                                        resiBitsLength = reqLength%8;                           
984                                        updateReqLength = 1;
985                                }
986                       
987                                resiBits = 0;
988                                if (resiBitsLength != 0) {
989                                        int kMod8 = k % 8;
990                                        int rightMovSteps = getRightMovingSteps(kMod8, resiBitsLength);
991                                        if (rightMovSteps > 0) {
992                                                int code = getRightMovingCode(kMod8, resiBitsLength);
993                                                resiBits = (tdps->residualMidBits[p] & code) >> rightMovSteps;
994                                        } else if (rightMovSteps < 0) {
995                                                int code1 = getLeftMovingCode(kMod8);
996                                                int code2 = getRightMovingCode(kMod8, resiBitsLength);
997                                                int leftMovSteps = -rightMovSteps;
998                                                rightMovSteps = 8 - leftMovSteps;
999                                                resiBits = (tdps->residualMidBits[p] & code1) << leftMovSteps;
1000                                                p++;
1001                                                resiBits = resiBits
1002                                                                | ((tdps->residualMidBits[p] & code2) >> rightMovSteps);
1003                                        } else // rightMovSteps == 0
1004                                        {
1005                                                int code = getRightMovingCode(kMod8, resiBitsLength);
1006                                                resiBits = (tdps->residualMidBits[p] & code);
1007                                                p++;
1008                                        }
1009                                        k += resiBitsLength;
1010                                }
1011
1012                                // recover the exact data
1013                                memset(curBytes, 0, 8);
1014                                leadingNum = leadNum[l++];
1015                                memcpy(curBytes, preBytes, leadingNum);
1016                                for (j = leadingNum; j < reqBytesLength; j++)
1017                                        curBytes[j] = tdps->exactMidBytes[curByteIndex++];
1018                                if (resiBitsLength != 0) {
1019                                        unsigned char resiByte = (unsigned char) (resiBits << (8 - resiBitsLength));
1020                                        curBytes[reqBytesLength] = resiByte;
1021                                }
1022
1023                                exactData = bytesToDouble(curBytes);
1024                                (*data)[index] = exactData + medianValue;
1025                                memcpy(preBytes,curBytes,8);
1026                        }
1027                }
1028
1029                /* Process Row-1 --> Row-r2-1 */
1030                for (ii = 1; ii < r2; ii++)
1031                {
1032                        /* Process Row-i data 0 */
1033                        index = kk*r23 + ii*r3;
1034                       
1035                        if(ii%blockSize==0)
1036                                II++;
1037                        JJ = 0;
1038                       
1039                        realPrecision = pwrErrBound[KK*R23+II*R3];                     
1040                        updateReqLength = 0;                                           
1041
1042                        type_ = type[index];
1043                        if (type_ != 0)
1044                        {
1045                                pred2D = (*data)[index-r3] + (*data)[index-r23] - (*data)[index-r23-r3];                               
1046                                (*data)[index] = pred2D + 2 * (type_ - exe_params->intvRadius) * realPrecision;
1047                        }
1048                        else
1049                        {
1050                                // compute resiBits
1051                                if(updateReqLength==0)
1052                                {
1053                                        computeReqLength_double(realPrecision, tdps->radExpo, &reqLength, &medianValue);
1054                                        reqBytesLength = reqLength/8;
1055                                        resiBitsLength = reqLength%8;                           
1056                                        updateReqLength = 1;
1057                                }
1058
1059                                resiBits = 0;
1060                                if (resiBitsLength != 0) {
1061                                        int kMod8 = k % 8;
1062                                        int rightMovSteps = getRightMovingSteps(kMod8, resiBitsLength);
1063                                        if (rightMovSteps > 0) {
1064                                                int code = getRightMovingCode(kMod8, resiBitsLength);
1065                                                resiBits = (tdps->residualMidBits[p] & code) >> rightMovSteps;
1066                                        } else if (rightMovSteps < 0) {
1067                                                int code1 = getLeftMovingCode(kMod8);
1068                                                int code2 = getRightMovingCode(kMod8, resiBitsLength);
1069                                                int leftMovSteps = -rightMovSteps;
1070                                                rightMovSteps = 8 - leftMovSteps;
1071                                                resiBits = (tdps->residualMidBits[p] & code1) << leftMovSteps;
1072                                                p++;
1073                                                resiBits = resiBits
1074                                                                | ((tdps->residualMidBits[p] & code2) >> rightMovSteps);
1075                                        } else // rightMovSteps == 0
1076                                        {
1077                                                int code = getRightMovingCode(kMod8, resiBitsLength);
1078                                                resiBits = (tdps->residualMidBits[p] & code);
1079                                                p++;
1080                                        }
1081                                        k += resiBitsLength;
1082                                }
1083
1084                                // recover the exact data
1085                                memset(curBytes, 0, 8);
1086                                leadingNum = leadNum[l++];
1087                                memcpy(curBytes, preBytes, leadingNum);
1088                                for (j = leadingNum; j < reqBytesLength; j++)
1089                                        curBytes[j] = tdps->exactMidBytes[curByteIndex++];
1090                                if (resiBitsLength != 0) {
1091                                        unsigned char resiByte = (unsigned char) (resiBits << (8 - resiBitsLength));
1092                                        curBytes[reqBytesLength] = resiByte;
1093                                }
1094
1095                                exactData = bytesToDouble(curBytes);
1096                                (*data)[index] = exactData + medianValue;
1097                                memcpy(preBytes,curBytes,8);
1098                        }
1099
1100                        /* Process Row-i data 1 --> data r3-1 */
1101                        for (jj = 1; jj < r3; jj++)
1102                        {
1103                                index = kk*r23 + ii*r3 + jj;
1104                                if(jj%blockSize==0)
1105                                        JJ++;
1106
1107                                realPrecision = pwrErrBound[KK*R23+II*R3+JJ];                   
1108                                updateReqLength = 0;                           
1109
1110                                type_ = type[index];
1111                                if (type_ != 0)
1112                                {
1113                                        pred3D = (*data)[index-1] + (*data)[index-r3] + (*data)[index-r23]
1114                                        - (*data)[index-r3-1] - (*data)[index-r23-r3] - (*data)[index-r23-1] + (*data)[index-r23-r3-1];                                 
1115                                        (*data)[index] = pred3D + 2 * (type_ - exe_params->intvRadius) * realPrecision;
1116                                }
1117                                else
1118                                {
1119                                        // compute resiBits
1120                                        if(updateReqLength==0)
1121                                        {
1122                                                computeReqLength_double(realPrecision, tdps->radExpo, &reqLength, &medianValue);
1123                                                reqBytesLength = reqLength/8;
1124                                                resiBitsLength = reqLength%8;                           
1125                                                updateReqLength = 1;
1126                                        }
1127                               
1128                                        resiBits = 0;
1129                                        if (resiBitsLength != 0) {
1130                                                int kMod8 = k % 8;
1131                                                int rightMovSteps = getRightMovingSteps(kMod8, resiBitsLength);
1132                                                if (rightMovSteps > 0) {
1133                                                        int code = getRightMovingCode(kMod8, resiBitsLength);
1134                                                        resiBits = (tdps->residualMidBits[p] & code) >> rightMovSteps;
1135                                                } else if (rightMovSteps < 0) {
1136                                                        int code1 = getLeftMovingCode(kMod8);
1137                                                        int code2 = getRightMovingCode(kMod8, resiBitsLength);
1138                                                        int leftMovSteps = -rightMovSteps;
1139                                                        rightMovSteps = 8 - leftMovSteps;
1140                                                        resiBits = (tdps->residualMidBits[p] & code1) << leftMovSteps;
1141                                                        p++;
1142                                                        resiBits = resiBits
1143                                                                        | ((tdps->residualMidBits[p] & code2) >> rightMovSteps);
1144                                                } else // rightMovSteps == 0
1145                                                {
1146                                                        int code = getRightMovingCode(kMod8, resiBitsLength);
1147                                                        resiBits = (tdps->residualMidBits[p] & code);
1148                                                        p++;
1149                                                }
1150                                                k += resiBitsLength;
1151                                        }
1152
1153                                        // recover the exact data
1154                                        memset(curBytes, 0, 8);
1155                                        leadingNum = leadNum[l++];
1156                                        memcpy(curBytes, preBytes, leadingNum);
1157                                        for (j = leadingNum; j < reqBytesLength; j++)
1158                                                curBytes[j] = tdps->exactMidBytes[curByteIndex++];
1159                                        if (resiBitsLength != 0) {
1160                                                unsigned char resiByte = (unsigned char) (resiBits << (8 - resiBitsLength));
1161                                                curBytes[reqBytesLength] = resiByte;
1162                                        }
1163
1164                                        exactData = bytesToDouble(curBytes);
1165                                        (*data)[index] = exactData + medianValue;
1166                                        memcpy(preBytes,curBytes,8);
1167                                }
1168                        }
1169                }
1170        }
1171
1172        free(pwrErrBound);
1173        free(leadNum);
1174        free(type);
1175        return;
1176}
1177
1178void decompressDataSeries_double_1D_pwrgroup(double** data, size_t dataSeriesLength, TightDataPointStorageD* tdps) 
1179{
1180        double *posGroups, *negGroups, *groups;
1181        double pos_01_group, neg_01_group;
1182        int *posFlags, *negFlags;
1183       
1184        updateQuantizationInfo(tdps->intervals);
1185       
1186        unsigned char* leadNum;
1187        double interval;// = (double)tdps->realPrecision*2;
1188       
1189        convertByteArray2IntArray_fast_2b(tdps->exactDataNum, tdps->leadNumArray, tdps->leadNumArray_size, &leadNum);
1190
1191        *data = (double*)malloc(sizeof(double)*dataSeriesLength);
1192
1193        int* type = (int*)malloc(dataSeriesLength*sizeof(int));
1194
1195        HuffmanTree* huffmanTree = createHuffmanTree(tdps->stateNum);
1196        decode_withTree(huffmanTree, tdps->typeArray, dataSeriesLength, type);
1197        SZ_ReleaseHuffman(huffmanTree); 
1198
1199        createRangeGroups_double(&posGroups, &negGroups, &posFlags, &negFlags);
1200       
1201        double realGroupPrecision;
1202        double realPrecision = tdps->realPrecision;
1203        char* groupID = decompressGroupIDArray(tdps->pwrErrBoundBytes, tdps->dataSeriesLength);
1204       
1205        //note that the groupID values here are [1,2,3,....,18] or [-1,-2,...,-18]
1206       
1207        double* groupErrorBounds = generateGroupErrBounds(confparams_dec->errorBoundMode, realPrecision, confparams_dec->pw_relBoundRatio);
1208        exe_params->intvRadius = generateGroupMaxIntervalCount(groupErrorBounds);
1209               
1210        size_t nbBins = (size_t)(1/confparams_dec->pw_relBoundRatio + 0.5);
1211        if(nbBins%2==1)
1212                nbBins++;
1213        exe_params->intvRadius = nbBins;
1214
1215       
1216
1217        unsigned char preBytes[8];
1218        unsigned char curBytes[8];
1219       
1220        memset(preBytes, 0, 8);
1221
1222        size_t curByteIndex = 0;
1223        int reqLength, reqBytesLength = 0, resiBitsLength = 0, resiBits; 
1224        unsigned char leadingNum;       
1225        double medianValue, exactData, curValue, predValue;
1226       
1227        medianValue = tdps->medianValue;
1228       
1229        size_t i, j, k = 0, p = 0, l = 0; // k is to track the location of residual_bit
1230                                                        // in resiMidBits, p is to track the
1231                                                        // byte_index of resiMidBits, l is for
1232                                                        // leadNum
1233                                                       
1234        int type_, updateReqLength = 0;
1235        char rawGrpID = 0, indexGrpID = 0;
1236        for (i = 0; i < dataSeriesLength; i++) 
1237        {
1238                rawGrpID = groupID[i];
1239               
1240                if(rawGrpID >= 2)
1241                {
1242                        groups = posGroups;
1243                        indexGrpID = rawGrpID - 2;
1244                }
1245                else if(rawGrpID <= -2)
1246                {
1247                        groups = negGroups;
1248                        indexGrpID = -rawGrpID - 2;             }
1249                else if(rawGrpID == 1)
1250                {
1251                        groups = &pos_01_group;
1252                        indexGrpID = 0;
1253                }
1254                else //rawGrpID == -1
1255                {
1256                        groups = &neg_01_group;
1257                        indexGrpID = 0;                 
1258                }
1259               
1260                type_ = type[i];
1261                switch (type_) {
1262                case 0:
1263                        // compute resiBits
1264                        if(updateReqLength==0)
1265                        {
1266                                computeReqLength_double(realPrecision, tdps->radExpo, &reqLength, &medianValue);
1267                                reqBytesLength = reqLength/8;
1268                                resiBitsLength = reqLength%8;   
1269                                updateReqLength = 1;   
1270                        }
1271                        resiBits = 0;
1272                        if (resiBitsLength != 0) {
1273                                int kMod8 = k % 8;
1274                                int rightMovSteps = getRightMovingSteps(kMod8, resiBitsLength);
1275                                if (rightMovSteps > 0) {
1276                                        int code = getRightMovingCode(kMod8, resiBitsLength);
1277                                        resiBits = (tdps->residualMidBits[p] & code) >> rightMovSteps;
1278                                } else if (rightMovSteps < 0) {
1279                                        int code1 = getLeftMovingCode(kMod8);
1280                                        int code2 = getRightMovingCode(kMod8, resiBitsLength);
1281                                        int leftMovSteps = -rightMovSteps;
1282                                        rightMovSteps = 8 - leftMovSteps;
1283                                        resiBits = (tdps->residualMidBits[p] & code1) << leftMovSteps;
1284                                        p++;
1285                                        resiBits = resiBits
1286                                                        | ((tdps->residualMidBits[p] & code2) >> rightMovSteps);
1287                                } else // rightMovSteps == 0
1288                                {
1289                                        int code = getRightMovingCode(kMod8, resiBitsLength);
1290                                        resiBits = (tdps->residualMidBits[p] & code);
1291                                        p++;
1292                                }
1293                                k += resiBitsLength;
1294                        }
1295
1296                        // recover the exact data       
1297                        memset(curBytes, 0, 8);
1298                        leadingNum = leadNum[l++];
1299                        memcpy(curBytes, preBytes, leadingNum);
1300                        for (j = leadingNum; j < reqBytesLength; j++)
1301                                curBytes[j] = tdps->exactMidBytes[curByteIndex++];
1302                        if (resiBitsLength != 0) {
1303                                unsigned char resiByte = (unsigned char) (resiBits << (8 - resiBitsLength));
1304                                curBytes[reqBytesLength] = resiByte;
1305                        }
1306                       
1307                        exactData = bytesToDouble(curBytes);
1308                        exactData = exactData + medianValue;
1309                        (*data)[i] = exactData;
1310                        memcpy(preBytes,curBytes,8);
1311                       
1312                        groups[indexGrpID] = exactData;
1313                       
1314                        break;
1315                default:
1316                        predValue = groups[indexGrpID]; //Here, groups[indexGrpID] is the previous value.
1317                        realGroupPrecision = groupErrorBounds[indexGrpID];
1318                        interval = realGroupPrecision*2;               
1319                       
1320                        curValue = predValue + (type_-exe_params->intvRadius)*interval;
1321                       
1322                        //groupNum = computeGroupNum_double(curValue);
1323                       
1324                        if((curValue>0&&rawGrpID<0)||(curValue<0&&rawGrpID>0))
1325                                curValue = 0;
1326                        //else
1327                        //{
1328                        //      realGrpID = fabs(rawGrpID)-2;
1329                        //      if(groupNum<realGrpID)
1330                        //              curValue = rawGrpID>0?pow(2,realGrpID):-pow(2,realGrpID);
1331                        //      else if(groupNum>realGrpID)
1332                        //              curValue = rawGrpID>0?pow(2,groupNum):-pow(2,groupNum);                         
1333                        //}     
1334                               
1335                        (*data)[i] = curValue;
1336                        groups[indexGrpID] = curValue;
1337                        break;         
1338                }
1339        }       
1340       
1341        free(leadNum);
1342        free(type);
1343       
1344        free(posGroups);
1345        free(negGroups);
1346        free(posFlags);
1347        free(negFlags);
1348        free(groupErrorBounds);
1349        free(groupID);
1350}
[9ee2ce3]1351
1352void decompressDataSeries_double_1D_pwr_pre_log(double** data, size_t dataSeriesLength, TightDataPointStorageD* tdps) {
1353
1354        decompressDataSeries_double_1D(data, dataSeriesLength, tdps);
1355        double threshold = tdps->minLogValue;
1356        if(tdps->pwrErrBoundBytes_size > 0){
1357                unsigned char * signs;
1358                sz_lossless_decompress(confparams_dec->losslessCompressor, tdps->pwrErrBoundBytes, tdps->pwrErrBoundBytes_size, &signs, dataSeriesLength);
1359
1360                for(size_t i=0; i<dataSeriesLength; i++){
1361                        if((*data)[i] < threshold) (*data)[i] = 0;
1362                        else (*data)[i] = exp2((*data)[i]);
1363                        if(signs[i]) (*data)[i] = -((*data)[i]);
1364                }
1365                free(signs);
1366        }
1367        else{
1368                for(size_t i=0; i<dataSeriesLength; i++){
1369                        if((*data)[i] < threshold) (*data)[i] = 0;
1370                        else (*data)[i] = exp2((*data)[i]);
1371                }
1372        }
1373
1374}
1375
1376void decompressDataSeries_double_2D_pwr_pre_log(double** data, size_t r1, size_t r2, TightDataPointStorageD* tdps) {
1377
1378        size_t dataSeriesLength = r1 * r2;
1379        decompressDataSeries_double_2D(data, r1, r2, tdps);
1380        double threshold = tdps->minLogValue;
1381        if(tdps->pwrErrBoundBytes_size > 0){
1382                unsigned char * signs;
1383                sz_lossless_decompress(confparams_dec->losslessCompressor, tdps->pwrErrBoundBytes, tdps->pwrErrBoundBytes_size, &signs, dataSeriesLength);
1384
1385                for(size_t i=0; i<dataSeriesLength; i++){
1386                        if((*data)[i] < threshold) (*data)[i] = 0;
1387                        else (*data)[i] = exp2((*data)[i]);
1388                        if(signs[i]) (*data)[i] = -((*data)[i]);
1389                }
1390                free(signs);
1391        }
1392        else{
1393                for(size_t i=0; i<dataSeriesLength; i++){
1394                        if((*data)[i] < threshold) (*data)[i] = 0;
1395                        else (*data)[i] = exp2((*data)[i]);
1396                }
1397        }
1398}
1399
1400void decompressDataSeries_double_3D_pwr_pre_log(double** data, size_t r1, size_t r2, size_t r3, TightDataPointStorageD* tdps) {
1401
1402        size_t dataSeriesLength = r1 * r2 * r3;
1403        decompressDataSeries_double_3D(data, r1, r2, r3, tdps);
1404        double threshold = tdps->minLogValue;
1405        if(tdps->pwrErrBoundBytes_size > 0){
1406                unsigned char * signs;
1407                sz_lossless_decompress(confparams_dec->losslessCompressor, tdps->pwrErrBoundBytes, tdps->pwrErrBoundBytes_size, &signs, dataSeriesLength);
1408
1409                for(size_t i=0; i<dataSeriesLength; i++){
1410                        if((*data)[i] < threshold) (*data)[i] = 0;
1411                        else (*data)[i] = exp2((*data)[i]);
1412                        if(signs[i]) (*data)[i] = -((*data)[i]);
1413                }
1414                free(signs);
1415        }
1416        else{
1417                for(size_t i=0; i<dataSeriesLength; i++){
1418                        if((*data)[i] < threshold) (*data)[i] = 0;
1419                        else (*data)[i] = exp2((*data)[i]);
1420                }
1421        }
1422}
1423
[2c47b73]1424#pragma GCC diagnostic pop
Note: See TracBrowser for help on using the repository browser.