Stata程序：最大公约数和最小公倍数

Stata连享会   主页 || 视频 || 推文 || 知乎

温馨提示： 定期 清理浏览器缓存，可以获得最佳浏览体验。

New！ lianxh 命令发布了：
随时搜索连享会推文、Stata 资源，安装命令如下：
  . ssc install lianxh
使用详情参见帮助文件 (有惊喜)：
  . help lianxh

✌ 课程详情： https://gitee.com/arlionn/Course   |   lianxh.cn

✨ 课程主页： https://gitee.com/arlionn/Course

⛳ Stata 系列推文：

• 全部 | Stata入门 | Stata教程 | Stata资源 | Stata命令
• 计量专题 | 论文写作 | 数据分享 | 专题课程
• 结果输出 | Stata绘图 | 数据处理 | Stata程序
• 回归分析 | 面板数据 | 交乘项-调节 | IV-GMM
• 内生性-因果推断 | 倍分法DID | 断点回归RDD | PSM-Matching | 合成控制法
• Probit-Logit | 时间序列 | 空间计量 | 分位数回归 | 生存分析 | SFA-DEA
• 文本分析-爬虫 | Python-R-Matlab | 机器学习
• Markdown | 工具软件 | 其它

作者： 冯超楠 (北京航空航天大学)
邮箱： fengcnhpy@126.com

目录

1. 引言

除了官方命令外，Stata 还允许用户自己编写 ado 程序文件，例如非官方命令 outreg2。ado 文件的编辑也相对简单，只需在 Do-file 编辑器中写入相应代码，然后保存为 .ado 格式。为了让大家对 ado 文件编辑有更好的认识，本文在推文「Stata 程序：10 分钟快乐编写 ado 文件」 的基础上，以计算最大公约数和最小公倍数实例来分享程序编写的具体思路和步骤。

同时，若是家中有正在学「公约数」和「公倍数」概念的小朋友，也可以让他们读读这篇推文，换个角度感受一下，或许小家伙会因此对「经济学」或「计量经济学」产生兴趣。

2. 算法原理

2.1 最大公约数

最大公约数是指两个或多个整数共有约数中最大的一个。在这里，我们将 $a$、$b$ 的最大公约数记为 $(a,b)$，同理，$a$、$b$、$c$ 的最大公约数记为 $(a,b,c)$，依次类推。求最大公约数有多种方法，常见的有：「质因数分解法」、「短除法」、「辗转相除法」、「更相减损法」。考虑到 Stata 循环语句的使用，本文选择「辗转相除法」计算最大公约数。其中，「辗转相除法」具体思路如下：

求 (319, 377)  

∵ 319 ÷ 377 = 0 (余 319)    
∴ (319, 377) = (377, 319)  

∵ 377 ÷ 319 = 1 (余 58)  
∴ (377, 319) = (319, 58)  

∵ 319 ÷ 58 = 5 (余 29)  
∴ (319, 58) = (58, 29)  

∵ 58 ÷ 29 = 2 (余 0)  
∴ (58, 29) = 29  

∴ (319, 377) = 29

对于用「辗转相除法」求 3 个数及以上的最大公约数，可以先求出其中 2 个数的最大公约数，再求这个最大公约数与第 3 个数的最大公约数，依次类推，直到最后一个数为止。并且，最后得到的最大公约数，就是所有数的最大公约数。

2.2 最小公倍数

最小公倍数是指两个或多个整数的公倍数中除 0 外最小的一个。在这里将 $a$、$b$ 的最小公倍数记为 $[a,b]$，同理，$a$、$b$、$c$ 的最小公倍数记为 $[a,b,c]$，依次类推。关于最小公倍数的求法有多种，例如分解质因数法和公式法。考虑到程序运行效率，本文主要采用公式法计算。其中，公式法具体思路如下：

由于最大公约数与最小公倍数间具有如下关系：

所以，求两个数的最小公倍数，就可以先求出它们的最大公约数，然后再用上述公式求出它们的最小公倍数。

求 [18, 20]

[18, 20] = 18 × 20 ÷ (18, 20) = 18 × 20 ÷ 2 = 180

对于求 3 个数及以上的最小公倍数，可以先求出其中 2 个数的最小公倍数，再求这个最小公倍数与第 3 个数的最小公倍数，以此类推，直到最后一个为止。最后得到的最小公倍数，就是所有数的最小公倍数。

3. Stata 实现

3.1 基础代码

3.1.1 公约数：1 个或 2 个自然数的最大公约数

