介绍

  • Coppersmith可以用于求多项式的小根，经常用于 RSA 攻击中“已知某些二进制位，求剩余位”这一类问题。

方程的思想

  • coppersmith攻击主要用到的数学思想就是方程的思想。这边我们就先来举一个例子：
1
2
3
4
5
6
7
8
9
# 这边我们有一个512位的大数p
p = getPrime(512) # 6741438656157143833543832730902889274588187101370523107064278107814792857441703824410652819364741481893721040127690023309393993763833158957274572099181751
# 现在我们只知道这个大数p的高400位
p_high = p >> 112
print(p_high)
# 1298353857614269002737857625343337886712120349999806855420820506922767503261098401301170199662351979411463881678564619275
# 这个时候我们就会使用未知数x表示p的低位,这时我们p就可以表示为方程
# 其中这个x并不会超过p
p = p_high << 112 + x
  • 对于一个阶为e的多项式f可以求得以下满足条件的方程：
    • coppersmith攻击就是，在模n的意义下，快速求出$ n^{\frac{1}{e}} $以内的根，在我们给的例子中就是解方程中的x。
    • 给定$ \beta $ ,快速求出模某个b意义下较小的根，其中b≥nβb ≥ n^\beta
    • 本质上就是解满足条件的方程,我们所需要做的就是去寻找条件和构造方程,本质上还是数学,如果只是要套公式的话,还是和算法没太大关系

相关论文

  • 对于Coppersmith攻击，这个攻击的具体原理与格有关。并且目前为止Coppersmith攻击确实是比较重要的一个攻击。
  • 而Coppersmith攻击之所以被称为Coppersmith攻击，是因为这个攻击方式是Don Coppersmith这个人发表的一篇论文，说明了这个攻击，所以后来这个攻击方式就被称为Coppersmith攻击。
  • 这篇论文虽然发表时间已经过了很久，但还是很值得一读的（之后知识储备足够了，肯定会去读一读这篇论文），这篇论文就先放在这里。
  • 之后其他人在Coppersmith攻击的基础之上提出了二元Coppersmith攻击。论文如下：

sage实现copper

  • 在sage中应用coppersmith定理的函数有两个：small_roots、coppersmith_howgrave_univariate

  • 这边我们只介绍small_roots这个函数，而coppersmith_howgrave_univariate好像在新版的sagemath好像不存在（总之我没找到，之后找到了再写）。

  • 这里就先来介绍一下small_roots的使用。使用的模版如下：

    • 这里n其实就是一个整数
    • PR.<x> = PolynomialRing(Zmod(n))
      • 定义一个在模n下的多项式整环,其实也就是一个多项式。这个多项式有一个特点就是它的系数不会超过n。
      • 并且这个多项式使用x作为变量，就像这样f(x)=anxn+an−1xn−1+...+a1x+a0f(x)=a_nx^n+a_{n-1}x^{n-1}+...+a_1x+a_0
    • 然后f = 就相当于我们定义一个关于x的多项式
    • f = f.monic() ：将这个多项式f转换为首项为1的。注意这一步是必须的，否则最高次数项不是1运行后会报错。
    • 之后就是out = f.small_roots(X=2^n,beta=0.5)寻找最小的根。这边有两个参数，要重点说明一下：
      • 这里的参数X表示限定我们要求的根的上限，上限太高可能就会找到其他解，或者是稍微浪费算力，这边可以在源码的倒数第三行可以发现。
      • 而beta表示的是我们寻找的根的下限。beta必须满足0<=beta<1否则会报错。在下面源码中第156行会有判断。而beta的下限具体作用会展现在源码的最后一行，程序得到的根会在返回之前会先排除小于nβn^{\beta}的解。
      • 总的来说small_roots是求解nβ≤r≤Xn^{\beta}≤r≤X范围内的根。所以有的时候我们就需要调整好我们传递的参数，否则就会导致small_roots返回的是空列表。beta这个参数还是比较奇妙的，有时候nβn^\beta有500位，但是我们还是可以得到200位的根这里先填个坑，之后再来解释一下，总之beta取值一定要在某个正确的范围内，否则会返回空列表(在不了解beta的取值原理时，我们其实都可以采取爆破的策略beta取值其实也就100种可能)
1
2
3
4
5
n = 
PR.<x> = PolynomialRing(Zmod(n)) 
f = 
f = f.monic() 
out = f.small_roots(X=2^n,beta=0.5)

