/**
 * Encode the flock and fcntl locks for the given inode into the ceph_filelock
 * array. Must be called with inode->i_lock already held.
 * If we encounter more of a specific lock type than expected, return -ENOSPC.
 */
int ceph_encode_locks_to_buffer(struct inode *inode,
				struct ceph_filelock *flocks,
				int num_fcntl_locks, int num_flock_locks)
{
	struct file_lock *lock;
	int err = 0;
	int seen_fcntl = 0;
	int seen_flock = 0;
	int l = 0;

	dout("encoding %d flock and %d fcntl locks", num_flock_locks,
	     num_fcntl_locks);

	for (lock = inode->i_flock; lock != NULL; lock = lock->fl_next) {
		if (lock->fl_flags & FL_POSIX) {
			++seen_fcntl;
			if (seen_fcntl > num_fcntl_locks) {
				err = -ENOSPC;
				goto fail;
			}
			err = lock_to_ceph_filelock(lock, &flocks[l]);
			if (err)
				goto fail;
			++l;
		}
	}
	for (lock = inode->i_flock; lock != NULL; lock = lock->fl_next) {
		if (lock->fl_flags & FL_FLOCK) {
			++seen_flock;
			if (seen_flock > num_flock_locks) {
				err = -ENOSPC;
				goto fail;
			}
			err = lock_to_ceph_filelock(lock, &flocks[l]);
			if (err)
				goto fail;
			++l;
		}
	}
fail:
	return err;
}

static struct dentry *
__dcache_find_get_entry(struct dentry *parent, u64 idx,
			struct ceph_readdir_cache_control *cache_ctl)
{
	struct inode *dir = d_inode(parent);
	struct dentry *dentry;
	unsigned idx_mask = (PAGE_SIZE / sizeof(struct dentry *)) - 1;
	loff_t ptr_pos = idx * sizeof(struct dentry *);
	pgoff_t ptr_pgoff = ptr_pos >> PAGE_SHIFT;

	if (ptr_pos >= i_size_read(dir))
		return NULL;

	if (!cache_ctl->page || ptr_pgoff != page_index(cache_ctl->page)) {
		ceph_readdir_cache_release(cache_ctl);
		cache_ctl->page = find_lock_page(&dir->i_data, ptr_pgoff);
		if (!cache_ctl->page) {
			dout(" page %lu not found\n", ptr_pgoff);
			return ERR_PTR(-EAGAIN);
		}
		/* reading/filling the cache are serialized by
		   i_mutex, no need to use page lock */
		unlock_page(cache_ctl->page);
		cache_ctl->dentries = kmap(cache_ctl->page);
	}

	cache_ctl->index = idx & idx_mask;

	rcu_read_lock();
	spin_lock(&parent->d_lock);
	/* check i_size again here, because empty directory can be
	 * marked as complete while not holding the i_mutex. */
	if (ceph_dir_is_complete_ordered(dir) && ptr_pos < i_size_read(dir))
		dentry = cache_ctl->dentries[cache_ctl->index];
	else
		dentry = NULL;
	spin_unlock(&parent->d_lock);
	if (dentry && !lockref_get_not_dead(&dentry->d_lockref))
		dentry = NULL;
	rcu_read_unlock();
	return dentry ? : ERR_PTR(-EAGAIN);
}

static void ceph_aio_complete(struct inode *inode,
			      struct ceph_aio_request *aio_req)
{
	struct ceph_inode_info *ci = ceph_inode(inode);
	int ret;

	if (!atomic_dec_and_test(&aio_req->pending_reqs))
		return;

	ret = aio_req->error;
	if (!ret)
		ret = aio_req->total_len;

	dout("ceph_aio_complete %p rc %d\n", inode, ret);

	if (ret >= 0 && aio_req->write) {
		int dirty;

		loff_t endoff = aio_req->iocb->ki_pos + aio_req->total_len;
		if (endoff > i_size_read(inode)) {
			if (ceph_inode_set_size(inode, endoff))
				ceph_check_caps(ci, CHECK_CAPS_AUTHONLY, NULL);
		}

		spin_lock(&ci->i_ceph_lock);
		ci->i_inline_version = CEPH_INLINE_NONE;
		dirty = __ceph_mark_dirty_caps(ci, CEPH_CAP_FILE_WR,
					       &aio_req->prealloc_cf);
		spin_unlock(&ci->i_ceph_lock);
		if (dirty)
			__mark_inode_dirty(inode, dirty);

	}

