%% $RCSfile: sdsprandsrc2.tlc,v $
%% $Revision: 1.9 $
%% $Date: 2001/04/25 18:24:31 $
%%
%% Copyright 1995-2001 The MathWorks, Inc.
%%
%% Abstract: Generate Uniform or Normal (Gaussian) Random Numbers 
%%

%implements "sdsprandsrc2" "C"

%% Function: BlockTypeSetup ================================================
%% Abstract:
%%
%function BlockTypeSetup(block, system) void

%% Render a DSP_InitializeSeed function once for use by all sdsprandsrc blocks.
%% Needed for both the Uniform and Gaussian cases.

%% First, cache the function prototype for DSP_InitializeSeed:
%%
%openfile InitSeedBuff
extern void DSP_InitializeSeed(uint32_T *urandSeed, real_T initSeed);
%closefile InitSeedBuff
%<DSPAddToFileHeader(InitSeedBuff)>

%% Next, cache the DSP_InitializeSeed function itself:
%%
%openfile InitSeedBuff
/* Function: DSP_InitializeSeed
 *  Bit-shift the given initial seed
 */
extern void DSP_InitializeSeed(uint32_T *urandSeed, real_T initSeed)
{
    const uint32_T maxseed = 2147483646;   /* 2^31-2 */
    const uint32_T seed0   = 1144108930;   /* Seed #6, starting from seed = 1 */
    const uint32_T bit16   = 32768;        /* 2^15   */

    *urandSeed = (uint32_T)initSeed;

    /* Interchange bits 1-15 and 17-31 */
    {
        int_T r = *urandSeed >> 16;
        int_T t = *urandSeed & bit16;
        *urandSeed = ((*urandSeed - (r << 16) - t) << 16) + t + r;
    }
    if (*urandSeed < 1) {
        *urandSeed = seed0;
    }
    if (*urandSeed > maxseed) {
        *urandSeed = maxseed;
    }
} /* end DSP_InitializeSeed */

%closefile InitSeedBuff

%<DSPAddToFile(InitSeedBuff)>

%endfunction  %% BlockTypeSetup


%% Function: BlockInstanceSetup ================================================
%% Abstract:
%%
%function BlockInstanceSetup(block, system) void

%assign src_type = block.SFcnParamSettings.SrcType

%if src_type == "Uniform"
    %<RenderUniformRandFcn()>
%else
    %<RenderUniformRandFcn()>
    %<RenderNormalRandFcn()>
%endif

%assign IS_COMPLEX = LibBlockOutputSignalIsComplex(0)
%assign INHERIT_ON = (block.SFcnParamSettings.InheritOn == "Yes")
%assign block = block + IS_COMPLEX + src_type + INHERIT_ON

%endfunction %% BlockInstanceSetup


%% Function: Start ================================================
%% Abstract:
%%      Initialize the real and/or imag seeds for all channels.
%%        Compute the first random seed
%%
%function Start(block, system) Output
    /* DSP Blockset Random Source (%<ParamSettings.FunctionName>) - %<Name> */
    /* Initialize the Random Seeds */
    %assign OUTPORT_NUM   = 0
    %assign numDims       = LibBlockOutputSignalNumDimensions(OUTPORT_NUM)
    %assign dims          = LibBlockOutputSignalDimensions(OUTPORT_NUM)
    %assign nchans        = DetermineChannels()
    %assign multipleChans = (nchans > 1)
    %if multipleChans
    {
    %endif
    %assign seedLen = LibGetNumberOfElements(InitSeed.Value)
    %%
    %% Determine datatype for the random seed:
    %%
    %assign urandDType = (IS_COMPLEX) ? "cuint32_T *" : "uint32_T *"
    %if multipleChans
        %<urandDType>urandSeed = (%<urandDType>)%<LibBlockDWorkAddr(RAND_SEED,"","",0)>;
    %endif
        %<DSP_RenderInitializeSeedLoop(nchans,multipleChans,IS_COMPLEX,seedLen)>  
    %if multipleChans
    }
    %endif

%endfunction