为简便起见，我们先计算 1 个或 2 个自然数的最大公约数。在 Stata 中，可以通过条件命令 if 来判断输入自然数的个数：

• 若为 1 个，则最大公约数是其本身；
• 若为 2 个，则通过 while 循环命令和 mod 取余命令实现「辗转相除法」计算最大公约数；

最后使用 dis 命令在屏幕上显示计算结果。详细代码如下：

* Greatest Common Denominator (GCD) of 2 integers
cap program drop gcd
program define gcd
    if "`2'" == "" {
        dis `1'
    }
    else {
        while `2' {
            local temp2 = `2'
            local 2 = mod(`1',`2')
            local 1 = `temp2'
        }
        dis `1'
    }
end

gcd 100
gcd 5 10

结果如下：

. gcd 100
100
. gcd 5 10
5

3.1.2 存储最大公约数的计算结果

考虑到计算最小公倍数时，需引用最大公约数的计算结果，因此我们要将该结果保存，而非仅在屏幕上显示。此时，可以使用 program define, rclass 命令，返回值将被存储到 r() 中。当然，还可以使用 eclass 和 sclass 命令，结果将分别保存在 e() 和 s() 中。详细代码如下：

* Greatest Common Denominator (GCD) of 2 integers
cap program drop gcd
program define gcd, rclass
    if "`2'" == "" {
        return scalar gcd = `1'
    }
    else {
        while `2' {
            local temp2 = `2'
            local 2 = mod(`1',`2')
            local 1 = `temp2'
        }
        return scalar gcd = `1'
    }
end

gcd 100
dis r(gcd)

gcd 5 10
dis r(gcd)

结果如下：

. gcd 100
. dis r(gcd)
100

. gcd 5 10
. dis r(gcd)
5

3.1.3 公倍数：1 个或 2 个自然数的最小公倍数

通过存储最大公约数的计算结果，我们可以调用 gcd 函数来计算最小公倍数。详细代码如下：

* Least Common Multiple (LCM) of 2 integers
cap program drop lcm
program define lcm, rclass
    if "`2'" == "" {
        return scalar lcm = `1'
    }
    else {
        gcd `1' `2'
        return scalar lcm = `1' * `2' / r(gcd)
    }
end

lcm 100
dis r(lcm)

lcm 5 10
dis r(lcm)

. lcm 100
. dis r(lcm)
100

. lcm 5 10
. dis r(lcm)
10

3.1.4 多个自然数的最大公约数和最小公倍数

在使用「辗转相除法」求几个数的最大公约数时，我们可以先求出其中任意 2 个数的最大公约数，再求这个最大公约数与第 3 个数的最大公约数，以此类推。因此，我们只需循环调用 gcd 函数，最小公倍数也是一样。详细代码如下：

* GCD of arbitrarily long list of integers
cap program drop gcdm
program define gcdm, rclass
    clear results
    foreach i of local 0 {
        gcd `i' `r(gcd)'
    }
    return scalar gcd = r(gcd)
	di "`r(gcd)'"
end

* LCM of arbitrarily long list of integers
cap program drop lcmm
program define lcmm, rclass
    clear results
    foreach i of local 0 {
        lcm `i' `r(lcm)'
    }
    return scalar lcm = r(lcm)
	di "`r(lcm)'"
end

gcdm 5 10 15
lcmm 5 10 15

. gcdm 5 10 15
5

. lcmm 5 10 15
30

3.2 增加 ado 文件功能

3.2.1 若用户输入的参数不是正整数，报错并提示

考虑到用户可能会输入负数或小数导致无法计算最大公约数和最小公倍数，因此我们需要在前文的基础上增加报错机制，即一旦用户输入内容无法用来计算最大公约数和最小公倍数，则报错并提示用户正确输入正整数。以计算最大公约数为例，详细代码如下：

* GCD of arbitrarily long list of integers

* Greatest Common Denominator (GCD) of long list of integers
capture program drop gcdm
program define gcdm,rclass
    clear results
    capture syntax anything(name=n)
    if _rc{  // 用户忘了输入数字
        dis as error "You must enter positive integers"
	exit
    }
 	tokenize "`n'"
	local j = 1
    while "``j''" !=""{     // `j' 表示房间号；``j'' 表示房间中的内容
	local k = ``j''       // 注意暂元的引用方式
	capture confirm integer number `k' // 用户输入了非整数或文字
    if _rc{
        dis as error "You must enter positive integers"
	exit
    }
    if `k'<= 0{  // 用户输入了负数
        dis as error "You must enter positive integers"
	exit
    }

	gcd `k' `r(gcd)'