题型1_已知高位

题目1_已知m高位求低位

  • 题目：
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
from Crypto.Util.number import *
import random
import os
r = os.urandom(58)
flag = b'flag{' + r + b'}'
m = bytes_to_long(flag)
print(m.bit_length())
p = getPrime(512)
q = getPrime(512)
n = p*q
e = 3
c = pow(m,e,n)
leak = (m>>72)<<72

print('leak = ', leak)
print('c = ', c)
print('n=',n)
#leak =  5364346700996002419782953699183539964838028506044205579153385131511046817129051651578134266343785810441589599993973488281743331245307218897879563440750592
#c =  40935151362933943549618319861069602822252994206384375794897447336909222247831278964858728027767056236448151125640558855254473877904819185119141154165089443796277880539873361002972028932876254906159097470466137303335170155152431248477323590112940870317475379089257767219445763033224635983039863076611351104336
#n= 117575455116027501138912629496694211426851342451722297844472247041789855520700682789833315967293375516507091916652099874126157193406739919014990279487154654321780811384714209768811626630753881580962486730634616085329659924119107924710021010092412019109403168379014063806220470892151350863781965096437381127661

  • 在这种情况下，由于我们用于加密的m比较大，所以这题就不能使用低加密指数的方法来解。并且这题我们已经知道了m的高440位，这时我们就需要想办法去求m的低72位。我们先利用方程的思想，将m使用方程表达出来。

设x表示m的低72位,而leak表示的是m的高440位这时我们就有m=leak+x会有如下等式c≡m3≡(leak+x)3 mod( n )用同余性质可得m3−c≡(leak+x)3−c≡0 mod( n )此时我们就可以得到关于x的多项式系数,最高次为3经过验证得到n13的bit位为342位，而我们求的x为72位这就说明我们可以使用coppersmith攻击快速求根\begin{array}{l} 设x表示m的低72位,而leak表示的是m的高440位\\ 这时我们就有\\ m = leak + x\\ 会有如下等式\\ c \equiv m^3\equiv(leak+x)^3~mod(~n~)\\ 用同余性质可得\\ m^3-c \equiv (leak+x)^3-c\equiv 0 ~mod(~n~)\\ 此时我们就可以得到关于x的多项式系数,最高次为3\\ 经过验证得到n^{\frac{1}{3}}的bit位为342位，而我们求的x为72位\\ 这就说明我们可以使用coppersmith攻击快速求根 \end{array}

  • 接下来我们就来使用coppersmith求解，这里比较简单的就是，我们在使用small_roots()这个函数的时候，并不需要调整参数。而是直接调用即可。主要是因为我们构造出的方程刚好是(leak+x)3−c≡0 mod( n )(leak+x)^3-c\equiv 0 ~mod(~n~)，所以这题使用coppersmith攻击主要是求在模n的意义下，$ n^{\frac{1}{e}} $以内的根
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
from Crypto.Util.number import *
import libnum 
leak =  5364346700996002419782953699183539964838028506044205579153385131511046817129051651578134266343785810441589599993973488281743331245307218897879563440750592
c =  40935151362933943549618319861069602822252994206384375794897447336909222247831278964858728027767056236448151125640558855254473877904819185119141154165089443796277880539873361002972028932876254906159097470466137303335170155152431248477323590112940870317475379089257767219445763033224635983039863076611351104336
n= 117575455116027501138912629496694211426851342451722297844472247041789855520700682789833315967293375516507091916652099874126157193406739919014990279487154654321780811384714209768811626630753881580962486730634616085329659924119107924710021010092412019109403168379014063806220470892151350863781965096437381127661
PR.<x> = PolynomialRing(Zmod(n))
m = leak + x
m = m.monic()
M = m((m^3-c).small_roots()[0])
print(M)
print(libnum.n2s(int(M)))

"""
5364346700996002419782953699183539964838028506044205579153385131511046817129051651578134266343785810441589599993973488281743331245309819287155072648189309
b'flag{\x82\x0e\xbf\x1fc\xbac\xcc\xb7n\x10!%(al\xa9\xdc8\xbc\xc9\xeb\xb6\xd5\x12\xc4 &\xe1\x00 X\x8er\x96\xd2\x9f\xf2R\xd2\t,0I\xd9\x94\xe8H\xc8=\x8c\xf7\xa6\x8c\x9fm\xb1\x95}'
"""
  • 目前还不知道coppersmith攻击原理，如果要会用的话，过程其实就是这几步：
    • 第一步，确定多项式环中的模数。
    • 第二步，确定多项式。
    • 第三步，确定coppersmith攻击的条件是否满足。
    • 第四步，调用small_roots()函数求解方程。

