Go: package main

import (
	"fmt"
	"math"
)

func main() {
	x := -1
	h := -1
	for {
		fmt.Scan(&x)
		fmt.Scan(&h)
		if x == 0 && h == 0 {
			return
		}

		f_x := float64(x) / 2.0

		sya := math.Sqrt(f_x*f_x + float64(h*h))
		fmt.Println(float64(x*x) + float64(x)*sya*2.0)

	}
}

C++: #define _USE_MATH_DEFINES
#define _CRT_SECURE_NO_DEPRECATE

#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <cmath>
#include <climits>
#include <cfloat>
#include <ctime>
#include <cassert>
#include <map>
#include <utility>
#include <set>
#include <iostream>
#include <memory>
#include <string>
#include <vector>
#include <algorithm>
#include <functional>
#include <sstream>
#include <complex>
#include <stack>
#include <queue>
#include <numeric>
#include <list>
#include <iomanip>
#include <fstream>


using namespace std;
const int MOD = 1000000009; // check!!!
const int INF = 100000000; //1e8

typedef long long ll;
typedef pair<int,int> Pii;
typedef pair<ll,ll> Pll;

#define FOR(i,n) for(int i = 0; i < (n); i++)
#define sz(c) ((int)(c).size())
#define ten(x) ((int)1e##x)
#define tenll(x) ((ll)1e##x)


const double EPS (1e-8);
const double PI=acos(-1.0);

#pragma region テ・ツ淞コテヲツ慊ャティツヲツ?ァツエツ?struct,typedef)

template<class T>
struct point{
    T x,y;
    point() : x(0), y(0) {}
    point(const T& x,const T& y) : x(x), y(y) {}
    point operator+(const point &a)const{ return point(x+a.x,y+a.y); }
    point& operator+=(const point &a){x += a.x; y += a.y; return *this;}
    point operator-(const point &a)const{ return point(x-a.x,y-a.y); }
    point& operator-=(const point &a){x -= a.x; y -= a.y; return *this;}
    point operator*(const double a)const{ return point(a*x,a*y); }
    point& operator*=(const double a){ x *= a; y *= a; return *this;}
    point operator/(const double a)const{ return point(x/a,y/a); }
    point& operator/=(const double a){ x /= a; y /= a; return *this;}
	
	//テヲツッツ氾ィツシツεァツ板ィ
	bool operator<(const point &a)const{return (x != a.x) ? (x < a.x) : (y < a.y);}
 
};

template<class T>
struct segment{
    point<T> a,b;
    segment() : a(point<T>()), b(point<T>()) {}
    segment(const point<T>& a,const point<T>& b) : a(a), b(b) {}
};

template<class T>
struct line{
    point<T> a,b;
    line() : a(point<T>()), b(point<T>()) {}
    line(const point<T>& a,const point<T>& b) : a(a), b(b) {}
};

template<class T>
struct circle{
	point<T> c;
	T r;
	circle(){}
	circle(const point<T> &c, T& r) : c(c),r(r) {}
};

typedef point<double> P;
typedef segment<double> S;
typedef line<double> L;
typedef circle<double> C;

#pragma endregion 

#pragma region テ・ツ淞コテヲツ慊ャティツヲツ?ァツエツ?テゥツ鳴「テヲツ閉ー)

template<class T> T SQ(T x){ return x*x; }
template<class T> T dist2(const point<T> &a,const point<T> &b){return SQ(a.x-b.x)+SQ(a.y-b.y);}
template<class T> T abs2(const point<T>& p){return SQ(p.x) + SQ(p.y);}
template<class T> T abs(const point<T>& p){return sqrt(abs2(p));}
template<class T> T dot(const point<T>& a,const point<T>& b)
{return a.x*b.x + a.y*b.y;}
template<class T> T cross(const point<T>& a,const point<T>& b)
{return a.x*b.y - a.y*b.x;}
template<class T> point<T> rot(const point<T>& a,const double theta){
        return point<T>(a.x*cos(theta)-a.y*sin(theta),
			a.x*sin(theta)+a.y*cos(theta));
}
template<class T>
double arg(const point<T> &a){	double t=atan2(a.y,a.x); return t<0?t+2*PI:t; }

enum{CCW=1,CW=-1,ON=0};
template<class T>
int ccw(const point<T> &a,point<T> b,point<T> c){
	b -= a; c -= a;
    double rdir=cross(b,c);
    double len = abs(b) * abs(c); //ティツェツ、テ・ツキツョテ・ツッツセテァツュツ?(テゥツ?淌・ツコツヲテ」ツ?古」ツつ?」ツ?ーテ」ツ??ヲツ卍づ」ツ?ッテ・ツ、ツ姪」ツ?凖」ツ?凝」ツ??、ツコツ古、ツケツ療」ツ?ォテァツスツョテ」ツ?催ヲツ渉崚」ツ?暗」ツつ?
    if(rdir> EPS * len) return CCW; //cテ」ツ?径bテ」ツつ暗」ツつ甘、ツクツ?テ・ツ渉催ヲツ卍づィツィツ暗・ツ堕ィテ」ツつ?
    if(rdir<-EPS * len) return CW; // cテ」ツ?径bテ」ツつ暗」ツつ甘、ツクツ?テヲツ卍づィツィツ暗・ツ堕ィテ」ツつ?
    return ON; // a,b,cテ」ツ?古、ツクツ?ァツ崢エテァツキツ堙、ツクツ?
}