local j = `j' + 1
	}
	return scalar gcdm = r(gcd)

end

gcdm
gcdm 10.5
gcdm -10
gcdm 100 200 50
dis r(gcdm)

. gcdm
You must enter positive integers

. gcdm 10.5
You must enter positive integers

. gcdm -10
You must enter positive integers

. gcdm 100 200 50

. dis r(gcdm)
50

为了检查用户输入内容是否均为正整数，我们需要对输入内容进行拆分，并逐个检查。在这里，首先使用 capture syntax anything(name=n) 解析命令行，将输入内容命名为 n，并结合 capture 命令检查用户是否输入内容。之后使用 tokenize 命令对字符串 n 进行拆分，默认以空格作为各个部分的分隔符，也可通过选项进行分隔符的指定。随后，使用 capture confirm integer number 或 if 命令检验用户输入了非整数、文字或负数，一旦发生用户输入内容符合以上几种无法计算最大公约数和最小公倍数的情况，则在屏幕上进行报错提示 "You must enter positive integers" 后直接运行 exit 命令。

3.2.2 产生新变量

通过前面的部分，我们已经大致了解了编写程序的思路和基本步骤。接着，我们试图将计算结果赋给指定变量名的新变量。以计算最大公约数为例，具体代码如下：

* Greatest Common Denominator (GCD) of 2 integers
capture program drop gcd
program define gcd,rclass
    if "`2'" == "" {
			return scalar gcd =`1'
	    }

    else {
        while `2' {
            local temp2 = `2'
            local 2 = mod(`1',`2')
            local 1 = `temp2'
         }

			return scalar gcd = `1'
	    }

end