%%%%%%%%%%%%%%%%%%%%% 
%% Function: Outputs =============================================================
%%
%function Outputs(block, system) Output   
    /* DSP Blockset Random Source (%<ParamSettings.FunctionName>) - %<Name> */
    /* Create the Random Numbers */
    %assign OUTPORT_NUM = 0
    %assign IS_UNIFORM  = CAST("Boolean",(src_type == "Uniform"))    
    %%
    %assign numDims       = LibBlockOutputSignalNumDimensions(OUTPORT_NUM)
    %assign dims          = LibBlockOutputSignalDimensions(OUTPORT_NUM)
    %assign nchans        = DetermineChannels()
    %if INHERIT_ON
      %assign frameSize   = dims[0]/nchans
    %else
      %assign frameSize   = (SFcnParamSettings.IsDiscrete=="Yes")?dims[0]:1
    %endif
    %%
    %assign multipleChans    = (nchans >1)
    %assign isScalar         = (!multipleChans && frameSize == 1)
    %assign isUniformScalar  = CAST("Boolean",0)
    %assign isGaussianScalar = CAST("Boolean",0)
    %%
    %% Determine datatype for the random seed:
    %assign urandDType = (IS_COMPLEX) ? "cuint32_T *" : "uint32_T *"
    %assign outDType   = (IS_COMPLEX) ? "creal_T *" : "real_T *"
    %%
    %if (IS_UNIFORM)
        %assign maxVal = Max.Value
        %assign maxLen = LibGetNumberOfElements(maxVal)
        %assign minVal = Min.Value
        %assign minLen = LibGetNumberOfElements(minVal)
        %assign isUniformScalar = CAST("Boolean",(minLen == 1 && maxLen == 1 && isScalar))
    %else
        %assign meanVal       = Mean.Value
        %assign meanLen       = LibGetNumberOfElements(meanVal)
        %assign varVal        = Variance.Value
        %assign varLen        = LibGetNumberOfElements(varVal)
        %assign isMeanComplex = CAST("Boolean",(Mean.ComplexSignal == "yes"))
        %assign meanDType     = (isMeanComplex) ? "creal_T *" : "real_T *"
        %assign isGaussianScalar = CAST("Boolean",(meanLen == 1 && varLen == 1 && isScalar))
    %endif
    %%
    %if isUniformScalar && nchans == 1
        %<DSP_ScalarUniformNumGen(IS_COMPLEX,OUTPORT_NUM)>\
    %elseif isGaussianScalar && nchans == 1
        %<DSP_ScalarGaussianNumGen(IS_COMPLEX,OUTPORT_NUM)>\
    %else
    {
        %<outDType>y           = (%<outDType>)%<LibBlockOutputSignalAddr(OUTPORT_NUM,"","",0)>;
        %<urandDType>urandSeed = (%<urandDType>)%<LibBlockDWorkAddr(RAND_SEED,"","",0)>;
        %%
        %if (IS_UNIFORM)
            real_T *pMin = (real_T *)%<LibBlockParameterAddr(Min,"","",0)>;
            real_T *pMax = (real_T *)%<LibBlockParameterAddr(Max,"","",0)>;
        %else
            %<meanDType>pMean = (%<meanDType>)%<LibBlockParameterAddr(Mean,"","",0)>;
            real_T     *pVar  = (real_T *)%<LibBlockParameterAddr(Variance,"","",0)>;
        %endif        
        %%
        %if multipleChans
            int_T i;
            for(i=0;i<%<nchans>;i++) {    
        %endif    
            %if (frameSize > 1) %%%<InputIsNonscalarFrame(frameSize,TID)>
                int_T j;
                for(j=0;j<%<frameSize>;j++ ) {
            %endif
                %%
                %<DSP_CallRandomNumberGenerator(IS_COMPLEX,IS_UNIFORM)>\
                %%
            %if (frameSize > 1) %%%<InputIsNonscalarFrame(frameSize,TID)>
                }
            %endif
            %%
        %if multipleChans
                urandSeed++;
                %if (IS_UNIFORM)
                    %if (maxLen > 1) 
                        pMax++;
                    %endif
                    %if (minLen > 1)
                        pMin++;
                    %endif 
                %else
                    %if (varLen  > 1) 
                        pVar++;
                    %endif
                    %if (meanLen > 1) 
                        pMean++;
                    %endif
                %endif            
            }
        %endif
    }
    %endif

%endfunction