#pragma endregion

#pragma region テ・ツ?・テ・ツ?コテ・ツ環?

std::ostream& operator<<(std::ostream& os, const P& point){return ( os << '(' << point.x << ',' << point.y << ')' );}
std::ostream& operator<<(std::ostream& os, const S& seg){return ( os << '{' << seg.a << ',' << seg.b << '}' );}
std::istream& operator>>(std::istream& is, P& point){return ( is >> point.x >> point.y );}

#pragma endregion

//テ、ツコツ、テァツつケテ」ツ?古・ツュツ佚・ツ慊ィテ」ツ?凖」ツつ凝」ツ?凝」ツ?ゥテ」ツ??」ツ??
bool intersected(const S &S1,const S &S2){
    if(max(S1.a.x,S1.b.x)+EPS<min(S2.a.x,S2.b.x)
    || max(S1.a.y,S1.b.y)+EPS<min(S2.a.y,S2.b.y)
    || max(S2.a.x,S2.b.x)+EPS<min(S1.a.x,S1.b.x)
    || max(S2.a.y,S2.b.y)+EPS<min(S1.a.y,S1.b.y)) return false;
    return ccw(S1.a,S1.b,S2.a)*ccw(S1.a,S1.b,S2.b)<=0
        && ccw(S2.a,S2.b,S1.a)*ccw(S2.a,S2.b,S1.b)<=0;
}

#pragma region ティツキツ敕ゥツ崢「

//テァツキツ堙・ツ按?」ツ?ィテァツつケテ」ツ?ョティツキツ敕ゥツ崢「(2テ、ツケツ?
double dist2(const S &S,const P &p){
    if(dot(S.b-S.a,p-S.a)<=0) return dist2(p,S.a);
    if(dot(S.a-S.b,p-S.b)<=0) return dist2(p,S.b);
    return (double)SQ(cross(S.b-S.a,p-S.a)) / dist2(S.a,S.b);
}

//テァツキツ堙・ツ按?」ツ?ィテァツつケテ」ツ?ョティツキツ敕ゥツ崢「
double dist(const S &S,const P &p){ return sqrt(dist2(S,p)); }

//テァツキツ堙・ツ按?」ツ?ィテァツキツ堙・ツ按?」ツ?ョティツキツ敕ゥツ崢「
double dist(const S &S1,const S &S2){
    if(intersected(S1,S2)) return 0;
    return sqrt(min(min(dist2(S1,S2.a),dist2(S1,S2.b)),
                    min(dist2(S2,S1.a),dist2(S2,S1.b))));
}

//テァツ崢エテァツキツ堙」ツ?ィテァツつケテ」ツ?ョティツキツ敕ゥツ崢「
double dist(const L &L,const P &p){
    return sqrt((double)SQ(cross(L.b-L.a,p-L.a)) / dist2(L.a,L.b));
}

//a -- bテ」ツ?ョテゥツ鳴禿」ツ?ョテ・ツ??」ツ?古・ツ?・テ」ツ?」テ」ツ?淌ヲツ卍づ」ツ?ョテ」ツ?‖ - bテゥツ鳴禿」ツ?ョティツキツ敕ゥツ崢「
//テ・ツ??」ツ?ョテゥツ鳴禿」ツ?ォ
double geodist(P a,P b,const C &c){
	double r=c.r;
	S s(a,b);
	if(dist(s,c.c)>r) 
		return abs(a-b);

	a-=c.c;
	b-=c.c;

	double L1=abs(a);
	double t1=arg(a)+acos(r/L1); if(t1<0) t1+=2*PI;
	double t2=arg(a)-acos(r/L1); if(t2<0) t2+=2*PI;

	double L2=abs(b);
	double t3=arg(b)+acos(r/L2); if(t3<0) t3+=2*PI;
	double t4=arg(b)-acos(r/L2); if(t4<0) t4+=2*PI;

	double theta=2*PI;
	theta=min(theta,min(abs(t1-t4),2*PI-abs(t1-t4)));
	theta=min(theta,min(abs(t2-t3),2*PI-abs(t2-t3)));
	return sqrt(L1*L1-r*r)+sqrt(L2*L2-r*r)+r*theta;
}

#pragma endregion

//テ・ツ、ツ堙ィツァツ津・ツスツ「テ」ツ?ョテゥツ敖「テァツゥツ?
double area(const vector<P>& v){
	double ans = 0;
	for (int i = 2; i < (int)v.size(); i++){
		ans += cross(v[i-1] - v[0],v[i] - v[0]);
	}
	return abs(ans) / 2;
}

