小球碰撞（校内小测）—— 浅谈利用扩展欧几里得算出一元二次方程

最新推荐文章于 2021-05-14 20:16:22 发布

PI_RE

最新推荐文章于 2021-05-14 20:16:22 发布

阅读量644

点赞数 3

本文链接：https://blog.csdn.net/weixin_43925239/article/details/87970780

版权

前言

在寒假集训中我学到了不少的新东西，也得到了许多感悟（~~当然还有很多天书一样的题~~）。

在集训的末尾，我们也是迎来了 ~~情理之中~~ 的检测考试（~~虽然还是一道都不会~~）。不过呢，在之后也是努力钻研了其中的几道题，于是乎，我也决定和大家一起分享一下，共同来探讨这一道 ~~堪称最简单的~~ 题（~~其实也就是因为这么久没写博客了来水一个而已~~）。

当然，通过某位不知名的da lao点拨，在这次及以后的博客中，我都不会像以前那样大片大片的讲解理论知识或者干巴巴的分析代码（我也相信很多人都像我一样当看到篇幅太多的时候就已经没有心情再看了），我决定只挑一些很重要的的来讲一下，其余的我相信大家的理解能力哈（手动滑稽），如果真的没看懂可以下去理解一下，第二天再来看（继续滑稽），或者通过评论来告诉我，我会竭尽所能一一解答的。

题目

题目描述

输入

输出

样例输入

样例1输入：
2 3 1


样例2输入：
232 320 34272

样例输出

样例1输出：
2.5


样例2输出：
1102656

我可以这么说，起码有80%的童鞋们像我一样，在看第一遍的时候绝对没有看懂题目是什么意思（da lao自动过滤）。

当然这也在情理之中嘛（我可能要被打诶），毕竟这其中包含了许多物理学和数学的高级知识（至少我是这么认为的）。

题目解析

好了不跟大家卖关子了，其实题目看起来这么多，真正有用的却在“物理知识补充”的最后一句话，让我来用大家听得懂的人话来解释一下，就是说给了你 a, b, p，让你求 $min \left \{ \frac{1}{2} * (a * va^{2} + b * vb^{2}) \right \}$ ，当然，这个va 和vb 可不是随便什么数都行，它必须满足 (a * va + b * vb) = p 。

更简单的来说，也就是求一个方程的解经过某种神秘的计算后使其结果更小，这样的话总能看懂了吧。

解析

在上面我们已经剖析开了这道题的实际意义，可是应该怎么解这么一个方程呢？（或许有人会通过一些很神奇的方法算出来，但那毕竟不是利用电脑求通解吧）在这里，我推荐一个很流弊厉害的方法，也正是这次的题目——拓展欧几里得方法。

扩展欧几里得

或许大家都知道欧几里得算法（也就是辗转相除法）他是通过不断取模得到一个最大公因数。可是有多少知道扩展欧几里得呢？其实这个是欧几里得的高级版，他不仅能够得到最大公因数，而且还可以保存一组特解。

在这里插入一个小知识，如果一个方程 ax + by = c 有解的话，那么就有一个条件： gcd (a, b) | c ，这个我们可以一起来推导一下：

我们可以令 (a, b) （就是a和b的最大公因数） = d，那么方程可以变成这样：k1 * d * x + k2 * d * y = c(k1 和k2 互质)

再合并同类项，就可以得到 d(k1 * x + k2 * y) = c。

再把d移过去，就变成了k1 * x + k2 * y = c / d。又因为k1、k2、x、y都是整数，那么必定等式右边也是整数，也就是说d | c。

当然，这个只是其中的一个推理，所有的扩欧还需要有 欧拉函数 的知识来作为支撑，那么因为篇幅原因（~~同时我可能讲不明白~~），所以这个具体的扩欧大家可以下去百度一下，我这里只最后甩出一个小的自己认为的推理：

有方程ax + by = c，(a, b) = 1，那么可以通过欧拉函数解出 $ax\equiv 1$ (mod b)的解x1（这个同余方程可以变成ax + by = 1）。

那么最开始的那个方程ax + by = c的解x和y就是x1 * c 和 y1 * c。

那么这个推理就可以跟上面的衔接上了，大家自己可以反复看一下这两个推理，希望对大家有用。

现给出扩欧实现代码：

int exgcd (int a, int b, int &x, int &y){    //x和y分别存的一组特解
    if (b == 0){
        x = 1;
        y = 0;
        return a;
    }
    int gcd = exgcd (b, a % b, y, x);
    y -= x * (a / b);
    return gcd;
}