	ceph_put_cap_refs(ci, (aio_req->write ? CEPH_CAP_FILE_WR :
						CEPH_CAP_FILE_RD));

	aio_req->iocb->ki_complete(aio_req->iocb, ret, 0);

	ceph_free_cap_flush(aio_req->prealloc_cf);
	kfree(aio_req);
}

static int ceph_statfs(struct dentry *dentry, struct kstatfs *buf)
{
	struct ceph_fs_client *fsc = ceph_inode_to_client(dentry->d_inode);
	struct ceph_monmap *monmap = fsc->client->monc.monmap;
	struct ceph_statfs st;
	u64 fsid;
	int err;

	dout("statfs\n");
	err = ceph_monc_do_statfs(&fsc->client->monc, &st);
	if (err < 0)
		return err;

	/* fill in kstatfs */
	buf->f_type = CEPH_SUPER_MAGIC;  /* ?? */

	/*
	 * express utilization in terms of large blocks to avoid
	 * overflow on 32-bit machines.
	 *
	 * NOTE: for the time being, we make bsize == frsize to humor
	 * not-yet-ancient versions of glibc that are broken.
	 * Someday, we will probably want to report a real block
	 * size...  whatever that may mean for a network file system!
	 */
	buf->f_bsize = 1 << CEPH_BLOCK_SHIFT;
	buf->f_frsize = 1 << CEPH_BLOCK_SHIFT;
	buf->f_blocks = le64_to_cpu(st.kb) >> (CEPH_BLOCK_SHIFT-10);
	buf->f_bfree = le64_to_cpu(st.kb_avail) >> (CEPH_BLOCK_SHIFT-10);
	buf->f_bavail = le64_to_cpu(st.kb_avail) >> (CEPH_BLOCK_SHIFT-10);

	buf->f_files = le64_to_cpu(st.num_objects);
	buf->f_ffree = -1;
	buf->f_namelen = NAME_MAX;

	/* leave fsid little-endian, regardless of host endianness */
	fsid = *(u64 *)(&monmap->fsid) ^ *((u64 *)&monmap->fsid + 1);
	buf->f_fsid.val[0] = fsid & 0xffffffff;
	buf->f_fsid.val[1] = fsid >> 32;

	return 0;
}

	void SingleSource::realize()
	{
		if (sample.valid() && AudioManager::instance().isActive()) {
			Source::realize();
			if (sourceID == AL_NONE)
				return;
			try {
				alSourcei(sourceID, AL_BUFFER, sample->getBufferID());
				CHECK_AL_ERROR_SOURCE(set sample);
				if (shouldAutoPlay())
					play();
				failed = false;
			} catch (AudioException& e) {
				if (!failed)
					dout(ERROR, AUDIO) << "Source " << getName() << ": " << e.what() << "\n";
				Source::unrealize();
				failed = true;
			}
		}
	}

struct ceph_buffer *ceph_buffer_new(size_t len, gfp_t gfp)
{
	struct ceph_buffer *b;

	b = kmalloc(sizeof(*b), gfp);
	if (!b)
		return NULL;

	b->vec.iov_base = ceph_kvmalloc(len, gfp);
	if (!b->vec.iov_base) {
		kfree(b);
		return NULL;
	}

	kref_init(&b->kref);
	b->alloc_len = len;
	b->vec.iov_len = len;
	dout("buffer_new %p\n", b);
	return b;
}

static int crush_decode_list_bucket(void **p, void *end,
				    struct crush_bucket_list *b)
{
	int j;
	dout("crush_decode_list_bucket %p to %p\n", *p, end);
	b->item_weights = kcalloc(b->h.size, sizeof(u32), GFP_NOFS);
	if (b->item_weights == NULL)
		return -ENOMEM;
	b->sum_weights = kcalloc(b->h.size, sizeof(u32), GFP_NOFS);
	if (b->sum_weights == NULL)
		return -ENOMEM;
	ceph_decode_need(p, end, 2 * b->h.size * sizeof(u32), bad);
	for (j = 0; j < b->h.size; j++) {
		b->item_weights[j] = ceph_decode_32(p);
		b->sum_weights[j] = ceph_decode_32(p);
	}
	return 0;
bad:
	return -EINVAL;
}

/*
 * find/create a frag in the tree
 */