//テ・ツ?クテ・ツ個?
vector<P> convex_hull(vector<P>& ps) {
  int n = ps.size(), k = 0;
  sort(ps.begin(), ps.end());
  vector<P> ch(2*n);
  for (int i = 0; i < n; ch[k++] = ps[i++]) // lower-hull
    while (k >= 2 && ccw(ch[k-2], ch[k-1], ps[i]) <= 0) --k;

  for (int i = n-2, t = k+1; i >= 0; ch[k++] = ps[i--]) // upper-hull
    while (k >= t && ccw(ch[k-2], ch[k-1], ps[i]) <= 0) --k;

  ch.resize(k-1);

  return ch;
}

#pragma region テ、ツコツ、テァツつケ

typedef pair<P,P> PD;
//テ・ツ??・ツ青古・ツ」ツォテ」ツ?ョテ、ツコツ、テァツつケ
PD GetContact(const C& p, const C &q){
	P a = q.c-p.c;
	//Point b = rot(a,PI / 2);
	P b(a.y,-a.x);
	double c = abs(a);
	b = b / c;

	double s = (p.r + q.r + c) / 2;
	double men = sqrt(s * (s - p.r) * (s - q.r) * (s - c));
	double l = 2 * men / c;
	double lp = sqrt(SQ(p.r) - SQ(l));
	
	a = a / c * lp;
	if(SQ(q.r) - SQ(p.r) - SQ(c) > 0){ //ティツ?ε」ツ?暗」ツ?ヲテ」ツ??」ツつ凝、ツクツ嘉ィツァツ津・ツスツ「テ」ツ?ョティツァツ津」ツ??0テ・ツコツヲテ、ツサツ・テ、ツクツ甘」ツ?ョテ・ツ?エテ・ツ青?
		a *= -1;
	}
	b = b * l;
	return PD(p.c + a + b,p.c + a - b);
}

// l : ax+by+c=0 テ」ツ?ォテ・ツ、ツ嘉ヲツ渉?
inline void calc_abc(const L &l,double &a,double &b,double &c){
    a=l.a.y-l.b.y;
    b=l.b.x-l.a.x;
    c=l.a.x*l.b.y-l.b.x*l.a.y;
}
 
//Pテ」ツ?凝」ツつ鵜テ」ツ?クテ」ツ?ョテ・ツ楪づァツキツ堙」ツ?ョティツカツウ
P perp_foot(const P &p,const L &l){
    double a,b,c;
    calc_abc(l,a,b,c);
    return p-P(a,b) * ((a*p.x+b*p.y+c)/(a*a+b*b));
}
 
//テ・ツ??」ツ?ィテァツ崢エテァツキツ堙」ツ?ョテ、ツコツ、テァツつケ
int intersect(const C &c,const L &l,P *p=NULL,P *q=NULL){
    P m=perp_foot(c.c,l);
    double d2=abs2(c.c-m);
    if(c.r*c.r+EPS<d2)   return 0;
    else if(abs(c.r*c.r-d2)<EPS){
        if(p&&q)    *p=*q=m;
        return 1;
    }
    else{
        if(p&&q){
            P v;
            double n0=abs2(l.a-m),n1=abs2(l.b-m);
            if(n0<n1)    v=l.b-m,n0=n1;
            else        v=l.a-m;
            v*=sqrt((c.r*c.r-d2)/n0);
            *p=m+v,*q=m-v;
        }
        return 2;
    }
}

//テ・ツ??」ツ?ィテァツキツ堙・ツ按?」ツ?ョテ、ツコツ、テァツつケ
int intersect(const C &c,const S &s,P *p=NULL,P *q=NULL){
    P pp,qq;
    int n=intersect(c,L(s.a,s.b),&pp,&qq);
    if(n==0)    return 0;
    else if(n==1){
        if(dot(pp-s.a,pp-s.b)<EPS){
            if(p&&q)    *p=*q=pp;
            return 1;
        }
        return 0;
    }
    else{
        bool bp=(dot(pp-s.a,pp-s.b)<EPS);
        bool bq=(dot(qq-s.a,qq-s.b)<EPS);
        if(bp&&bq){
            if(p&&q)    *p=pp,*q=qq;
            return 2;
        }
        if(bp||bq){
            if(p&&q)    *p=(bp?pp:qq);
            return 1;
        }
        return 0;
    }
}

//テァツ崢エテァツキツ堙」ツ?ョテ、ツコツ、テァツつケ
P intersect(const L &L1,const L &L2){
	double a1=cross(L1.b-L1.a,L2.b-L2.a);
	double a2=cross(L1.b-L1.a,L1.b-L2.a);
	if(abs(a1)<EPS) return L1.a;
	return L2.a+(L2.b-L2.a) * (a2/a1);
}

#pragma endregion


int main(){
	int x,h;
	while(cin>>x>>h,x&&h){
		double l = sqrt(2) / 2 * x;
		double m = sqrt(SQ(l)+SQ(h));
		double n = sqrt(SQ(m)-SQ(x/2.0));
		double ans = 2*x*n+x*x;
		printf("%.15f\n",ans);
	}

	return 0;
}