那么有了这个代码，我们就能够求出gcd和一组特解了。

我们令通过扩欧求出的特解为va1和vb1，那么原方程的特解va0就等于va1 * (p / d)，vb0 * (p / d)。而方程的通解就为 $\begin{Bmatrix} va = va1 +\frac{mb}{d} * t\\ vb = vb1 - \frac{ma}{d} * t \end{Bmatrix}$ （实在不知道怎么把右半边的括号删去，大家凑合着看吧）

那么E就等于 $\frac{ma(va0 + \frac{mb}{d}t)^{2} + mb (vb0 - \frac{ma}{d}t)^{2}}{2}$ $\frac{ma(va0 + \frac{mb}{d}t)^{2} + mb (vb0 - \frac{ma}{d}t)^{2}}{2}$ 。然后最后再把这两个平方拆开，就能够得到它： $\frac{\frac{ma*mb(ma + mb)}{d^{2}}t^{2} + \frac{2ma*mb(va0 - vb0)}{d}t + ma*va0^{2} + mb * vb0^{2}}{2}$ 。

可是的话，得到这一大串有什么用呢？或许会有很多人跟我一样在最开始的时候有这种疑问。But，我们在仔细看看，这里面的字母，是不是除了t之外，我们全都已经得到了？也就是说，除了t以外，其他都可以看成常量了对吧。

再进一步说，这个E，也就是一个关于t的函数！，它随着t的变化而变化。

从这个式子中，我们似乎得到了所有的信息，那么也应该继续往下走了。

二次函数

当然，大家千万别忘了数学补充的那句黑体重要的话：

二次函数的对称轴为 $-\frac{b}{2a}$ 。

这句话其实非常的重要！有思路比较快的童鞋可能已经想到了。没想到的也没关系，我们一起来看：

上面已经说过，E是关于t的函数，同时，题目也要我们求出E的最小值，这么说来，只需要t最小，那么E也就会随之最小了吧。

参考代码

#include <cstdio>
#include <cmath>
#include <cstring>
#include <cstdlib>
#include <algorithm>
#define min(a, b) a < b ? a : b
#define ld double
#define LL long long
 
void read (LL &x) {
    x = 0;
    LL f = 1;
    char c = getchar ();
    while (c < '0' || c > '9'){
        if (c == '-') f = -1;
        c = getchar ();
    }
    while (c >= '0' && c <= '9'){
        x = (x << 1) + (x << 3) + c - 48;
        c = getchar ();
    }
    x *= f;
}
 
void print (LL x){
    if (x < 0){
        putchar ('-');
        x = ~(x) + 1;
    }
    if (x / 10) print (x / 10);
    putchar (x % 10 + 48);
}
 
LL ma, mb, p, va0, vb0, x, y, d;
 
ld js (){
    ld a = (1.0 * ma * mb * (ma + mb)) / (1.0 *(d * d));
    ld b = (2.0 * ma * mb * (va0 - vb0)) / (1.0 * d);
    return -b / (2 * a);
}
 
ld js1 (LL t){
    ld a = (1.0 * ma * mb * (ma + mb)) / (1.0 * d * d) * (1.0 * t * t),
       b = (2.0 * ma * mb * (va0 - vb0)) / (1.0 * d) * (1.0 * t),
       c = (1.0 * ma * va0 * va0 + 1.0 * mb * vb0 * vb0);
    return 0.5 * (a + b + c);
}
 
LL exgcd (LL a, LL b, LL &x, LL &y){
    if (!b){
        x = 1, y = 0;
        return a;
    }
    LL k = exgcd (b, a % b, y, x);
    y -= x * (a / b);
    return k;
}
 
int main (){
    read (ma); read (mb); read (p);
    d = exgcd (ma, mb, x, y);
    if (p % d){
        print (-1);
        putchar (10);
        return 0;
    }
    va0 = x * (p / d), vb0 = y * (p / d);
    LL t = floor (js ()), t1 = t + 1;
    ld E = min (js1 (t), js1 (t1));
    /*if ( (E) == E) print ( (E));
    else printf ("%.1lf", E);*/
    if (floor (E) == E)
        printf ("%.0lf", E);
    else printf ("%.1lf", E);
    putchar (10);
}

这确确实实是一种计算方式，而且看起来也并没有问题，思路也是严格按照上面的推理写成的。

可是经过一番调试（加上老师的数据）后，我发现了一个问题，而且这个问题出在js1上。

