#include <stdio.h>  
#include <math.h>  
#include <stdlib.h>  
#include <time.h>  
  
// 定义π常量  
#define PI 3.14159265358979323846  
  
// 椭圆参数  
static double a = 2.0, b = 1.0;  
  
// 目标函数 f(t) = (a*cos(t)-x)*(a*sin(t)) + (b*sin(t)-y)*b*cos(t)// 这是距离平方对t的导数的一半，求其零点即可找到极值点  
double f(double t, double x, double y) {  
    return (a * cos(t) - x) * (a * sin(t)) + (b * sin(t) - y) * b * cos(t);  
}  
  
// 目标函数的导数 f'(t) = (a² + b²) * (cos²(t) - sin²(t)) - a*x*cos(t) + b*y*sin(t)double df(double t, double x, double y) {  
    double cos_t = cos(t);  
    double sin_t = sin(t);  
  
    return (a * a + b * b) * (cos_t * cos_t - sin_t * sin_t)  
           - a * x * cos_t  
           + b * y * sin_t;  
}  
  
// 点到椭圆的距离  
double distance(double t, double x, double y) {  
    double ellipse_x = a * cos(t);  
    double ellipse_y = b * sin(t);  
    return sqrt((ellipse_x - x) * (ellipse_x - x) +  
                (ellipse_y - y) * (ellipse_y - y));  
}  
  
// 牛顿下山法求解，增加初始角度参数  
// 成功返回1，失败返回0  
int newton_downhill(double x, double y, double initial_t, double* t_result,  
                    double e1, double e2, int max_iter) {  
    double t0 = initial_t; // 使用传入的初始角度  
    double lambda = 1.0; // 下山因子  
    double t, ft, dft;  
  
    for (int i = 0; i < max_iter; i++) {  
        ft = f(t0, x, y);  
        dft = df(t0, x, y);  
  
        // 检查导数是否过小，避免除以零  
        if (fabs(dft) < e2) {  
            // 导数接近零，可能已找到极值点或遇到平台  
            if (fabs(ft) < e1) {  
                *t_result = t0;  
                return 1; // 收敛  
            }  
            return 0; // 失败  
        }  
  
        // 牛顿下山法迭代  
        t = t0 - lambda * ft / dft;  
  
        // 处理角度周期性，确保t在[0, 2*PI)范围内  
        while (t >= 2 * PI) t -= 2 * PI;  
        while (t < 0) t += 2 * PI;  
  
        // 检查收敛条件  
        if (fabs(f(t, x, y)) < e1 && fabs(t - t0) < e2) {  
            *t_result = t;  
            return 1; // 成功  
        }  
  
        // 下山条件：新点的函数值绝对值更小  
        if (fabs(f(t, x, y)) < fabs(ft)) {  
            t0 = t;  
            lambda = 1.0; // 恢复下山因子  
        } else {  
            lambda *= 0.5; // 减小下山因子  
            if (lambda < 1e-10) { // 防止lambda过小导致死循环  
                return 0;  
            }  
        }  
    }  
  
    return 0; // 达到最大迭代次数  
}  
  
// 多初始点尝试（提高成功率）  
double find_nearest_point(double x, double y, double* best_t) {  
    double initial_angles[] = {0.0, PI/4, PI/2, 3*PI/4, PI, 5*PI/4, 3*PI/2, 7*PI/4};  
    int num_angles = 8;  
    double min_distance = 1e10;  
    double t_result;  
    int success_count = 0;  
  
    for (int i = 0; i < num_angles; i++) {  
        // 将初始角度传递给牛顿下山法  
        if (newton_downhill(x, y, initial_angles[i], &t_result, 1e-9, 1e-9, 100)) {  
            success_count++;  
            double dist = distance(t_result, x, y);  
            if (dist < min_distance) {  
                min_distance = dist;  
                *best_t = t_result;  
            }  
        }  
    }  
  
    printf("成功收敛次数: %d/%d\n", success_count, num_angles);  
    return min_distance;  
}  
  
// 验证算法  
void test_algorithm() {  
    srand((unsigned int)time(NULL));  
    int num_tests = 20;  
    int success_count = 0;  
  
    printf("=== 牛顿下山法求解点到椭圆最近距离验证 (修正版) ===\n");  
    printf("椭圆方程: x²/%.1f² + y²/%.1f² = 1\n\n", a, b);  
  
    for (int i = 0; i < num_tests; i++) {  
        // 生成随机测试点  
        double test_x = (rand() % 200 - 100) / 10.0; // [-10, 10]  
        double test_y = (rand() % 200 - 100) / 10.0; // [-10, 10]  
  
        double best_t;  
        double min_dist = find_nearest_point(test_x, test_y, &best_t);  
  
        if (min_dist < 1e9) { // 使用一个比1e10小的值来判断成功  
            success_count++;  
            double nearest_x = a * cos(best_t);  
            double nearest_y = b * sin(best_t);  
  
            printf("测试 %2d: 点(%.2f, %.2f) -> 最近点(%.4f, %.4f), 距离=%.6f\n",  
                   i+1, test_x, test_y, nearest_x, nearest_y, min_dist);  
        } else {  
            printf("测试 %2d: 点(%.2f, %.2f) -> 计算失败\n", i+1, test_x, test_y);  
        }  
    }  
  
    printf("\n=== 统计结果 ===\n");  
    printf("总测试次数: %d\n", num_tests);  
    printf("成功次数: %d\n", success_count);  
    printf("成功率: %.2f%%\n", (success_count * 100.0) / num_tests);  
}  
  
int main() {  
    test_algorithm();  
    return 0;  
}