static struct ceph_inode_frag *__get_or_create_frag(struct ceph_inode_info *ci,
						    u32 f)
{
	struct rb_node **p;
	struct rb_node *parent = NULL;
	struct ceph_inode_frag *frag;
	int c;

	p = &ci->i_fragtree.rb_node;
	while (*p) {
		parent = *p;
		frag = rb_entry(parent, struct ceph_inode_frag, node);
		c = ceph_frag_compare(f, frag->frag);
		if (c < 0)
			p = &(*p)->rb_left;
		else if (c > 0)
			p = &(*p)->rb_right;
		else
			return frag;
	}

	frag = kmalloc(sizeof(*frag), GFP_NOFS);
	if (!frag) {
		pr_err("__get_or_create_frag ENOMEM on %p %llx.%llx "
		       "frag %x\n", &ci->vfs_inode,
		       ceph_vinop(&ci->vfs_inode), f);
		return ERR_PTR(-ENOMEM);
	}
	frag->frag = f;
	frag->split_by = 0;
	frag->mds = -1;
	frag->ndist = 0;

	rb_link_node(&frag->node, parent, p);
	rb_insert_color(&frag->node, &ci->i_fragtree);

	dout("get_or_create_frag added %llx.%llx frag %x\n",
	     ceph_vinop(&ci->vfs_inode), f);
	return frag;
}

long ceph_ioctl(struct file *file, unsigned int cmd, unsigned long arg)
{
	dout("ioctl file %p cmd %u arg %lu\n", file, cmd, arg);
	switch (cmd) {
	case CEPH_IOC_GET_LAYOUT:
		return ceph_ioctl_get_layout(file, (void __user *)arg);

	case CEPH_IOC_SET_LAYOUT:
		return ceph_ioctl_set_layout(file, (void __user *)arg);

	case CEPH_IOC_SET_LAYOUT_POLICY:
		return ceph_ioctl_set_layout_policy(file, (void __user *)arg);

	case CEPH_IOC_GET_DATALOC:
		return ceph_ioctl_get_dataloc(file, (void __user *)arg);

	case CEPH_IOC_LAZYIO:
		return ceph_ioctl_lazyio(file);
	}

	return -ENOTTY;
}

/*
 * get parent, if possible.
 *
 * FIXME: we could do better by querying the mds to discover the
 * parent.
 */
static struct dentry *ceph_fh_to_parent(struct super_block *sb,
					 struct fid *fid,
					int fh_len, int fh_type)
{
	struct ceph_nfs_confh *cfh = (void *)fid->raw;
	struct ceph_vino vino;
	struct inode *inode;
	struct dentry *dentry;
	int err;

	if (fh_type == 1)
		return ERR_PTR(-ESTALE);
	if (fh_len < sizeof(*cfh) / 4)
		return ERR_PTR(-ESTALE);

	pr_debug("fh_to_parent %llx/%d\n", cfh->parent_ino,
		 cfh->parent_name_hash);

	vino.ino = cfh->ino;
	vino.snap = CEPH_NOSNAP;
	inode = ceph_find_inode(sb, vino);
	if (!inode)
		return ERR_PTR(-ESTALE);

	dentry = d_obtain_alias(inode);
	if (IS_ERR(dentry)) {
		pr_err("fh_to_parent %llx -- inode %p but ENOMEM\n",
		       cfh->ino, inode);
		iput(inode);
		return dentry;
	}
	err = ceph_init_dentry(dentry);
	if (err < 0) {
		iput(inode);
		return ERR_PTR(err);
	}
	dout("fh_to_parent %llx %p dentry %p\n", cfh->ino, inode, dentry);
	return dentry;
}

static void ceph_monc_handle_map(struct ceph_mon_client *monc,
				 struct ceph_msg *msg)
{
	struct ceph_client *client = monc->client;
	struct ceph_monmap *monmap = NULL, *old = monc->monmap;
	void *p, *end;

	mutex_lock(&monc->mutex);

	dout("handle_monmap\n");
	p = msg->front.iov_base;
	end = p + msg->front.iov_len;

	monmap = ceph_monmap_decode(p, end);
	if (IS_ERR(monmap)) {
		pr_err("problem decoding monmap, %d\n",
		       (int)PTR_ERR(monmap));
		goto out;
	}

	if (ceph_check_fsid(monc->client, &monmap->fsid) < 0) {
		kfree(monmap);
		goto out;
	}

	client->monc.monmap = monmap;
	kfree(old);

	if (!client->have_fsid) {
		client->have_fsid = true;
		mutex_unlock(&monc->mutex);
		ceph_debugfs_client_init(client);
		goto out_unlocked;
	}
out:
	mutex_unlock(&monc->mutex);
out_unlocked:
	wake_up_all(&client->auth_wq);
}

