verify/miscellaneous/integral.test.cpp

• View this file on GitHub
• Last update: 2025-07-19 21:22:35+09:00
• Problem: https://judge.yosupo.jp/problem/aplusb

Depends on

• src/contest/template.hpp
• src/miscellaneous/integrate.hpp

Code

#define PROBLEM "https://judge.yosupo.jp/problem/aplusb"
#include "src/contest/template.hpp"
#include "src/miscellaneous/integrate.hpp"

int main(){
    for(int r=1;r<=20;r++){
        double out=quad(-r,r,[&](db x){
            return quad(-r,r,[&](db y){
                return quad(-r,r,[&](db z){
                    return x*x+y*y+z*z<r*r;
                },100);
            },100);
        },100);
        double ans=4*PI*r*r*r/3;
        auto err=abs(out-ans)/max(ans,1.0);
        assert(err<1e-4);
    }
    int a,b;
    cin >> a >> b;
    cout << a+b << "\n";
}

#line 1 "verify/miscellaneous/integral.test.cpp"
#define PROBLEM "https://judge.yosupo.jp/problem/aplusb"
#line 2 "src/contest/template.hpp"
#include<bits/stdc++.h>
#include<ext/pb_ds/assoc_container.hpp>
#include<ext/pb_ds/tree_policy.hpp>
 
using namespace std;
using namespace __gnu_pbds;
 
using ll = long long;
using db = long double;
using vi = vector<int>;
using vl = vector<ll>;
using vd = vector<db>;
using pii = pair<int,int>;
using pll = pair<ll,ll>;
using pdd = pair<db,db>;
const int INF=0x3fffffff;
const ll LINF=0x1fffffffffffffff;
const db DINF=numeric_limits<db>::infinity();
const db EPS=1e-9;
const db PI=acos(db(-1));
 
template<class T>
using ordered_set = tree<T,null_type,less<T>,rb_tree_tag,tree_order_statistics_node_update>;
 
mt19937 rng(chrono::steady_clock::now().time_since_epoch().count());
mt19937_64 rng64(chrono::steady_clock::now().time_since_epoch().count());
#line 2 "src/miscellaneous/integrate.hpp"

/**
 * Author: Teetat T.
 * Date: 2025-07-19
 * Description: find definite integral with Simpson's method, error proportional to $dx^4$.
 */

template<class F>
db quad(db a,db b,const F &f,int n){
    db res=0;
    db dx=(b-a)/n;
    db fl=0,fr=f(a);
    for(int i=0;i<n;i++){
        db l=a+dx*i,r=l+dx;
        fl=fr;
        fr=f(r);
        db fm=f((l+r)/2);
        res+=fl+4*fm+fr;
    }
    return res*dx/6;
}
#line 4 "verify/miscellaneous/integral.test.cpp"

int main(){
    for(int r=1;r<=20;r++){
        double out=quad(-r,r,[&](db x){
            return quad(-r,r,[&](db y){
                return quad(-r,r,[&](db z){
                    return x*x+y*y+z*z<r*r;
                },100);
            },100);
        },100);
        double ans=4*PI*r*r*r/3;
        auto err=abs(out-ans)/max(ans,1.0);
        assert(err<1e-4);
    }
    int a,b;
    cin >> a >> b;
    cout << a+b << "\n";
}

Back to top page