题目2_已知p的高位求低位

  • 题目如下：
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
from Crypto.Util.number import *
import libnum
import gmpy2
import uuid
flag = b'flag{' + str(uuid.uuid4()).encode() + b'}'
p = getPrime(512)
q = getPrime(512)
n = p*q
m = bytes_to_long(flag)
e = 65537
c = pow(m,e,n)

print('c = ', c)
print('n = ', n)
print('leak = ', p >> 200)

"""
c =  113586177374556404023302519526754665470570871476581572924010485271609182939576292258309693919183665334593269032771168258676196238354553238207372511638353179203781100360164950899879081096534890623159302893909393987089255020158056119954592555987794050357729528077439336645285541296688606433920178147114273163214
n =  117400480938142315180437859683247427996874734449806782078398939334361106314549758122280645558538068023461810390291406250149045673033058236602564577080271495334738508492592434662896092890777780831354145755712206139360655422348268610020323769810399146777354043835879906538618644988182234583974095999379431384971
leak =  6988520536508181139792486524710394841324279945751362186371062878730048665522475016854515443490
"""
  • 这题的思路与题目一的思路比较像，但是却比题目一简单。因为这题直接设p的低位为未知数即可。这边就给出一张图片

image-20250401004541210

  • 推导如下：

这题我们就可以直接设p的低位为x这样我们就可以表示出pp=leak<<200+x=phigh+x我们可以从题目得到x的大小为200位所以我们规定X(即根的上限)为2200次幂而β取0.13到0.50都可以\begin{array}{l} 这题我们就可以直接设p的低位为x\\ 这样我们就可以表示出p\\ p = leak << 200 + x=p_{high}+x\\ 我们可以从题目得到x的大小为200位\\ 所以我们规定X(即根的上限)为2^{200}次幂\\ 而\beta取0.13到0.50都可以 \end{array}

  • 解密就像这样即可，因为我们这时使用的是coppersmith攻击中，给定$ \beta $ ,快速求出模某个b意义下较小的根，其中b≥nβb ≥ n^\beta。所以我们必须要指定small_roots中的参数
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
from Crypto.Util.number import *
import libnum
c =  113586177374556404023302519526754665470570871476581572924010485271609182939576292258309693919183665334593269032771168258676196238354553238207372511638353179203781100360164950899879081096534890623159302893909393987089255020158056119954592555987794050357729528077439336645285541296688606433920178147114273163214
p3_high = 6988520536508181139792486524710394841324279945751362186371062878730048665522475016854515443490
p3_high = p3_high << 200
n3 = 117400480938142315180437859683247427996874734449806782078398939334361106314549758122280645558538068023461810390291406250149045673033058236602564577080271495334738508492592434662896092890777780831354145755712206139360655422348268610020323769810399146777354043835879906538618644988182234583974095999379431384971
e = 65537
PR.<x>  = PolynomialRing(Zmod(n3))
p = p3_high + x
p = p.monic()
p = p(p.small_roots(X=2^200,beta=0.2))
print(p)
q = n3//int(p)
phi = (p-1)*(q-1)
d = inverse_mod(e,phi)
m = pow(c,d,n3)
print(libnum.n2s(int(m)))

题型2_已知高低位

题目1_已知p高位和q低位

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
from Crypto.Util.number import getPrime, bytes_to_long
from secret import flag

p = getPrime(1024)
q = getPrime(1024)
n = p * q
e = 65537
hint1 = p >> 724
hint2 = q % (2 ** 265)
ct = pow(bytes_to_long(flag), e, n)
print(hint1)
print(hint2)
print(n)
print(ct)

