rem Normal Distribution .... Rev 3.1
      rem A J Tooth // December 2006
      
      *FLOAT64
      
      rem Input section
      input "What value for the Mean ",mu
      input "What value for the Standard Deviation ",sig
      input "What left-hand value for the Random Variable ",Xl
      input "What right-hand value for the Random Variable ",Xr
      
      rem Estimate integral
      P_Xl=fn_normal(mu,sig,Xl)
      P_Xr=fn_normal(mu,sig,Xr)
      
      PA_Xl=fn_norm(mu,sig,Xl)
      PA_Xr=fn_norm(mu,sig,Xr)
      
      print:print "Probability that X is less than ";Xl;" is: ";P_Xl
      print:print "Probability that X is less than ";Xr;" is: ";P_Xr
      
      print:print "ASM_Probability that X is less than ";Xl;" is: ";PA_Xl
      print:print "ASM_Probability that X is less than ";Xr;" is: ";PA_Xr
      
      end
      rem End of Program ===================================================
      rem ==================================================================
      
      rem Normal Distribution Function
      def fn_normal(Mu,Sig,X)
      local cum, stp, var, Res, iRes
      
      stp=0.0001
      cum=0.0
      var=-6.0
      
      rem Estimate integral
      repeat
        cum+=stp*exp(-(var*var)/2)
        var+=stp
      until var>((X-Mu)/Sig)
      cum=cum/sqr(2*pi)
      
      rem Round off to 4 decimal places
      Res=cum*10000
      iRes=1.0*int(Res)
      if Res-iRes<0.5 then Res=iRes else Res=iRes+1.0
      =Res/10000
      rem ==================================================================
      
      rem Calls Normal Assembler function
      def fn_norm(Mu,Sig,X)
      local A#,B#,C#,Ans#
      private stp%, normal%, cum%, var%, k%, lam%, itmp%, ftmp%, res%, Z%
      
      if Z%=0 then
        dim normal% 2000, cum% 7, var% 7, stp% 7, k% 3, lam% 7, itmp% 3, ftmp% 7, res% 7
        
        |stp%=0.0001 : rem Step value
        |cum%=0.0
        |var%=-6.0
        
        rem Dual-pass assembly, in case of labels
        for pass&=0 to 2 step 2
          proc_norm(pass&)
        next pass&
        
        Z%=1
      else
        rem Reset cumulative variables to original values
        |cum%=0.0
        |var%=-6.0
      endif
      
      Ans#=0.0
      A#=1.0*Mu
      B#=1.0*Sig
      C#=1.0*X
      
      call normal%,A#,B#,C#,Ans#
      
      =Ans#
      rem ============================================================
      
      rem Assembly Language Routine for the Normal Distribution
      def proc_norm(opt&)
      P%=normal%
      
      [opt opt&
      
      finit
      mov esi,[ebp+12]
      fld qword [esi]          ;X
      mov esi,[ebp+2]
      fld qword [esi]          ;Mu
      fsubp st1,st0
      mov esi,[ebp+7]
      fld qword [esi]          ;Sig
      fdivp st1,st0            ;(X-Mu)/Sig
      fstp qword [ftmp%]
      
      .rept
      fld qword [var%]
      fmul st0,st0
      fld1
      fadd st0,st0
      fdivp st1,st0
      fchs
      
      fldl2e                   ;st0=Log2(e)
      fmulp st1,st0
      fist dword [k%]          ;Int(st0) which = k
      fild dword [k%]
      fsubp st1,st0
      fstp qword [lam%]        ;st0 - Int(st0)
      
      mov ecx,[k%]             ;Is k<0, k=0, k>0 ??
      cmp ecx,0
      jl blow
      je zero
      
      mov eax,2                ;When k>0
      dec cl
      shl eax,cl               ;2^k
      mov [itmp%],eax
      fild dword [itmp%]
      jmp over                 ;When k>0
      
      .blow                    ;When k<0
      neg ecx
      mov eax,2
      dec cl
      shl eax,cl
      mov [itmp%],eax
      fild dword [itmp%]
      fld1
      fdivrp st1,st0           ;2^(-|k|)
      jmp over                 ;When k<0
      
      .zero
      fld1                     ;When k=0
      
      .over
      fld qword [lam%]
      f2xm1                    ;2^st0 - 1
      fld1
      faddp st1,st0
      fmulp st1,st0
      fld qword [stp%]
      fmulp st1,st0
      fld qword [cum%]
      faddp st1,st0
      fstp qword [cum%]
      
      fld qword [var%]
      fld qword [stp%]
      faddp st1,st0
      fst qword [var%]         ;var=var+stp
      fld qword [ftmp%]
      
      fsubrp st1,st0
      ftst                     ;Test whether st1>st0
      fstsw ax
      fstp qword [res%]        ;Discard st0
      and ah,1
      
      cmp ah,0
      je near rept
      
      fldpi
      fadd st0,st0
      fsqrt
      fld qword [cum%]
      fdivrp st1,st0
      
      fild dword [i10000%]     ;Round And truncate result
      fmulp st1,st0
      fistp dword [itmp%]
      fild dword [itmp%]       ;Accurate t0 4 decimal places
      fild dword [i10000%]
      fdivp st1,st0            ;Round And truncate result
      
      mov esi,[ebp+17]
      fstp qword [esi]
      
      jmp fin
      .i10000% DD 10000
      .fin ret
      ]
      endproc
      rem ============================================================