1 /* This file, checksum.c, contains the checksum-related routines in the */
2 /* FITSIO library. */
3
4 /* The FITSIO software was written by William Pence at the High Energy */
5 /* Astrophysic Science Archive Research Center (HEASARC) at the NASA */
6 /* Goddard Space Flight Center. */
7
8 #include <string.h>
9 #include <stdlib.h>
10 #include "fitsio2.h"
11 /*------------------------------------------------------------------------*/
ffcsum(fitsfile * fptr,long nrec,unsigned long * sum,int * status)12 int ffcsum(fitsfile *fptr, /* I - FITS file pointer */
13 long nrec, /* I - number of 2880-byte blocks to sum */
14 unsigned long *sum, /* IO - accumulated checksum */
15 int *status) /* IO - error status */
16 /*
17 Calculate a 32-bit 1's complement checksum of the FITS 2880-byte blocks.
18 This routine is based on the C algorithm developed by Rob
19 Seaman at NOAO that was presented at the 1994 ADASS conference,
20 published in the Astronomical Society of the Pacific Conference Series.
21 This uses a 32-bit 1's complement checksum in which the overflow bits
22 are permuted back into the sum and therefore all bit positions are
23 sampled evenly.
24 */
25 {
26 long ii, jj;
27 unsigned short sbuf[1440];
28 unsigned long hi, lo, hicarry, locarry;
29
30 if (*status > 0)
31 return(*status);
32 /*
33 Sum the specified number of FITS 2880-byte records. This assumes that
34 the FITSIO file pointer points to the start of the records to be summed.
35 Read each FITS block as 1440 short values (do byte swapping if needed).
36 */
37 for (jj = 0; jj < nrec; jj++)
38 {
39 ffgbyt(fptr, 2880, sbuf, status);
40
41 #if BYTESWAPPED
42
43 ffswap2( (short *)sbuf, 1440); /* reverse order of bytes in each value */
44
45 #endif
46
47 hi = (*sum >> 16);
48 lo = *sum & 0xFFFF;
49
50 for (ii = 0; ii < 1440; ii += 2)
51 {
52 hi += sbuf[ii];
53 lo += sbuf[ii+1];
54 }
55
56 hicarry = hi >> 16; /* fold carry bits in */
57 locarry = lo >> 16;
58
59 while (hicarry | locarry)
60 {
61 hi = (hi & 0xFFFF) + locarry;
62 lo = (lo & 0xFFFF) + hicarry;
63 hicarry = hi >> 16;
64 locarry = lo >> 16;
65 }
66
67 *sum = (hi << 16) + lo;
68 }
69 return(*status);
70 }
71 /*-------------------------------------------------------------------------*/
ffesum(unsigned long sum,int complm,char * ascii)72 void ffesum(unsigned long sum, /* I - accumulated checksum */
73 int complm, /* I - = 1 to encode complement of the sum */
74 char *ascii) /* O - 16-char ASCII encoded checksum */
75 /*
76 encode the 32 bit checksum by converting every
77 2 bits of each byte into an ASCII character (32 bit word encoded
78 as 16 character string). Only ASCII letters and digits are used
79 to encode the values (no ASCII punctuation characters).
80
81 If complm=TRUE, then the complement of the sum will be encoded.
82
83 This routine is based on the C algorithm developed by Rob
84 Seaman at NOAO that was presented at the 1994 ADASS conference,
85 published in the Astronomical Society of the Pacific Conference Series.
86 */
87 {
88 unsigned int exclude[13] = { 0x3a, 0x3b, 0x3c, 0x3d, 0x3e, 0x3f, 0x40,
89 0x5b, 0x5c, 0x5d, 0x5e, 0x5f, 0x60 };
90 unsigned long mask[4] = { 0xff000000, 0xff0000, 0xff00, 0xff };
91
92 int offset = 0x30; /* ASCII 0 (zero) */
93
94 unsigned long value;
95 int byte, quotient, remainder, ch[4], check, ii, jj, kk;
96 char asc[32];
97
98 if (complm)
99 value = 0xFFFFFFFF - sum; /* complement each bit of the value */
100 else
101 value = sum;
102
103 for (ii = 0; ii < 4; ii++)
104 {
105 byte = (value & mask[ii]) >> (24 - (8 * ii));
106 quotient = byte / 4 + offset;
107 remainder = byte % 4;
108 for (jj = 0; jj < 4; jj++)
109 ch[jj] = quotient;
110
111 ch[0] += remainder;
112
113 for (check = 1; check;) /* avoid ASCII punctuation */
114 for (check = 0, kk = 0; kk < 13; kk++)
115 for (jj = 0; jj < 4; jj += 2)
116 if ((unsigned char) ch[jj] == exclude[kk] ||
117 (unsigned char) ch[jj+1] == exclude[kk])
118 {
119 ch[jj]++;
120 ch[jj+1]--;
121 check++;
122 }
123
124 for (jj = 0; jj < 4; jj++) /* assign the bytes */
125 asc[4*jj+ii] = ch[jj];
126 }
127
128 for (ii = 0; ii < 16; ii++) /* shift the bytes 1 to the right */
129 ascii[ii] = asc[(ii+15)%16];
130
131 ascii[16] = '\0';
132 }
133 /*-------------------------------------------------------------------------*/
ffdsum(char * ascii,int complm,unsigned long * sum)134 unsigned long ffdsum(char *ascii, /* I - 16-char ASCII encoded checksum */
135 int complm, /* I - =1 to decode complement of the */
136 unsigned long *sum) /* O - 32-bit checksum */
137 /*
138 decode the 16-char ASCII encoded checksum into an unsigned 32-bit long.
139 If complm=TRUE, then the complement of the sum will be decoded.
140
141 This routine is based on the C algorithm developed by Rob
142 Seaman at NOAO that was presented at the 1994 ADASS conference,
143 published in the Astronomical Society of the Pacific Conference Series.
144 */
145 {
146 char cbuf[16];
147 unsigned long hi = 0, lo = 0, hicarry, locarry;
148 int ii;
149
150 /* remove the permuted FITS byte alignment and the ASCII 0 offset */
151 for (ii = 0; ii < 16; ii++)
152 {
153 cbuf[ii] = ascii[(ii+1)%16];
154 cbuf[ii] -= 0x30;
155 }
156
157 for (ii = 0; ii < 16; ii += 4)
158 {
159 hi += (cbuf[ii] << 8) + cbuf[ii+1];
160 lo += (cbuf[ii+2] << 8) + cbuf[ii+3];
161 }
162
163 hicarry = hi >> 16;
164 locarry = lo >> 16;
165 while (hicarry || locarry)
166 {
167 hi = (hi & 0xFFFF) + locarry;
168 lo = (lo & 0xFFFF) + hicarry;
169 hicarry = hi >> 16;
170 locarry = lo >> 16;
171 }
172
173 *sum = (hi << 16) + lo;
174 if (complm)
175 *sum = 0xFFFFFFFF - *sum; /* complement each bit of the value */
176
177 return(*sum);
178 }
179 /*------------------------------------------------------------------------*/
ffpcks(fitsfile * fptr,int * status)180 int ffpcks(fitsfile *fptr, /* I - FITS file pointer */
181 int *status) /* IO - error status */
182 /*
183 Create or update the checksum keywords in the CHDU. These keywords
184 provide a checksum verification of the FITS HDU based on the ASCII
185 coded 1's complement checksum algorithm developed by Rob Seaman at NOAO.
186 */
187 {
188 char datestr[20], checksum[FLEN_VALUE], datasum[FLEN_VALUE];
189 char comm[FLEN_COMMENT], chkcomm[FLEN_COMMENT], datacomm[FLEN_COMMENT];
190 int tstatus;
191 long nrec;
192 LONGLONG headstart, datastart, dataend;
193 unsigned long dsum, olddsum, sum;
194 double tdouble;
195
196 if (*status > 0) /* inherit input status value if > 0 */
197 return(*status);
198
199 /* generate current date string and construct the keyword comments */
200 ffgstm(datestr, NULL, status);
201 strcpy(chkcomm, "HDU checksum updated ");
202 strcat(chkcomm, datestr);
203 strcpy(datacomm, "data unit checksum updated ");
204 strcat(datacomm, datestr);
205
206 /* write the CHECKSUM keyword if it does not exist */
207 tstatus = *status;
208 if (ffgkys(fptr, "CHECKSUM", checksum, comm, status) == KEY_NO_EXIST)
209 {
210 *status = tstatus;
211 strcpy(checksum, "0000000000000000");
212 ffpkys(fptr, "CHECKSUM", checksum, chkcomm, status);
213 }
214
215 /* write the DATASUM keyword if it does not exist */
216 tstatus = *status;
217 if (ffgkys(fptr, "DATASUM", datasum, comm, status) == KEY_NO_EXIST)
218 {
219 *status = tstatus;
220 olddsum = 0;
221 ffpkys(fptr, "DATASUM", " 0", datacomm, status);
222
223 /* set the CHECKSUM keyword as undefined, if it isn't already */
224 if (strcmp(checksum, "0000000000000000") )
225 {
226 strcpy(checksum, "0000000000000000");
227 ffmkys(fptr, "CHECKSUM", checksum, chkcomm, status);
228 }
229 }
230 else
231 {
232 /* decode the datasum into an unsigned long variable */
233
234 /* olddsum = strtoul(datasum, 0, 10); doesn't work on SUN OS */
235
236 tdouble = atof(datasum);
237 olddsum = (unsigned long) tdouble;
238 }
239
240 /* close header: rewrite END keyword and following blank fill */
241 /* and re-read the required keywords to determine the structure */
242 if (ffrdef(fptr, status) > 0)
243 return(*status);
244
245 if ((fptr->Fptr)->heapsize > 0)
246 ffuptf(fptr, status); /* update the variable length TFORM values */
247
248 /* write the correct data fill values, if they are not already correct */
249 if (ffpdfl(fptr, status) > 0)
250 return(*status);
251
252 /* calc size of data unit, in FITS 2880-byte blocks */
253 if (ffghadll(fptr, &headstart, &datastart, &dataend, status) > 0)
254 return(*status);
255
256 nrec = (long) ((dataend - datastart) / 2880);
257 dsum = 0;
258
259 if (nrec > 0)
260 {
261 /* accumulate the 32-bit 1's complement checksum */
262 ffmbyt(fptr, datastart, REPORT_EOF, status);
263 if (ffcsum(fptr, nrec, &dsum, status) > 0)
264 return(*status);
265 }
266
267 if (dsum != olddsum)
268 {
269 /* update the DATASUM keyword with the correct value */
270 snprintf(datasum, FLEN_VALUE, "%lu", dsum);
271 ffmkys(fptr, "DATASUM", datasum, datacomm, status);
272
273 /* set the CHECKSUM keyword as undefined, if it isn't already */
274 if (strcmp(checksum, "0000000000000000") )
275 {
276 strcpy(checksum, "0000000000000000");
277 ffmkys(fptr, "CHECKSUM", checksum, chkcomm, status);
278 }
279 }
280
281 if (strcmp(checksum, "0000000000000000") )
282 {
283 /* check if CHECKSUM is still OK; move to the start of the header */
284 ffmbyt(fptr, headstart, REPORT_EOF, status);
285
286 /* accumulate the header checksum into the previous data checksum */
287 nrec = (long) ((datastart - headstart) / 2880);
288 sum = dsum;
289 if (ffcsum(fptr, nrec, &sum, status) > 0)
290 return(*status);
291
292 if (sum == 0 || sum == 0xFFFFFFFF)
293 return(*status); /* CHECKSUM is correct */
294
295 /* Zero the CHECKSUM and recompute the new value */
296 ffmkys(fptr, "CHECKSUM", "0000000000000000", chkcomm, status);
297 }
298
299 /* move to the start of the header */
300 ffmbyt(fptr, headstart, REPORT_EOF, status);
301
302 /* accumulate the header checksum into the previous data checksum */
303 nrec = (long) ((datastart - headstart) / 2880);
304 sum = dsum;
305 if (ffcsum(fptr, nrec, &sum, status) > 0)
306 return(*status);
307
308 /* encode the COMPLEMENT of the checksum into a 16-character string */
309 ffesum(sum, TRUE, checksum);
310
311 /* update the CHECKSUM keyword value with the new string */
312 ffmkys(fptr, "CHECKSUM", checksum, "&", status);
313
314 return(*status);
315 }
316 /*------------------------------------------------------------------------*/
ffupck(fitsfile * fptr,int * status)317 int ffupck(fitsfile *fptr, /* I - FITS file pointer */
318 int *status) /* IO - error status */
319 /*
320 Update the CHECKSUM keyword value. This assumes that the DATASUM
321 keyword exists and has the correct value.
322 */
323 {
324 char datestr[20], chkcomm[FLEN_COMMENT], comm[FLEN_COMMENT];
325 char checksum[FLEN_VALUE], datasum[FLEN_VALUE];
326 int tstatus;
327 long nrec;
328 LONGLONG headstart, datastart, dataend;
329 unsigned long sum, dsum;
330 double tdouble;
331
332 if (*status > 0) /* inherit input status value if > 0 */
333 return(*status);
334
335 /* generate current date string and construct the keyword comments */
336 ffgstm(datestr, NULL, status);
337 strcpy(chkcomm, "HDU checksum updated ");
338 strcat(chkcomm, datestr);
339
340 /* get the DATASUM keyword and convert it to a unsigned long */
341 if (ffgkys(fptr, "DATASUM", datasum, comm, status) == KEY_NO_EXIST)
342 {
343 ffpmsg("DATASUM keyword not found (ffupck");
344 return(*status);
345 }
346
347 tdouble = atof(datasum); /* read as a double as a workaround */
348 dsum = (unsigned long) tdouble;
349
350 /* get size of the HDU */
351 if (ffghadll(fptr, &headstart, &datastart, &dataend, status) > 0)
352 return(*status);
353
354 /* get the checksum keyword, if it exists */
355 tstatus = *status;
356 if (ffgkys(fptr, "CHECKSUM", checksum, comm, status) == KEY_NO_EXIST)
357 {
358 *status = tstatus;
359 strcpy(checksum, "0000000000000000");
360 ffpkys(fptr, "CHECKSUM", checksum, chkcomm, status);
361 }
362 else
363 {
364 /* check if CHECKSUM is still OK */
365 /* rewrite END keyword and following blank fill */
366 if (ffwend(fptr, status) > 0)
367 return(*status);
368
369 /* move to the start of the header */
370 ffmbyt(fptr, headstart, REPORT_EOF, status);
371
372 /* accumulate the header checksum into the previous data checksum */
373 nrec = (long) ((datastart - headstart) / 2880);
374 sum = dsum;
375 if (ffcsum(fptr, nrec, &sum, status) > 0)
376 return(*status);
377
378 if (sum == 0 || sum == 0xFFFFFFFF)
379 return(*status); /* CHECKSUM is already correct */
380
381 /* Zero the CHECKSUM and recompute the new value */
382 ffmkys(fptr, "CHECKSUM", "0000000000000000", chkcomm, status);
383 }
384
385 /* move to the start of the header */
386 ffmbyt(fptr, headstart, REPORT_EOF, status);
387
388 /* accumulate the header checksum into the previous data checksum */
389 nrec = (long) ((datastart - headstart) / 2880);
390 sum = dsum;
391 if (ffcsum(fptr, nrec, &sum, status) > 0)
392 return(*status);
393
394 /* encode the COMPLEMENT of the checksum into a 16-character string */
395 ffesum(sum, TRUE, checksum);
396
397 /* update the CHECKSUM keyword value with the new string */
398 ffmkys(fptr, "CHECKSUM", checksum, "&", status);
399
400 return(*status);
401 }
402 /*------------------------------------------------------------------------*/
ffvcks(fitsfile * fptr,int * datastatus,int * hdustatus,int * status)403 int ffvcks(fitsfile *fptr, /* I - FITS file pointer */
404 int *datastatus, /* O - data checksum status */
405 int *hdustatus, /* O - hdu checksum status */
406 /* 1 verification is correct */
407 /* 0 checksum keyword is not present */
408 /* -1 verification not correct */
409 int *status) /* IO - error status */
410 /*
411 Verify the HDU by comparing the value of the computed checksums against
412 the values of the DATASUM and CHECKSUM keywords if they are present.
413 */
414 {
415 int tstatus;
416 double tdouble;
417 unsigned long datasum, hdusum, olddatasum;
418 char chksum[FLEN_VALUE], comm[FLEN_COMMENT];
419
420 if (*status > 0) /* inherit input status value if > 0 */
421 return(*status);
422
423 *datastatus = -1;
424 *hdustatus = -1;
425
426 tstatus = *status;
427 if (ffgkys(fptr, "CHECKSUM", chksum, comm, status) == KEY_NO_EXIST)
428 {
429 *hdustatus = 0; /* CHECKSUM keyword does not exist */
430 *status = tstatus;
431 }
432 if (chksum[0] == '\0')
433 *hdustatus = 0; /* all blank checksum means it is undefined */
434
435 if (ffgkys(fptr, "DATASUM", chksum, comm, status) == KEY_NO_EXIST)
436 {
437 *datastatus = 0; /* DATASUM keyword does not exist */
438 *status = tstatus;
439 }
440 if (chksum[0] == '\0')
441 *datastatus = 0; /* all blank checksum means it is undefined */
442
443 if ( *status > 0 || (!(*hdustatus) && !(*datastatus)) )
444 return(*status); /* return if neither keywords exist */
445
446 /* convert string to unsigned long */
447
448 /* olddatasum = strtoul(chksum, 0, 10); doesn't work w/ gcc on SUN OS */
449 /* sscanf(chksum, "%u", &olddatasum); doesn't work w/ cc on VAX/VMS */
450
451 tdouble = atof(chksum); /* read as a double as a workaround */
452 olddatasum = (unsigned long) tdouble;
453
454 /* calculate the data checksum and the HDU checksum */
455 if (ffgcks(fptr, &datasum, &hdusum, status) > 0)
456 return(*status);
457
458 if (*datastatus)
459 if (datasum == olddatasum)
460 *datastatus = 1;
461
462 if (*hdustatus)
463 if (hdusum == 0 || hdusum == 0xFFFFFFFF)
464 *hdustatus = 1;
465
466 return(*status);
467 }
468 /*------------------------------------------------------------------------*/
ffgcks(fitsfile * fptr,unsigned long * datasum,unsigned long * hdusum,int * status)469 int ffgcks(fitsfile *fptr, /* I - FITS file pointer */
470 unsigned long *datasum, /* O - data checksum */
471 unsigned long *hdusum, /* O - hdu checksum */
472 int *status) /* IO - error status */
473
474 /* calculate the checksums of the data unit and the total HDU */
475 {
476 long nrec;
477 LONGLONG headstart, datastart, dataend;
478
479 if (*status > 0) /* inherit input status value if > 0 */
480 return(*status);
481
482 /* get size of the HDU */
483 if (ffghadll(fptr, &headstart, &datastart, &dataend, status) > 0)
484 return(*status);
485
486 nrec = (long) ((dataend - datastart) / 2880);
487
488 *datasum = 0;
489
490 if (nrec > 0)
491 {
492 /* accumulate the 32-bit 1's complement checksum */
493 ffmbyt(fptr, datastart, REPORT_EOF, status);
494 if (ffcsum(fptr, nrec, datasum, status) > 0)
495 return(*status);
496 }
497
498 /* move to the start of the header and calc. size of header */
499 ffmbyt(fptr, headstart, REPORT_EOF, status);
500 nrec = (long) ((datastart - headstart) / 2880);
501
502 /* accumulate the header checksum into the previous data checksum */
503 *hdusum = *datasum;
504 ffcsum(fptr, nrec, hdusum, status);
505
506 return(*status);
507 }
508
509