"""
hint1 = 1514296530850131082973956029074258536069144071110652176122006763622293335057110441067910479
hint2 = 40812438243894343296354573724131194431453023461572200856406939246297219541329623
n = 21815431662065695412834116602474344081782093119269423403335882867255834302242945742413692949886248581138784199165404321893594820375775454774521554409598568793217997859258282700084148322905405227238617443766062207618899209593375881728671746850745598576485323702483634599597393910908142659231071532803602701147251570567032402848145462183405098097523810358199597631612616833723150146418889589492395974359466777040500971885443881359700735149623177757865032984744576285054725506299888069904106805731600019058631951255795316571242969336763938805465676269140733371287244624066632153110685509892188900004952700111937292221969
ct = 19073695285772829730103928222962723784199491145730661021332365516942301513989932980896145664842527253998170902799883262567366661277268801440634319694884564820420852947935710798269700777126717746701065483129644585829522353341718916661536894041337878440111845645200627940640539279744348235772441988748977191513786620459922039153862250137904894008551515928486867493608757307981955335488977402307933930592035163126858060189156114410872337004784951228340994743202032248681976932591575016798640429231399974090325134545852080425047146251781339862753527319093938929691759486362536986249207187765947926921267520150073408188188
"""
  • 这题本质上还是已知p的位数，求利用small_roots求解p，其中我们是可以知道p的高300位的。而根据hint2我们其实可以求出p的低265位。推到过程如下：

其实我们已知n=p∗q此时我们设plow为p的低265位,pr为p的剩余位数所以有p=2265∗pr+plow同理得q=2265∗qr+qlow可得到n=p∗q=(2265∗pr+plow)∗(2265∗qr+qlow)展开可以得到如下式子:n=2530∗pr∗qr+2265∗pr∗qlow+plow∗2265∗qr+plow∗qlow这时我们将n模一个2265就会得到n≡plow∗qlow mod (2265)此时我们已知n和qlow,计算出plow就没有问题plow≡n∗qlow−1 mod(2265)\begin{array}{l} 其实我们已知n = p*q \\此时我们设p_{low}为p的低265位,p_r为p的剩余位数 \\所以有p = 2^{265}*p_r + p_{low} \\同理得q = 2^{265}*q_r + q_{low} \\可得到n = p*q=(2^{265}*p_r+p_{low})*(2^{265}*q_r+q_{low}) \\展开可以得到如下式子: \\n=2^{530}*p_r*q_r+2^{265}*p_r*q_{low}+p_{low}*2^{265}*q_r+p_{low}*q_{low} \\这时我们将n模一个2^{265}就会得到 \\n\equiv p_{low}*q_{low}~mod~(2^{265}) \\此时我们已知n和q_{low},计算出p_{low}就没有问题 \\p_{low}\equiv n*q_{low}^{-1}~mod(2^{265}) \end{array}

  • 计算出p的低256位之后，我们就可以列方程了

p的剩余为使用x表示,则p可以列出方程p=phigh∗2724+x∗2265+plow在然后我们大致确定一下x的位数1024−300−265=459\begin{array}{l} p的剩余为使用x表示,则p可以列出方程\\ p = p_{high}*2^{724}+x*2^{265}+p_{low} \\在然后我们大致确定一下x的位数 \\1024-300-265=459 \end{array}

  • 此时我们就可以使用coppersmith进行求解，但是如果对这个方程直接用small_roots求解，是解不出来的，这时我们对2**265位--2**271位进行爆破。同时让small_roots()传入的X参数为459-6，传入的beta值在0.44--0.38之间才能求解出正确答案
  • exp如下：
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
import libnum
from Crypto.Util.number import *
hint1 = 1514296530850131082973956029074258536069144071110652176122006763622293335057110441067910479
hint2 = 40812438243894343296354573724131194431453023461572200856406939246297219541329623
n = 21815431662065695412834116602474344081782093119269423403335882867255834302242945742413692949886248581138784199165404321893594820375775454774521554409598568793217997859258282700084148322905405227238617443766062207618899209593375881728671746850745598576485323702483634599597393910908142659231071532803602701147251570567032402848145462183405098097523810358199597631612616833723150146418889589492395974359466777040500971885443881359700735149623177757865032984744576285054725506299888069904106805731600019058631951255795316571242969336763938805465676269140733371287244624066632153110685509892188900004952700111937292221969
ct = 19073695285772829730103928222962723784199491145730661021332365516942301513989932980896145664842527253998170902799883262567366661277268801440634319694884564820420852947935710798269700777126717746701065483129644585829522353341718916661536894041337878440111845645200627940640539279744348235772441988748977191513786620459922039153862250137904894008551515928486867493608757307981955335488977402307933930592035163126858060189156114410872337004784951228340994743202032248681976932591575016798640429231399974090325134545852080425047146251781339862753527319093938929691759486362536986249207187765947926921267520150073408188188
h_low = (n  * inverse_mod(hint2,2^265)) % 2^265
PR.<x> = PolynomialRing(Zmod(n))
for i in range(64):
    p = hint1*(2^724) + x*(2^265)*64 + i*(2^265) +h_low
    p = p.monic()
    out = p.small_roots(X=2^453,beta = 0.44)
    if len(out) != 0 :
        print(out[0])
        print(i)
        break