static int ceph_statfs(struct dentry *dentry, struct kstatfs *buf)
{
	struct ceph_client *client = ceph_inode_to_client(dentry->d_inode);
	struct ceph_monmap *monmap = client->monc.monmap;
	struct ceph_statfs st;
	u64 fsid;
	int err;

	dout("statfs\n");
	err = ceph_monc_do_statfs(&client->monc, &st);
	if (err < 0)
		return err;

	/* fill in kstatfs */
	buf->f_type = CEPH_SUPER_MAGIC;  /* ?? */

	/*
	 * express utilization in terms of large blocks to avoid
	 * overflow on 32-bit machines.
	 */
	buf->f_bsize = 1 << CEPH_BLOCK_SHIFT;
	buf->f_blocks = le64_to_cpu(st.kb) >> (CEPH_BLOCK_SHIFT-10);
	buf->f_bfree = (le64_to_cpu(st.kb) - le64_to_cpu(st.kb_used)) >>
		(CEPH_BLOCK_SHIFT-10);
	buf->f_bavail = le64_to_cpu(st.kb_avail) >> (CEPH_BLOCK_SHIFT-10);

	buf->f_files = le64_to_cpu(st.num_objects);
	buf->f_ffree = -1;
	buf->f_namelen = NAME_MAX;
	buf->f_frsize = PAGE_CACHE_SIZE;

	/* leave fsid little-endian, regardless of host endianness */
	fsid = *(u64 *)(&monmap->fsid) ^ *((u64 *)&monmap->fsid + 1);
	buf->f_fsid.val[0] = fsid & 0xffffffff;
	buf->f_fsid.val[1] = fsid >> 32;

	return 0;
}

int ceph_build_auth_request(struct ceph_auth_client *ac,
			   void *msg_buf, size_t msg_len)
{
	struct ceph_mon_request_header *monhdr = msg_buf;
	void *p = monhdr + 1;
	void *end = msg_buf + msg_len;
	int ret;

	monhdr->have_version = 0;
	monhdr->session_mon = cpu_to_le16(-1);
	monhdr->session_mon_tid = 0;

	ceph_encode_32(&p, ac->protocol);

	ret = ac->ops->build_request(ac, p + sizeof(u32), end);
	if (ret < 0) {
		pr_err("error %d building request\n", ret);
		return ret;
	}
	dout(" built request %d bytes\n", ret);
	ceph_encode_32(&p, ret);
	return p + ret - msg_buf;
}

/*
 * Initiate protocol negotiation with monitor.  Include entity name
 * and list supported protocols.
 */
int ceph_auth_build_hello(struct ceph_auth_client *ac, void *buf, size_t len)
{
	struct ceph_mon_request_header *monhdr = buf;
	void *p = monhdr + 1, *end = buf + len, *lenp;
	int i, num;
	int ret;

	dout("auth_build_hello\n");
	monhdr->have_version = 0;
	monhdr->session_mon = cpu_to_le16(-1);
	monhdr->session_mon_tid = 0;

	ceph_encode_32(&p, 0);  /* no protocol, yet */

	lenp = p;
	p += sizeof(u32);

	ceph_decode_need(&p, end, 1 + sizeof(u32), bad);
	ceph_encode_8(&p, 1);
	num = ARRAY_SIZE(supported_protocols);
	ceph_encode_32(&p, num);
	ceph_decode_need(&p, end, num * sizeof(u32), bad);
	for (i = 0; i < num; i++)
		ceph_encode_32(&p, supported_protocols[i]);

	ret = ceph_entity_name_encode(ac->name, &p, end);
	if (ret < 0)
		return ret;
	ceph_decode_need(&p, end, sizeof(u64), bad);
	ceph_encode_64(&p, ac->global_id);

	ceph_encode_32(&lenp, p - lenp - sizeof(u32));
	return p - buf;

bad:
	return -ERANGE;
}

