binaries/src/fasta34/mmgetaa.c

   1 /* mmgetaa.c - functions for mmap()ed access to libraries */
   2
   3 /* copyright (c) 1999,2000 William R. Pearson */
   4
   5 /* version 0 September, 1999 */
   6
   7 /*
   8   This is one of two alternative files that can be used to
   9   read a database.  The two files are nmgetaa.c, and mmgetaa.c
  10   (nxgetaa.c has been retired).
  11
  12   nmgetlib.c and mmgetaa.c are used together. nmgetlib.c provides
  13   the same functions as nxgetaa.c if memory mapping is not used,
  14   mmgetaa.c provides the database reading functions if memory
  15   mapping is used. The decision to use memory mapping is made on
  16   a file-by-file basis.
  17 */
  18
  19 /* $Name: fa_34_26_5 $ - $Id: mmgetaa.c,v 1.41 2006/04/12 18:00:02 wrp Exp $ */
  20
  21 #include <stdio.h>
  22 #include <stdlib.h>
  23 #include <unistd.h>
  24 #include <string.h>
  25 #include <errno.h>
  26
  27 #include <sys/types.h>
  28 #include <sys/stat.h>
  29 #include <sys/mman.h>
  30 #include <fcntl.h>
  31
  32 #define MAXLINE 512
  33 #define EOSEQ 0
  34
  35 #define XTERNAL
  36 #include "uascii.h"
  37 /* #include "upam.h" */
  38 #undef XTERNAL
  39
  40 #ifdef SUPERFAMNUM
  41 extern int nsfnum;      /* number of superfamily numbers */
  42 extern int sfnum[10];   /* superfamily number from types 0 and 5 */
  43 extern int nsfnum_n;
  44 extern int sfnum_n[10];
  45 static char tline[MAXLINE];
  46 #endif
  47
  48 #define GCGBIN 6
  49
  50 #ifndef MAP_FILE
  51 #define MAP_FILE 0
  52 #endif
  53
  54 #include "defs.h"
  55 #include "mm_file.h"
  56
  57 extern MM_OFF bl2_long8_cvt(int64_t);
  58 extern int bl2_uint4_cvt(int);
  59
  60
  61 long crck(char *, int);
  62 extern void src_int4_read(FILE *fd,  int *val);
  63 extern void src_long4_read(FILE *fd,  long  *valp);
  64 extern void src_long8_read(FILE *fd,  int64_t *val);
  65
  66 /* load_mmap() loads the d_pos[] and s_pos[] arrays for rapid access */
  67
  68 struct lmf_str *
  69 load_mmap(FILE *libi,   /* fd for already open ".xin" file */
  70           char *sname,  /* name of sequence database file */
  71           int lib_type, /* 0-Fasta, 5-vms_pir, 6-gcg_binary */
  72           int ldnaseq,  /* 1 for DNA, 0 for protein */
  73           struct lmf_str *m_fd)
  74 {
  75   char format[4];
  76   int i, lib_aa;
  77   MM_OFF f_size;
  78   long lf_size;
  79   struct stat statbuf;
  80   int max_cnt;
  81   MM_OFF *d_pos_arr, *s_pos_arr;
  82   int mm_flag, mm64_flag;
  83   int *tmp_pos_arr;
  84
  85   /* first check that the necessary indices are up-to-date */
  86   /* read the offsets in ".xin" file */
  87   if (fread(format,1,4,libi)==0) {
  88     fprintf(stderr," cannot read .xin format\n");
  89     return NULL;
  90   }
  91
  92   mm64_flag = (format[2]==1);   /* 4 bytes or 8 bytes for long? */
  93
  94 #ifndef BIG_LIB64
  95   if (mm64_flag) {return NULL;}
  96 #endif
  97
  98   if (format[3]!=lib_type) {
  99     fprintf(stderr," cannot read format %d != lib_type %d\n",
 100             format[3],lib_type);
 101     return NULL;
 102   }
 103
 104   src_int4_read(libi,&lib_aa);
 105   if (lib_aa == ldnaseq) { /* database residue mismatch */
 106     fprintf(stderr," residue type mismatch %s != %s (.xin) in %s\n",
 107             (lib_aa ? "DNA" : "prot."),(ldnaseq ? "prot." : "DNA"),
 108             sname);
 109     return NULL;
 110   }
 111
 112   /* everything looks good, allocate an lmf_str */
 113
 114   m_fd->lib_aa = lib_aa;
 115
 116   /* get get file size from index */
 117   if (mm64_flag) src_long8_read(libi,&f_size);
 118   else {
 119     src_long4_read(libi,&lf_size);
 120     f_size = lf_size;
 121   }
 122
 123   /* now, start to open mmap()ed file */
 124   mm_flag=((m_fd->mmap_fd=open(sname,O_RDONLY))>=0);
 125   if (!mm_flag) {
 126     fprintf(stderr," cannot open %s for mmap()", sname);
 127     perror("...");
 128     return NULL;        /* file did not open */
 129   }
 130
 131   /* fstat the library file and get size */
 132   if(fstat(m_fd->mmap_fd, &statbuf) < 0) {
 133     fprintf(stderr," cannot stat %s for mmap()", sname);
 134     perror("...");
 135     m_fd->mm_flg = 0;
 136     goto finish;
 137   }
 138
 139   /* check for identical sizes - if different, do not mmap */
 140   if (f_size != statbuf.st_size) {
 141     fprintf(stderr," %s file size (%lld) and expected size (%ld) don't match\n",
 142             sname,statbuf.st_size,f_size);
 143     mm_flag = 0;
 144     goto finish;
 145   }
 146
 147   /* the index file and library file are open and the sizes match */
 148   /* allocate the m_file struct and  map the file */
 149
 150   m_fd->st_size = statbuf.st_size;
 151   if((m_fd->mmap_base =
 152       mmap(NULL, m_fd->st_size, PROT_READ,
 153            MAP_FILE | MAP_SHARED, m_fd->mmap_fd, 0)) == (char *) -1) {
 154     mm_flag = 0;
 155 #ifdef DEBUG
 156     fprintf(stderr," cannot mmap %s", sname);
 157     perror("...");
 158 #endif
 159   }
 160  finish:
 161   close(m_fd->mmap_fd);
 162   if (!mm_flag) { return NULL; }
 163
 164   /* now finish reading the index file */
 165   src_int4_read(libi,&max_cnt);
 166
 167   if (mm64_flag) {
 168     src_long8_read(libi,&m_fd->tot_len);
 169   }
 170   else {
 171     src_long4_read(libi,&lf_size);
 172     m_fd->tot_len = lf_size;
 173   }
 174   src_long4_read(libi,&lf_size);
 175   m_fd->max_len = lf_size;
 176
 177 #ifdef DEBUG
 178   fprintf(stderr,
 179           "%s\tformat: %c%c%d %d; max_cnt: %d; tot_len: %lld max_len: %ld\n",
 180           sname,format[0],format[1],format[2],format[3],
 181           max_cnt,m_fd->tot_len,m_fd->max_len);
 182 #endif
 183
 184   /* allocate array of description pointers */
 185   if (!mm64_flag) {
 186     if ((tmp_pos_arr=(int *)calloc(max_cnt+1,sizeof(int)))==NULL) {
 187       fprintf(stderr," cannot allocate %d for tmp_pos array\n",
 188               max_cnt+1);
 189     }
 190   }
 191
 192   if ((d_pos_arr=(MM_OFF *)calloc(max_cnt+1, sizeof(MM_OFF)))==NULL) {
 193     fprintf(stderr," cannot allocate %d for desc. array\n",max_cnt+1);
 194     exit(1);
 195   }
 196
 197   /* read m_fd->d_pos[max_cnt+1] */
 198   if (mm64_flag) {
 199     if (fread(d_pos_arr,sizeof(MM_OFF),max_cnt+1,libi)!=
 200         max_cnt+1) {
 201       fprintf(stderr," error reading desc. offsets: %s\n",sname);
 202       return NULL;
 203     }
 204   }
 205   else {
 206     if (fread(tmp_pos_arr,sizeof(int),max_cnt+1,libi)!=
 207         max_cnt+1) {
 208       fprintf(stderr," error reading desc. offsets: %s\n",sname);
 209       return NULL;
 210     }
 211 #ifdef DEBUG
 212     fprintf(stderr,"d_pos_crc: %ld\n",
 213             crck((char *)tmp_pos_arr,sizeof(int)*(max_cnt+1)));
 214 #endif
 215   }
 216
 217
 218 #ifndef IS_BIG_ENDIAN
 219   if (mm64_flag)
 220     for (i=0; i<=max_cnt; i++) {
 221       d_pos_arr[i] = bl2_long8_cvt(d_pos_arr[i]);
 222     }
 223   else
 224     for (i=0; i<=max_cnt; i++) {
 225       d_pos_arr[i] = bl2_uint4_cvt(tmp_pos_arr[i]);
 226     }
 227 #else
 228   if (!mm64_flag) {
 229     for (i=0; i<=max_cnt; i++) {
 230       d_pos_arr[i] = tmp_pos_arr[i];
 231     }
 232   }
 233 #endif
 234
 235 #ifdef DEBUG
 236   for (i=0; i<max_cnt-1; i++) {
 237     if (d_pos_arr[i+1] <= d_pos_arr[i] )
 238       fprintf(stderr," ** dpos_error [%d]\t%ld\t%ld\n",
 239               i,d_pos_arr[i],d_pos_arr[i+1]);
 240   }
 241 #endif
 242
 243   /* allocate array of sequence pointers */
 244   if ((s_pos_arr=(MM_OFF *)calloc(max_cnt+1,sizeof(MM_OFF)))==NULL) {
 245     fprintf(stderr," cannot allocate %d for seq. array\n",max_cnt+1);
 246     exit(1);
 247   }
 248
 249   /* read m_fd->s_pos[max_cnt+1] */
 250   if (mm64_flag) {
 251     if (fread(s_pos_arr,sizeof(long),max_cnt+1,libi)!=
 252         max_cnt+1) {
 253       fprintf(stderr," error reading seq offsets: %s\n",sname);
 254       return NULL;
 255     }
 256   }
 257   else {
 258     if (fread(tmp_pos_arr,sizeof(int),max_cnt+1,libi)!=
 259         max_cnt+1) {
 260       fprintf(stderr," error reading seq offsets: %s\n",sname);
 261       return NULL;
 262     }
 263 #ifdef DEBUG
 264     fprintf(stderr,"s_pos_crc: %ld\n",
 265             crck((char *)tmp_pos_arr,sizeof(int)*(max_cnt+1)));
 266 #endif
 267   }
 268
 269 #ifndef IS_BIG_ENDIAN
 270   if (mm64_flag)
 271     for (i=0; i<=max_cnt; i++)
 272       s_pos_arr[i] = bl2_long8_cvt(s_pos_arr[i]);
 273   else
 274     for (i=0; i<=max_cnt; i++)
 275       s_pos_arr[i] = (long)bl2_uint4_cvt(tmp_pos_arr[i]);
 276 #else
 277   if (!mm64_flag)
 278     for (i=0; i<=max_cnt; i++)
 279       s_pos_arr[i] = (long)tmp_pos_arr[i];
 280 #endif
 281
 282 #ifdef DEBUG
 283   for (i=1; i<max_cnt-1; i++) {
 284     if (s_pos_arr[i+1]<s_pos_arr[i])
 285       fprintf(stderr," ** spos_error [%d]\t%ld\t%ld\n",
 286               i,s_pos_arr[i],s_pos_arr[i]);
 287   }
 288 #endif
 289
 290   if (!mm64_flag) free(tmp_pos_arr);
 291
 292   m_fd->max_cnt = max_cnt;
 293   m_fd->d_pos_arr = d_pos_arr;
 294   m_fd->s_pos_arr = s_pos_arr;
 295   m_fd->lpos = 0;
 296
 297   /*   check_mmap(m_fd,-2); */
 298
 299   return m_fd;
 300 }
 301
 302 char *mgets (char *s, int n, struct lmf_str *m_fd)
 303 {
 304   char *cs, *mfp;
 305
 306   mfp = m_fd->mmap_addr;
 307   cs = s;
 308
 309   while (--n > 0 && (*mfp != (char)EOF))
 310     if ((*cs++ = *mfp++) == '\n') break;
 311   *cs = '\0';
 312
 313   m_fd->mmap_addr = mfp;
 314   return (*mfp == (char)EOF && cs == s) ? NULL : s;
 315 }
 316
 317 int
 318 agetlibm(unsigned char *seq,
 319          int maxs,
 320          char *libstr,
 321          int n_libstr,
 322          fseek_t *libpos,
 323          int *lcont,
 324          struct lmf_str *m_fd,
 325          long *l_off)
 326 {
 327   register unsigned char *cp, *seqp;
 328   register int *ap;
 329   char *desc;
 330   int lpos;             /* entry number in library */
 331   long l;
 332   unsigned char *seqm, *seqm1;
 333   char *bp;
 334   static long seq_len;
 335   static unsigned char *cp_max;
 336 #ifdef SUPERFAMNUM
 337   char *bp1, *bpa, *tp;
 338   int i;
 339 #endif
 340
 341   *l_off = 1;
 342
 343   lpos = m_fd->lpos;
 344
 345   seqp = seq;
 346   seqm = &seq[maxs-9];
 347   seqm1 = seqm-1;
 348
 349   ap = m_fd->sascii;
 350
 351   if (*lcont==0) {
 352     if (lpos >= m_fd->max_cnt) return (-1);
 353     seq_len = m_fd->d_pos_arr[lpos+1] - m_fd->s_pos_arr[lpos];
 354     if (seq_len < 0 || (seq_len > m_fd->max_len && seq_len > (m_fd->max_len*5)/4)) {
 355       fprintf(stderr," ** sequence over-run: %ld at %d\n",seq_len,lpos);
 356       return(-1);
 357     }
 358     *libpos = (fseek_t)lpos;
 359
 360     desc = m_fd->mmap_base+m_fd->d_pos_arr[lpos]+1;
 361     strncpy(libstr,desc,n_libstr-1);
 362     libstr[n_libstr-1]='\0';
 363     if ((bp=strchr(libstr,'\r'))!=NULL) *bp='\0';
 364     if ((bp=strchr(libstr,'\n'))!=NULL) *bp='\0';
 365     if (n_libstr > MAX_UID) {
 366       bp = libstr;
 367       while (*bp++) if ( *bp=='\001' || *bp=='\t') *bp=' ';
 368     }
 369
 370     for (bp = desc; *bp && (*bp != '\n'); *bp++ )
 371       if (*bp == '@' && !strncmp(bp+1,"C:",2)) sscanf(bp+3,"%ld",l_off);
 372
 373 #ifdef SUPERFAMNUM
 374     sfnum[0]=nsfnum=0;
 375     strncpy(tline,desc,sizeof(tline));
 376     tline[MAXLINE-1]='\0';
 377     if ((bp=strchr(tline,'\n'))!=NULL) *bp='\0';
 378     if ((bp=strchr(tline,' ')) && (bp=strchr(bp+1,SFCHAR))) {
 379       if ((bpa = strchr(bp+1,'\001'))!=NULL) *bpa = '\0';
 380       if ((bp1=strchr(bp+1,SFCHAR))==NULL) {
 381         fprintf(stderr," second %c missing: %s\n",SFCHAR,tline);
 382       }
 383       else {
 384         *bp1 = '\0';
 385         i = 0;
 386         if ((tp = strtok(bp+1," \t"))!=NULL) {
 387           sfnum[i++] = atoi(tp);
 388           while ((tp = strtok((char *)NULL," \t")) != (char *)NULL) {
 389             sfnum[i++] = atoi(tp);
 390             if (i>=9) break;
 391           }
 392         }
 393         sfnum[nsfnum=i]= 0;
 394         if (nsfnum>1) sf_sort(sfnum,nsfnum);
 395         else {
 396           if (nsfnum<1) fprintf(stderr," found | but no sfnum: %s\n",libstr);
 397         }
 398       }
 399     }
 400     else {
 401       sfnum[0] = nsfnum = 0;
 402       }
 403 #endif
 404
 405     m_fd->mmap_addr = m_fd->mmap_base+m_fd->s_pos_arr[lpos];
 406     cp_max = (unsigned char *)(m_fd->mmap_addr+seq_len);
 407   }
 408
 409   for (cp=(unsigned char *)m_fd->mmap_addr; seqp<seqm1; ) {
 410     if ((*seqp++=ap[*cp++])<NA &&
 411         (*seqp++=ap[*cp++])<NA &&
 412         (*seqp++=ap[*cp++])<NA &&
 413         (*seqp++=ap[*cp++])<NA &&
 414         (*seqp++=ap[*cp++])<NA &&
 415         (*seqp++=ap[*cp++])<NA &&
 416         (*seqp++=ap[*cp++])<NA &&
 417         (*seqp++=ap[*cp++])<NA &&
 418         (*seqp++=ap[*cp++])<NA &&
 419         (*seqp++=ap[*cp++])<NA) continue;
 420     --seqp;
 421     if (cp >= cp_max) break;
 422   }
 423   m_fd->mmap_addr = (char *)cp;
 424
 425   if (seqp>=seqm1) (*lcont)++;
 426   else {
 427     *lcont=0;
 428     lpos++;
 429     m_fd->lpos = lpos;
 430   }
 431   *seqp = EOSEQ;
 432   /*   if ((int)(seqp-seq)==0) return 1; */
 433   return (int)(seqp-seq);
 434 }
 435
 436 void
 437 aranlibm(char *str,
 438          int cnt,
 439          fseek_t libpos,
 440          char *libstr,
 441          struct lmf_str *m_fd)
 442 {
 443   char *bp;
 444   int llen;
 445   int lpos;
 446
 447   lpos = (int) libpos;
 448
 449   llen = m_fd->s_pos_arr[lpos]-m_fd->d_pos_arr[lpos];
 450   if (llen >= cnt) llen = cnt-1;
 451
 452   strncpy(str,m_fd->mmap_base+m_fd->d_pos_arr[lpos]+1,llen);
 453   str[llen]='\0';
 454   if ((bp = strchr(str,'\r'))!=NULL) *bp='\0';
 455   if ((bp = strchr(str,'\n'))!=NULL) *bp='\0';
 456   bp = str;
 457   while (*bp++) if ( *bp=='\001' || *bp=='\t') *bp=' ';
 458   m_fd->lpos = lpos;
 459 }
 460
 461 /* there is no vgetlibm() because vgetlibm() and agetlibm() are
 462    identical - the difference in the two file formats relates to the
 463    location of the sequence, which is already available in spos_arr[].
 464
 465    however vranlibm must accomodate both type 5 and 6 files;
 466    type 6 has extra stuff after the seq_id.
 467 */
 468
 469 void
 470 vranlibm(char *str,
 471          int cnt,
 472          fseek_t libpos,
 473          char *libstr,
 474          struct lmf_str *m_fd)
 475 {
 476   char *bp, *mp;
 477   int llen;
 478   int lpos;
 479
 480   lpos = (int)libpos;
 481
 482   llen = m_fd->s_pos_arr[lpos]-m_fd->d_pos_arr[lpos];
 483
 484   mp = m_fd->mmap_base+m_fd->d_pos_arr[lpos];
 485
 486   strncpy(str,mp+4,20);
 487   str[20]='\0';
 488   if ((bp=strchr(str,' '))!=NULL) *(bp+1) = '\0';
 489   else if ((bp=strchr(str,'\n'))!=NULL) *bp = ' ';
 490   bp = strchr(mp,'\n');
 491
 492   llen -= (bp-mp)-5;
 493   if (llen >  cnt-strlen(str)) llen = cnt-strlen(str)-1;
 494
 495   strncat(str,bp+1,llen);
 496   if ((bp = strchr(str,'\n'))!=NULL) *bp='\0';
 497   str[cnt-1]='\0';
 498   m_fd->lpos = lpos;
 499 }
 500
 501 void
 502 close_mmap(struct lmf_str *m_fd) {
 503   free(m_fd->s_pos_arr);
 504   free(m_fd->d_pos_arr);
 505   if (m_fd->mm_flg) {
 506     munmap(m_fd->mmap_base,m_fd->st_size);
 507     free(m_fd);
 508   }
 509   m_fd->mm_flg=0;
 510 }
 511
 512 #ifndef min
 513 #define min(x,y) ((x) > (y) ? (y) : (x))
 514 #endif
 515
 516 static int gcg_bton[4]={2,4,1,3};
 517
 518 int
 519 gcg_getlibm(unsigned char *seq,
 520             int maxs,
 521             char *libstr,
 522             int n_libstr,
 523             fseek_t *libpos,
 524             int *lcont,
 525             struct lmf_str *m_fd,
 526             long *l_off)
 527 {
 528   char dummy[20];
 529   char gcg_date[6];
 530   char gcg_type[10];
 531   register unsigned char *cp, *seqp, stmp;
 532   register int *ap, lpos;
 533   unsigned char *seqm, *seqm1;
 534   long r_block, b_block, r_fact, r16_block;
 535
 536   *l_off = 1;
 537
 538   seqp = seq;
 539   seqm = &seq[maxs-9];
 540   seqm1 = seqm-1;
 541
 542   ap = m_fd->sascii;
 543   lpos = m_fd->lpos;
 544
 545   if (*lcont==0) {
 546     if (lpos >= m_fd->max_cnt) return (-1);
 547     sscanf(m_fd->mmap_base+m_fd->d_pos_arr[lpos]+4,"%s %s %s %s %ld\n",
 548            libstr,gcg_date,gcg_type,dummy,&(m_fd->gcg_len));
 549
 550     m_fd->gcg_binary = (gcg_type[0]=='2');
 551
 552     libstr[12]='\0';
 553     *libpos = lpos;
 554     m_fd->mmap_addr = m_fd->mmap_base+m_fd->s_pos_arr[lpos];
 555   }
 556
 557   r_block = b_block = min((size_t)(seqm-seqp),m_fd->gcg_len);
 558   if (m_fd->gcg_binary) {
 559     r_block = (r_block+3)/4;
 560   }
 561
 562   cp=(unsigned char *)m_fd->mmap_addr;
 563   if (!m_fd->gcg_binary) {
 564     r_fact = 1;
 565     r16_block = r_block/16;
 566     while (r16_block-- > 0) {
 567       *seqp++ = ap[*cp++];
 568       *seqp++ = ap[*cp++];
 569       *seqp++ = ap[*cp++];
 570       *seqp++ = ap[*cp++];
 571       *seqp++ = ap[*cp++];
 572       *seqp++ = ap[*cp++];
 573       *seqp++ = ap[*cp++];
 574       *seqp++ = ap[*cp++];
 575       *seqp++ = ap[*cp++];
 576       *seqp++ = ap[*cp++];
 577       *seqp++ = ap[*cp++];
 578       *seqp++ = ap[*cp++];
 579       *seqp++ = ap[*cp++];
 580       *seqp++ = ap[*cp++];
 581       *seqp++ = ap[*cp++];
 582       *seqp++ = ap[*cp++];
 583     }
 584     while (seqp<seq+r_block) *seqp++ = ap[*cp++];
 585   }
 586   else if (m_fd->gcg_binary) {
 587     r_fact = 4;
 588     r16_block = r_block/8;
 589     while(r16_block-- > 0) {
 590       stmp = *cp++;
 591       *seqp++ = gcg_bton[(stmp>>6) &3];
 592       *seqp++ = gcg_bton[(stmp>>4) &3];
 593       *seqp++ = gcg_bton[(stmp>>2) &3];
 594       *seqp++ = gcg_bton[(stmp) &3];
 595       stmp = *cp++;
 596       *seqp++ = gcg_bton[(stmp>>6) &3];
 597       *seqp++ = gcg_bton[(stmp>>4) &3];
 598       *seqp++ = gcg_bton[(stmp>>2) &3];
 599       *seqp++ = gcg_bton[(stmp) &3];
 600       stmp = *cp++;
 601       *seqp++ = gcg_bton[(stmp>>6) &3];
 602       *seqp++ = gcg_bton[(stmp>>4) &3];
 603       *seqp++ = gcg_bton[(stmp>>2) &3];
 604       *seqp++ = gcg_bton[(stmp) &3];
 605       stmp = *cp++;
 606       *seqp++ = gcg_bton[(stmp>>6) &3];
 607       *seqp++ = gcg_bton[(stmp>>4) &3];
 608       *seqp++ = gcg_bton[(stmp>>2) &3];
 609       *seqp++ = gcg_bton[(stmp) &3];
 610       stmp = *cp++;
 611       *seqp++ = gcg_bton[(stmp>>6) &3];
 612       *seqp++ = gcg_bton[(stmp>>4) &3];
 613       *seqp++ = gcg_bton[(stmp>>2) &3];
 614       *seqp++ = gcg_bton[(stmp) &3];
 615       stmp = *cp++;
 616       *seqp++ = gcg_bton[(stmp>>6) &3];
 617       *seqp++ = gcg_bton[(stmp>>4) &3];
 618       *seqp++ = gcg_bton[(stmp>>2) &3];
 619       *seqp++ = gcg_bton[(stmp) &3];
 620       stmp = *cp++;
 621       *seqp++ = gcg_bton[(stmp>>6) &3];
 622       *seqp++ = gcg_bton[(stmp>>4) &3];
 623       *seqp++ = gcg_bton[(stmp>>2) &3];
 624       *seqp++ = gcg_bton[(stmp) &3];
 625       stmp = *cp++;
 626       *seqp++ = gcg_bton[(stmp>>6) &3];
 627       *seqp++ = gcg_bton[(stmp>>4) &3];
 628       *seqp++ = gcg_bton[(stmp>>2) &3];
 629       *seqp++ = gcg_bton[(stmp) &3];
 630     }
 631
 632     while (seqp < seq+4*r_block) {
 633       stmp = *cp++;
 634       *seqp++ = gcg_bton[(stmp>>6) &3];
 635       *seqp++ = gcg_bton[(stmp>>4) &3];
 636       *seqp++ = gcg_bton[(stmp>>2) &3];
 637       *seqp++ = gcg_bton[(stmp) &3];
 638     }
 639   }
 640   if (r_fact * r_block >= m_fd->gcg_len) {
 641     *lcont = 0;
 642     m_fd->lpos++;
 643   }
 644   else {
 645     if (m_fd->gcg_binary) b_block = 4*r_block;
 646     m_fd->gcg_len -= b_block;
 647     (*lcont)++;
 648   }
 649
 650   seq[b_block] = EOSEQ;
 651   /*   if (b_block==0) return 1; else */
 652   return b_block;
 653 }
 654
 655 void lget_ann_m(struct lmf_str *lm_fd, char *libstr, int n_libstr);
 656
 657 int
 658 lgetlibm(unsigned char *seq,
 659          int maxs,
 660          char *libstr,
 661          int n_libstr,
 662          fseek_t *libpos,
 663          int *lcont,
 664          struct lmf_str *m_fd,
 665          long *l_off)
 666 {
 667   register unsigned char *cp, *seqp;
 668   register int *ap, lpos;
 669   unsigned char *seqm, *seqm1;
 670
 671   *l_off = 1;
 672
 673   seqp = seq;
 674   seqm = &seq[maxs-11];
 675   seqm1 = seqm-1;
 676
 677   lpos = m_fd->lpos;
 678   ap = m_fd->sascii;
 679
 680   if (*lcont==0) {
 681     if (lpos >= m_fd->max_cnt) return (-1);
 682
 683     if (n_libstr <= 21) {
 684       strncpy(libstr,m_fd->mmap_base+m_fd->d_pos_arr[lpos]+12,12);
 685       libstr[12]='\0';
 686     }
 687     else {
 688       lget_ann_m(m_fd,libstr,n_libstr);
 689     }
 690     *libpos = lpos;
 691
 692     m_fd->mmap_addr = m_fd->mmap_base+m_fd->s_pos_arr[lpos];
 693     cp = (unsigned char *)m_fd->mmap_addr;
 694   }
 695   else cp = (unsigned char *)m_fd->mmap_addr;
 696
 697   while (seqp<seqm1) {
 698     if (*cp=='/' && *(cp-1)=='\n') break;
 699     if ((*seqp++=ap[*cp++])<NA &&
 700         (*seqp++=ap[*cp++])<NA &&
 701         (*seqp++=ap[*cp++])<NA &&
 702         (*seqp++=ap[*cp++])<NA &&
 703         (*seqp++=ap[*cp++])<NA &&
 704         (*seqp++=ap[*cp++])<NA &&
 705         (*seqp++=ap[*cp++])<NA &&
 706         (*seqp++=ap[*cp++])<NA &&
 707         (*seqp++=ap[*cp++])<NA &&
 708         (*seqp++=ap[*cp++])<NA &&
 709         (*seqp++=ap[*cp++])<NA) continue;
 710     --seqp;
 711     if (*cp=='\n' && *(cp+1)==' ') cp += 11;
 712   }
 713
 714   if (seqp>=seqm1) {
 715     (*lcont)++;
 716     m_fd->mmap_addr = (char *)cp;
 717   }
 718   else {
 719     *lcont=0;
 720     m_fd->lpos++;
 721   }
 722
 723   *seqp = EOSEQ;
 724   return (int)(seqp-seq);
 725 }
 726
 727 void
 728 lget_ann_m(struct lmf_str *lm_fd, char *libstr, int n_libstr) {
 729   char *bp, *bp_gid, locus[120], desc[120], acc[120], ver[120];
 730
 731   /* copy in locus from lm_fd->lline */
 732   strncpy(locus,&lm_fd->mmap_addr[12],sizeof(locus));
 733   if ((bp=strchr(locus,' '))!=NULL) *(bp+1) = '\0';
 734
 735   /* get description */
 736   mgets(desc,sizeof(desc),lm_fd);
 737   while (desc[0]!='D' || desc[1]!='E' || strncmp(desc,"DEFINITION",10))
 738     mgets(desc,sizeof(desc),lm_fd);
 739   if ((bp = strchr(&desc[12],'\n'))!=NULL) *bp='\0';
 740
 741   /* get accession */
 742   mgets(acc,sizeof(acc),lm_fd);
 743   while (acc[0]!='A' || acc[1]!='C' || strncmp(acc,"ACCESSION",9)) {
 744     mgets(acc,sizeof(acc),lm_fd);
 745     if (acc[0]=='O' && acc[1]=='R' && strncmp(acc,"ORIGIN",6)==0)
 746       break;
 747   }
 748   if ((bp = strchr(&acc[12],'\n'))!=NULL) *bp='\0';
 749   if ((bp = strchr(&acc[12],' '))!=NULL) *bp='\0';
 750
 751   /* get version */
 752   mgets(ver,sizeof(ver),lm_fd);
 753   while (ver[0]!='V' || ver[1]!='E' || strncmp(ver,"VERSION",7)) {
 754     mgets(ver,sizeof(ver),lm_fd);
 755     if (ver[0]=='O' && ver[1]=='R' && strncmp(ver,"ORIGIN",6)==0)
 756       break;
 757   }
 758   if ((bp = strchr(&ver[12],'\n'))!=NULL) *bp='\0';
 759
 760       /* extract gi:123456 from version line */
 761   bp_gid = strchr(&ver[12],':');
 762   if (bp_gid != NULL) {
 763     if ((bp=strchr(bp_gid+1,' '))!=NULL) *bp='\0';
 764     bp_gid++;
 765   }
 766   if ((bp = strchr(&ver[12],' '))!=NULL) *bp='\0';
 767
 768       /* build up FASTA header line */
 769   if (bp_gid != NULL) {
 770     strncpy(libstr,"gi|",n_libstr-1);
 771     strncat(libstr,bp_gid,n_libstr-4);
 772     strncat(libstr,"|gb|",n_libstr-20);
 773   }
 774   else {libstr[0]='\0';}
 775
 776   /* if we have a version number, use it, otherwise accession,
 777          otherwise locus/description */
 778
 779   if (ver[0]=='V') {
 780     strncat(libstr,&ver[12],n_libstr-1-strlen(libstr));
 781     strncat(libstr,"|",n_libstr-1-strlen(libstr));
 782   }
 783   else if (acc[0]=='A') {
 784     strncat(libstr,&acc[12],n_libstr-1-strlen(libstr));
 785     strncat(libstr," ",n_libstr-1-strlen(libstr));
 786   }
 787
 788   strncat(libstr,locus,n_libstr-1-strlen(libstr));
 789   strncat(libstr,&desc[11],n_libstr-1-strlen(libstr));
 790   libstr[n_libstr-1]='\0';
 791 }
 792
 793 void
 794 lranlibm(char *str,
 795          int cnt,
 796          fseek_t seek,
 797          char *libstr,
 798          struct lmf_str *m_fd)
 799 {
 800   char *bp, *llp;
 801   char acc[MAXLINE], desc[MAXLINE];
 802
 803   llp = m_fd->mmap_addr = m_fd->mmap_base + m_fd->d_pos_arr[seek];
 804
 805   lget_ann_m(m_fd,str,cnt);
 806
 807   str[cnt-1]='\0';
 808
 809   m_fd->lpos = seek;
 810 }
 811
 812 static int check_status=0;
 813
 814 void
 815 check_mmap(struct lmf_str *m_fd,long ntt) {
 816
 817   int i, seq_len, ok_stat;
 818
 819   ok_stat = 1;
 820   if ( ++check_status > 5) return;
 821
 822   fprintf(stderr," ** checking %s %ld**\n", m_fd->lb_name,ntt);
 823   for (i=0; i<m_fd->max_cnt; i++) {
 824     seq_len = m_fd->d_pos_arr[i+1] - m_fd->s_pos_arr[i];
 825     if (seq_len < 0 || (seq_len > m_fd->max_len && seq_len > (m_fd->max_len*5)/4)) {
 826       fprintf(stderr,"%d:\t%ld\t%ld\t%ld\n",
 827               i,m_fd->d_pos_arr[i],m_fd->s_pos_arr[i],
 828               m_fd->d_pos_arr[i+1]-m_fd->s_pos_arr[i]);
 829       ok_stat=0;
 830     }
 831   }
 832   if (ok_stat) {
 833     if (check_status) fprintf(stderr," ** check_mmap OK %s %ld**\n",
 834                               m_fd->lb_name,ntt);
 835   }
 836 }
 837
 838 #ifdef DEBUG
 839 /*  C H K 3  --  Compute a type-3 Kermit block check.  */
 840 /*
 841  Calculate the 16-bit CRC of a null-terminated string using a byte-oriented
 842  tableless algorithm invented by Andy Lowry (Columbia University).  The
 843  magic number 010201 is derived from the CRC-CCITT polynomial x^16+x^12+x^5+1.
 844  Note - this function could be adapted for strings containing imbedded 0's
 845  by including a length argument.
 846 */
 847 long
 848 crck(s,n)
 849     char *s; int n;
 850 {
 851     unsigned int c, q;
 852     long crc = 0;
 853
 854     while (n-->0) {
 855         c = *s++;
 856         /* if (parity)*/
 857         c &= 0177;
 858         q = (crc ^ c) & 017;            /* Low-order nibble */
 859         crc = (crc >> 4) ^ (q * 010201);
 860         q = (crc ^ (c >> 4)) & 017;     /* High order nibble */
 861         crc = (crc >> 4) ^ (q * 010201);
 862     }
 863     return(crc);
 864 }
 865 #endif