The VMS SharkOpenVMS Notes: Fun with Floats

  1. The information and software presented on this web site are intended for educational use only by OpenVMS application developers and OpenVMS system attendants.
  2. The information and software presented on this web site are provided free of charge.
  3. The information and software presented on this web site are presented to you as-is. I will not be held responsible in any way if the information and software presented on this web site damages your computer system, business or organization (sounds like the legal warning from a Microsoft shrink-wrap seal, eh?)
  4. Is this text too small? You have two options:
    1. hold down the CTRL key while rolling the mouse wheel (zoom-in, zoom-out)
    2. use your keyboard like so:
      • hit: CTRL with "-" key to zoom smaller
      • hit: CTRL with "+" key to zoom larger
      • hit: CTRL with zero key to reset zoom

Menu

DFT Demo Bug Overview

Food-for-thought. I have been rereading technical books on DFT-FFT (Discreet Fourier Transform - Fast Fourier Transform) and most contain example programs written in BASIC (probably since they are readable by most people). About half of the books amusingly contain algebra directly coded into BASIC without any regard to the short comings of computer data types. For example, the following example contains a FOR-NEXT loop based upon floats like so:

    FOR I = 0 TO 2*PI STEP PI/8
I'm certain everyone reading this already knows the number-of-loops will change based upon how 2*PI and PI/8 are represented internally. Surprisingly, the VMS BASIC version of the program works properly when variables are changed from DOUBLE to XFLOAT but this could be just a fluke.

Obviously the author should have used a FOR-NEXT loop based upon integers which would manipulate some fraction of PI (unless it was his desire to keep his examples very uncomplicated)

Example Program written for GWBASIC (IBM Compatible PCs)

1 REM =============================================================
2 REM  title   : DFT1_0.bas (fig 1.3 on Page 8)
3 REM  book    : Understanding the FFT (c) Anders E. Zonst
4 REM  language: GWBASIC for Windows PC
5 REM =============================================================
10 REM *** DFT1.0 - GENERATE SQUARE WAVE ***
12 INPUT "NUMBER OF TERMS";N
20 PI = 3.14159265358#
30 FOR I = 0 TO 2*PI STEP PI/8
32 Y=0
40 FOR J=1 TO N STEP 2: Y=Y+SIN(J*I)/J: NEXT
50 PRINT Y
60 NEXT I
70 END

Example Program ported to VMS BASIC (for Alpha and Itanium Servers)

Notes about floats:
  1. SINGLE, DOUBLE, and GFLOAT represent traditional VMS floating point data types
  2. SFLOAT, TFLOAT, and XFLOAT represent new IEEE floating point data types supported on Alpha and Itanium but not VAX
  3. the following examples were compiled with Alpha BASIC V1.6-000 on OpenVMS-8.3
    • a line-by-line comparison of the next two displays prove that DOUBLE is no better than TFLOAT. In fact, it appears that DOUBLE only has 15 decimal digits of precision rather than the advertised 16. Could this be a bug? Perhaps.
  4. a line-by-line comparison of the middle two displays prove that XFLOAT is better than TFLOAT 
10 %title "dft1_0"
   declare string constant k_program = "dft1_0"
!=============================================================
! title    : DFT1_0.bas (fig 1.3 on Page 8)
! book     : Understanding the FFT (c) Anders E. Zonst
! language : GWBASIC for Windows PC
!=============================================================
!10 REM       *** DFT1.0 - GENERATE SQUARE WAVE  ***
!12 INPUT "NUMBER OF TERMS";N
!20 PI = 3.14159265358#
!30 FOR I = 0 TO 2*PI STEP PI/8
!32 Y=0
!40 FOR J=1 TO N STEP 2: Y=Y+SIN(J*I)/J: NEXT
!50 PRINT Y
!60 NEXT I
!70 END
!=============================================================
! language : HP BASIC for OpenVMS
! port     : Neil Rieck
! notes    : change lexical %hack to refine results
!-------------------------------------------------------------
! SINGLE ( 32-bit)  .29 * 10^-38   to 1.7  * 10^38    6 digits
! DOUBLE ( 64-bit)  .29 * 10^-38   to 1.7  * 10^38   16 digits
! GFLOAT ( 64-bit)  .56 * 10^-308  to  .90 * 10^308  15 digits
! SFLOAT ( 32-bit) 1.18 * 10^-38   to 3.40 * 10^38    6 digits
! TFLOAT ( 64-bit) 2.23 * 10^-308  to 1.80 * 10^308  15 digits
! XFLOAT (128-bit) 6.48 * 10^-4966 to 1.19 * 10^4932 33 digits
!=============================================================
   option type=explicit         ! required by HP BASIC for OpenVMS
   option angle=radians         ! be sure of trig units
   %let  %hack=2%
   %if   %hack=0% %then
     option size= (real double, integer long)
     declare double i, j, y, n, z
     print "-i-real type: double"
   %end %if
   %if   %hack=1% %then
     option size= (real tfloat, integer long)
     declare tfloat i, j, y, n, z
     print "-i-real type: tfloat"
   %end %if
   %if   %hack=2% %then
     option size= (real xfloat, integer long)
     declare xfloat i, j, y, n, z
     print "-i-real type: xfloat"
   %end %if