不知道为什么，每次调试时都可以看出来里面的局部变量a、b、c出问题了。而且据我所推断，可能是爆空间了？？！但是因为调试器的查看界面的原因（~~加上自身……~~），我也并不能很确定到底发生了些什么，所以如果有da lao看到了（并且有兴趣），可以帮本蒟调试一下，可以私信，当然也可以评论（在此谢谢dalao了）。这是那组特殊数据：24 52 400000000。

不过在换一种方式以后，竟然就奇妙的ac了（同学们都解释不出来原因是什么。。。），所以还是那句诚恳地话，如果有dalao可以在评论回复一下，在此感激不尽！！！

#include <cstdio>
#include <cmath>
#include <cstring>
#include <cstdlib>
#include <algorithm>
#define min(a, b) a < b ? a : b
#define ld double
#define LL long long
 
void read (LL &x) {
    x = 0;
    LL f = 1;
    char c = getchar ();
    while (c < '0' || c > '9'){
        if (c == '-') f = -1;
        c = getchar ();
    }
    while (c >= '0' && c <= '9'){
        x = (x << 1) + (x << 3) + c - 48;
        c = getchar ();
    }
    x *= f;
}
 
void print (LL x){
    if (x < 0){
        putchar ('-');
        x = ~(x) + 1;
    }
    if (x / 10) print (x / 10);
    putchar (x % 10 + 48);
}
 
LL ma, mb, p, va0, vb0, x, y, d;
 
ld js (){    //这里有两种方法，一种是把括号去掉化简了的，不过都是对的
    /*ld a = (1.0 * ma * mb * (ma + mb)) / (1.0 *(d * d));
    ld b = (2.0 * ma * mb * (va0 - vb0)) / (1.0 * d);
    return -b / (2 * a);*/
    return 1.0 * (vb0 - va0) * d / (1.0 * (ma + mb));
}
 
ld js1 (LL t){
    ld a = (1.0 * ma * mb * (ma + mb)) / (1.0 * d * d) * (1.0 * t * t),
       b = (2.0 * ma * mb * (va0 - vb0)) / (1.0 * d) * (1.0 * t),
       c = (1.0 * ma * va0 * va0 + 1.0 * mb * vb0 * vb0);
    return 0.5 * (a + b + c);
}
 
LL exgcd (LL a, LL b, LL &x, LL &y){
    if (!b){
        x = 1, y = 0;
        return a;
    }
    LL k = exgcd (b, a % b, y, x);
    y -= x * (a / b);
    return k;
}
 
int main (){
    read (ma); read (mb); read (p);
    d = exgcd (ma, mb, x, y);
    if (p % d){
        print (-1);
        putchar (10);
        return 0;
    }
    va0 = x * (p / d), vb0 = y * (p / d);
    LL t = floor (js ()), t1 = t + 1;
    ld E1 = (1ll * ma * (va0 + t * mb / d) * (va0 + t * mb / d) + mb * (vb0 - t * ma / d) * (vb0 - t * ma / d)) / 2.0,
       E2 = (1ll * ma * (va0 + t1 * mb / d) * (va0 + t1 * mb / d) + mb * (vb0 - t1 * ma / d) * (vb0 - t1 * ma / d)) / 2.0;
    if (E1 < E2){
        if (floor (E1) == E1) printf ("%.0lf", E1);
        else printf ("%.1lf", E1);
    }
    else{
        if (floor (E2) == E2) printf ("%.0lf", E2);
        else printf ("%.1lf", E2);
    }
    /*ld E = min (js1 (t), js1 (t1));
    if ( (E) == E) print ( (E));
    else printf ("%.1lf", E);
    if (floor (E) == E)
        printf ("%.0lf", E);
    else printf ("%.1lf", E);*/
    putchar (10);
}

后记

在真正搞懂以后好像确实不是特别难对吧（只是融合了许多只是而已）。

所以可以看出来，只要认真搞懂，我相信在座的大家（当然也包括我啦哈哈）都能成为那“遥不可及”的dalao哒！！

PI_RE

关注

3
点赞
踩
0

收藏

觉得还不错? 一键收藏
1
评论
小球碰撞（校内小测）—— 浅谈利用扩展欧几里得算出一元二次方程

目录前言题目题目描述输入输出样例输入样例输出题目解析解析扩展欧几里得二次函数参考代码后记前言在寒假集训中我学到了不少的新东西，也得到了许多感悟（当然还有很多天书一样的题）。在集训的末尾，我们也是迎来了情理之中的检测考试（虽然还是一道都不会）。不过呢，在之后也是努力钻研了其中的几道题，于是乎，我也决定和大家一起分享一下，共同来探讨这...
复制链接

扫一扫