#include <stdlib.h>  
main () {  
    /* Are we little or big endian?  Originally from Harbison and Steele. */  
    union  
    {  
long l;  
char c[sizeof(long)];  
    } u;  
    u.l = 1;  
    exit (u.c[sizeof(long)-1] == 1);  
    /* Return 0 (success) if we’re little-endian, or 1 (fail) if big-endian */  
}

/*  
 * Explore IEEE single- and double-precision floats.  
 *  
 * Invoked with an argument, it reprints that number as a double, hex  
 * and binary.  Without argument, it reads numbers from stdin.  The  
 * argument is either a double or an integer (hex, with leading ’0x’,  
 * eg 0x3f800000) as argument.  
 *  
 * $Id$  
 */  
 
#include <config.h>  
 
#include <stdio.h>  
#include <stdlib.h>  
#include <ctype.h>  
 
#if HAVE_ERRNO_H  
#include <errno.h>  
#else  
extern int errno;  
#endif  
 
/* bytesex.h defines BIGENDIAN=1 (true) or 0 */  
#include "bytesex.h"  
 
#ifndef DOUBLEPRECISION  
#define DOUBLEPRECISION 1 /* double-precision by default */  
#endif  
 
#if DOUBLEPRECISION  
#define EXPONENTWIDTH 11  
#define MANTISSAWIDTH 20 /* actually just the part of the  
   mantissa in the MSB */  
#define ZEROMANTISSA(F) ((F).ieee.m == 0 && (F).ieee.m2 == 0)  
#define NINTS 2 /* number of ints in a double */  
typedef double Float;  
#else  
#define EXPONENTWIDTH 8  
#define MANTISSAWIDTH 23  
#define ZEROMANTISSA(F) ((F).ieee.m == 0)  
#define NINTS 1 /* number of ints in a float */  
typedef float Float;  
#endif  
/* Create an integer with EXPONENTWIDTH set bits.  Relies on  
   zero-filling when left-shifting. */  
#define FULLEXPONENT ~(~0 << EXPONENTWIDTH)  
 
/* Define the structure we’ll use to do the manipulations */  
typedef union {  
    Float f;  
    unsigned int i[NINTS];  
    struct {  
#if DOUBLEPRECISION  
#if BIGENDIAN  
unsigned int s  : 1;  
unsigned int e  : 11;  
unsigned int m  : 20;  
unsigned int m2 : 32;  
#else  
unsigned int m2 : 32;  
unsigned int m  : 20;  
unsigned int e  : 11;  
unsigned int s  : 1;  
#endif  
#else  /* DOUBLEPRECISION */  
#if BIGENDIAN  
unsigned int s : 1;  
unsigned int e : 8;  
unsigned int m : 23;  
#else  
unsigned int m : 23;  
unsigned int e : 8;  
unsigned int s : 1;  
#endif  
#endif  
    } ieee;  
} ieeefloat;  
 
/* Function prototypes */  
/* Utility functions */  
int pisnan (const Float darg); /* true if argument is a NaN */  
int pisinf (const Float darg); /* true if argument is infinite */  
int pisfinite (const Float darg); /* true if argument is finite */  
int pisnormal (const Float darg); /* true if argument is normal */  
Float pnan (const char *tagp); /* return a NaN */  
Float pinfinity (Float darg); /* return +ve or -ve infinity */  
 
/* Functions which do the work of this demo program */  
char *pbits (unsigned int i, unsigned int n);  
int parse_int (const char *s, unsigned int i[NINTS]);  
int reprint_number (char *s, Float farg);  
 
 
 
/* Following are utility functions, which are portable */  
 
/* Returns non-zero if argument is a NaN -- all bits in exponent set,  
   and mantissa non-zero. */  
int pisnan (const Float darg)  
{  
    ieeefloat d;  
    d.f = darg;  
#if DOUBLEPRECISION  
    return (d.ieee.e == FULLEXPONENT && !(d.ieee.m == 0 && d.ieee.m2 == 0));  
#else  
    return (d.ieee.e == FULLEXPONENT && d.ieee.m != 0);  
#endif  
}  
 
