source: thirdparty/SZ/sz/src/szd_double_pwr.c @ 2c47b73

Revision 2c47b73, 39.5 KB checked in by Hal Finkel <hfinkel@…>, 6 years ago (diff)

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