* Greatest Common Denominator (GCD) of 2 integers
capture program drop gcdm
program define gcdm,rclass
    clear results
    capture syntax anything(name=n)[,GENerate(string) replace]
    if _rc{  // 用户忘了输入数字
        dis as error "You must enter positive integers"
	exit
    }
	if "`generate'" != ""{
	    if "`replace'" != "" {
		    capture drop `generate'
		}
         gen `generate' = .

     }

 	tokenize "`n'"
	local j = 1
    while "``j''" !=""{     // `j' 表示房间号；``j'' 表示房间中的内容
	local k = ``j''       // 注意暂元的引用方式
	capture confirm integer number `k' // 用户输入了非整数或文字
    if _rc{
        dis as error "You must enter positive integers"
	exit
    }
    if `k'<= 0{  // 用户输入了负数
        dis as error "You must enter positive integers"
	exit
    }

	gcd `k' `r(gcd)'
	local j = `j' + 1

	}
	if "`generate'" != ""{
    replace `generate'= r(gcd)
	}
	return scalar gcdm = r(gcd)

//	di r(gcdm)

//	di "`r(gcdm)'"

end

clear
set obs 3
gen x = 1
gen mm=1
gcdm 100 200 50 30 20 60 80 60, gen(mm) replace
dis r(gcdm)
list if _n<=1

gcdm 100 200 50,gen(mm) replace
dis r(gcdm)
list if _n<=1

. gcdm 100 200 50 30 20 60 80 60, gen(mm) replace
(3 missing values generated)
(3 real changes made)

. dis r(gcdm)
10

. list if _n<=1

     +--------+
     | x   mm |
     |--------|
  1. | 1   10 |
     +--------+

. 
. gcdm 100 200 50,gen(mm) replace
(3 missing values generated)
(3 real changes made)

. dis r(gcdm)
50

. list if _n<=1

     +--------+
     | x   mm |
     |--------|
  1. | 1   50 |
     +--------+

capture syntax anything(name=n) 后选项 [,GENerate(string) replace] 的添加，允许用户使用命令 generate 生成新变量并存储计算结果，若使用 generate replace 命令则会替换原有的变量。

4. 帮助文档

在编写完一个命令之后，无论是自己还是他人使用，我们都希望能有帮助文档对命令的用法进行说明。由于本文只是为大家简要介绍如何编写 ado 文件，因此编写的帮助命令也是最基础的。以 gcdm 为例，在 Do-file 编辑器中写完帮助文件后，保存为 .sthlp 格式即可。

{smcl}
{* 9 Dec 2020}{...}
{hline}
help for {hi:gcdm}
{hline}

{title:Title}

{p 4 4 2}
{bf:gcdm} —— Compute the Greatest Common Denominator (GCD) of Positive Integers

{title:Syntax}

{p 4 4 2}
{cmdab:gcdm} {anything},
[
{cmdab:g:enerate}{cmd:(}varname{cmd:)} {cmd:replace}
]

{title:Description}

{p 4 4 2}
{cmd:gcdm} makes it easy for users to compute the greatest common denominator (GCD) of positive integers
and the result is stored in r(gcdm).

{title:Options}
{phang}
{cmdab:g:enerate(}{it:varname}{cmd:)} Specifies the name for the greatest common denominator (GCD) variable.
{p_end}
{phang}
{cmd:replace} Specifies that an existing variable is to be overwritten with parameter draws.
{p_end}

{title:Examples}

{p 4 4 2} *- Calculate the greatest common denominator (GCD) of 5 10 15 30. {p_end}
{p 4 4 2}{inp:.} {stata `"gcdm 5 10 15 30"'}{p_end}
{p 4 4 2}{inp:.} {stata `"dis r(gcdm)"'}{p_end}

{p 4 4 2} *- Calculate the greatest common denominator (GCD) of 5 10 15 30 and generate a new variable to store it {p_end}
{p 4 4 2}{inp:.} {stata `"clear"'}{p_end}
{p 4 4 2}{inp:.} {stata `"set obs 3"'}{p_end}
{p 4 4 2}{inp:.} {stata `"gen x = 1"'}{p_end}
{p 4 4 2}{inp:.} {stata `"gcdm 5 10 15 30,gen(gcd_result) replace"'}{p_end}

{title:Author}

{p 4 4 2}
{cmd:Chaonan,Feng}{break}
School of Economics and Management, Beihang University.{break}
E-mail: {browse "mailto:fengcnhpy@126.com":fengcnhpy@126.com}. {break}

5. 命令发布

最后，我们可以将自己编写的命令发布到相关平台上供用户使用，例如 ssc、github、gitee 等，详细过程见「Stata 程序：10 分钟快乐编写 ado 文件」。

6. 参考资料

• Computing greatest common denominators and least common multiples in Stata -Link-
• Manipulating locals via macro lists -Link-

7. 相关推文

Note：产生如下推文列表的 Stata 命令为：
  lianxh 程序, m
安装最新版 lianxh 命令：
  ssc install lianxh, replace

• 专题：Stata命令

  • moremata程序包手动安装方法
• 专题：Stata程序

  • Stata程序：10 分钟快乐编写 ado 文件
  • Stata程序：Monte-Carlo-模拟之产生符合特定分布的随机数
  • Stata程序：我的程序多久能跑完？
  • Stata程序：暂元-(local)-和全局暂元-(global)
  • Stata程序：切割文件路径和文件名
  • Stata程序：在我的程序中接纳另一个程序的所有选项
  • Stata程序：是否有类似-Python-中的-zip()-函数
  • Stata 程序：数值求解极值的基本思想
• 专题：合成控制法

  • Stata：合成控制法程序分享
• 专题：内生性-因果推断

  • Robins - Causal Inference: What if - 数据和程序

相关课程

免费公开课

• 直击面板数据模型 - 连玉君，时长：1小时40分钟，课程主页
• Stata 33 讲 - 连玉君, 每讲 15 分钟.
• Stata 小白的取经之路 - 龙志能，时长：2 小时，课程主页
• 部分直播课 课程资料下载 (PPT，dofiles等)

最新课程-直播课

• Note: 部分课程的资料，PPT 等可以前往 连享会-直播课 主页查看，下载。

⛳ 课程主页

⛳ 课程主页

关于我们

• Stata连享会 由中山大学连玉君老师团队创办，定期分享实证分析经验。
• 连享会-主页 和 知乎专栏，400+ 推文，实证分析不再抓狂。直播间 有很多视频课程，可以随时观看。
• 公众号关键词搜索/回复 功能已经上线。大家可以在公众号左下角点击键盘图标，输入简要关键词，以便快速呈现历史推文，获取工具软件和数据下载。常见关键词：课程, 直播, 视频, 客服, 模型设定, 研究设计, stata, plus, 绘图, 编程, 面板, 论文重现, 可视化, RDD, DID, PSM, 合成控制法 等

连享会小程序：扫一扫，看推文，看视频……

扫码加入连享会微信群，提问交流更方便

✏ 连享会学习群-常见问题解答汇总：
✨ https://gitee.com/arlionn/WD

New！ lianxh 命令发布了：
随时搜索连享会推文、Stata 资源，安装命令如下：
  . ssc install lianxh
使用详情参见帮助文件 (有惊喜)：
  . help lianxh