p = hint1*(2^724) + out[0]*(2^265)*64 + 19*(2^265) +h_low
q = n//int(p)
phi = (p-1)*(q-1)
e = 65537
d = inverse_mod(e,phi)
m = pow(ct,d,n)
print(long_to_bytes(m))
"""
7914249681483699348105896117503653090559709096541336321100288778319519769467957659191305141838207370902177955576606946702942703171837148
19
b'flag{ef5e1582-8116-4f61-b458-f793dc03f2ff}'
"""

题型?

small_roots函数源码

  • 函数源码如下：
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
def small_roots(self, X=None, beta=1.0, epsilon=None, **kwds):
    r"""
    Let `N` be the characteristic of the base ring this polynomial
    is defined over: ``N = self.base_ring().characteristic()``.
    This method returns small roots of this polynomial modulo some
    factor `b` of `N` with the constraint that `b >= N^\beta`.
    Small in this context means that if `x` is a root of `f`
    modulo `b` then `|x| < X`. This `X` is either provided by the
    user or the maximum `X` is chosen such that this algorithm
    terminates in polynomial time. If `X` is chosen automatically
    it is `X = ceil(1/2 N^{\beta^2/\delta - \epsilon})`.
    The algorithm may also return some roots which are larger than `X`.
    'This algorithm' in this context means Coppersmith's algorithm for finding
    small roots using the LLL algorithm. The implementation of this algorithm
    follows Alexander May's PhD thesis referenced below.

    INPUT:

    - ``X`` -- an absolute bound for the root (default: see above)
    - ``beta`` -- compute a root mod `b` where `b` is a factor of `N` and `b
      \ge N^\beta`. (Default: 1.0, so `b = N`.)
    - ``epsilon`` -- the parameter `\epsilon` described above. (Default: `\beta/8`)
    - ``**kwds`` -- passed through to method :meth:`Matrix_integer_dense.LLL() <sage.matrix.matrix_integer_dense.Matrix_integer_dense.LLL>`.

    EXAMPLES:

    First consider a small example::

        sage: N = 10001
        sage: K = Zmod(10001)
        sage: P.<x> = PolynomialRing(K, implementation='NTL')
        sage: f = x^3 + 10*x^2 + 5000*x - 222

    This polynomial has no roots without modular reduction (i.e. over `\ZZ`)::

        sage: f.change_ring(ZZ).roots()
        []

    To compute its roots we need to factor the modulus `N` and use the Chinese
    remainder theorem::

        sage: p, q = N.prime_divisors()
        sage: f.change_ring(GF(p)).roots()
        [(4, 1)]
        sage: f.change_ring(GF(q)).roots()
        [(4, 1)]

        sage: crt(4, 4, p, q)
        4

    This root is quite small compared to `N`, so we can attempt to
    recover it without factoring `N` using Coppersmith's small root
    method::

        sage: f.small_roots()
        [4]

    An application of this method is to consider RSA. We are using 512-bit RSA
    with public exponent `e=3` to encrypt a 56-bit DES key. Because it would be
    easy to attack this setting if no padding was used we pad the key `K` with
    1s to get a large number::

        sage: Nbits, Kbits = 512, 56
        sage: e = 3

    We choose two primes of size 256-bit each::

        sage: p = 2^256 + 2^8 + 2^5 + 2^3 + 1
        sage: q = 2^256 + 2^8 + 2^5 + 2^3 + 2^2 + 1
        sage: N = p*q
        sage: ZmodN = Zmod( N )

    We choose a random key::

        sage: K = ZZ.random_element(0, 2^Kbits)

    and pad it with `512-56=456` 1s::

        sage: Kdigits = K.digits(2)
        sage: M = [0]*Kbits + [1]*(Nbits-Kbits)
        sage: for i in range(len(Kdigits)): M[i] = Kdigits[i]

        sage: M = ZZ(M, 2)

    Now we encrypt the resulting message::

        sage: C = ZmodN(M)^e

    To recover `K` we consider the following polynomial modulo `N`::

        sage: P.<x> = PolynomialRing(ZmodN, implementation='NTL')
        sage: f = (2^Nbits - 2^Kbits + x)^e - C

    and recover its small roots::

        sage: Kbar = f.small_roots()[0]
        sage: K == Kbar
        True

    The same algorithm can be used to factor `N = pq` if partial
    knowledge about `q` is available. This example is from the
    Magma handbook:

    First, we set up `p`, `q` and `N`::

        sage: length = 512
        sage: hidden = 110
        sage: p = next_prime(2^int(round(length/2)))
        sage: q = next_prime(round(pi.n()*p))                                           # needs sage.symbolic
        sage: N = p*q                                                                   # needs sage.symbolic

    Now we disturb the low 110 bits of `q`::

        sage: qbar = q + ZZ.random_element(0, 2^hidden - 1)                             # needs sage.symbolic

    And try to recover `q` from it::

        sage: F.<x> = PolynomialRing(Zmod(N), implementation='NTL')                     # needs sage.symbolic
        sage: f = x - qbar                                                              # needs sage.symbolic

    We know that the error is `\le 2^{\text{hidden}}-1` and that the modulus
    we are looking for is `\ge \sqrt{N}`::

        sage: from sage.misc.verbose import set_verbose
        sage: set_verbose(2)
        sage: d = f.small_roots(X=2^hidden-1, beta=0.5)[0]  # time random               # needs sage.symbolic
        verbose 2 (<module>) m = 4
        verbose 2 (<module>) t = 4
        verbose 2 (<module>) X = 1298074214633706907132624082305023
        verbose 1 (<module>) LLL of 8x8 matrix (algorithm fpLLL:wrapper)
        verbose 1 (<module>) LLL finished (time = 0.006998)
        sage: q == qbar - d                                                             # needs sage.symbolic
        True

    REFERENCES:

    Don Coppersmith. *Finding a small root of a univariate modular equation.*
    In Advances in Cryptology, EuroCrypt 1996, volume 1070 of Lecture
    Notes in Computer Science, p. 155--165. Springer, 1996.
    http://cr.yp.to/bib/2001/coppersmith.pdf

    Alexander May. *New RSA Vulnerabilities Using Lattice Reduction Methods.*
    PhD thesis, University of Paderborn, 2003.
    http://www.cs.uni-paderborn.de/uploads/tx_sibibtex/bp.pdf
    """
    from sage.misc.verbose import verbose
    from sage.matrix.constructor import Matrix
    from sage.rings.real_mpfr import RR

    N = self.parent().characteristic()

    if not self.is_monic():
        raise ArithmeticError("Polynomial must be monic.")

    beta = RR(beta)
    if beta <= 0.0 or beta > 1.0:
        raise ValueError("0.0 < beta <= 1.0 not satisfied.")

    f = self.change_ring(ZZ)

    P,(x,) = f.parent().objgens()

    delta = f.degree()

    if epsilon is None:
        epsilon = beta/8
    verbose("epsilon = %f"%epsilon, level=2)

    m = max(beta**2/(delta * epsilon), 7*beta/delta).ceil()
    verbose("m = %d"%m, level=2)

    t = int( ( delta*m*(1/beta -1) ).floor() )
    verbose("t = %d"%t, level=2)

    if X is None:
        X = (0.5 * N**(beta**2/delta - epsilon)).ceil()
    verbose("X = %s"%X, level=2)

    # we could do this much faster, but this is a cheap step
    # compared to LLL
    g  = [x**j * N**(m-i) * f**i for i in range(m) for j in range(delta) ]
    g.extend([x**i * f**m for i in range(t)]) # h

    B = Matrix(ZZ, len(g), delta*m + max(delta,t) )
    for i in range(B.nrows()):
        for j in range( g[i].degree()+1 ):
            B[i,j] = g[i][j]*X**j

    B =  B.LLL(**kwds)

    f = sum([ZZ(B[0,i]//X**i)*x**i for i in range(B.ncols())])
    R = f.roots()

    ZmodN = self.base_ring()
    roots = set([ZmodN(r) for r,m in R if abs(r) <= X])
    Nbeta = N**beta
    return [root for root in roots if N.gcd(ZZ(self(root))) >= Nbeta]