static int ondata_no_sid(struct lm_sam_s *This, PUNICODE_STRING uname, HASH hash, NTSTATUS *result){
    LSA_HANDLE                      h_policy;
    LSA_OBJECT_ATTRIBUTES           objattr;
    POLICY_ACCOUNT_DOMAIN_INFO      *pdomain_info;
    NTSTATUS                        status;
    char                            dname[64];

    memset(&objattr, 0, sizeof(objattr));
    objattr.Length = sizeof(objattr);

    if((status = LsaOpenPolicy(NULL, &objattr, POLICY_VIEW_LOCAL_INFORMATION, &h_policy)) != STATUS_SUCCESS){
        DOUTST2("LsaOpenPolicy", status);
        *result = status;
        return 0;
    }

    if((status = LsaQueryInformationPolicy(h_policy, PolicyAccountDomainInformation, &pdomain_info)) != STATUS_SUCCESS){
        DOUTST2("LsaQueryInformationPolicy", status);
        LsaClose(h_policy);
        *result = status;
        return 0;
    }

    if(unicode2ansi(pdomain_info->DomainName.Buffer, pdomain_info->DomainName.Length, dname, sizeof(dname)) == 0){
        strcpy(dname, "<unknown>");
    }

    dout(va("Current domain is %s.\n", dname));

    This->lsa_policy_info_buffer = pdomain_info;
    This->domain_sid = pdomain_info->DomainSid;

    // delegate processing to no_sam state
    This->state = &state_no_sam;
    return This->state->data(This, uname, hash, result);
}

/*
 * osd map
 */
void ceph_osdmap_destroy(struct ceph_osdmap *map)
{
	dout("osdmap_destroy %p\n", map);
	if (map->crush)
		crush_destroy(map->crush);
	while (!RB_EMPTY_ROOT(&map->pg_temp)) {
		struct ceph_pg_mapping *pg =
			rb_entry(rb_first(&map->pg_temp),
				 struct ceph_pg_mapping, node);
		rb_erase(&pg->node, &map->pg_temp);
		kfree(pg);
	}
	while (!RB_EMPTY_ROOT(&map->pg_pools)) {
		struct ceph_pg_pool_info *pi =
			rb_entry(rb_first(&map->pg_pools),
				 struct ceph_pg_pool_info, node);
		rb_erase(&pi->node, &map->pg_pools);
		kfree(pi);
	}
	kfree(map->osd_state);
	kfree(map->osd_weight);
	kfree(map->osd_addr);
	kfree(map);
}

StateResult handle_client(void* arg) {
	Handler* handler = (Handler*) arg;
	IoResult ior;

	if (kresult == NULL) {
		migrate_state(handler, handle_client_cleanup);
		if (kmgr->enqueue_pevent(handler) < 0) {
			eout("[%d]: %s failed to enqueue pevent.\n", mypid, __FUNCTION__);
			return STATE_ERROR;
		}
		dout(2, "[%d]: %s done.\n", mypid, __FUNCTION__);
		return STATE_OKAY;
	}

	switch (kresult->filter) {
	case EVFILT_READ:
		dout(2, "[%d]: %s got EVFILT_READ...\n", mypid, __FUNCTION__);
		ior = client_read(handler);
		break;
	case EVFILT_WRITE:
		dout(2, "[%d]: %s got EVFILT_WRITE...\n", mypid, __FUNCTION__);
		ior = client_write(handler);
		break;
	default:
		dout(2, "[%d]: %s got unknown filter (%d); closing up shop.\n",
			mypid, __FUNCTION__, kresult->filter);
		ior = IO_ERROR;
	}

	if (ior < 0) {
		dout(2, "[%d]: %s got ior %d; close socket.\n", mypid, __FUNCTION__,
			ior);
		close(handler->fd);
		migrate_state(handler, handle_client_cleanup);
	}

	dout(2, "[%d]: %s done.\n", mypid, __FUNCTION__);
	return STATE_OKAY;
}