%% Function: RenderUniformRandFcn ==============================================
%% Abstract:
%%     Render the DSP_UniformRand function and prototype.
%%     Only render the content ONCE.  Any additional calls
%%     are simply ignored.
%%
%function RenderUniformRandFcn() void

%assign database_entry = "sdsprandsrc_uniformrand_fcn"
%assign model_cache    = "::CompiledModel." + database_entry

%% Check info so that we do not define this function more than once:
%if !EXISTS("%<model_cache>")

    %% Retain definition to prevent multiple identical defines:
    %%
    %assign %<database_entry> = 1
    %assign ::CompiledModel = ::CompiledModel + %<database_entry>
    %undef %<database_entry>  %%Remove from block scope

    %% First, cache the function prototype for DSP_UniformRand:
    %%
    %% Render a DSP_UniformRand function once for use by all sdsprandsrc blocks.

    %openfile DSP_RandBuff
    /* DSP Blockset Random Source block Uniform random number generator */
    extern real_T DSP_UniformRand(uint32_T *seed);
    %closefile DSP_RandBuff
    %<DSPAddToFileHeader(DSP_RandBuff)>

    %% Next, cache the DSP_UniformRand function itself:
    %%
    %openfile DSP_RandBuff
    /*
     * DSP Blockset Random Source block
     * Uniform random number generator
     * Generates random number in range (0,1)
     */
    extern real_T DSP_UniformRand(uint32_T *seed)     /* pointer to a running seed */
    {
        const uint32_T IA = 16807;                /* magic multiplier = 7^5    */
        const uint32_T IM = 2147483647;            /* modulus = 2^31-1            */  
        const uint32_T IQ = 127773;                /* IM div IA                */
        const uint32_T IR = 2836;                /* IM modulo IA                */
        const real_T   S  = 4.656612875245797e-10;     /* reciprocal of 2^31-1    */

        uint32_T hi   = *seed / IQ;
        uint32_T lo   = *seed % IQ;
        int32_T  test = IA * lo - IR * hi;    /* never overflows */

        *seed = ((test < 0) ? (unsigned int)(test + IM) : (unsigned int)test);

        return( (real_T) ((*seed) * S) );

    } /* end DSP_UniformRand */

    %closefile DSP_RandBuff
    %<DSPAddToFile(DSP_RandBuff)>

%endif
%endfunction %% RenderUniformRandFcn


%% Function: RenderNormalRandFcn ==============================================
%% Abstract:
%%     Render the DSP_NormalRand function and prototype.
%%     Only render the content ONCE.  Any additional calls
%%     are simply ignored.
%%
%function RenderNormalRandFcn() void

%assign database_entry = "sdsprandsrc_normrand_fcn"
%assign model_cache    = "::CompiledModel." + database_entry

%% Check info so that we do not define this function more than once:
%if !EXISTS("%<model_cache>")

    %% Retain definition to prevent multiple identical defines:
    %%
    %assign %<database_entry> = 1
    %assign ::CompiledModel = ::CompiledModel + %<database_entry>
    %undef %<database_entry>  %%Remove from block scope

    %% First, cache the function prototype for DSP_NormalRand:
    %%
    %openfile DSP_RandBuff
    /* DSP Blockset Random Source block Gaussian random number generator */
    extern real_T DSP_NormalRand(uint32_T *seed);
    %closefile DSP_RandBuff
    %<DSPAddToFileHeader(DSP_RandBuff)>

    %% Next, cache the DSP_NormalRand function itself:
    %%
    %openfile DSP_RandBuff
    /* Function: DSP_NormalRand 
     *  Normal (Gaussian) random number generator 
     */
    extern real_T DSP_NormalRand(unsigned int *seed)
    {
        real_T sr, si, t;
        do {
            sr = 2.0 * DSP_UniformRand(seed) - 1.0;
            si = 2.0 * DSP_UniformRand(seed) - 1.0;
            t  = sr * sr + si * si;
        } while (t > 1.0);

        return(sr * sqrt((-2.0 * log(t)) / t));
    } /* end DSP_NormalRand */

    %closefile DSP_RandBuff
    %<DSPAddToFile(DSP_RandBuff)>

%endif
%endfunction  %% RenderNormalRandFcn