11 margin #0, 132
13 print "-i-program: "+ k_program +" (generate square wave)"
   print "display float precision"
   print "======================="
   print       "              1234567890123456789012345678901234"
   print       "pi (ref) =  3.14159265358979323846264338327950288419716939937510"
   print using "pi       = ##.##################################";pi
   print using "pi*2     = ##.##################################";pi*2.0
   print using "pi/8     = ##.##################################";pi/8.0
   print using "sin(1.0) = ##.##################################";sin(1.0)
   print using "sin(pi)  = ##.##################################";sin(pi)
19 input "NUMBER OF TERMS";n
20 z = 0
25 print       "Line  Y 1234567890123456789012345678901234  I 1234567890123456789012345678901234"
30 for i = 0 to 2*PI step PI/8
32 y = 0
40 for j=1 to n step 2 \ y=y+sin(j*i)/j \ next j
50 print using "#### ##.##################################";z;y;
51 print using    " ##.##################################";i
55 z = z + 1
60 next i
70 end

Sample Run #1 (VMS double): 16 lines are output from the loop

 $ run  DFT1_0
-i-real type: double (supposed to be 16 digits of precision
-i-program: dft1_0 (generate square wave)
display float precision
=======================
              1234567890123456789012345678901234
pi (ref) =  3.14159265358979323846264338327950288419716939937510 (yellow indicates 16 digits)
pi       =  3.1415926535897900000000000000000000 (bold indicates 15 digits of precision)
pi*2     =  6.2831853071795900000000000000000000 (bold indicates 15 digits of precision)
pi/8     =  0.3926990816987240000000000000000000 (bold indicates 15 digits of precision)
sin(1.0) =  0.8414709848078970000000000000000000 (bold indicates 15 digits of precision)
sin(pi)  =  0.0000000000000001224646799147350000 (bold indicates 15 digits of precision)
NUMBER OF TERMS? 1000
Line  Y 1234567890123456789012345678901234  I 1234567890123456789012345678901234
   0  0.0000000000000000000000000000000000  0.0000000000000000000000000000000000
   1  0.7867047098266320000000000000000000  0.3926990816987240000000000000000000
   2  0.7846910587375420000000000000000000  0.7853981633974480000000000000000000
   3  0.7859393587706950000000000000000000  1.1780972450961700000000000000000000
   4  0.7848981638974470000000000000000000  1.5707963267949000000000000000000000 (red indicates diff from TFLOAT)
   5  0.7859393587706950000000000000000000  1.9634954084936200000000000000000000
   6  0.7846910587375420000000000000000000  2.3561944901923400000000000000000000
   7  0.7867047098266320000000000000000000  2.7488935718910700000000000000000000
   8  0.0000000000000602504127728754000000  3.1415926535897900000000000000000000
   9 -0.7867047098266310000000000000000000  3.5342917352885200000000000000000000
  10 -0.7846910587375420000000000000000000  3.9269908169872400000000000000000000
  11 -0.7859393587706940000000000000000000  4.3196898986859700000000000000000000
  12 -0.7848981638974470000000000000000000  4.7123889803846900000000000000000000
  13 -0.7859393587706930000000000000000000  5.1050880620834200000000000000000000
  14 -0.7846910587375440000000000000000000  5.4977871437821400000000000000000000
  15 -0.7867047098266310000000000000000000  5.8904862254808600000000000000000000 <<<--- oops, loop ended early
$

Sample Run #2 (IEEE tfloat): 16 lines are output from the loop

$ run  DFT1_0
-i-real type: tfloat (15 digits of precision)
-i-program: dft1_0 (generate square wave)
display float precision
=======================
              1234567890123456789012345678901234
pi (ref) =  3.14159265358979323846264338327950288419716939937510 (yellow indicates 15 digits)
pi       =  3.1415926535897900000000000000000000 (bold indicates 15 digits of precision)
pi*2     =  6.2831853071795900000000000000000000 (bold indicates 15 digits of precision)
pi/8     =  0.3926990816987240000000000000000000 (bold indicates 15 digits of precision)
sin(1.0) =  0.8414709848078970000000000000000000 (bold indicates 15 digits of precision)
sin(pi)  =  0.0000000000000001224646799147350000 (bold indicates 15 digits of precision)
NUMBER OF TERMS? 1000
Line  Y 1234567890123456789012345678901234  I 1234567890123456789012345678901234
   0  0.0000000000000000000000000000000000  0.0000000000000000000000000000000000
   1  0.7867047098266320000000000000000000  0.3926990816987240000000000000000000
   2  0.7846910587375420000000000000000000  0.7853981633974480000000000000000000
   3  0.7859393587706950000000000000000000  1.1780972450961700000000000000000000
   4  0.7848981638974460000000000000000000  1.5707963267949000000000000000000000 (red indicates diff from DOUBLE)
   5  0.7859393587706950000000000000000000  1.9634954084936200000000000000000000
   6  0.7846910587375420000000000000000000  2.3561944901923400000000000000000000
   7  0.7867047098266310000000000000000000  2.7488935718910700000000000000000000
   8  0.0000000000000617440257325753000000  3.1415926535897900000000000000000000
   9 -0.7867047098266310000000000000000000  3.5342917352885200000000000000000000
  10 -0.7846910587375420000000000000000000  3.9269908169872400000000000000000000
  11 -0.7859393587706940000000000000000000  4.3196898986859700000000000000000000
  12 -0.7848981638974460000000000000000000  4.7123889803846900000000000000000000
  13 -0.7859393587706940000000000000000000  5.1050880620834100000000000000000000
  14 -0.7846910587375410000000000000000000  5.4977871437821400000000000000000000
  15 -0.7867047098266350000000000000000000  5.8904862254808600000000000000000000 <<<--- oops, loop ended early
$

Sample Run #3 (IEEE xfloat): 17 lines are output from the loop

$ run  DFT1_0
-i-real type: xfloat (33 digits of precision)
-i-program: dft1_0 (generate square wave)
display float precision
=======================
              1234567890123456789012345678901234
pi (ref) =  3.14159265358979323846264338327950288419716939937510 (yellow indicates 33 digits)
pi       =  3.1415926535897932384626433832795000 (bold indicates 33 digits of precision)
pi*2     =  6.2831853071795864769252867665590000 (bold indicates 33 digits of precision)
pi/8     =  0.3926990816987241548078304229099380 (bold indicates 33 digits of precision)
sin(1.0) =  0.8414709848078965066525023216302990 (bold indicates 33 digits of precision)
sin(pi)  =  0.0000000000000000000000000000000009 (???)
NUMBER OF TERMS? 1000
Line  Y 1234567890123456789012345678901234  I 1234567890123456789012345678901234
   0  0.0000000000000000000000000000000000  0.0000000000000000000000000000000000
   1  0.7867047098266324129536091420549900  0.3926990816987241548078304229099380
   2  0.7846910587375418025179337058983590  0.7853981633974483096156608458198760
   3  0.7859393587706949526402299139422830  1.1780972450961724644234912687298100
   4  0.7848981638974458096461601533451350  1.5707963267948966192313216916397500
   5  0.7859393587706949526402299139422840  1.9634954084936207740391521145496900
   6  0.7846910587375418025179337058983610  2.3561944901923449288469825374596300
   7  0.7867047098266324129536091420549900  2.7488935718910690836548129603695600
   8  0.0000000000000000000000000000002319  3.1415926535897932384626433832795000
   9 -0.7867047098266324129536091420549900  3.5342917352885173932704738061894400
  10 -0.7846910587375418025179337058983620  3.9269908169872415480783042290993800
  11 -0.7859393587706949526402299139422840  4.3196898986859657028861346520093200
  12 -0.7848981638974458096461601533451350  4.7123889803846898576939650749192500
  13 -0.7859393587706949526402299139422820  5.1050880620834140125017954978291900
  14 -0.7846910587375418025179337058983460  5.4977871437821381673096259207391300
  15 -0.7867047098266324129536091420549890  5.8904862254808623221174563436490700
  16 -0.0000000000000000000000000000012410  6.2831853071795864769252867665590000 <<<--- better
$

Example Program repaired (17 output lines guaranteed)

10 %title "dft1_0_alt"
   declare string constant k_program = "dft1_0_alt"
!=============================================================
! title    : DFT1_0.bas (fig 1.3 on Page 8)
! book     : Understanding the FFT (c) Anders E. Zonst
! language : GWBASIC for Windows PC
!=============================================================
!10 REM       *** DFT1.0 - GENERATE SQUARE WAVE  ***
!12 INPUT "NUMBER OF TERMS";N
!20 PI = 3.14159265358#
!30 FOR I = 0 TO 2*PI STEP PI/8
!32 Y=0
!40 FOR J=1 TO N STEP 2: Y=Y+SIN(J*I)/J: NEXT
!50 PRINT Y
!60 NEXT I
!70 END
!=============================================================
! language : HP BASIC for OpenVMS
! port     : Neil Rieck
! notes    : change lexical %hack to refine results
!-------------------------------------------------------------
! SINGLE ( 32-bit)  .29 * 10^-38   to 1.7  * 10^38    6 digits
! DOUBLE ( 64-bit)  .29 * 10^-38   to 1.7  * 10^38   16 digits
! GFLOAT ( 64-bit)  .56 * 10^-308  to  .90 * 10^308  15 digits
! SFLOAT ( 32-bit) 1.18 * 10^-38   to 3.40 * 10^38    6 digits
! TFLOAT ( 64-bit) 2.23 * 10^-308  to 1.80 * 10^308  15 digits
! XFLOAT (128-bit) 6.48 * 10^-4966 to 1.19 * 10^4932 33 digits
!=============================================================
   option type=explicit         ! required by HP BASIC for OpenVMS
   option angle=radians         ! be sure of trig units
   %let  %hack=2%
   %if   %hack=0% %then
     option size= (real double, integer long)
     declare double i, j, y, n, z
     print "-i-real type: double"
   %end %if
   %if   %hack=1% %then
     option size= (real tfloat, integer long)
     declare tfloat i, j, y, n, z
     print "-i-real type: tfloat"
   %end %if
   %if   %hack=2% %then
     option size= (real xfloat, integer long)
     declare xfloat i, j, y, n, z
     print "-i-real type: xfloat"
   %end %if
11 margin #0, 132
13 print "-i-program: "+ k_program +" (generate square wave)"
   print "display float precision"
   print "======================="
   print       "              1234567890123456789012345678901234"
   print       "pi (ref) =  3.14159265358979323846264338327950288419716939937510"
   print using "pi       = ##.##################################";pi
   print using "pi*2     = ##.##################################";pi*2.0
   print using "pi/8     = ##.##################################";pi/8.0
   print using "sin(1.0) = ##.##################################";sin(1.0)
   print using "sin(pi)  = ##.##################################";sin(pi)
19 input "NUMBER OF TERMS";n
25 print       "Line  Y 1234567890123456789012345678901234  I 1234567890123456789012345678901234"
30 for z = 0 to 16
31 i = z * pi / 8 
32 y = 0
40 for j=1 to n step 2 \ y=y+sin(j*i)/j \ next j
50 print using "#### ##.##################################";z;y;
51 print using    " ##.##################################";i
60 next z
70 end

Sample Run #3a (IEEE xfloat): 17 lines guaranteed

$ run  DFT1_0_alt
-i-real type: xfloat
-i-program: dft1_0_alt (generate square wave)
display float precision
=======================
              1234567890123456789012345678901234
pi (ref) =  3.14159265358979323846264338327950288419716939937510
pi       =  3.1415926535897932384626433832795000
pi*2     =  6.2831853071795864769252867665590000
pi/8     =  0.3926990816987241548078304229099380
sin(1.0) =  0.8414709848078965066525023216302990
sin(pi)  =  0.0000000000000000000000000000000009
NUMBER OF TERMS? 1000
Line  Y 1234567890123456789012345678901234  I 1234567890123456789012345678901234
   0  0.0000000000000000000000000000000000  0.0000000000000000000000000000000000
   1  0.7867047098266324129536091420549900  0.3926990816987241548078304229099380
   2  0.7846910587375418025179337058983590  0.7853981633974483096156608458198760
   3  0.7859393587706949526402299139422830  1.1780972450961724644234912687298100
   4  0.7848981638974458096461601533451350  1.5707963267948966192313216916397500
   5  0.7859393587706949526402299139422840  1.9634954084936207740391521145496900
   6  0.7846910587375418025179337058983580  2.3561944901923449288469825374596300
   7  0.7867047098266324129536091420549900  2.7488935718910690836548129603695600
   8  0.0000000000000000000000000000004285  3.1415926535897932384626433832795000
   9 -0.7867047098266324129536091420549910  3.5342917352885173932704738061894400
  10 -0.7846910587375418025179337058983580  3.9269908169872415480783042290993800
  11 -0.7859393587706949526402299139422840  4.3196898986859657028861346520093200
  12 -0.7848981638974458096461601533451350  4.7123889803846898576939650749192500
  13 -0.7859393587706949526402299139422820  5.1050880620834140125017954978291900
  14 -0.7846910587375418025179337058983460  5.4977871437821381673096259207391300
  15 -0.7867047098266324129536091420549910  5.8904862254808623221174563436490700
  16 -0.0000000000000000000000000000008570  6.2831853071795864769252867665590000
$

Comparisons between run 3 and run 3a

 Comparison (16 iterations - one wave):
line Y from DFT1_0 Y from DFT1_0_alt Comments
0000 0.000000000000000000000000000000000 0.000000000000000000000000000000000 same
0008 0.000000000000000000000000000000232 0.000000000000000000000000000000429 a little higher
0016 -0.000000000000000000000000000001241 -0.000000000000000000000000000000857 a little lower
Comparison (160 iterations - ten waves):
line Y from DFT1_0 Y from DFT1_0_alt Comments
0160 -0.000000000000000000000000000096384 -0.000000000000000000000000000007041 ten time lower
Comparison (1600 iterations - one hundred waves):
line Y from DFT1_0 Y from DFT1_0_alt Comments
1600 -0.000000000000000000000000009314737 -0.000000000000000000000000000093896 ten times lower

Float Precision Test

1000    !=====================================================================
        ! title  : Float_Precision_Test.bas
        ! author : NSR
        ! created: 2011-01-30
        ! notes  : this test was done with Alpha BASIC V1.6-000 on OpenVMS-8.3
        !=====================================================================
        option type=explicit                            !
        option size=(real xfloat,integer long)          !
        declare xfloat xf, &
                tfloat tf, &
                sfloat sf, &
                gfloat gf, &
                double d,  &
                single s,  &
                string yada$
2000    main:
        print "title      : Float Precision Test"
        print "environment: Alpha BASIC V1.6-000 on OpenVMS-8.3"
        yada$ = "0.04444444444444444444444444444444444"
        gosub convert_n_display
        yada$ = "0.05555555555555555555555555555555555"
        gosub convert_n_display
        yada$ = "0.01234567891234567891234567891234567"
        gosub convert_n_display
        goto fini
        !
        convert_n_display:
        xf = real(yada$,XFLOAT)
        tf = xf
        sf = xf
        gf = xf
        d  = xf
        s  = xf
        print       "test data = "; yada$
        print using "xfloat    = #.###################################";xf
        print using "tfloat    = #.###################################";tf
        print using "sfloat    = #.###################################";sf
        print using "gfloat    = #.###################################";gf
        print using "double    = #.###################################";d
        print using "single    = #.###################################";s
        print       "               0000000001111111111222222222233333"
        print       "               1234567890123456789012345678901234 precision"
        print
        return
        !
32000   fini:
        end

Float Precision Test: Output 

$ run  FLOAT_PRECISION_DEMO

title      : Float Precision Test
environment: Alpha BASIC V1.6-000 on OpenVMS-8.3
float precision test - HP BASIC for OpenVMS Alpha
test data = 0.04444444444444444444444444444444444
xfloat    = 0.04444444444444444444444444444444440 33 digits
tfloat    = 0.04444444444444440000000000000000000 15 digits
sfloat    = 0.04444440000000000000000000000000000  6 digits
gfloat    = 0.04444444444444440000000000000000000 15 digits
double    = 0.04444444444444440000000000000000000 15 digits (shouldn't this be 16?)
single    = 0.04444440000000000000000000000000000  6 digits
               0000000001111111111222222222233333
               1234567890123456789012345678901234 precision

test data = 0.05555555555555555555555555555555555
xfloat    = 0.05555555555555555555555555555555560
tfloat    = 0.05555555555555560000000000000000000
sfloat    = 0.05555560000000000000000000000000000
gfloat    = 0.05555555555555560000000000000000000
double    = 0.05555555555555560000000000000000000
single    = 0.05555560000000000000000000000000000
               0000000001111111111222222222233333
               1234567890123456789012345678901234 precision

test data = 0.01234567891234567891234567891234567
xfloat    = 0.01234567891234567891234567891234570
tfloat    = 0.01234567891234570000000000000000000
sfloat    = 0.01234570000000000000000000000000000
gfloat    = 0.01234567891234570000000000000000000
double    = 0.01234567891234570000000000000000000
single    = 0.01234570000000000000000000000000000
               0000000001111111111222222222233333
               1234567890123456789012345678901234 precision
$

External Links

Back to OpenVMS
Back to Home
Neil Rieck
Kitchener - Waterloo - Cambridge, Ontario, Canada.