/* returns non-zero if argument is infinity -- all bits in exponent  
   set, and mantissa zero. */  
int pisinf (const Float darg)  
{  
    ieeefloat d;  
    d.f = darg;  
#if DOUBLEPRECISION  
    return (d.ieee.e == FULLEXPONENT && d.ieee.m == 0 && d.ieee.m2 == 0);  
#else  
    return (d.ieee.e == FULLEXPONENT && d.ieee.m == 0);  
#endif  
}  
 
/* returns non-zero if argument is finite (normal, subnormal or zero)  
   -- not all bits in exponent set. */  
int pisfinite (const Float darg)  
{  
    ieeefloat d;  
    d.f = darg;  
    return (d.ieee.e != FULLEXPONENT);  
}  
 
/* Returns non-zero if argument is normal -- exponent is not all-ones  
 * nor all-zeros.  
 *  
 * Note that this refers to the argument as a _double_, and might not  
 * give the expected results if it’s applied to a single which is  
 * promoted to a double when given as a parameter.  
 */  
int pisnormal (const Float darg)  
{  
    ieeefloat d;  
    d.f = darg;  
    return (d.ieee.e != FULLEXPONENT && d.ieee.e != 0);  
}  
 
/* returns a NaN.  Argument is ignored */  
Float pnan (const char *tagp)  
{  
    ieeefloat d;  
    d.ieee.s  = 0;  
    d.ieee.e  = FULLEXPONENT;  
    /* Set leading bit in mantissa -- quiet NaN. */  
    d.ieee.m = (1<<(MANTISSAWIDTH-1));  
#if DOUBLEPRECISION  
    d.ieee.m2 = 0;  
#endif  
    return d.f;  
}  
 
/* returns positive or negative infinity, depending on the sign of darg */  
Float pinfinity (Float darg)  
{  
    ieeefloat d;  
    d.f = darg; /* copies sign */  
    d.ieee.e = FULLEXPONENT;  
    d.ieee.m = 0;  
#if DOUBLEPRECISION  
    d.ieee.m2 = 0;  
#endif  
    return d.f;  
}  
 
/*  
 * pbits: given an integer i (assumed 32 bits),  
 * return the rightmost n (<=32) bits of the integer expressed in binary.  
 */  
char *pbits (unsigned int i, unsigned int n)  
{  
    static char s[33];  
    char *p;  
    /* printf ("pbits(%x/%d)", i, n); */  
    s[n] = ’\0’;                        /* terminate the string */  
    for (p=s+n-1; p>=s; p--, i>>=1)  
        *p = ((i&1) ? ’1’ : ’0’);  
    return s;  
}  
 
#if DOUBLEPRECISION  
/* Take a string of (up to) 16 hex digits, remove leading ’0x’, and  
 * any spaces, and convert the first 8 into i[0], and the second 8 into  
 * i[1].  If the string is shorter than 16 digits, pad the remainder  
 * with ’0’.  
 *  
 * Only need this for converting 16-digit strings to d-p floats.  
 *  
 * Return non-zero on error.  
 */  
