Abstract:
There is provided a system and method for statistical timing analysis of an electrical circuit. The system includes at least one parameter input, a statistical static timing analyzer, and at least one output. The at least one parameter input is for receiving parameters of the electrical circuit. At least one of the parameters has at least one of a non-Gaussian probability distribution and a non-linear delay effect. The statistical static timing analyzer is for calculating at least one of a signal arrival time and a signal required time for the electrical circuit using the at least one parameter. The at least one output is for outputting the at least one of the signal arrival time and the signal required time.

Description:
CROSS-REFERENCE TO RELATED APPLICATIONS  
       [0001]     This is a non-provisional application claiming the benefit of provisional application Ser. No. 60/620,039, entitled “System and Method for Accommodating Non-Gaussian and Non-linear Sources of Variation in Statistical Static Timing Analysis”, filed on Oct. 19, 2004, which is incorporated by reference herein. 
     
    
     FIELD OF THE INVENTION  
       [0002]     The present invention generally relates to statistical static timing analysis and, more particularly, to a system and method for accommodating non-Gaussian and nonlinear sources of variation in statistical static timing analysis.  
       BACKGROUND OF THE INVENTION  
       [0003]     It is commonly recognized that electrical characteristics of transistors and interconnects are not the same for different chips and even for the same chip at different time moments or chip locations. Variation of electrical characteristics can be due to variation of process parameters, changing of environmental conditions and even chips aging (Hot Carriers Injections, Negative Bias Temperature Instability, electromigration, and so forth). Variation of electrical characteristics results in variations of gate timing characteristics. The traditional conservative way to handle these variations is to consider so-called process corners at which the gates have the worst combinations of delays. Then chips are designed so that they can properly function at all process corners assuming that as a result they will function at any other combination of gate delays.  
         [0004]     Unfortunately, with decreasing transistor size and interconnect width, the variation of electrical characteristics is getting proportionally larger. Therefore, the approach to design for process corners, which used to work well, now results in too conservative and non-optimal designs because most design efforts and chip resources are spent to make chips function at very low-probability combinations of electrical characteristics.  
         [0005]     An alternative approach to designing chips as described above is to consider actual statistical characteristics of process parameter variations and use them to compute statistical characteristics of a designed circuit. For digital circuits, this approach is known as statistical timing analysis. There are several varieties of statistical timing analysis.  
         [0006]     One of the most useful for circuit analysis and optimization is parameterized statistical static timing analysis (STA). According to this technique, gate delays and signal arrival times are represented as functions of process parameters. 
 
 A=f ( X   1   , X   2 , . . . )  (1) 
 
         [0007]     All the parameters are assumed independent. This assumption significantly simplifies the analysis but does not limit its applicability because independence can be obtained by a principal component analysis technique. Using this representation, the parameterized STA computes a statistical approximation of the circuit timing characteristics (arrival and required arrival times, delay, timing slack) as functions of the same parameters.  
         [0008]     The process parameters can be considered as either globally correlated or fully random. Both types of parameters exhibit variation due to the manufacturing process. The correlated parameters have the same values for all gates and wires of the analyzed circuit while the fully random parameters vary for each gate independently. Usually, it is convenient to combine all the fully random parameters into one term.  
         [0009]     The parameterized statistical STA can be either path-based or block-based. Path-based statistical STA analyzes each signal propagation path separately and computes the probability distribution for circuit delay as the probabilistic maximum of all paths delays. Usually this requires enumeration of all signal propagation paths and integration in the space of parameters variations, which is an inefficient computational procedure.  
         [0010]     A more efficient technique of parameterized STA is so-called block-based statistical STA. This technique is very similar to traditional deterministic STA. It computes signal arrival times (or signal required arrival times) as functions of process parameters for each circuit node in their topological order similarly to propagating arrival times by a deterministic STA. It is convenient for both implementation and application. That kind of statistical timing analyzer can be easily implemented on the base of existing deterministic timing analyzer. This type of timing analysis lends itself to incremental operation, whereby after a change of the circuit is made, timing can be queried efficiently.  
         [0011]     For computational efficiency, parameterized statistical block-based STA assumes that all parameters variations have normal Gaussian distributions and that gate and wire delays depend on parameters linearly. Without loss of generality, Gaussian distributions with 0 mean and unit standard deviation are assumed. This assumption allows representation of gate delays in first-order canonical form, as described by C. Visweswariah in U.S. patent application Ser. No. 10/666,353, entitled “System and Method for Statistical Timing Analysis of Digital Circuits”, filed on Sep. 19, 2003, incorporated herein by reference in its entirety (hereinafter referred to as the “System and Method for Statistical Timing Analysis of Digital Circuits patent”):  
             A   =       a   0     +       ∑     i   =   1     n     ⁢         a   i     ·   Δ     ⁢           ⁢     X   i         +       a     n   +   1       ⁢   Δ   ⁢           ⁢     R   a                 (   2   )             
 
 where: a o  is a mean value; ΔX i  is a variation of parameter X i , 
 
         [0012]     ΔX i =X i −x o  where x 0  is the mean value of X i ; a i  is the sensitivity of the gate delay to parameter variation ΔX i ; ΔR a  is a random variable responsible for uncorrelated variation of the gate delay; and a n+1  is a sensitivity of the gate delay to uncorrelated variation ΔR a . The following United States Patent Applications, commonly assigned to the assignee herein, are incorporated by reference herein in their entireties: U.S. patent application Ser. No. 10/665,092, entitled “System and Method for Incremental Statistical Timing Analysis of Digital Circuits”, filed on Sep. 18, 2003; U.S. patent application Ser. No. 10/666,353, entitled “System and Method for Statistical Timing Analysis of Digital Circuits”, filed on Sep. 19, 2003; and U.S. patent application Ser. No. 10/666,470, entitled “System and Method for Probabilistic Criticality Prediction of Digital Circuits”, filed on Sep. 19, 2003.  
         [0013]     Sensitivity coefficients express dependency of a gate delay on a process parameter. The last coefficient expresses a fully random part of delay variation. Using unit Gaussian but not arbitrary Gaussian distributions simplifies the analysis but is not restrictive because parameter distributions can always be normalized, e.g., by accumulating their non-zero means into delay means and properly adjusting sensitivity coefficients.  
         [0014]     Block-based STA computes arrival times at each circuit node using two basic operations: incrementing arrival time by a gate delay and computing worst-case arrival time. The incrementing of the arrival time by a gate delay corresponds to propagating signals from a gate input to its output. This operation is performed by summation of the arrival time at the gate input and gate delay. The computing of the worst-case arrival time selects the worst-case signal from the signals arriving at a gate&#39;s inputs. The worst signal can be the latest or the earliest one depending on the type of timing analysis. The computing of the worst-case arrival time is performed by computing the maximum or the minimum of several arrival times. Herein, only the case of computing the latest arrival time is discussed. However, it is to be appreciated that the proposed approach can be easily extended to computation of the earliest arrival time using obvious symmetry between minimum and maximum functions, while maintaining the spirit of the present invention.  
         [0015]     In deterministic timing analysis, these basic operations are trivial. In the statistical case, the situation is more complex because expressions, and not simply numbers, are operated upon.  
         [0016]     Using the first-order canonical form of gate delays and applying some approximation, it is possible to compute signal arrival times at circuit nodes and circuit delay in the same linear canonical form. Incrementing an arrival time by a gate delay given in the first-order canonical form Equation 2 is known. Incrementing the arrival time is performed by summing the corresponding sensitivity coefficients, as described in the System and Method for Statistical Timing Analysis of Digital Circuits patent and, thus, it is not considered herein.  
         [0017]     The computation of the maximum of the two arrival times represented in the first-order form Equation 2 is significantly more difficult because the maximum is a non-linear function. Therefore, parameterized statistical STA computes an approximate first-order representation of the maximum of the two arrival times.  
         [0018]     Linear and Gaussian assumptions are very convenient for statistical parameterized STA because it is possible to use approximate analytical formulae for computing canonical forms of arrival times. Analytical formulae make statistical timing analysis fast and computationally efficient, which is important for implementing a statistical approach in circuit synthesis and optimization. Unfortunately, process parameters may have probability distributions, which are significantly different from Gaussian. It is impossible to approximate certain asymmetric distributions by Gaussian ones with a reasonable error. For example, via resistance distributions are asymmetric. Another problem arises from the fact that some parameters can affect delay in a non-linear way. Usually a linear approximation is justified by assuming that process parameter variation and the corresponding delay variation is sufficiently small. However, with the reduction of transistor sizes, process variation is getting higher. Therefore, a linear approximation of delay dependence on variation is not always sufficiently accurate. For example, delay of a gate depends on transistor channel length non-linearly.  
         [0019]     Accordingly, it would be desirable and highly advantageous to have a system and method for accommodating non-Gaussian and nonlinear sources of variation in statistical static timing analysis.  
       SUMMARY OF THE INVENTION  
       [0020]     These and other drawbacks and disadvantages of the prior art are addressed by the present invention, which is directed to a system and method for accommodating non-Gaussian and non-linear sources of variation in statistical static timing analysis.  
         [0021]     According to an aspect of the present invention, there is provided a system for statistical timing analysis of an electrical circuit. At least one parameter input is for receiving parameters of the electrical circuit. At least one of the parameters has at least one of a non-Gaussian probability distribution and a non-linear delay effect. A statistical static timing analyzer is for calculating at least one of a signal arrival time and a signal required time for the electrical circuit using the at least one parameter. At least one output is for outputting the at least one of the signal arrival time and the signal required time.  
         [0022]     According to another aspect of the present invention, there is provided a method for statistical timing analysis of an electrical circuit. Parameters of the electrical circuit are received. At least one of the parameters has at least one of a non-Gaussian probability distribution and a non-linear delay effect. A statistical static timing analysis is performed to calculate at least one of a signal arrival time and a signal required time for the electrical circuit using the at least one parameter. The at least one of the signal arrival time and the signal required time are output.  
         [0023]     These and other aspects, features and advantages of the present invention will become apparent from the following detailed description of exemplary embodiments, which is to be read in connection with the accompanying drawings. 
     
    
     BRIEF DESCRIPTION OF THE DRAWINGS  
       [0024]     The present invention may be better understood in accordance with the following exemplary figures, in which:  
         [0025]      FIG. 1  is a block diagram illustrating a parameterized statistical static timing analysis tool, according to an illustrative embodiment of the present invention;  
         [0026]      FIG. 2  is a flow diagram illustrating a method for computing parameters of the first-order canonical form C appr  approximating the maximum of two arrival times in first-order canonical form with linear Gaussian parameters, according to an illustrative embodiment of the present invention;  
         [0027]      FIG. 3  is a flow diagram illustrating a method for computing a tightness probability with non-Gaussian parameters, according to an illustrative embodiment of the present invention;  
         [0028]      FIG. 4  is a flow diagram illustrating a method for computing a mean value of the maximum of two canonical forms with non-Gaussian parameters, according to an illustrative embodiment of the present invention;  
         [0029]      FIG. 5  is a flow diagram illustrating a method for computing a tightness probability with non-linear and/or non-Gaussian parameters, according to an illustrative embodiment of the present invention;  
         [0030]      FIG. 6  is a flow diagram illustrating a method for computing the mean value of the maximum of two canonical forms with non-linear and/or non-Gaussian parameters, according to an illustrative embodiment of the present invention; and  
         [0031]      FIG. 7  is a diagram illustrating a grid  700  for numerical integration over a region of non-Gaussian and/or non-linear parameters, according to an illustrative embodiment of the present invention.  
     
    
     DETAILED DESCRIPTION OF PREFERRED EMBODIMENTS  
       [0032]     The present invention is directed to a system and method for accommodating non-Gaussian and nonlinear sources of variation in statistical static timing analysis.  
         [0033]     Digital systems typically work on the basis of the circuit doing some work or computation in each clock cycle or “tick.” In a certain clock cycle, a given signal may not switch, may switch once, or may switch many times, depending on the inputs applied to the circuit. Within each clock cycle, for each and every signal in the system, we are very interested in knowing two things to ensure correct timing: namely the early mode arrival time (or “early arrival time” for short) and the late mode arrival time (or “late arrival time” for short).  
         [0034]     The early mode arrival time is the earliest time at which the signal could possibly switch (i.e., change from the stable logical state at which it was during the previous clock cycle). It is preferred to have the previous cycle “settled down” and recorded correct logic values before the logic values changed or switched during the current clock cycle. The earliest time at which the signal is ALLOWED to switch and still let the circuit work properly is referred to as the “early required time”.  
         [0035]     The late mode arrival time is the latest time at which the signal will stop switching and become stable. It would be beneficial to make sure that the circuit has completed its function in time to meet the clock cycle constraint. The time at which the signal has to stop switching and become stable at its final stable logical state to allow correct functioning of the circuit is called the late required time.  
         [0036]     Timers compute arrival times by “forward propagation” of timing values through a timing graph, and required times by “backward propagation” of timing values through a timing graph.  
         [0037]     The present description illustrates the principles of the present invention. It will thus be appreciated that those skilled in the art will be able to devise various arrangements that, although not explicitly described or shown herein, embody the principles of the invention and are included within its spirit and scope.  
         [0038]     All examples and conditional language recited herein are intended for pedagogical purposes to aid the reader in understanding the principles of the invention and the concepts contributed by the inventor to furthering the art, and are to be construed as being without limitation to such specifically recited examples and conditions.  
         [0039]     Moreover, all statements herein reciting principles, aspects, and embodiments of the invention, as well as specific examples thereof, are intended to encompass both structural and functional equivalents thereof. Additionally, it is intended that such equivalents include both currently known equivalents as well as equivalents developed in the future, i.e., any elements developed that perform the same function, regardless of structure.  
         [0040]     Thus, for example, it will be appreciated by those skilled in the art that the block diagrams presented herein represent conceptual views of illustrative circuitry embodying the principles of the invention. Similarly, it will be appreciated that any flow charts, flow diagrams, state transition diagrams, pseudocode, and the like represent various processes which may be substantially represented in computer readable media and so executed by a computer or processor, whether or not such computer or processor is explicitly shown.  
         [0041]     The functions of the various elements shown in the figures may be provided through the use of dedicated hardware as well as hardware capable of executing software in association with appropriate software. When provided by a processor, the functions may be provided by a single dedicated processor, by a single shared processor, or by a plurality of individual processors, some of which may be shared. Moreover, explicit use of the term “processor” or “controller” should not be construed to refer exclusively to hardware capable of executing software, and may implicitly include, without limitation, digital signal processor (“DSP”) hardware, read-only memory (“ROM”) for storing software, random access memory (“RAM”), and non-volatile storage.  
         [0042]     Other hardware, conventional and/or custom, may also be included. Similarly, any switches shown in the figures are conceptual only. Their function may be carried out through the operation of program logic, through dedicated logic, through the interaction of program control and dedicated logic, or even manually, the particular technique being selectable by the implementer as more specifically understood from the context.  
         [0043]     In the claims hereof, any element expressed as a means for performing a specified function is intended to encompass any way of performing that function including, for example, a) a combination of circuit elements that performs that function or b) software in any form, including, therefore, firmware, microcode or the like, combined with appropriate circuitry for executing that software to perform the function. The invention as defined by such claims resides in the fact that the functionalities provided by the various recited means are combined and brought together in the manner which the claims call for. Applicant thus regards any means that can provide those functionalities as equivalent to those shown herein.  
         [0044]      FIG. 1  is a block diagram illustrating a parameterized statistical static timing analysis tool  102 , according to an illustrative embodiment of the present invention. The first input to the analysis tool  102  is the circuit netlist  100  representing the structure of the circuit to be analyzed. The second input is a set of timing assertions  110 . These typically include arrival times at the primary inputs, required arrival times at the primary outputs, information about the phases of the clock, and details of external loads that are driven by the primary outputs. The assertions can be in the form of deterministic numbers or independent probability distributions or correlated probability distributions. The third input is a set of parameterized delay models  120 . These allow the timer to determine the delay of a gate or wire as a function not only of traditional delay-model variables (like input slew or rise/fall time, and output load) but also as a function of the sources of variation. For example, a first-order linear model may be employed, like the one shown below:  
         delay   =       a   0     +       ∑     i   =   1     n     ⁢       a   i     ⁢   Δ   ⁢           ⁢     x   i         +       a     n   +   1       ⁢   Δ   ⁢           ⁢   R         ,         
 where the delay consists of a deterministic (constant) portion a 0 , a correlated (or global) portion  
         ∑     i   =   1     n     ⁢       a   i     ⁢   Δ   ⁢           ⁢     x   i             
 and an independent (or local) portion a n+1 ΔR. The number of sources of variation is n, and a i ,i= 1 , . . . , n are the sensitivities of the delay to the sources of variation x i ,i=1, . . . , n and a n+1  is the sensitivity to an independent random source of variation R. The notation Δx i  denotes the deviation of x i  from its mean or nominal value, and ΔR denotes the deviation of R from its mean or nominal value. It is to be understood that the delay models can be stored in a pre-characterization step or calculated on the fly as required. The format in which they are stored could include analytical delay equations or table models. The next input is information about the statistics of the sources of variation  130 . This input typically has a list of the sources of variation with a mean value and standard deviation for each source of variation. Any correlations between the sources of variation are specified here. 
 
         [0045]     The probabilistic or statistical static timing module  140  accepts all of these inputs and produces a novel parameterized statistical timing report  150 . This report typically includes arrival times, required arrival times, slacks and slews at each node of the circuit. These timing quantities are not single numbers, but rather they are probability distributions. The information in the timing report can take many forms, including one or more of: mean value and variance for each timing quantity, a parameterized representation of the distribution of each timing quantity, a graphical representation of the distribution of each timing quantity, and a correlation report between these various timing quantities. Various automatic audits such as checking for excessive sensitivities can be built into the timing report. The excessive sensitivities are of interest since they are preferred to be reduced in order to improve the robustness of the circuit. It is to be understood that the circuit being timed can be very large and can consist of millions of gates and wires. All the information mentioned above could be extremely voluminous, and so it is ordinary practice to provide options to selectively output the required information, or even calculate or graphically display all this information on-demand.  
         [0046]     A novel method and apparatus are provided for handling process parameters with non-Gaussian probability distributions and/or affecting gate delay non-linearly. The most obvious way to handle such complications as arbitrary non-Gaussian distributions and non-linear dependences is to apply a purely numerical approach. However, this automatically results in a loss of computational efficiency. On the other hand, the necessity to handle arbitrary distributions and arbitrary non-linear effects is preferable to obtain accurate results. Herein, a combined technique is described which processes linear parameters having Gaussian distributions analytically similar to the prior art parameterized statistical timing analysis. Non-linear and non-Gaussian parameters are processed numerically which allows for avoidance of any restrictions on probability distributions or types of non-linearity. Such an approach is efficient for cases when linear Gaussian approximation is accurate enough for most parameters and only a few parameters demonstrate complex non-linear and/or non-Gaussian behavior. Experiments have shown that it is the most common practical case.  
         [0047]     The present invention extends the prior art parameterized statistical block-based STA proposed in the System and Method for Statistical Timing Analysis of Digital Circuits patent. The present invention uses a generalized version of the linear canonical form for representing gate delays. Moreover, the present invention computes an approximation of all arrival times in the same generalized canonical form, matching the exact mean and standard deviation in a manner completely similar to the original statistical STA. Therefore, the present invention is very convenient for implementation in a statistical timing analyzer. The present invention was implemented by adding several functions, which handle non-linear and non-Gaussian effects without any significant changes of the existing code.  
         [0048]     A description will now be given regarding first-order approximation of the maximum of arrival times depending on Gaussian parameters.  
         [0049]     Assume that arrival times A and B have the first-order representation:  
               A   =       a   0     +       ∑     i   =   1     n     ⁢         a   i     ·   Δ     ⁢           ⁢     X   i         +         a     n   +   1       ·   Δ     ⁢           ⁢     R   a           ⁢     
     ⁢     B   =       b   0     +       ∑     i   =   1     n     ⁢         b   i     ·   Δ     ⁢           ⁢     X   i         +         b     n   +   1       ·   Δ     ⁢           ⁢     R   b                   (   3   )             
 
         [0050]     All parameters variations have Gaussian probability distributions.  
         [0051]     According to the System and Method for Statistical Timing Analysis of Digital Circuits patent, C=max(A,B) can be approximately expressed in the same first-order form:  
               C   appr     =       c   0     +       ∑     i   =   1     n     ⁢         c   i     ·   Δ     ⁢           ⁢     X   i         +         c     n   +   1       ·   Δ     ⁢           ⁢     R   c                 (   4   )             
 
 having the same mean and variance values as actual C=max (A, B). 
 
         [0052]     The values of c 0 , c i  and c n+1  are computed by the following method of  FIG. 2 .  FIG. 2  is a flow diagram illustrating a method  200  for computing parameters of the first-order canonical form C appr  approximating the maximum of two arrival times in first-order canonical form with linear Gaussian parameters, according to an illustrative embodiment of the present invention.  
         [0053]     A start block  202  passes control to a function block  205 . The function block  205  computes variances and covariance of arrival times A and B as follows:  
                   σ   A   2     =       ∑     i   =   1       n   +   1       ⁢     a   i   2         ,       σ   B   2     =       ∑     i   =   1       n   +   1       ⁢     b   i   2         ,     r   =       ∑     i   =   1     n     ⁢       a   i     ⁢     b   i             ⁢     
     ⁢     θ   =         σ   A   2     +     σ   B   2     -     2   ⁢   r                   (   5   )             
 
         [0054]     and passes control to a function block  210 . The function block  210  computes the tightness probability T A  that arrival time A is larger than B, i.e., T A =P(A&gt;B) as follows:  
               T   A     =       ∫       B   -   A     &lt;   0       ⁢         p   ⁡     (       Δ   ⁢           ⁢   X     ,     Δ   ⁢           ⁢     R   a       ,     Δ   ⁢           ⁢     R   b         )       ·           ⁢     ⅆ   Δ       ⁢           ⁢     X   ·     ⅆ   Δ       ⁢           ⁢       R   a     ·     ⅆ   Δ       ⁢           ⁢     R   b                 (   6   )             
 
 where: p(ΔX, ΔR a , ΔR b ) is joint probability density function. For parameters, having Gaussian distribution tightness probability T A  can be computed by the following closed form analytical expressions:  
                 T   A     =       Φ   ⁡     (         a   0     -     b   0       θ     )       ⁢           ⁢   where   ⁢     :         ⁢     
     ⁢         Φ   ⁡     (   y   )       =       ∫     -   ∞     y     ⁢       ϕ   ⁡     (   x   )       ⁢           ⁢     ⅆ   x           ,       ϕ   ⁡     (   x   )       =       1       2   ⁢   π         ⁢     exp   ⁡     (     -       x   2     2       )                     (   7   )             
 
 and passes control to a function block  215 . The function block  215  computes the mean c 0  and variance σ c   2  for C=max(A,B) as follows:  
                 c   0     =         a   0     ⁢     T   A       +       b   0     ⁡     (     1   -     T   A       )       +     θϕ   ⁡     (         a   0     -     b   0       θ     )           ⁢     
     ⁢       σ   C   2     =         (       σ   A   2     +     a   0   2       )     ⁢     T   A       +       (       σ   B   2     +     b   0   2       )     ⁢     (     1   -     T   A       )       +       (       a   0     +     b   0       )     ⁢     θϕ   ⁡     (         a   0     -     b   0       θ     )         -     c   0   2                 (   8   )             
 
 where T A , θ, and φ(x) are defined in the previous step; and passes control to a function block  220 . The function block  220  computes the coefficients c i  (1≦i≦n) as follows: 
 
 c   i   =T   A   a   i +(1 −T   A ) b   i   (9) 
 
 and passes control to a function block  225 . The function block  225  computes the coefficient c n+1  to make the total variance of C appr  the same as the variance of C=max(A,B), and passes control to an end block  230 . 
 
         [0055]     A description will now be given regarding the inventive first-order approximation of the maximum of arrival times depending both on Gaussian and non-Gaussian parameters.  
         [0056]     Assume that arrival times A and B have the first-order representation with parameter variations having both Gaussian and non-Gaussian probability distributions. Gaussian and non-Gaussian parameters in the first-order representations are separated and the following formulae are obtained:  
               A   =       a   0     +       ∑     i   =   1       n   G       ⁢         a     G   ,   i       ·   Δ     ⁢           ⁢     X     G   ,   i           +       ∑     i   =   1       n   N       ⁢         a     N   ,   i       ·   Δ     ⁢           ⁢     X     N   ,   i           +         a     n   +   1       ·   Δ     ⁢           ⁢     R   a           ⁢     
     ⁢     B   =       b   0     +       ∑     i   =   1       n   G       ⁢         b     G   ,   i       ·   Δ     ⁢           ⁢     X     G   ,   i           +       ∑     i   =   1       n   N       ⁢         b     N   ,   i       ·   Δ     ⁢           ⁢     X     N   ,   i           +         b     n   +   1       ·   Δ     ⁢           ⁢     R   b                   (   10   )             
 
 where ΔX G,i  are parameters with Gaussian probability distribution; ΔX N,i  are parameters with non-Gaussian probability distribution; and ΔR a  and ΔR b  are parameters responsible for uncorrelated variation of the arrival times. Uncorrelated variations are assumed to be Gaussian. 
 
         [0057]     The incrementing of the arrival time given in the form of Equation 10 by a gate delay given in the same form is obvious. Similar to the pure Gaussian case as described in the System and Method for Statistical Timing Analysis of Digital Circuits patent, the correspondent coefficients of their first-order forms are summed. A much more difficult problem is to derive an approximate first-order representation for C=max (A, B):  
               C   appr     =       c   0     +       ∑     i   =   1       n   G       ⁢         c     G   ,   i       ·   Δ     ⁢           ⁢     X     G   ,   i           +       ∑     i   =   1       n   N       ⁢         c     N   ,   i       ·   Δ     ⁢           ⁢     X     N   ,   i           +         c     n   +   1       ·   Δ     ⁢           ⁢     R   c                 (   11   )             
 
 having the same mean and variance values as the exact distribution of C=max(A,B). This problem is solved by computing proper values of c 0 , c G,i , c N,i  and c n+1 . 
 
         [0058]     The same idea is used as in the case with pure Gaussian distributions as described in the System and Method for Statistical Timing Analysis of Digital Circuits patent. First, the tightness probability T A  is computed, which is the probability that arrival time A is more than arrival time B. Then, coefficients c G,i  and c N,i  are expressed as linear combinations of coefficients a G,i  and, b G,i  and a N,i  and b N,i  correspondingly: 
 
 c   G,i   =T   A   ·a   G,i +(1 −T   A )· b   G,i   c   N,i   =T   A   ·a   N,i +(1 −T   A )· b   N,i   (12) 
 
         [0059]     Coefficients c 0  and c n+1  are computed so that the mean and variance values for C appr  are the same as for actual C=max (A, B).  
         [0060]     Implementation of this approach requires computing: tightness probabilities T A ; and mean and variance of C=max (A, B).  
         [0061]     These values are computed by separating integration over the subspaces of parameters with Gaussian and non-Gaussian distributions. The integration over the subspace of Gaussian parameters can be done analytically. The integration over the subspace of non-Gaussian parameters is done numerically.  
         [0062]     A description will now be given regarding tightness probability computation.  
         [0063]     For arrival times A and B, having first-order representations as per Equation 10, the tightness probability T A =P(A&gt;B)=P(B−A&lt;0) is defined as follows:  
               T   A     =       ∫       B   -   A     &lt;   U       ⁢         p   ⁡     (       Δ   ⁢           ⁢     X   N       ,     Δ   ⁢           ⁢     X   G       ,     Δ   ⁢           ⁢     R   a       ,     Δ   ⁢           ⁢     R   b         )       ·           ⁢     ⅆ   Δ       ⁢           ⁢       X   N     ·     ⅆ   Δ       ⁢           ⁢       X   G     ·     ⅆ   Δ       ⁢           ⁢       R   a     ·     ⅆ   Δ       ⁢           ⁢     R   b                 (   13   )             
 
 where p(ΔX N , ΔX G , ΔR a , ΔR b ) is a joint probability density function of all the parameter variations. Using the independence of the parameters, this joint probability density function can be decomposed into a product of Gaussian and non-Gaussian density functions as follows: 
 
 p (Δ X   N   , ΔX   G   , ΔR   a   , ΔR   b )= p (Δ X   N )· p (Δ X   G   , ΔR   a   , ΔR   b )  (14) 
 
         [0064]     By substituting expressions in Equation 10 for A and B in the inequality B−A&lt;0, the following inequality is obtained:  
                 b   0     -     a   0     +       ∑     i   =   1       n   N       ⁢           ⁢       (       b     N   ,   i       -     a     N   ,   i         )     ⁢   Δ   ⁢           ⁢     X     N   .   i           +       ∑     i   =   1       n   G       ⁢           ⁢       (       b     G   ,   i       -     a     G   ,   i         )     ⁢   Δ   ⁢           ⁢     X     G   ,   i           +       b     n   +   1       ⁢   Δ   ⁢           ⁢     R   b       -       a     n   +   1       ⁢   Δ   ⁢           ⁢     R   a         &lt;   0           (   15   )             
 
         [0065]     After rearranging terms, the following is obtained:  
                   ∑     i   =   1       n   G       ⁢           ⁢       (       b     G   ,   i       -     a     G   ,   i         )     ⁢   Δ   ⁢           ⁢     X     G   ,   i           +       b     n   +   1       ⁢   Δ   ⁢           ⁢     R   b       -       a     n   +   1       ⁢   Δ   ⁢           ⁢     R   a         &lt;       a   0     -     b   0     +       ∑     i   =   1       n   N       ⁢           ⁢       (       a     N   ,   i       -     b     N   ,   i         )     ⁢   Δ   ⁢           ⁢     X     N   ,   i                     (   16   )             
 
         [0066]     Now this form of the inequality B−A&lt;0 will be used to specify the integration region in Equation 13. Using the representation of probability density function as a product of Gaussian and non-Gaussian parts per Equation 14, the tightness probability of Equation 13 is represented as follows:  
               T   A     =       ∫     -   ∞     ∞     ⁢       p   ⁡     (     Δ   ⁢           ⁢     X   N       )       ⁢     Q   ⁡     (     Δ   ⁢           ⁢     X   N       )       ⁢           ⁢     ⅆ   Δ     ⁢           ⁢     X   N                 (   17   )             
 
         [0067]     where Q(ΔX N ) is the following integral:  
               Q   ⁡     (     Δ   ⁢           ⁢     X   N       )       =       ∫           ∑     i   =   1       n   G       ⁢           ⁢       (       b     G   ,   i       -     a     G   ,   i         )     ⁢   Δ   ⁢           ⁢     X     G   ,   i           -       a     n   +   1       ⁢   Δ   ⁢           ⁢     R   a       +       b     n   +   1       ⁢   Δ   ⁢           ⁢     R   b         &lt;       a   0     -     b   0     +       ∑     i   =   1       n   N       ⁢           ⁢       (       a     N   ,   i       -     b     N   .   i         )     ⁢   Δ   ⁢           ⁢     X     N   ,   i                         ⁢       p   (           ⁢       Δ   ⁢           ⁢     X   G       ,     Δ   ⁢           ⁢     R   a       ,     Δ   ⁢           ⁢     R   b         )     ⁢     ⅆ   Δ     ⁢           ⁢     X   G     ⁢     ⅆ   Δ     ⁢           ⁢     R   a     ⁢     ⅆ   Δ     ⁢           ⁢     R   b                 (   18   )             
 
         [0068]     The integration in Equation 18 is performed in the subspace of parameters with Gaussian distributions only. The integration is performed over the region obtained by the cross section of the region B−A&lt;0 of the whole space with the hyper-plane ΔX N =const, where const represents a vector of values that do not vary during the integration. In other words, Q(ΔX N ) is a conditional probability P(B−A&lt;0|ΔX N =const), i.e., the probability that B−A&lt;0 at the condition of fixed values ΔX N .  
         [0069]     It can be seen that Q(ΔX N ) can be considered as the tightness probability of variables A G  and B G  defined by the following first-order forms:  
                 A   G     =       a     G   ,   0       +       ∑     i   =   1       n   G       ⁢           ⁢       a     G   ,   i       ⁢   Δ   ⁢           ⁢     X     G   ,   i           +       a     n   +   1       ⁢   Δ   ⁢           ⁢     R   a           ⁢     
     ⁢       B   G     =       b     G   ,   0       +       ∑     i   =   1       n   N       ⁢           ⁢       b     G   ,   i       ⁢   Δ   ⁢           ⁢     X     G   ,   i           +       b     n   +   1       ⁢   Δ   ⁢           ⁢     R   b                   (   20   )             
 
 where a G,0  and b G,0  are mean values of A G  and B G :  
                 a     G   ,   0       =       a   0     +       ∑     i   =   1       n   N       ⁢           ⁢       a     N   ,   i       ⁢   Δ   ⁢           ⁢     X     N   ,   i               ⁢     
     ⁢       b     G   ,   0       =       b   0     +       ∑     i   =   1       n   N       ⁢           ⁢       b     N   ,   i       ⁢   Δ   ⁢           ⁢     X     N   ,   i                       (   20   )             
 
         [0070]     Here, ΔX N  is considered not as a vector of random variables but as a vector of parameters having constant values during tightness probability computation. In fact, A G  and B G  are exactly A and B, where ΔX N  are considered as fixed parameters and not random variable.  
         [0071]     Obviously A G  and B G  depend only on parameter variations having Gaussian distributions. Therefore, using the approach of procedure  200  of  FIG. 2  for computing the tightness probability for the Gaussian case, Q(ΔX N ) can be computed analytically in closed form per Equation 7.  
         [0072]     Substituting the analytical expression of Q(ΔX N ) into Equation 17, the tightness probability T A  is computed by numerical integration over the subspace of non-Gaussian parameter variations ΔX N . The numerical integration can be performed using any known method for multidimensional numerical integration. For example, for the simplest numerical integration approach, the tightness probability can be computed by the following method described herein below with respect to  FIG. 3 .  
         [0073]      FIG. 3  is a flow diagram illustrating a method  300  for computing a tightness probability with non-Gaussian parameters, according to an illustrative embodiment of the present invention.  
         [0074]     A start block  302  passes control to a function block  305  that defines a finite region U of non-Gaussian parameter variations where their joint probability density function has non-negligible values, and passes control to a function block  310 . The function block  310  builds an orthogonal grid (uniform or non-uniform) in region U, and passes control to a function block  315 . The number of grid cells controls the accuracy of the procedure: the more the grid cells, the higher the accuracy and the higher the computational cost. The function block  315  sets the tightness probability T A =0, sets a counter of grid cells k=1, and passes control to a function block  320 . The function block  320  computes the average values of non-Gaussian parameters ΔX k,N  for the k-th grid cell, and passes control to a function block  325 . The function block  325  computes the average value of the joint probability  10  density function of non-Gaussian parameters p(ΔX k,N ) for the k-th grid cell, and passes control to a function block  330 . The function block  330  computes variances and covariance of the linear Gaussian part of arrival times A and B as follows:  
         σ     G   ,   A     2     =         ∑     i   =   1       n   G       ⁢           ⁢     a     G   ,   i     2       +     a     n   +   1     2           
         σ     G   ,   B     2     =         ∑     i   =   1       n   G       ⁢           ⁢     b     G   ,   i     2       +     b     n   +   1     2           
         r   G     =       ∑     i   =   1       n   G       ⁢           ⁢       a     G   ,   i       ⁢     b     G   ,   i               
         θ   G     =         σ     G   ,   A     2     +     σ     G   ,   B     2     -     2   ⁢           ⁢     r   G               
 
 and passes control to a function block  335 . The function block  335  computes the tightness probability T A,k  that arrival time A is larger than B, i.e., T A =P(A&gt;B) when non-Gaussian parameters are set to their average values ΔX k,N  in the k-th cell as follows:  
         T     A   ,   k       =         Q   k     ⁡     (     Δ   ⁢           ⁢     X     k   ,   N         )       =     Φ   (         a   0     +       ∑     i   =   1       n   N       ⁢           ⁢       a     N   ,   i       ⁢   Δ   ⁢           ⁢     X     N   ,   i           -     b   0     -       ∑     i   =   1       n   N       ⁢           ⁢       b     N   ,   i       ⁢   Δ   ⁢           ⁢     X     N   ,   i               θ   G       )           
       where:       
           Φ   ⁡     (   y   )       =       ∫     -   ∞     y     ⁢       ϕ   ⁡     (   x   )       ⁢           ⁢     ⅆ   x           ,       ϕ   ⁡     (   x   )       =       1       2   ⁢           ⁢   π         ⁢     exp   ⁡     (       x   2     2     )               
 
 and passes control to a function block  340 . The function block  340  computes the volume ΔV k  of the k-th cell as a product of its dimensions, and passes control to a function block  345 . The function block  345  increments the tightness probability TA by the contribution of the k-th cell as follows: 
 
 T   A   =T   A   +p (Δ X   k,N ) T   A,k   ΔV   k  
 
 and passes control to a function block  350 . The function block  350  increments the value of k by 1 (i.e., k=k+1), and passes control to a decision block  355 . The decision block  355  determines whether the current value of k is less than the number of cells. If the current value of k is less than the number of cells, then control is passed back to function block  320 . Otherwise, if the current value of k is not less than the number of cells, then control is passed to an end block  360 . 
 
         [0075]     The complexity of the method of  FIG. 3  depends only on the number of non-Gaussian parameters. Therefore, it is efficient for the case when we have many Gaussian and only a few non-Gaussian parameters, which is a frequent case in practice of statistical timing analysis.  
         [0076]     A description will now be given regarding mean computation.  
         [0077]     For arrival times A and B having first-order representations (Equation 10) the mean value c 0  of arrival time C=max(A,B) is defined as follows:  
               c   0     =       ∫     -   ∞     ∞     ⁢       max   ⁡     (     A   ,   B     )       ⁢     p   ⁡     (       Δ   ⁢           ⁢     X   N       ,     Δ   ⁢           ⁢     X   G       ,     Δ   ⁢           ⁢     R   a       ,     Δ   ⁢           ⁢     R   b         )       ⁢           ⁢     ⅆ   Δ     ⁢           ⁢     X   N     ⁢     ⅆ   Δ     ⁢           ⁢     X   G     ⁢     ⅆ   Δ     ⁢           ⁢     R   a     ⁢     ⅆ   Δ     ⁢           ⁢     R   b                 (   21   )             
 
         [0078]     p(ΔX N , ΔX G , ΔR a , ΔR b ) is a joint probability density function of all the parameter variations. Using the independence of the parameters, this joint probability density function can be decomposed into a product of Gaussian and non-Gaussian density functions per Equation 14. Using this representation, the mean value of C=max(A,B) can be expressed in the following way:  
               c   0     =       ∫     -   ∞     ∞     ⁢       p   ⁡     (     Δ   ⁢           ⁢     X   N       )       ⁢     (       ∫     -   ∞     ∞     ⁢       max   ⁡     (     A   ,   B     )       ⁢     p   ⁡     (       Δ   ⁢           ⁢     X   G       ,     Δ   ⁢           ⁢     R   a       ,     Δ   ⁢           ⁢     R   b         )       ⁢           ⁢     ⅆ   Δ     ⁢           ⁢     X   G     ⁢           ⁢     ⅆ   Δ     ⁢           ⁢     R   a     ⁢     ⅆ   Δ     ⁢           ⁢     R   b         )     ⁢     ⅆ   Δ     ⁢           ⁢     X   N                 (   22   )             
 
         [0079]     The inner integral is a function of variations with Gaussian distributions only. Denoting it as D(ΔX N ), this equation is transformed into:  
               c   0     =       ∫     -   ∞     ∞     ⁢       D   ⁡     (     Δ   ⁢           ⁢     X   N       )       ⁢     p   ⁡     (     Δ   ⁢           ⁢     X   N       )       ⁢           ⁢     ⅆ   Δ     ⁢           ⁢     X   N                 (   23   )             
 
 where D(ΔX N ) is the following integral:  
               D   ⁡     (     Δ   ⁢           ⁢     X   N       )       =       ∫     -   ∞     ∞     ⁢       max   ⁡     (     A   ,   B     )       ⁢     p   ⁡     (       Δ   ⁢           ⁢     X   G       ,     Δ   ⁢           ⁢     R   a       ,     Δ   ⁢           ⁢     R   b         )       ⁢           ⁢     ⅆ   Δ     ⁢           ⁢     X   G     ⁢           ⁢     ⅆ   Δ     ⁢           ⁢     R   a     ⁢     ⅆ   Δ     ⁢           ⁢     R   b                 (   24   )             
 
         [0080]     Arrival times A and B can be considered as first-order representations A G  and B G  defined by Equation 19 with means defined by Equation 20. The means are functions of parameter variations with non-Gaussian distributions. In Equation 24, similar to the tightness probability computations, ΔX N  are considered to be parameters having a fixed value during the integration process. Then again, A G  and B G  can be considered as functions depending on parameter variations with only Gaussian distributions.  
         [0081]     Integration in Equation 24 does not involve parameter variations with non-Gaussian distributions. Those variations are only external parameters assumed to be constant during the integration. That formula can be rewritten using A G  and B G  instead of A and B as follows:  
               D   ⁡     (     Δ   ⁢           ⁢     X   N       )       =       ∫     -   ∞     ∞     ⁢       max   ⁡     (       A   G     ,     B   G       )       ⁢     p   ⁡     (       Δ   ⁢           ⁢     X   G       ,     Δ   ⁢           ⁢     R   a       ,     Δ   ⁢           ⁢     R   b         )       ⁢           ⁢     ⅆ   Δ     ⁢           ⁢     X   G     ⁢           ⁢     ⅆ   Δ     ⁢           ⁢     R   a     ⁢     ⅆ   Δ     ⁢           ⁢     R   b                 (   25   )             
 
         [0082]     Now it can be seen that it is an expression for the mean of A G  and B G  depending on parameter variations with only Gaussian distributions. Therefore, the results of the procedure  200  of  FIG. 2  are applied and the analytical expression of Equation 8 is used for computing D(ΔX N ). Using that analytical expression for D(ΔX N ), the mean value for max(A,B) is computed by Equation 23 by numerical integration over the subspace of non-Gaussian parameter variations.  
         [0083]     For example, in the simplest numerical integration approach, the mean value of C=max(A,B) is computed by the following method described herein below with respect to  FIG. 4 .  
         [0084]      FIG. 4  is a flow diagram illustrating a method  400  for computing a mean value of the maximum of two canonical forms with non-Gaussian parameters, according to an illustrative embodiment of the present invention.  
         [0085]     A start block  402  passes control to a function block  405 . The function block  405  defines a finite region U of non-Gaussian parameter variations where their joint probability density function has non-negligible values, and passes control to a function block  410 . The function block  410  builds an orthogonal grid (uniform or non-uniform) in region U, and passes control to a function block  415 . The function block  415  sets the mean value c 0 =0, sets a counter of grid cells k=1, and passes control to a function block  420 . The function block  420  computes the average values of non-Gaussian parameters ΔX k,N  for the k-th grid cell, and passes control to a function block  425 . The function block  425  computes the average value of the joint probability density function of non-Gaussian parameters p(ΔX k,N ) for the k-th grid cell, and passes control to a function block  430 . The function block  430  computes variances and covariance of the linear Gaussian part of arrival times A and B as follows  
         σ     G   ,   A     2     =         ∑     i   =   1       n   G       ⁢           ⁢     a     G   ,   i     2       +     a     n   +   1     2           
         σ     G   ,   B     2     =         ∑     i   =   1       n   G       ⁢           ⁢     b     G   ,   i     2       +     b     n   +   1     2           
         r   G     =       ∑     i   =   1       n   G       ⁢       a     G   ,   i       ⁢     b     G   ,   i               
         θ   G     =         σ     G   ,   A     2     +     σ     G   ,   B     2     -     2   ⁢     r   G               
 
 and passes control to a function block  435 . The function block  435  computes the tightness probability T A,k  that arrival time A is larger than B, i.e., T A =P(A&gt;B) when non-Gaussian parameter variations are set to their average values ΔX k,N  in the k-th cell as follows:  
         T     A   ,   k       =         Q   k     ⁡     (     Δ   ⁢           ⁢     X     k   ,   N         )       =     Φ   (         a   0     +       ∑     i   =   1       n   N       ⁢           ⁢       a     N   ,   i       ⁢   Δ   ⁢           ⁢     X     N   ,   i           -     b   0     -       ∑     i   =   1       n   N       ⁢           ⁢       b     N   ,   i       ⁢   Δ   ⁢           ⁢     X     N   ,   i               θ   G       )           
       where:       
           Φ   ⁡     (   y   )       =       ∫     -   ∞     ∞     ⁢       ϕ   ⁡     (   x   )       ⁢           ⁢     ⅆ   x           ,       ϕ   ⁡     (   x   )       =       1       2   ⁢           ⁢   π         ⁢     exp   ⁡     (     -       x   2     2       )               
 
 and control passes to function block  440 . The function block  440  computes the mean value of max(A,B) when non-Gaussian parameter variations are set to their average values ΔX k,N  in the k-th cell as follows:  
         c     0   ,   k       =         D   k     ⁡     (     Δ   ⁢           ⁢     X     k   ,   N         )       =         (       a   0     +       ∑     i   =   1       n   N       ⁢           ⁢       a     N   ,   i       ⁢   Δ   ⁢           ⁢     X     N   ,   i             )     ⁢     T     A   ,   k         +       (       b   0     +       ∑     i   =   1       n   N       ⁢       b     N   ,   i       ⁢   Δ   ⁢           ⁢     X     N   ,   i             )     ⁢     (     1   -     T     A   ,   k         )       +       θ   G     ⁢     ϕ   (         a   0     +       ∑     i   =   1       n   N       ⁢           ⁢       a     N   ,   i       ⁢   Δ   ⁢           ⁢     X     N   ,   i           -     b   0     -       ∑     i   =   1       n   N       ⁢           ⁢       b     N   ,   i       ⁢   Δ   ⁢           ⁢     X     N   ,   i               θ   G       )               
 
 and passes control to a function block  445 . The function block  445  computes the volume ΔV k  of the k-th cell as the product of its dimensions, and passes control to a function block  450 . The function block  450  increments the mean value c 0  by the contribution of the k-th cell as follows: 
 
 c   0   =c   0   +p (Δ X   k,N ) c   0,k   ΔV   k  
 
 and passes control to a function block  455 . The function block  455  increments the value of k by 1 (i.e., k=k+1), and passes control to a decision block  460 . The decision block  460  determines whether the current value of k is less than the number of cells. If the current value of k is less than the number of cells, then control is passed back to function block  420 . Otherwise, if the current value of k is not less than the number of cells, then control is passed to an end block  465 . 
 
         [0086]     The method of mean computation is very similar to the method of tightness probability computation. The major difference is only integration limits of the inner integral which does not depend on the values of non-Gaussian parameters and the integrand function of the inner integral which depends on non-Gaussian parameters.  
         [0087]     A description will now be given regarding variance computation.  
         [0088]     For variance computation, the same technique as for mean computation is used. To begin, the second moment of C=max(A,B) is computed by applying all the transformations described above with respect to the mean computation.  
         [0089]     For arrival times A and B having first-order representations (Equation 10), the second moment of arrival time C=max(A,B) is defined as:  
                 C   ^     2     =       ∫     -   ∞     ∞     ⁢         (     max   ⁡     (     A   ,   B     )       )     2     ⁢     p   ⁡     (       Δ   ⁢           ⁢     X   N       ,     Δ   ⁢           ⁢     X   G       ,     Δ   ⁢           ⁢     R   a       ,     Δ   ⁢           ⁢     R   b         )       ⁢           ⁢     ⅆ   Δ     ⁢           ⁢     X   N     ⁢     ⅆ   Δ     ⁢           ⁢     X   G     ⁢     ⅆ   Δ     ⁢           ⁢     R   a     ⁢     ⅆ   Δ     ⁢           ⁢     R   b                 (   26   )             
 
 where p(ΔX G , ΔX N , ΔR a , ΔR b ) is a joint probability density function for all the parameter variations. The joint probability density function is a product of Gaussian and non-Gaussian density functions per Equation 14. Using this representation, Equation 26 is rewritten as follows:  
                 C   ^     2     =       ∫     -   ∞     ∞     ⁢         p   ⁡     (     Δ   ⁢           ⁢     X   N       )       ·     (       ∫     -   ∞     ∞     ⁢         (     max   ⁡     (     A   ,   B     )       )     2     ⁢     p   ⁡     (       Δ   ⁢           ⁢     X   G       ,     Δ   ⁢           ⁢     R   a       ,     Δ   ⁢           ⁢     R   b         )       ⁢     ⅆ   Δ     ⁢           ⁢     X   G     ⁢     ⅆ   Δ     ⁢           ⁢     R   a     ⁢     ⅆ   Δ     ⁢           ⁢     R   b         )       ⁢     ⅆ   Δ     ⁢           ⁢     X   N                 (   27   )             
 
         [0090]     The inner integral is calculated over a Gaussian sub-space only. The inner integral depends on parameters having non-Gaussian distributions. This equation is rewritten as follows:  
                   C   ^     2     =       ∫     -   ∞     ∞     ⁢       F   ⁡     (     Δ   ⁢           ⁢     X   N       )       ⁢     p   ⁡     (     Δ   ⁢           ⁢     X   N       )       ⁢           ⁢     ⅆ   Δ     ⁢           ⁢     X   N           ⁢     
     ⁢     where:             (   28   )                 F   ⁡     (     Δ   ⁢           ⁢     X   N       )       =       ∫     -   ∞     ∞     ⁢         (     max   ⁡     (     A   ,   B     )       )     2     ⁢     p   ⁡     (       Δ   ⁢           ⁢     X   G       ,     Δ   ⁢           ⁢     R   a       ,     Δ   ⁢           ⁢     R   b         )       ⁢     ⅆ   Δ     ⁢           ⁢     X   G     ⁢     ⅆ   Δ     ⁢           ⁢     R   a     ⁢     ⅆ   Δ     ⁢           ⁢     R   b                 (   29   )             
 
         [0091]     Arrival times A and B can be considered as first-order forms A G  and B G  defined by Equation 19. Their means are defined by Equation 20 as a function of parameter variations with non-Gaussian distributions. A G  and B G  depend on parameter variations with only Gaussian distributions.  
         [0092]     Integration in Equation 29 does not involve parameter variations with non-Gaussian distributions. Non-Gaussian variations are only external parameters assumed to be constant during the integration. This formula is rewritten using A G  and B G  instead of A and B as follows:  
               F   ⁡     (     Δ   ⁢           ⁢     X   N       )       =       ∫     -   ∞     ∞     ⁢         (     max   ⁡     (       A   G     ,     B   G       )       )     2     ⁢     p   ⁡     (       Δ   ⁢           ⁢     X   G       ,     Δ   ⁢           ⁢     R   a       ,     Δ   ⁢           ⁢     R   b         )       ⁢           ⁢     ⅆ   Δ     ⁢           ⁢     X   G     ⁢     ⅆ   Δ     ⁢           ⁢     R   a     ⁢     ⅆ   Δ     ⁢           ⁢     R   b                 (   30   )             
 
         [0093]     It is exactly a formula for the second moment of max(A G ,B G ) where A G  and B G  depend only on Gaussian parameters. It is known that the second moment can be expressed in terms of mean and variance. Using the analytical formulae for mean and variance per Equation 8, an analytical expression for F(ΔX N ) is obtained.  
               F   ⁡     (     Δ   ⁢           ⁢     X   N       )       =         (       σ     G   ,   A     2     +     a     G   ,   0     2       )     ⁢     T     G   ,   A         +       (       σ     G   ,   B     2     +     b     G   ,   0     2       )     ⁢     (     1   -     T     G   ,   A         )       +       (       a     G   ,   0       +     b     G   ,   0         )     ⁢     θ   G     ⁢     ϕ   ⁡     (         a     G   ,   0       -     b     G   ,   0           θ   G       )                   (   31   )             
 
 where: 
 
         [0094]     A G  and B G  are defined by equation 19, a G,0 , b G,0 , σ G,A   2  and σ G,B   2  are means and variances of A G  and B G , T G,A  and T G,B  are tightness probabilities of A G  and B G , and θ G  is a θ parameter for A G  and B G .  
         [0095]     Using this formula for F(ΔX N ), the value of the second moment of C=max (A, B) is computed by Equation 28 performing numerical integration over the space of parameter variations with non-Gaussian distributions.  
         [0096]     Using the second moment and mean, the variance is computed by subtracting the squared mean from the second moment. The method for computing the second moment is similar to the method of  FIG. 4  for computing the mean value. The only difference is using analytical Equation 31 for variance F(ΔX N ) instead of the formula for mean.  
         [0097]     A description will now be given regarding the first-order approximation of the maximum of arrival times depending both on linear and non-linear parameters.  
         [0098]     The same ideas are used here as were applied for handling non-Gaussian parameters hereinabove. The first-order representation of arrival times A and B with non-linear parameters is further generalized as follows:  
               A   =       a   0     +       ∑     i   =   1       n   L       ⁢         a     L   ,   i       ·   Δ     ⁢           ⁢     X     L   ,   i           +       f   A     ⁡     (     Δ   ⁢           ⁢     X   NL       )       +         a     n   +   1       ·   Δ     ⁢           ⁢     R   a           ⁢     
     ⁢     B   =       b   0     +       ∑     i   =   1       n   L       ⁢         b     L   ,   i       ·   Δ     ⁢           ⁢     X     L   ,   i           +       f   B     ⁡     (     Δ   ⁢           ⁢     X   NL       )       +         b     n   +   1       ·   Δ     ⁢           ⁢     R   b                   (   32   )             
 
 where: ΔX L,i  are linear parameters with Gaussian probability distribution; ΔX NL,i  are non linear parameters with arbitrary probability distribution; f A  and f B  are arbitrary functions describing dependence of arrival times on non-linear parameters (these functions can be represented, for example, by tables using a discretization approach); and ΔR a  and ΔR b  are parameters responsible for uncorrelated variation of the arrival times. They are assumed linear and Gaussian. 
 
         [0099]     Incrementing an arrival time given in the form of Equation 12 by a gate delay given in the same form is obvious. Method for Statistical Timing Analysis of Digital Circuits patent, corresponding coefficients of the first-order forms representing arrival time and gate delay are summed. Summation of non-linear terms is not very difficult if they are represented by tables.  
         [0100]     A much more difficult problem is to derive an approximate first-order representation for C=max (A, B):  
               C   appr     =       c   0     +       ∑     i   =   1       n   L       ⁢         c     L   ,   i       ·   Δ     ⁢           ⁢     X     L   ,   i           +       f   c     ⁡     (     Δ   ⁢           ⁢     X   NL       )       +         c     n   +   1       ·   Δ     ⁢           ⁢     R   c                 (   33   )             
 
 having the same mean and variance values as the actual C=max(A,B). It is solved by computing proper values of c 0 , c L,i , c n +1 and a proper function f C . 
 
         [0101]     First, the tightness probability T A  is computed, which is the probability that arrival time A is more than arrival time B. Then, coefficients c L,i  are expressed as linear combinations of coefficients a L,I  and b L,i  as follows: 
 
 c   L,i   =T   A   ·a   L,i +(1 −T   A )·b L,i   ( 34 ) 
 
         [0102]     Similarly, function f c (ΔX NL ) is expressed as a linear combination of functions f A  and f B  as follows: 
 
 f   C (Δ X   NL )= T   A   f   A (Δ X   NL )+(1 −T   A ) f   B (Δ X   NL )  (35) 
 
         [0103]     After that, coefficients c 0  and c n+1  are computed so that the mean and variance values for C appr  are exactly the same as for the actual C=max(A,B).  
         [0104]     Implementation of this approach requires computing the following: tightness probability T A ; and mean and variance of C=max(A,B). These values are computed by separating the integration over the subspace of linear Gaussian parameters ΔX L  and over the subspace of the other parameters. The integration over the subspace of Gaussian parameters can be done analytically. The integration over the subspace of non-linear and/or non-Gaussian parameters is done numerically.  
         [0105]     A description will now be given regarding the tightness probability computation.  
         [0106]     For arrival times A and B, having first-order representations per Equation 32, the tightness probability T A =P(A&gt;B)=P(B−A&lt;0) is defined as:  
               T   A     =       ∫       B   -   A     &lt;   0       ⁢         p   ⁡     (       Δ   ⁢           ⁢     X   NL       ,     Δ   ⁢           ⁢     X   L       ,     Δ   ⁢           ⁢     R   a       ,     Δ   ⁢           ⁢     R   b         )       ·           ⁢     ⅆ   Δ       ⁢           ⁢       X   NL     ·     ⅆ   Δ       ⁢           ⁢       X   L     ·     ⅆ   Δ       ⁢           ⁢       R   a     ·     ⅆ   Δ       ⁢           ⁢     R   b                 (   36   )             
 
 where p(ΔX NL , ΔX L , ΔR a , ΔR b ) is the joint probability density function of all the parameter variations. Using parameter independence, this joint probability density function can be decomposed into a product of two joint density functions 
 
 p (ΔX NL   , ΔX   L   , ΔR   a   , ΔR   b )= p (Δ X   NL )· p (Δ X   L   , ΔR   a   , ΔR   b )  (37) 
 
         [0107]     One of the two joint density functions is a joint probability density function of linear Gaussian parameters, including parameters responsible for uncorrelated variations. The other one is a joint probability density function of all the other parameters.  
         [0108]     Substituting expressions per Equation 32 for A and B in inequality B−A&lt;0 and rearranging terms, the following inequality is obtained:  
                   ∑     i   =   1       n   L       ⁢       (       b     L   ,   i       -     a     L   ,   i         )     ⁢   Δ   ⁢           ⁢     X     L   ,   i           +       b     n   +   1       ⁢   Δ   ⁢           ⁢     R   b       -       a     n   +   1       ⁢   Δ   ⁢           ⁢     R   a         &lt;       a   0     -     b   0     +       f   A     ⁡     (     Δ   ⁢           ⁢     X   NL       )       -       f   B     ⁡     (     Δ   ⁢           ⁢     X   NL       )                 (   38   )             
 
         [0109]     This form of the inequality B−A&lt;0 is used to define the integration region in Equation 36. Using the representation of probability density function as a product of two parts per Equation 37, the tightness probability per Equation 36 is represented as follows:  
               T   A     =       ∫     -   ∞     ∞     ⁢       p   ⁡     (     Δ   ⁢           ⁢     X   NL       )       ⁢     Q   ⁡     (     Δ   ⁢           ⁢     X   NL       )       ⁢           ⁢     ⅆ   Δ     ⁢           ⁢     X   NL                 (   39   )             
 
 where Q(ΔX NL ) is the following integral:  
               Q   ⁡     (     Δ   ⁢           ⁢     X   NL       )       =       ∫           ∑     i   =   1       n   L       ⁢       (       b     L   ,   i       -     a     L   ,   i         )     ⁢   Δ   ⁢           ⁢     X     L   ,   i           -       a     n   +   1       ⁢   Δ   ⁢           ⁢     R   a       +       b     n   +   1       ⁢   Δ   ⁢           ⁢     R   b         &lt;       a   0     -     b   0     +       f   A     ⁡     (     Δ   ⁢           ⁢     X   NL       )       -       f   B     ⁡     (     Δ   ⁢           ⁢     X   NL       )             ⁢       p   ⁡     (       Δ   ⁢           ⁢     X   L       ,     Δ   ⁢           ⁢     R   a       ,     Δ   ⁢           ⁢     R   b         )       ⁢           ⁢     ⅆ   Δ     ⁢           ⁢     X   L     ⁢     ⅆ   Δ     ⁢           ⁢     R   a     ⁢     ⅆ   Δ     ⁢           ⁢     R   b                 (   40   )             
 
         [0110]     The integration here is performed over the subspace of linear Gaussian parameters. It is performed over the region constructed by the cross section of the region B−A&lt;0 of the whole variation space with the hyper-plane ΔX NL =const. “const” means a vector of values that does not vary during the integration. In other words, Q(ΔX NL ) is a conditional probability P(B−A&lt;0|ΔX NL =const), i.e., the probability that B−A&lt;0 at the condition of fixed values ΔX NL .  
         [0111]     It can be seen that Q(ΔX NL ) can be considered as the tightness probability of variables A L  and B L :  
                 A   L     =       a     L   ,   0       +       ∑     i   =   1       n   L       ⁢       a     L   ,   i       ⁢   Δ   ⁢           ⁢     X     L   ,   i           +       a     n   +   1       ⁢   Δ   ⁢           ⁢     R   a           ⁢     
     ⁢       B   L     =       b     L   ,   0       +       ∑     i   =   1       n   L       ⁢       b     L   ,   i       ⁢   Δ   ⁢           ⁢     X     L   ,   i           +       b     n   +   1       ⁢   Δ   ⁢           ⁢     R   b                   (   41   )             
 
 where a L,0  and b L,0  are mean values of A L  and B L : 
 
 a   L,0   =a   0   +f   A (Δ X   NL )  b   L,0   =b   0   +f   B (Δ X   NL )  (42) 
 
         [0112]     Here ΔX NL  is considered not as a vector of random variables but as a vector of parameters having constant values during tightness probability computation. In fact, A L  and B L  are exactly A and B where ΔX NL  are considered as fixed parameters rather than random variables.  
         [0113]     A L  and B L  depend only on linear parameter variations having Gaussian distributions. Therefore, Q(ΔX NL ) is computed analytically using Equation 7.  
         [0114]     Substituting the analytical expression for Q(ΔX NL ) into Equation 39, the tightness probability T A  is computed by numerical integration over the subspace of non-linear and/or non-Gaussian parameters variations ΔX NL . The numerical integration can be performed using any known method for multidimensional numerical integration. For example, in the simple case of a numerical integration approach, the tightness probability can be computed by the following generalization of the method of  FIG. 3  described herein below with respect to  FIG. 5 .  
         [0115]      FIG. 5  is a flow diagram illustrating a method  500  for computing a tightness probability with non-linear and/or non-Gaussian parameters, according to an illustrative embodiment of the present invention.  
         [0116]     A start block  502  passes control to a function block  505 . The function block  505  defines a finite region U of non-linear and/or non-Gaussian parameter variations where their joint probability density function has non-negligible values, and passes control to a function block  510 . The function block  510  builds an orthogonal grid (uniform or non-uniform) in region U, and passes control to a function block  515 . The function block  515  sets the tightness probability T A =0, sets a counter of grid cells k=1, and passes control to a function block  520 . The function block  520  computes the average values of non-linear and/or non-Gaussian parameters ΔX k,NL  for k-th grid cell, and passes control to a function block  525 . The function block  525  computes the average value of the joint probability density function of non-linear and/or non-Gaussian parameters p(ΔX k,NL ) for the k-th grid cell, and passes control to a function block  530 . The function block  530  computes the variances and covariance of the linear Gaussian part of arrival times A and B as follows:  
         σ     L   ,   A     2     =         ∑     i   =   1       n   L       ⁢     a     L   ,   i     2       +     a     n   +   1     2           
         σ     L   ,   B     2     =         ∑     i   =   1       n   L       ⁢     b     L   ,   i     2       +     b     n   +   1     2           
         r   L     =       ∑     i   =   1       n   L       ⁢       a     L   ,   i       ⁢     b     G   ,   i               
         θ   L     =         σ     L   ,   A     2     +     σ     L   ,   B     2     -     2   ⁢     r   L               
 
 and passes control to a function block  535 . The function block  535  computes the tightness probability T A,k  that arrival time A is larger than B, i.e., T A =P(A&gt;B) when non-linear and/or non-Gaussian parameter variations are set equal to their average values ΔX k,NL  in the k-th cell as follows:  
         T     A   ,   k       =         Q   k     ⁡     (     Δ   ⁢           ⁢     X     k   ,   NL         )       =     Φ   ⁡     (         a   0     +       f   A     ⁡     (     Δ   ⁢           ⁢     X   NL       )       -     b   0     -       f   B     ⁡     (     Δ   ⁢           ⁢     X   NL       )           θ   L       )             
       where   ⁢     :         
           Φ   ⁡     (   y   )       =       ∫     -   ∞     y     ⁢       ϕ   ⁡     (   x   )       ⁢           ⁢     ⅆ   x           ,       ϕ   ⁡     (   x   )       =       1       2   ⁢           ⁢   π         ⁢     exp   ⁡     (     -       x   2     2       )               
 
 and passes control to a function block  540 . The function block  540  computes the volume ΔV k  of the k-th cell as the product of its dimensions, and passes control to a function block  545 . The function block  545  increments the tightness probability TA by the contribution of the k-th cell as follows: 
 
 T   A   =T   A   +p (Δ X   k,NL ) T   A,k   ΔV   k  
 
 and passes control to a function block  550 . The function block  550  increments the value of k by 1 (i.e., k=k+1), and passes control to a decision block  555 . The decision block  555  determines whether the current value of k is less than the number of cells. If the current value of k is less than the number of cells, then control is passed back to function block  520 . Otherwise, if the current value of k is not less than the number of cells, then control is passed to an end block  560 . 
 
         [0117]     The complexity of the method of  FIG. 5  depends only on the number of the parameters that are non-linear and/or non-Gaussian. Therefore, the method is very efficient when we have many Gaussian linear parameters and only a few others, which is frequently the case in statistical timing analysis.  
         [0118]     A description will now be given regarding mean computation.  
         [0119]     For arrival times A and B having first-order representations (Equation 32) the mean value c 0  of arrival time C=max(A,B) is defined as follows:  
               c   0     =       ∫     -   ∞     ∞     ⁢       max   ⁡     (     A   ,   B     )       ⁢     p   ⁡     (       Δ   ⁢           ⁢     X   NL       ,     Δ   ⁢           ⁢     X   L       ,     Δ   ⁢           ⁢     R   a       ,     Δ   ⁢           ⁢     R   b         )       ⁢           ⁢     ⅆ   Δ     ⁢           ⁢     X   NL     ⁢     ⅆ   Δ     ⁢           ⁢     X   L     ⁢     ⅆ   Δ     ⁢           ⁢     R   a     ⁢     ⅆ   Δ     ⁢           ⁢     R   b                 (   43   )             
 
 where p(ΔX L , ΔX NL , ΔR a , ΔR b ) is the joint probability density function of all the parameter variations. Using parameter independence, the joint probability density function can be decomposed into a product of two joint probability density functions per Equation 37. Using this representation, the mean value of C=max(A,B) is expressed as follows:  
               c   0     =       ∫     -   ∞     ∞     ⁢       p   ⁡     (     Δ   ⁢           ⁢     X   NL       )       ⁢     (       ∫     -   ∞     ∞     ⁢       max   ⁡     (     A   ,   B     )       ⁢     p   ⁡     (       Δ   ⁢           ⁢     X   L       ,     Δ   ⁢           ⁢     R   a       ,     Δ   ⁢           ⁢     R   b         )       ⁢           ⁢     ⅆ   Δ     ⁢           ⁢     X   L     ⁢           ⁢     ⅆ   Δ     ⁢           ⁢     R   a     ⁢     ⅆ   Δ     ⁢           ⁢     R   b         )     ⁢     ⅆ   Δ     ⁢           ⁢     X   NL                 (   44   )             
 
         [0120]     The inner integral is a function of non-linear and/or non-Gaussian parameters. This equation is transformed into:  
               c   0     =       ∫     -   ∞     ∞     ⁢       D   ⁡     (     Δ   ⁢           ⁢     X   NL       )       ⁢     p   ⁡     (     Δ   ⁢           ⁢     X   NL       )       ⁢           ⁢     ⅆ   Δ     ⁢           ⁢     X   NL                 (   45   )             
 
 where D(ΔX NL ) is the following integral:  
               D   ⁡     (     Δ   ⁢           ⁢     X   NL       )       =       ∫     -   ∞     ∞     ⁢       max   ⁡     (     A   ,   B     )       ⁢     p   ⁡     (       Δ   ⁢           ⁢     X   L       ,     Δ   ⁢           ⁢     R   a       ,     Δ   ⁢           ⁢     R   b         )       ⁢           ⁢     ⅆ   Δ     ⁢           ⁢     X   L     ⁢     ⅆ   Δ     ⁢           ⁢     R   a     ⁢     ⅆ   Δ     ⁢           ⁢     R   b                 (   46   )             
 
         [0121]     Arrival times A and B can be considered to be first-order representations A L  and B L  defined by Equation 41 with means defined by Equation 42. The means are functions of non-linear and/or non-Gaussian parameter variations. Therefore, A L  and B L  can be considered as functions of linear parameters with Gaussian distributions.  
         [0122]     In Equation 46, similarly to the tightness probability computations, ΔX NL  are considered as parameters having fixed values during the integration process. Integration in Equation 46 involves only linear Gaussian parameters. This formula is rewritten using A L  and B L  instead of A and B as follows:  
               D   ⁡     (     Δ   ⁢           ⁢     X   NL       )       =       ∫     -   ∞     ∞     ⁢       max   ⁡     (       A   L     ,     B   L       )       ⁢     p   ⁡     (       Δ   ⁢           ⁢     X   L       ,     Δ   ⁢           ⁢     R   a       ,     Δ   ⁢           ⁢     R   b         )       ⁢           ⁢     ⅆ   Δ     ⁢           ⁢     X   L     ⁢     ⅆ   Δ     ⁢           ⁢     R   a     ⁢     ⅆ   Δ     ⁢           ⁢     R   b                 (   47   )             
 
         [0123]     Now using the fact that A L  and B L  depend only on linear parameter variations with Gaussian distributions, the analytical expression of Equation 8 can be used for computing D(ΔX NL ). Using that analytical expression, the mean value of max(A,D) is computed by Equation 45, performing numerical integration over the space of non-linear and/or non-Gaussian parameters.  
         [0124]     For the simplest numerical integration approach, the tightness probability can be computed by the following generalization of the method of  FIG. 4  described herein after with respect to  FIG. 6 .  
         [0125]      FIG. 6  is a flow diagram illustrating a method for computing the mean value of the maximum of two canonical forms with non-linear and/or non-Gaussian parameters, according to an illustrative embodiment of the present invention.  
         [0126]     A start block  602  passes control to a function block  605 . The function block  605  defines a finite region U of non-linear and/or non-Gaussian parameter variations where their joint probability density function has non-negligible values, and passes control to a function block  610 . The function block  610  builds an orthogonal grid (uniform or non-uniform) in region U, and passes control to a function block  615 . The function block  615  sets the mean value c 0 =0, sets a counter of grid cells k=1, and passes control to a function block  620 . The function block  620  computes the average values of non-linear and/or non-Gaussian parameters ΔX k,NL  for the k-th grid cell, and passes control to a function block  625 . The function block  625  computes the average value of the joint probability density function of non linear and/or non-Gaussian parameters p(ΔX k,NL ), and passes control to a function block  630 . The function block  630  computes the variances and covariance of the linear Gaussian part of arrival times A and B as follows:  
               σ     L   ,   A     2     =         ∑     i   =   1       n   L       ⁢           ⁢     a     L   ,   i     2       +     a     n   +   1     2                     σ     L   ,   B     2     =         ∑     i   =   1       n   L       ⁢           ⁢     b     L   ,   i     2       +     b     n   =   1     2                     r   L     =       ∑     i   =   1       n   L       ⁢           ⁢       a     L   ,   i       ⁢     b     L   ,   i                         θ   L     =         σ     L   ,   A     2     +     σ     L   ,   B     2     -     2   ⁢     r   L                     
 
 and passes control to a function block  635 . The function block  635  computes the tightness probability T A,k  that arrival time A is larger than B, i.e., T A =P(A&gt;B) when non-Gaussian parameter variations are set equal to their average values ΔX k,N  in the k-th cell as follows:  
         T     A   ,   k       =         Q   k     ⁡     (     Δ   ⁢           ⁢     X     k   ,   NL         )       =     Φ   ⁡     (         a   0     +       f   A     ⁡     (     Δ   ⁢           ⁢     X   NL       )       -     b   0     -       f   B     ⁡     (     Δ   ⁢           ⁢     X   NL       )           θ   L       )             
       where   ⁢     :         
           Φ   ⁡     (   y   )       =       ∫     -   ∞     y     ⁢       ϕ   ⁡     (   x   )       ⁢           ⁢     ⅆ   x           ,       ϕ   ⁡     (   x   )       =       1       2   ⁢           ⁢   π         ⁢     exp   ⁡     (     -       x   2     2       )               
 
 and passes control to a function block  640 . The function block  640  computes the mean value of max(A,B) with non-Gaussian parameter variations set equal to their average values ΔX k,NL  in the k-th cell as follows:  
         c     0   ,   k       =         D   k     ⁡     (     Δ   ⁢           ⁢     X     k   ,   NL         )       =         (       a   0     +       f   A     ⁡     (     Δ   ⁢           ⁢     X   NL       )         )     ⁢     T     A   ,   k         +       (       b   0     +       f   B     ⁡     (     Δ   ⁢           ⁢     X   NL       )         )     ⁢     (     1   -     T     A   ,   k         )       +       θ   L     ⁢     ϕ   ⁡     (         a   0     +       f   A     ⁡     (     Δ   ⁢           ⁢     X   NL       )       -     b   0     -       f   B     ⁡     (     Δ   ⁢           ⁢     X   NL       )           θ   L       )                 
 
 and passes control to a function block  645 . The function block  645  computes the volume ΔV k  of the k-th cell as the product of its dimensions, and passes control to a function block  650 . The function block  650  increments the mean value c 0  by the contribution of the k-th cell as follows: 
 
 c   o   =c   o   +p (Δ X   k,NL ) c   0,k   ΔV   k  
 
 and passes control to a function block  655 . The function block  655  increments the value of k by 1 (i.e., k=k+1), and passes control to a decision block  660 . The decision block  660  determines whether the current value of k is less than the number of cells. If the current value of k is less than the number of cells, then control is passed back to function block  620 . Otherwise, if the current value of k is not less than the number of cells, then control is passed to an end block  665 . 
 
         [0127]     The method for mean computation is very similar to the method for the tightness probability computation. The major difference is only the integration limits of the inner integrand which do not depend on the values of non-Gaussian parameters and the integrant function of the inner integral which depends on non-Gaussian parameters.  
         [0128]     A description will now be given regarding variance computation.  
         [0129]     For variance computation, the same technique is used as that for mean computation. Initially, the second moment of C=max (A, B) is computed by applying all the transformations described above with respect to mean computation.  
         [0130]     For arrival times A and B having first-order representations per Equation 32, the second moment of arrival  
                 C   ^     2     =       ∫     -   ∞     ∞     ⁢         (     max   ⁢           ⁢     (     A   ,   B     )       )     2     ⁢     p   ⁡     (       Δ   ⁢           ⁢     X   NL       ,     Δ   ⁢           ⁢     X   L       ,     Δ   ⁢           ⁢     R   a       ,     Δ   ⁢           ⁢     R   b         )       ⁢           ⁢     ⅆ   Δ     ⁢           ⁢     X   NL     ⁢     ⅆ   Δ     ⁢           ⁢     X   L     ⁢     ⅆ   Δ     ⁢           ⁢     R   a     ⁢     ⅆ   Δ     ⁢           ⁢     R   b                 (   48   )             
 
 where p(ΔX NL , ΔX L , ΔR a , ΔR b ) is the joint probability density function for all the parameter variations. Using parameter independence, the joint probability density function is decomposed into a product of two joint probability density functions per Equation 37. Using this representation, Equation 48 is rewritten as follows:  
                 C   ^     2     =       ∫     -   ∞     ∞     ⁢         p   ⁡     (     Δ   ⁢           ⁢     X   NL       )       ·     (       ∫     -   ∞     ∞     ⁢         (     max   ⁡     (     A   ,   B     )       )     2     ⁢     p   ⁡     (       Δ   ⁢           ⁢     X   L       ,     Δ   ⁢           ⁢     R   a       ,     Δ   ⁢           ⁢     R   b         )       ⁢           ⁢     ⅆ   Δ     ⁢           ⁢     X   L     ⁢           ⁢     ⅆ   Δ     ⁢           ⁢     R   a     ⁢     ⅆ   Δ     ⁢           ⁢     R   b         )       ⁢     ⅆ   Δ     ⁢           ⁢     X   NL                 (   49   )             
 
         [0131]     The inner integral is performed over the space of linear Gaussian parameters. This equation is rewritten as follows:  
                 C   ^     2     =       ∫     -   ∞     ∞     ⁢       F   ⁡     (     Δ   ⁢           ⁢     X   NL       )       ⁢     p   ⁡     (     Δ   ⁢           ⁢     X   NL       )       ⁢     ⅆ   Δ     ⁢           ⁢     X   NL                 (   50   )             
 
 where:  
               F   ⁡     (     Δ   ⁢           ⁢     X   NL       )       =       ∫     -   ∞     ∞     ⁢         (     max   ⁡     (     A   ,   B     )       )     2     ⁢     p   ⁡     (       Δ   ⁢           ⁢     X   L       ,     Δ   ⁢           ⁢     R   a       ,     Δ   ⁢           ⁢     R   b         )       ⁢           ⁢     ⅆ   Δ     ⁢           ⁢     X   L     ⁢     ⅆ   Δ     ⁢           ⁢     R   a     ⁢     ⅆ   Δ     ⁢           ⁢     R   b                 (   51   )             
 
         [0132]     Using Equation 41, arrival times A and B can be considered as first-order forms A L  and B L  of linear Gaussian parameters. Their means are defined by Equation 42 as a function of non-linear and/or non-Gaussian parameters.  
         [0133]     Integration in Equation 51 involves only linear parameter variations with Gaussian distributions. Non-linear and non-Gaussian parameters are only external parameters assumed to be constant during the integration. This formula is rewritten using A L  and B L  instead of A and B:  
               F   ⁡     (     Δ   ⁢           ⁢     X   NL       )       =       ∫     -   ∞     ∞     ⁢         (     max   ⁡     (       A   L     ,     B   L       )       )     2     ⁢     p   ⁡     (       Δ   ⁢           ⁢     X   L       ,     Δ   ⁢           ⁢     R   a       ,     Δ   ⁢           ⁢     R   b         )       ⁢           ⁢     ⅆ   Δ     ⁢           ⁢     X   L     ⁢     ⅆ   Δ     ⁢           ⁢     R   a     ⁢     ⅆ   Δ     ⁢           ⁢     R   b                 (   52   )             
 
         [0134]     It is exactly a formula for the second moment of max(A L ,B L ), where A L  and B L  depend only on linear Gaussian parameters. Using the expression of the second moment through mean and variance and applying analytical formulas for mean and variance per Equation 8, the analytical expression for F(ΔX NL ) is obtained as follows:  
               F   ⁢     (     Δ   ⁢           ⁢     X   NL       )       =         (       σ     L   ,   A     2     +     a     L   ,   0     2       )     ⁢     T     L   ,   A         +       (       σ     L   ,   B     2     +     b     L   ,   0     2       )     ⁢     (     1   -     T     L   ,   A         )       +       (       a     L   ,   0       +     b     L   ,   0         )     ⁢     θ   L     ⁢     ϕ   ⁡     (         a     L   ,   0       -     b     L   ,   0           θ   L       )                   (   53   )             
 
 where: 
 
         [0135]     A L  and B L  defined by Equation 41, a L,0 , b L,0 , σ L,A   2  and σ L,B   2  are the means and variances of A L  and B L , T L,A  and T L,B  are the tightness probabilities of A L  and B L , and θL is a θ parameter for A L  and B L .  
         [0136]     Using this formula for F(ΔX NL ), the value of the second moment of C=max(A,B) is computed by Equation 50, performing numerical integration over the space of parameter variations with non-Gaussian distributions.  
         [0137]     Using the second moment and mean, the variance is computed by subtracting the squared mean from the second moment. The method for computing the second moment is similar to the method of  FIG. 6  for computing the mean value.  
         [0138]      FIG. 7  is a diagram illustrating a grid  700  for numerical integration over a region of non-Gaussian and/or non-linear parameters, according to an illustrative embodiment of the present invention. The integration region has two dimensions corresponding to process parameters variations ΔX 1  and ΔX 2 . This example shows a grid consisting of 9 cells  710 . However, for reasonable accuracy, the numerical integration usually uses a finer grid with a larger number of cells.  
         [0139]     While the description of the present invention provided herein focused on the max operation, it is to be appreciated that one of ordinary skill in the art can apply the teaching to a min operation. Moreover, while the description of the present invention provided herein focused on the computation of the latest arrival time, it is to be appreciated that one of ordinary skill in the art can apply the teaching to a computation of the earliest arrival time. Further, while the description of the present invention provided herein focused on combinational circuits, it is to be appreciated that one of ordinary skill in the art can apply the teaching to sequential circuits. Additionally, while the description of the present invention provided herein focused on the case when a gate delay is a separable function, it is to be appreciated that one of ordinary skill in the art can apply the teaching to a case when a gate delay is a non-separable function.  
         [0140]     These and other features, advantages, and variations of the present invention may be readily ascertained by one of ordinary skill in the pertinent art based on the teachings herein. It is to be understood that the teachings of the present invention may be implemented in various forms of hardware, software, firmware, special purpose processors, or combinations thereof.  
         [0141]     Most preferably, the teachings of the present invention are implemented as a combination of hardware and software. Moreover, the software is preferably implemented as an application program tangibly embodied on a program storage unit. The application program may be uploaded to, and executed by, a machine comprising any suitable architecture. Preferably, the machine is implemented on a computer platform having hardware such as one or more central processing units (“CPU”), a random access memory (“RAM”), and input/output (“I/O”) interfaces. The computer platform may also include an operating system and microinstruction code. The various processes and functions described herein may be either part of the microinstruction code or part of the application program, or any combination thereof, which may be executed by a CPU. In addition, various other peripheral units may be connected to the computer platform such as an additional data storage unit and a printing unit.  
         [0142]     It is to be further understood that, because some of the constituent system components and methods depicted in the accompanying drawings are preferably implemented in software, the actual connections between the system components or the process function blocks may differ depending upon the manner in which the present invention is programmed. Given the teachings herein, one of ordinary skill in the pertinent art will be able to contemplate these and similar implementations or configurations of the present invention.  
         [0143]     Although the illustrative embodiments have been described herein with reference to the accompanying drawings, it is to be understood that the present invention is not limited to those precise embodiments, and that various changes and modifications may be effected therein by one of ordinary skill in the pertinent art without departing from the scope or spirit of the present invention. All such changes and modifications are intended to be included within the scope of the present invention as set forth in the appended claims.