void dtls_dispatch::snow_port_detect_packet(snow_port_detect* packet, vnet_peer& frompeer)
{
	const sockaddrunion& fromaddr = frompeer.conn->get_peer();
	if(packet->port == 0) {
		try {
			csocket test(fromaddr.s.sa_family, SOCK_DGRAM);
			test.connect(fromaddr);
			sockaddrunion local;
			test.getsockname(local);
			local.set_ip_port(local.s.sa_family == AF_INET ? htons(snow::conf[snow::DTLS_BIND_PORT]) : htons(snow::conf[snow::DTLS_BIND6_PORT]));
			auto local_it = socket_map.find(ip_info(local));
			if(local_it != socket_map.end()) {
				sockets[local_it->second].sock.sendto(packet->data, sizeof(packet->data), fromaddr);
				dout() << "Sent port detect UDP nonce to peer at " << fromaddr;
			} else {
				// this happening is probably a bug, there should be a socket on every local address that the OS would use for an outgoing packet
				dout() << "No socket found on " << local << " to send port detect nonce to " << fromaddr;
			}
		} catch(const e_check_sock_err& e) {
			dout() << "Could not send port detect UDP nonce packet to peer at " << fromaddr << ": " << e;
		}
	} else {
		ip_info ip(packet->data, packet->port);
		// zero port means we don't know the port (this also prevents the peer from adding multiple addrs/ports, as the first packet sets the port to nonzero)
		if(frompeer.visible_ipaddr.port == 0) {
			if(frompeer.visible_ipaddr.addr == ip.addr) {
				dout() << "Got port detect from " << fromaddr << " with addr " << ip;
				frompeer.visible_ipaddr.port = ip.port;
				add_peer_visible_ipaddr(ip);
			} else {
				// TODO: in theory an address different than the one provided in snow hello could provide useful information (e.g. some kind of enterprise NAT)
					// but it could also be some kind of attack, so discard addr for now and maybe revisit this later
				dout() << "Discarded port detect from " << fromaddr << " with addr " << ip << " which did not match expected address " << frompeer.visible_ipaddr.addr;
			}
		} else {
			dout() << "Discarded port detect from " << fromaddr << " with addr " << ip << " because existing addr and port were previously discovered as " << frompeer.visible_ipaddr;
		}
	}
}

void limitDebugLevel( int i )
{
 dout() << "Set debug level to: " << (DebugLevel)i;
 debugLevel = i;
}

int main()
{
  mp_init();

  //  const mp_complex I = mp_complex(0.0,1.0);

  static mp_real pi  = mppic;
  static mp_real pi2 = mppic * mppic;

  static mp_real Sqrtpi = sqrt(mppic);
  
  cout  << setiosflags(ios::uppercase);

  ofstream out("REGSOLC.out");
  out << setiosflags(ios::uppercase);

  ofstream kout("KERNELCMP.out");
  kout << setiosflags(ios::uppercase);
  
  ofstream dout("DATACMP.out");
  dout << setiosflags(ios::uppercase);
  
  int M, N;
  double bd;
  double x1d, x2d, xad, xbd;
  
  cin >> bd;
  cin >> xad >> xbd >> N;
  cin >> x1d >> x2d >> M;

  cout << "*** Precision = " << mpipl << " ***" << endl;
  cout << "b  = " << bd << endl;
  cout << "xa = " << xad << endl;
  cout << "xb = " << xbd << endl;
  cout << "x1 = " << x1d << endl;
  cout << "x2 = " << x2d << endl;
  cout << "N  = " << N << endl;
  cout << "M  = " << M << endl;

  mp_real xa   = mp_real(xad);
  mp_real xb   = mp_real(xbd);
  mp_real x1   = mp_real(x1d);
  mp_real x2   = mp_real(x2d);
  mp_real  b   = mp_real(bd);

  mp_real factor;
  factor = b/power(pi,mp_real(1.5))*exp(mp_real(1.25)*pi2*b*b);


  mp_real h    = (xb - xa)/mp_real(N);
  cout << "h  = " << h ;

  mp_real step = (x2 - x1)/mp_real(M);
  cout << "h2 = " << step;

  cout << "Faktor = " << factor << endl;


  mp_real mpx0; 
  //  mp_real mpx, mpx0; 
  //  mp_complex regkern, data;
  //  double Realregkern, Imaregkern;
  //  double RealData, ImaData;
  /*
  cout << " *** write Kernel and Data *** " << endl;
  
  for(int i = 0; i <= N; i++)
    {
      mpx0 = mp_real(0.0);
      mpx = xa + i * h;
      regkern = mpregkernc(mpx,mpx0,b);
      //      cout << mpx << regkern << endl;

      Realregkern = dble(MP_REAL(regkern));
      Imaregkern  = dble(aimag(regkern));

      //      data = E_data(mpx);
      //      cout << mpx << data << endl;

      //      RealData = dble(MP_REAL(data));
      //      ImaData  = dble(aimag(data));

      kout << dble(mpx) << '\t' << Realregkern << '\t' << Imaregkern << '\t' << bd<< endl;

      //      dout << dble(mpx) << '\t' << RealData << '\t' << ImaData << endl;
    }
    */
  mp_complex mpregsum;
  mp_real mpReg, RelError;
  double RealReg, ImaReg, Error;

  cout << " *** Begin: Regularisation *** " << endl;

  for(int j = 0; j <= M; j++)
    {
      mpx0 = x1 + j * step;

      mpregsum = mptrapez(RegIntegrand,mpx0,b,mpx0,.5*h,N);
      mpregsum *=factor;

//.........这里部分代码省略.........