int parse_int (const char *s, unsigned int i[2])  
{  
    char digits0[9], digits1[9], *p;  
    int ndigits = 0;  
 
    if (!(s[0] == ’0’ && s[1] == ’x’))  
return 1;  
    p = digits0;  
    s += 2;  
    while (ndigits <= 16)  
    {  
while (isspace(*s)) /* skip blanks */  
    s++;  
if (*s == ’\0’) /* pad remainder */  
    *p++ = ’0’;  
else  
    *p++ = *s++;  
ndigits++;  
if (ndigits == 8) /* done 8 digits - switch to second 8*/  
    p = digits1;  
    }  
    digits0[8] = digits1[8] = ’\0’; /* terminate the strings */  
 
    i[0] = strtoul (digits0, (char**)NULL, 16); /* ...and convert them */  
    i[1] = strtoul (digits1, (char**)NULL, 16);  
 
    /* printf ("parse_int: digits=%s %s --> %08x %08x\n",  
       digits0, digits1, i[0], i[1]); */  
    return 0;  
}  
#else  
/* same, but simpler for single precision */  
int parse_int (const char *s, unsigned int i[1])  
{  
    char digits0[9], *p;  
    int ndigits;  
 
    if (!(s[0] == ’0’ && s[1] == ’x’))  
return 1;  
    p = digits0;  
    s += 2;  
    for (ndigits=0; ndigits<=8; ndigits++)  
    {  
while (isspace(*s)) /* skip blanks */  
    s++;  
if (*s == ’\0’) /* pad remainder */  
    *p++ = ’0’;  
else  
    *p++ = *s++;  
    }  
    digits0[8] = ’\0’; /* terminate the strings */  
 
    *i = strtoul (digits0, (char**)NULL, 16); /* ...and convert them */  
 
    /* printf ("parse_int: digits=%s --> %08x\n", digits0, *i); */  
    return 0;  
}  
#endif  
 
 
/*  
 * Convert string to number and display in various formats.  
 * If s is NULL, then use the number farg instead.  If the input  
 * string cannot be interpreted, return non-zero, or zero on success.  
 */  
int reprint_number (char *s, Float farg)  
{  
    ieeefloat F;  
    if (s == NULL)  
F.f = farg;  
    else  
    {  
if (s[0]==’0’ && s[1]==’x’) { /* it’s an integer */  
    if (parse_int (s, F.i)) {  
printf ("can’t parse integer <%s>\n", s);  
return 1;  
    }  
} else {  
            char *endptr;  
            errno = 0;  
    F.f = strtod (s, &endptr);  
            if (errno != 0)  
                return 1;  
            for (; isspace(*endptr); endptr++)  
                ;  
            if (*endptr != ’\0’)  
                return 1;  
        }  
    }  
 
#if DOUBLEPRECISION  
    /* 49 bits of precision is 49*log_10(2)=14.75 dec.digits of precision */  
printf ("double  %24.16e\n   hex  %08x %08x\n",  
    F.f, F.i[0], F.i[1]);  
#else  
printf (" float  %12.7e\n   hex  %08x\n", F.f, F.i[0]);  
#endif  
    printf ("binary  %s ", pbits (F.ieee.s, 1));  
    printf ("%s ", pbits (F.ieee.e, EXPONENTWIDTH));  
    printf ("%s ", pbits (F.ieee.m, MANTISSAWIDTH));  
#if DOUBLEPRECISION  
    printf ("%s", pbits (F.ieee.m2, 32));  
#endif  
    printf ("\n  type  ");  
 
    /* now print out what type of number it is */  
    if (F.ieee.e == 0)  
        if (ZEROMANTISSA(F))  
            printf ("%s zero",  
                    (F.ieee.s ? "Negative" : "Positive"));  
        else  
            printf ("Subnormal");  
    else if (F.ieee.e == FULLEXPONENT)  
        if (ZEROMANTISSA(F))  
            printf ("%s infinity",  
                    (F.ieee.s ? "Negative" : "Positive"));  
        else  
            printf ("%s NaN",  
                    (F.ieee.m & (1<<(MANTISSAWIDTH-1))  
     ? "Quiet"  
     : "Signalling"));  
    else  
        printf ("Normal");  
 
    /* test the pis??? routines */  
    printf (" (%c%c%c%c)\n",  
    (pisnan    (F.f) ? ’N’ : ’n’),  
    (pisinf    (F.f) ? ’I’ : ’i’),  
    (pisfinite (F.f) ? ’F’ : ’f’),  
    (pisnormal (F.f) ? ’L’ : ’l’));  
 
    return 0;  
}  
 
int main (int argc, char **argv)  
{  
    Float f = 0;  
 
    printf ("%s-precision %s endian...\n",  
    (DOUBLEPRECISION ? "Double" : "Single"),  
    (BIGENDIAN ? "big" : "little"));  
 
    if (argc > 1)  
        reprint_number (argv[1], 0);  
    else  
    {  
        char line[80];  
        while (fgets (line, sizeof(line), stdin) != NULL)  
    if (line[0] == ’=’) {  
switch (line[1]) {  
  case ’n’:  
    reprint_number (NULL, pnan(""));  
    break;  
  case ’i’:  
    reprint_number (NULL, pinfinity(+1.0));  
    break;  
  case ’z’:  
    reprint_number (NULL, 0);  
    break;  
  default:  
    printf ("eh?\n");  
    break;  
}  
            } else {  
if (reprint_number (line, 0))  
                    printf("Enter float, 0x... integer, =n, =i, =z, "  
                           "or ^D to exit\n");  
            }  
    }  
    exit (0);  
}

#include <stdio.h>  
 
/* bytesex.h defines BIGENDIAN=1 (true) or 0 */  
#include "bytesex.h"  
 
/*  
 * pbits: given an integer i, put the rightmost n bits of the integer  
 * expressed in binary, into the string s.  
 */  
void pbits (unsigned int i, unsigned int n, char *s)  
{  
    char *p;  
    for (p=s+n-1; p>=s; p--, i>>=1)  
        *p = ((i&1) ? ’1’ : ’0’);  
    return;  
}  
 
/* Given a float, return a (static) string holding its bit pattern */  
char *print_float (float f)  
{  
    static char b[35];  
    union {  
        float f;  
        struct {  
#if BIGENDIAN  
            unsigned int s : 1;  
            unsigned int e : 8;  
            unsigned int m : 23;  
#else  
            unsigned int m : 23;  
            unsigned int e : 8;  
            unsigned int s : 1;  
#endif  
        } ieee;  
    } ieeefloat;  
 
    ieeefloat.f = f;  
 
    pbits (ieeefloat.ieee.s, 1,  b);  
    b[1] = ’ ’;  
    pbits (ieeefloat.ieee.e, 8,  &b[2]);  
    b[10] = ’ ’;  
    pbits (ieeefloat.ieee.m, 23, &b[11]);  
    b[34] = ’\0’;  
 
    return b;  
}  
 
 
main ()  
{  
    float f, a, b, c;  
    double fd, ad, bd, cd;  
 
    printf ("\n  You lose accuracy adding small numbers to large ones...\n");  
    /* 2^23+1 works OK... */  
    f = 8388608.0;  
    printf ("2^23 = %e = %s\n", f, print_float (f));  
    f = f + 1.0;  
    printf ("  +1 = %e = %s\n", f, print_float (f));  
    /* ...but 2^24+1=2^24 - catastrophic loss of precision */  
    f = 16777216.0;  
    printf ("2^24 = %e = %s\n", f, print_float (f));  
    f = f + 1.0;  
    printf ("  +1 = %e = %s\n", f, print_float (f));  
 
    printf ("\n\n  Display the dread effects of catastrophic cancellation\n\  
  by calculating a*b-a*c, in several different ways.  These are equivalent\n\  
  arithmetically, but not computationally.\n");  
    /* Now show the effect of catastrophic cancellation by calculating  
     * ab-ac, when a, b and c are such that the terms nearly cancel.  
     * To generate numbers which best demonstrate the effect,  
     * set a = (1+da/2^12), b=(1+db/2^12), etc.  Thus  
     * a * b = (1 + db/2^12 + da/2^12 + da.db/2^24).  
     * Put da=1.0.  Now pick db=1-, so that final term is just below 1/2^24,  
     * then dc=1+, so that final term is just above 1/2^24.  
     * Thus a*b rounds down, and a*c rounds up.  
     * Both are quite accurate, but large errors are revealed when  
     * they’re subtracted.  We can ameliorate this by rewriting the expression.  
     */  
 
    /* First, do the calculation in double, to get relative errors. */  
    /* The cancellation is ignorable in double */  
    ad = 1.0 + ((double)1/4096);  
    bd = 1.0 + ((double)0.9/4096);  
    cd = 1.0 + ((double)1.1/4096);  
    fd = ad*bd-ad*cd;  
 
    a = 1.0 + ((float)1/4096);  /* a=1.000244; */  
    b = 1.0 + ((float)0.9/4096); /* b=1.000220; */  
    c = 1.0 + ((float)1.1/4096); /* c=1.000269; */  
 
    /* first method - naive */  
    f = a*b-a*c;  
    printf ("a=%e  b=%e  c=%e\n", a, b, c);  
    printf ("a*b-a*c       = %e  (error=%e)\n", f, ((double)f/fd-1.0));  
 
    /* pre-subtract the nearly-equal b and c */  
    f = a*(b-c);  
    printf ("a*(b-c)       = %e  (%e)\n", f, ((double)f/fd-1.0));  
 
    /* rewrite the expression, to calculate a * ((b-1) - (c-1)).  
       Thus b-1 and c-1 have full accuracy */  
    b = ((float)0.9/4096);      /* = (above b) -1 */  
    c = ((float)1.1/4096);  
    f = a*(b-c);  
 
    printf ("a((b-1)-(c-1))= %e  (%e)\n", f, ((double)f/fd-1.0));  
    /* Can’t do the same trick with a.  If we calculate (a-1)*() + 1*(),  
       we don’t get any improvement.  We’re not carelessly discarding  
       accuracy, now - we can’t keep any more than this. */  
 
    printf ("\n  ...and further illustrate what’s happening by showing\n\  
  b and b-1.  Note the extra accuracy in the latter.\n");  
    /* Display b and (b-1).  Note extra accuracy in latter. */  
    b = 1.0 + ((float)0.9/4096);/* = (above b) */  
    printf ("1+1/10000 = %14.7e = %s\n", b, print_float (b));  
    b = ((float)0.9/4096);      /* = (above b) -1 */  
    printf ("  1/10000 = %14.7e = %s\n", b, print_float (b));;  
 
    printf ("\n\n  Display Infinity and NaN\n");  
    /* Display NaNs and Infinity */  
    /* Don’t just write a=1.0/0.0, since compiler can warn about this */  
    a = 0.0;                    /* and log(0) */  
    a = 1/a;  
    b = 1/a;  
    c = a*b;  
    printf ("a = 1/0.0 = %14e = %s\n", a, print_float (a));  
    printf ("b = 1/a   = %14e = %s\n", b, print_float (b));  
    printf ("c = a*b   = %14e = %s\n", c, print_float (c));  
}  

      program silly  
      implicit none  
      integer n, i  
      real determinant  
      parameter (n=2)  
      real twobytwo(2,2) /4*-1/  
      do i = 1, 100000  
         call mkidentity (twobytwo, n)  
      enddo  
      print *, determinant (twobytwo)  
      end

      subroutine mkidentity (matrix, dim)  
      implicit none  
      integer dim  
      real matrix (dim,dim)  
      integer m,n  
      do m = 1, dim  
         do n = 1, dim  
            if (m.eq.n) then  
               matrix(m,n) = 1.  
            else  
               matrix(m,n) = 0.  
            endif  
         enddo  
      enddo  
      return  
      end

      real function determinant(m)  
      implicit none  
      real m(2,2)  
      determinant = m(1,1) * m(2,2) - m(1,2) * m(2,1)  
      return  
      end

#include <stdio.h>  
 
int bananas (int i)  
{  
    char buf[80];  
    sprintf (buf, "banana number %d\n", i);  
    if (i != 0)  
        fputs (buf, stdout);  
    else  
    {  
        int *zp = 0;  
        printf ("banana split: %d\n", *zp);  
    }  
    return 0;  
}  
 
int main (void)  
{  
    printf ("Entering the bananas function: 1\n");  
    bananas (1);  
    printf ("No bananas left!\n");  
    bananas (0);  
    printf ("We made it out!\n");  
}

#include <stdio.h>  
 
int main (void)  
{  
    float m[3][3];  
    int i, j, n;  
 
    /*  
     * Note the C declaration of this Fortran function.  Although  
     * we’re passing a matrix, we declare the first argument just  
     * float* .  This is because a dynamically-dimensioned array can’t  
     * be expressed in C.  
     */  
    extern void addarr (float *m, int *n);  
 
    n = 3;                      /* dimension of matrix */  
    for (i=0; i<3; i++)  
        for (j=0; j<3; j++)  
            m[i][j] = (float)(10*i+j);  
 
    printf ("Before:\n");  
    for (i=0; i<3; i++)  
    {  
        printf ("m[%d][j]: ", i);  
        for (j=0; j<3; j++)  
            printf ("%6.0f", m[i][j]);  
        printf ("\n");  
    }  
 
    /* Pass the array, coerced to (float*). */  
    /* &m[0][0] (address of the first element of the array) would also work. */  
    addarr ((float*)m, &n);  
 
    printf ("\nAfter:\n");  
    for (i=0; i<3; i++)  
    {  
        printf ("m[%d][j]: ", i);  
        for (j=0; j<3; j++)  
            printf ("%6.0f", m[i][j]);  
        printf ("\n");  
    }  
}  
 
/* cadd2 adds twice its second argument to its first. */  
void cadd2 (float *a1, float a2)  
{  
    a2 *= 2;            /* double a2 - this doesn’t change a2 in caller */  
    *a1 += a2;          /* add a2 to what a1 points to */  
}

c   Given a square array m(n,n), call cadd2 on each element  
      subroutine addarr (m, n)  
      real m(n,n)  
      integer n  
      integer i,j  
      real addval  
 
      addval = m(n,n)  
      do j=1,n  
         do i=1,n  
            call cadd2 (m(i,j), %val(addval))  
         enddo  
      enddo  
 
      return  
      end

# read into maple with   > read ‘maple.ms‘:  
 
# Define a 2-d gaussian  
gi:=amp*exp(-(aparam^2/sa^2+bparam^2/sb^2)/2):  
aparam:=cos(theta)*(xc-x0)+sin(theta)*(yc-y0):  
bparam:=-sin(theta)*(xc-x0)+cos(theta)*(yc-y0):  
 
# Declare arrays.  We expressed the gaussian in terms of intelligible  
# identifiers, but we want these to be output as array references when  
# we produce the Fortran at the end  
xa := array(1..5):  
ca := array(1..1):  
 
x0 := xa[1]:  
y0 := xa[2]:  
sa := xa[3]:  
sb := xa[4]:  
theta := xa[5]:  
amp := ca[1]:  
 
 
# aadf are the variables we’ll differentiate with respect to  
aadf:=[x0,y0,sa,sb,theta]:  
# ...and array of differentials.  
dyda:=array(1..6):  
for i to 5 do  
        dyda[i]:=diff(gi,aadf[i])  
od:  
dyda[6] := gi:  
 
# Output the result as Fortran  
fortran (dyda, filename=‘gaussianab.f‘, optimized);

; Common IDL functions  
 
; Read the contents of an array from a file.  
;  
; The file is in ‘Fortran format’: it holds the contents of the array  
; as a list of numbers, with the left-most array index changing  
; fastest.  For simplicity, the function does no error checking, and  
; assumes that the array is the right length.  
pro read_file, a, fn  
        openr, iu, fn, /get_lun    ; open the file to read, allocate unit  
        readf, iu, a  
        free_lun, iu               ; close the file and release the unit.  
end  
 
 
; Read columns of numbers into a 2-d array.  
;  
; Given a file with  n  columns and  m  rows, this takes an  n x m  
; array and fills it with the contents of the file.  
;  
; It ignores lines beginning #  
pro read_columns, a, fn  
        dims = size (a)            ; get array dimensions  
        line = ’’  
        print, "cols=",dims[1]," rows=",dims[2]  
        openr, iu, fn, /get_lun  
        linevals = fltarr(dims[1])  
        i = 0  
        while not eof(iu) do begin  
                readf, iu, line  
                if (strmid(line,0,1) ne ’#’) then begin  
                        reads,line,linevals  
                        if (n_elements(linevals) ne dims[1]) then begin  
                                print,’line ’,i,": expected’,dims[1], $  
                                        ’ got’,n_elements(linevals)  
                                return  
                        endif  
                        if (i ge dims[2]) then message,’too many lines in file’  
                        a(*,i) = linevals  
                        i = i+1  
                endif  
        end  
        if (i ne dims[2]) then message,’too few lines in file’  
        close, iu  
end